MOM6
Data Types | Functions/Subroutines | Variables
mom_ice_shelf Module Reference

Detailed Description

Implements the thermodynamic aspects of ocean / ice-shelf interactions, along with a crude placeholder for a later implementation of full ice shelf dynamics, all using the MOM framework and coding style.

section_ICE_SHELF

This module implements the thermodynamic aspects of ocean/ice-shelf inter-actions using the MOM framework and coding style.

Derived from code by Chris Little, early 2010.

The ice-sheet dynamics subroutines do the following: initialize_shelf_mass - Initializes the ice shelf mass distribution.

  • Initializes h_shelf, h_mask, area_shelf_h
  • CURRENTLY: initializes mass_shelf as well, but this is unnecessary, as mass_shelf is initialized based on h_shelf and density_ice immediately afterwards. Possibly subroutine should be renamed update_shelf_mass - updates ice shelf mass via netCDF file USER_update_shelf_mass (TODO). solo_time_step - called only in ice-only mode. shelf_calc_flux - after melt rate & fluxes are calculated, ice dynamics are done. currently mass_shelf is updated immediately after ice_shelf_advect in fully dynamic mode.

NOTES: be aware that hmask(:,:) has a number of functions; it is used for front advancement, for subroutines in the velocity solve, and for thickness boundary conditions (this last one may be removed). in other words, interfering with its updates will have implications you might not expect.

Overall issues: Many variables need better documentation and units and the subgrid on which they are discretized.

ICE_SHELF equations

The three fundamental equations are: Heat flux

\[ \qquad \rho_w C_{pw} \gamma_T (T_w - T_b) = \rho_i \dot{m} L_f \]

Salt flux

\[ \qquad \rho_w \gamma_s (S_w - S_b) = \rho_i \dot{m} S_b \]

Freezing temperature

\[ \qquad T_b = a S_b + b + c P \]

where ....

References

Asay-Davis, Xylar S., Stephen L. Cornford, Benjamin K. Galton-Fenzi, Rupert M. Gladstone, G. Hilmar Gudmundsson, David M. Holland, Paul R. Holland, and Daniel F. Martin. Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: MISMIP v. 3 (MISMIP+), ISOMIP v. 2 (ISOMIP+) and MISOMIP v. 1 (MISOMIP1). Geoscientific Model Development 9, no. 7 (2016): 2471.

Goldberg, D. N., et al. Investigation of land ice-ocean interaction with a fully coupled ice-ocean model: 1. Model description and behavior. Journal of Geophysical Research: Earth Surface 117.F2 (2012).

Goldberg, D. N., et al. Investigation of land ice-ocean interaction with a fully coupled ice-ocean model: 2. Sensitivity to external forcings. Journal of Geophysical Research: Earth Surface 117.F2 (2012).

Holland, David M., and Adrian Jenkins. Modeling thermodynamic ice-ocean interactions at the base of an ice shelf. Journal of Physical Oceanography 29.8 (1999): 1787-1800.

Data Types

type  ice_shelf_cs
 Control structure that contains ice shelf parameters and diagnostics handles. More...
 

Functions/Subroutines

subroutine, public shelf_calc_flux (state, fluxes, Time, time_step, CS, forces)
 Calculates fluxes between the ocean and ice-shelf using the three-equations formulation (optional to use just two equations). See ICE_SHELF equations. More...
 
subroutine change_thickness_using_melt (ISS, G, time_step, fluxes, rho_ice, debug)
 Changes the thickness (mass) of the ice shelf based on sub-ice-shelf melting. More...
 
subroutine, public add_shelf_forces (G, US, CS, forces, do_shelf_area)
 This subroutine adds the mechanical forcing fields and perhaps shelf areas, based on the ice state in ice_shelf_CS. More...
 
subroutine add_shelf_pressure (G, US, CS, fluxes)
 This subroutine adds the ice shelf pressure to the fluxes type. More...
 
subroutine, public add_shelf_flux (G, US, CS, state, fluxes)
 Updates surface fluxes that are influenced by sub-ice-shelf melting. More...
 
subroutine, public initialize_ice_shelf (param_file, ocn_grid, Time, CS, diag, forces, fluxes, Time_in, solo_ice_sheet_in)
 Initializes shelf model data, parameters and diagnostics. More...
 
subroutine initialize_shelf_mass (G, param_file, CS, ISS, new_sim)
 Initializes shelf mass based on three options (file, zero and user) More...
 
subroutine update_shelf_mass (G, CS, ISS, Time)
 Updates the ice shelf mass using data from a file. More...
 
subroutine, public ice_shelf_save_restart (CS, Time, directory, time_stamped, filename_suffix)
 Save the ice shelf restart file. More...
 
subroutine, public ice_shelf_end (CS)
 Deallocates all memory associated with this module. More...
 
subroutine, public solo_time_step (CS, time_step, nsteps, Time, min_time_step_in)
 This routine is for stepping a stand-alone ice shelf model without an ocean. More...
 

Variables

integer id_clock_shelf
 CPU Clock for the ice shelf code. More...
 
integer id_clock_pass
 CPU Clock for group pass calls. More...
 

Function/Subroutine Documentation

◆ add_shelf_flux()

subroutine, public mom_ice_shelf::add_shelf_flux ( type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(ice_shelf_cs), pointer  CS,
type(surface), intent(inout)  state,
type(forcing), intent(inout)  fluxes 
)

Updates surface fluxes that are influenced by sub-ice-shelf melting.

Parameters
[in,out]gThe ocean's grid structure.
[in]usA dimensional unit scaling type
csThis module's control structure.
[in,out]stateSurface ocean state
[in,out]fluxesA structure of surface fluxes that may be used/updated.

Definition at line 870 of file MOM_ice_shelf.F90.

870  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
871  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
872  type(ice_shelf_CS), pointer :: CS !< This module's control structure.
873  type(surface), intent(inout) :: state!< Surface ocean state
874  type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be used/updated.
875 
876  ! local variables
877  real :: Irho0 !< The inverse of the mean density [m3 kg-1].
878  real :: frac_area !< The fractional area covered by the ice shelf [nondim].
879  real :: shelf_mass0 !< Total ice shelf mass at previous time (Time-dt).
880  real :: shelf_mass1 !< Total ice shelf mass at current time (Time).
881  real :: delta_mass_shelf!< Change in ice shelf mass over one time step [kg s-1]
882  real :: taux2, tauy2 !< The squared surface stresses [Pa].
883  real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [Pa].
884  real :: asu1, asu2 !< Ocean areas covered by ice shelves at neighboring u-
885  real :: asv1, asv2 !< and v-points [m2].
886  real :: fraz !< refreezing rate [kg m-2 s-1]
887  real :: mean_melt_flux !< spatial mean melt flux [kg s-1] or [kg m-2 s-1] at various points in the code.
888  real :: sponge_area !< total area of sponge region [m2]
889  real :: t0 !< The previous time (Time-dt) [s].
890  type(time_type) :: Time0!< The previous time (Time-dt)
891  real, dimension(SZDI_(G),SZDJ_(G)) :: last_mass_shelf !< Ice shelf mass
892  !! at at previous time (Time-dt) [kg m-2]
893  real, dimension(SZDI_(G),SZDJ_(G)) :: last_h_shelf !< Ice shelf thickness [Z ~> m]
894  !! at at previous time (Time-dt)
895  real, dimension(SZDI_(G),SZDJ_(G)) :: last_hmask !< Ice shelf mask
896  !! at at previous time (Time-dt)
897  real, dimension(SZDI_(G),SZDJ_(G)) :: last_area_shelf_h !< Ice shelf area [m2]
898  !! at at previous time (Time-dt)
899  type(ice_shelf_state), pointer :: ISS => null() !< A structure with elements that describe
900  !! the ice-shelf state
901 
902  real :: kv_rho_ice ! The viscosity of ice divided by its density [m5 kg-1 s-1]
903  real, parameter :: rho_fw = 1000.0 ! fresh water density
904  character(len=160) :: mesg ! The text of an error message
905  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
906  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
907  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
908 
909  if ((cs%grid%isc /= g%isc) .or. (cs%grid%iec /= g%iec) .or. &
910  (cs%grid%jsc /= g%jsc) .or. (cs%grid%jec /= g%jec)) &
911  call mom_error(fatal,"add_shelf_flux: Incompatible ocean and ice shelf grids.")
912 
913  iss => cs%ISS
914 
915  call add_shelf_pressure(g, us, cs, fluxes)
916 
917  ! Determine ustar and the square magnitude of the velocity in the
918  ! bottom boundary layer. Together these give the TKE source and
919  ! vertical decay scale.
920 
921  if (cs%debug) then
922  if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then
923  call uvchksum("tau[xy]_shelf", state%taux_shelf, state%tauy_shelf, &
924  g%HI, haloshift=0)
925  endif
926  endif
927 
928  if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then
929  call pass_vector(state%taux_shelf, state%tauy_shelf, g%domain, to_all, cgrid_ne)
930  ! GMM: melting is computed using ustar_shelf (and not ustar), which has already
931  ! been passed, I so believe we do not need to update fluxes%ustar.
932 ! Irho0 = 1.0 / CS%Rho0
933 ! do j=js,je ; do i=is,ie ; if (fluxes%frac_shelf_h(i,j) > 0.0) then
934  ! ### THIS SHOULD BE AN AREA WEIGHTED AVERAGE OF THE ustar_shelf POINTS.
935  ! taux2 = 0.0 ; tauy2 = 0.0
936  ! asu1 = (ISS%area_shelf_h(i-1,j) + ISS%area_shelf_h(i,j))
937  ! asu2 = (ISS%area_shelf_h(i,j) + ISS%area_shelf_h(i+1,j))
938  ! asv1 = (ISS%area_shelf_h(i,j-1) + ISS%area_shelf_h(i,j))
939  ! asv2 = (ISS%area_shelf_h(i,j) + ISS%area_shelf_h(i,j+1))
940  ! if ((asu1 + asu2 > 0.0) .and. associated(state%taux_shelf)) &
941  ! taux2 = (asu1 * state%taux_shelf(I-1,j)**2 + &
942  ! asu2 * state%taux_shelf(I,j)**2 ) / (asu1 + asu2)
943  ! if ((asv1 + asv2 > 0.0) .and. associated(state%tauy_shelf)) &
944  ! tauy2 = (asv1 * state%tauy_shelf(i,J-1)**2 + &
945  ! asv2 * state%tauy_shelf(i,J)**2 ) / (asv1 + asv2)
946 
947  ! fluxes%ustar(i,j) = MAX(CS%ustar_bg, US%m_to_Z*US%T_to_s*sqrt(Irho0 * sqrt(taux2 + tauy2)))
948 ! endif ; enddo ; enddo
949  endif
950 
951  if (cs%active_shelf_dynamics .or. cs%override_shelf_movement) then
952  do j=jsd,jed ; do i=isd,ied
953  if (g%areaT(i,j) > 0.0) &
954  fluxes%frac_shelf_h(i,j) = iss%area_shelf_h(i,j) * us%m_to_L**2*g%IareaT(i,j)
955  enddo ; enddo
956  endif
957 
958  do j=js,je ; do i=is,ie ; if (iss%area_shelf_h(i,j) > 0.0) then
959  frac_area = fluxes%frac_shelf_h(i,j) !### Should this be 1-frac_shelf_h?
960  if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0
961  if (associated(fluxes%sw_vis_dir)) fluxes%sw_vis_dir(i,j) = 0.0
962  if (associated(fluxes%sw_vis_dif)) fluxes%sw_vis_dif(i,j) = 0.0
963  if (associated(fluxes%sw_nir_dir)) fluxes%sw_nir_dir(i,j) = 0.0
964  if (associated(fluxes%sw_nir_dif)) fluxes%sw_nir_dif(i,j) = 0.0
965  if (associated(fluxes%lw)) fluxes%lw(i,j) = 0.0
966  if (associated(fluxes%latent)) fluxes%latent(i,j) = 0.0
967  if (associated(fluxes%evap)) fluxes%evap(i,j) = 0.0
968  if (associated(fluxes%lprec)) then
969  if (iss%water_flux(i,j) > 0.0) then
970  fluxes%lprec(i,j) = us%kg_m3_to_R*us%m_to_Z*us%T_to_s*frac_area*iss%water_flux(i,j)*cs%flux_factor
971  else
972  fluxes%lprec(i,j) = 0.0
973  fluxes%evap(i,j) = us%kg_m3_to_R*us%m_to_Z*us%T_to_s*frac_area*iss%water_flux(i,j)*cs%flux_factor
974  endif
975  endif
976 
977  if (associated(fluxes%sens)) &
978  fluxes%sens(i,j) = -frac_area*iss%tflux_ocn(i,j)*cs%flux_factor
979  if (associated(fluxes%salt_flux)) &
980  fluxes%salt_flux(i,j) = frac_area * iss%salt_flux(i,j)*cs%flux_factor
981  endif ; enddo ; enddo
982 
983  ! keep sea level constant by removing mass in the sponge
984  ! region (via virtual precip, vprec). Apply additional
985  ! salt/heat fluxes so that the resultant surface buoyancy
986  ! forcing is ~ 0.
987  ! This is needed for some of the ISOMIP+ experiments.
988 
989  if (cs%constant_sea_level) then
990  !### This code has lots of problems with hard coded constants and the use of
991  !### of non-reproducing sums. It needs to be refactored. -RWH
992 
993  if (.not. associated(fluxes%salt_flux)) allocate(fluxes%salt_flux(ie,je))
994  if (.not. associated(fluxes%vprec)) allocate(fluxes%vprec(ie,je))
995  fluxes%salt_flux(:,:) = 0.0 ; fluxes%vprec(:,:) = 0.0
996 
997  mean_melt_flux = 0.0; sponge_area = 0.0
998  do j=js,je ; do i=is,ie
999  frac_area = fluxes%frac_shelf_h(i,j)
1000  if (frac_area > 0.0) &
1001  mean_melt_flux = mean_melt_flux + (iss%water_flux(i,j)) * iss%area_shelf_h(i,j)
1002 
1003  !### These hard-coded limits need to be corrected. They are inappropriate here.
1004  if (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0) then
1005  sponge_area = sponge_area + us%L_to_m**2*g%areaT(i,j)
1006  endif
1007  enddo ; enddo
1008 
1009  ! take into account changes in mass (or thickness) when imposing ice shelf mass
1010  if (cs%override_shelf_movement .and. cs%mass_from_file) then
1011  t0 = time_type_to_real(cs%Time) - cs%time_step
1012 
1013  ! just compute changes in mass after first time step
1014  if (t0>0.0) then
1015  time0 = real_to_time(t0)
1016  last_hmask(:,:) = iss%hmask(:,:) ; last_area_shelf_h(:,:) = iss%area_shelf_h(:,:)
1017  call time_interp_external(cs%id_read_mass, time0, last_mass_shelf)
1018  last_h_shelf(:,:) = last_mass_shelf(:,:) / cs%rho_ice
1019 
1020  ! apply calving
1021  if (cs%min_thickness_simple_calve > 0.0) then
1022  call ice_shelf_min_thickness_calve(g, last_h_shelf, last_area_shelf_h, last_hmask, &
1023  cs%min_thickness_simple_calve)
1024  ! convert to mass again
1025  last_mass_shelf(:,:) = last_h_shelf(:,:) * cs%rho_ice
1026  endif
1027 
1028  shelf_mass0 = 0.0; shelf_mass1 = 0.0
1029  ! get total ice shelf mass at (Time-dt) and (Time), in kg
1030  do j=js,je ; do i=is,ie
1031  ! just floating shelf (0.1 is a threshold for min ocean thickness)
1032  if (((1.0/cs%density_ocean_avg)*state%ocean_mass(i,j) > 0.1) .and. &
1033  (iss%area_shelf_h(i,j) > 0.0)) then
1034  shelf_mass0 = shelf_mass0 + (last_mass_shelf(i,j) * iss%area_shelf_h(i,j))
1035  shelf_mass1 = shelf_mass1 + (iss%mass_shelf(i,j) * iss%area_shelf_h(i,j))
1036  endif
1037  enddo ; enddo
1038  call sum_across_pes(shelf_mass0); call sum_across_pes(shelf_mass1)
1039  delta_mass_shelf = (shelf_mass1 - shelf_mass0)/cs%time_step
1040 ! delta_mass_shelf = (shelf_mass1 - shelf_mass0)* &
1041 ! (rho_fw/CS%density_ice)/CS%time_step
1042 ! write(mesg,*)'delta_mass_shelf = ',delta_mass_shelf
1043 ! call MOM_mesg(mesg,5)
1044  else! first time step
1045  delta_mass_shelf = 0.0
1046  endif
1047  else ! ice shelf mass does not change
1048  delta_mass_shelf = 0.0
1049  endif
1050 
1051  call sum_across_pes(mean_melt_flux)
1052  call sum_across_pes(sponge_area)
1053 
1054  ! average total melt flux over sponge area
1055  mean_melt_flux = (mean_melt_flux+delta_mass_shelf) / sponge_area !kg/(m^2 s)
1056 
1057  ! apply fluxes
1058  do j=js,je ; do i=is,ie
1059  ! Note the following is hard coded for ISOMIP
1060  if (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0) then
1061  fluxes%vprec(i,j) = -mean_melt_flux * cs%density_ice/1000. ! evap is negative
1062  fluxes%sens(i,j) = fluxes%vprec(i,j) * cs%Cp * cs%T0 ! W /m^2
1063  ! Rescale fluxes%vprec to the proper units.
1064  fluxes%vprec(i,j) = us%kg_m3_to_R*us%m_to_Z*us%T_to_s * fluxes%vprec(i,j)
1065  fluxes%salt_flux(i,j) = fluxes%vprec(i,j) * cs%S0*1.0e-3 ! kg (salt)/(m^2 s)
1066  endif
1067  enddo ; enddo
1068 
1069  if (cs%debug) then
1070  write(mesg,*) 'Mean melt flux (kg/(m^2 s)), dt = ', mean_melt_flux, cs%time_step
1071  call mom_mesg(mesg)
1072  call mom_forcing_chksum("After constant sea level", fluxes, g, cs%US, haloshift=0)
1073  endif
1074 
1075  endif !constant_sea_level
1076 

References add_shelf_pressure(), mom_error_handler::mom_error(), mom_error_handler::mom_mesg(), mom_time_manager::real_to_time(), and mom_domains::to_all.

Referenced by shelf_calc_flux().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ add_shelf_forces()

subroutine, public mom_ice_shelf::add_shelf_forces ( type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(ice_shelf_cs), pointer  CS,
type(mech_forcing), intent(inout)  forces,
logical, intent(in), optional  do_shelf_area 
)

This subroutine adds the mechanical forcing fields and perhaps shelf areas, based on the ice state in ice_shelf_CS.

Parameters
[in,out]gThe ocean's grid structure.
[in]usA dimensional unit scaling type
csThis module's control structure.
[in,out]forcesA structure with the driving mechanical forces
[in]do_shelf_areaIf true find the shelf-covered areas.

Definition at line 760 of file MOM_ice_shelf.F90.

760  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
761  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
762  type(ice_shelf_CS), pointer :: CS !< This module's control structure.
763  type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
764  logical, optional, intent(in) :: do_shelf_area !< If true find the shelf-covered areas.
765 
766  real :: kv_rho_ice ! The viscosity of ice divided by its density [m5 kg-1 s-1].
767  real :: press_ice ! The pressure of the ice shelf per unit area of ocean (not ice) [Pa].
768  logical :: find_area ! If true find the shelf areas at u & v points.
769  type(ice_shelf_state), pointer :: ISS => null() ! A structure with elements that describe
770  ! the ice-shelf state
771 
772  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
773  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
774  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
775 
776  if ((cs%grid%isc /= g%isc) .or. (cs%grid%iec /= g%iec) .or. &
777  (cs%grid%jsc /= g%jsc) .or. (cs%grid%jec /= g%jec)) &
778  call mom_error(fatal,"add_shelf_forces: Incompatible ocean and ice shelf grids.")
779 
780  iss => cs%ISS
781 
782  find_area = .true. ; if (present(do_shelf_area)) find_area = do_shelf_area
783 
784  if (find_area) then
785  ! The frac_shelf is set over the widest possible area. Could it be smaller?
786  do j=jsd,jed ; do i=isd,ied-1
787  forces%frac_shelf_u(i,j) = 0.0
788  if ((g%areaT(i,j) + g%areaT(i+1,j) > 0.0)) & ! .and. (G%areaCu(I,j) > 0.0)) &
789  forces%frac_shelf_u(i,j) = ((iss%area_shelf_h(i,j) + iss%area_shelf_h(i+1,j)) / &
790  (us%L_to_m**2*g%areaT(i,j) + us%L_to_m**2*g%areaT(i+1,j)))
791  enddo ; enddo
792  do j=jsd,jed-1 ; do i=isd,ied
793  forces%frac_shelf_v(i,j) = 0.0
794  if ((g%areaT(i,j) + g%areaT(i,j+1) > 0.0)) & ! .and. (G%areaCv(i,J) > 0.0)) &
795  forces%frac_shelf_v(i,j) = ((iss%area_shelf_h(i,j) + iss%area_shelf_h(i,j+1)) / &
796  (us%L_to_m**2*g%areaT(i,j) + us%L_to_m**2*g%areaT(i,j+1)))
797  enddo ; enddo
798  call pass_vector(forces%frac_shelf_u, forces%frac_shelf_v, g%domain, to_all, cgrid_ne)
799  endif
800 
801  !### Consider working over a smaller array range.
802  do j=jsd,jed ; do i=isd,ied
803  press_ice = (iss%area_shelf_h(i,j) * us%m_to_L**2*g%IareaT(i,j)) * (cs%g_Earth * iss%mass_shelf(i,j))
804  if (associated(forces%p_surf)) then
805  if (.not.forces%accumulate_p_surf) forces%p_surf(i,j) = 0.0
806  forces%p_surf(i,j) = forces%p_surf(i,j) + press_ice
807  endif
808  if (associated(forces%p_surf_full)) then
809  if (.not.forces%accumulate_p_surf) forces%p_surf_full(i,j) = 0.0
810  forces%p_surf_full(i,j) = forces%p_surf_full(i,j) + press_ice
811  endif
812  enddo ; enddo
813 
814  ! For various reasons, forces%rigidity_ice_[uv] is always updated here. Note
815  ! that it may have been zeroed out where IOB is translated to forces and
816  ! contributions from icebergs and the sea-ice pack added subsequently.
817  !### THE RIGIDITY SHOULD ALSO INCORPORATE AREAL-COVERAGE INFORMATION.
818  kv_rho_ice = cs%kv_ice / cs%density_ice
819  do j=js,je ; do i=is-1,ie
820  if (.not.forces%accumulate_rigidity) forces%rigidity_ice_u(i,j) = 0.0
821  forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + &
822  kv_rho_ice * min(iss%mass_shelf(i,j), iss%mass_shelf(i+1,j))
823  enddo ; enddo
824  do j=js-1,je ; do i=is,ie
825  if (.not.forces%accumulate_rigidity) forces%rigidity_ice_v(i,j) = 0.0
826  forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + &
827  kv_rho_ice * min(iss%mass_shelf(i,j), iss%mass_shelf(i,j+1))
828  enddo ; enddo
829 
830  if (cs%debug) then
831  call uvchksum("rigidity_ice_[uv]", forces%rigidity_ice_u, forces%rigidity_ice_v, &
832  g%HI, symmetric=.true.)
833  call uvchksum("frac_shelf_[uv]", forces%frac_shelf_u, forces%frac_shelf_v, &
834  g%HI, symmetric=.true.)
835  endif
836 

References mom_error_handler::mom_error(), and mom_domains::to_all.

Referenced by initialize_ice_shelf(), mom_main(), shelf_calc_flux(), mom_ocean_model_mct::update_ocean_model(), and mom_ocean_model_nuopc::update_ocean_model().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ add_shelf_pressure()

subroutine mom_ice_shelf::add_shelf_pressure ( type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(ice_shelf_cs), intent(in)  CS,
type(forcing), intent(inout)  fluxes 
)
private

This subroutine adds the ice shelf pressure to the fluxes type.

Parameters
[in,out]gThe ocean's grid structure.
[in]usA dimensional unit scaling type
[in]csThis module's control structure.
[in,out]fluxesA structure of surface fluxes that may be updated.

Definition at line 841 of file MOM_ice_shelf.F90.

841  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
842  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
843  type(ice_shelf_CS), intent(in) :: CS !< This module's control structure.
844  type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be updated.
845 
846  real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [Pa].
847  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
848  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
849 
850  if ((cs%grid%isc /= g%isc) .or. (cs%grid%iec /= g%iec) .or. &
851  (cs%grid%jsc /= g%jsc) .or. (cs%grid%jec /= g%jec)) &
852  call mom_error(fatal,"add_shelf_pressure: Incompatible ocean and ice shelf grids.")
853 
854  do j=js,je ; do i=is,ie
855  press_ice = (cs%ISS%area_shelf_h(i,j) * us%m_to_L**2*g%IareaT(i,j)) * (cs%g_Earth * cs%ISS%mass_shelf(i,j))
856  if (associated(fluxes%p_surf)) then
857  if (.not.fluxes%accumulate_p_surf) fluxes%p_surf(i,j) = 0.0
858  fluxes%p_surf(i,j) = fluxes%p_surf(i,j) + press_ice
859  endif
860  if (associated(fluxes%p_surf_full)) then
861  if (.not.fluxes%accumulate_p_surf) fluxes%p_surf_full(i,j) = 0.0
862  fluxes%p_surf_full(i,j) = fluxes%p_surf_full(i,j) + press_ice
863  endif
864  enddo ; enddo
865 

References mom_error_handler::mom_error().

Referenced by add_shelf_flux(), and initialize_ice_shelf().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ change_thickness_using_melt()

subroutine mom_ice_shelf::change_thickness_using_melt ( type(ice_shelf_state), intent(inout)  ISS,
type(ocean_grid_type), intent(inout)  G,
real, intent(in)  time_step,
type(forcing), intent(inout)  fluxes,
real, intent(in)  rho_ice,
logical, intent(in), optional  debug 
)
private

Changes the thickness (mass) of the ice shelf based on sub-ice-shelf melting.

Parameters
[in,out]gThe ocean's grid structure.
[in,out]issA structure with elements that describe the ice-shelf state
[in]time_stepThe time step for this update [s].
[in,out]fluxesstructure containing pointers to any possible thermodynamic or mass-flux forcing fields.
[in]rho_iceThe density of ice-shelf ice [kg m-2 Z-1 ~> kg m-3].
[in]debugIf present and true, write chksums

Definition at line 707 of file MOM_ice_shelf.F90.

707  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
708  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
709  !! the ice-shelf state
710  real, intent(in) :: time_step !< The time step for this update [s].
711  type(forcing), intent(inout) :: fluxes !< structure containing pointers to any possible
712  !! thermodynamic or mass-flux forcing fields.
713  real, intent(in) :: rho_ice !< The density of ice-shelf ice [kg m-2 Z-1 ~> kg m-3].
714  logical, optional, intent(in) :: debug !< If present and true, write chksums
715 
716  ! locals
717  real :: I_rho_ice
718  integer :: i, j
719 
720  i_rho_ice = 1.0 / rho_ice
721 
722  do j=g%jsc,g%jec ; do i=g%isc,g%iec
723  if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
724  ! first, zero out fluxes applied during previous time step
725  if (associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0
726  if (associated(fluxes%sens)) fluxes%sens(i,j) = 0.0
727  if (associated(fluxes%frac_shelf_h)) fluxes%frac_shelf_h(i,j) = 0.0
728  if (associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = 0.0
729 
730  if (iss%water_flux(i,j) / rho_ice * time_step < iss%h_shelf(i,j)) then
731  iss%h_shelf(i,j) = iss%h_shelf(i,j) - iss%water_flux(i,j) / rho_ice * time_step
732  else
733  ! the ice is about to melt away, so set thickness, area, and mask to zero
734  ! NOTE: this is not mass conservative should maybe scale salt & heat flux for this cell
735  iss%h_shelf(i,j) = 0.0
736  iss%hmask(i,j) = 0.0
737  iss%area_shelf_h(i,j) = 0.0
738  endif
739  endif
740  enddo ; enddo
741 
742  call pass_var(iss%area_shelf_h, g%domain)
743  call pass_var(iss%h_shelf, g%domain)
744  call pass_var(iss%hmask, g%domain)
745 
746  !### combine this with the loops above.
747  do j=g%jsd,g%jed ; do i=g%isd,g%ied
748  if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
749  iss%mass_shelf(i,j) = iss%h_shelf(i,j)*rho_ice
750  endif
751  enddo ; enddo
752 
753  call pass_var(iss%mass_shelf, g%domain)
754 

Referenced by shelf_calc_flux().

Here is the caller graph for this function:

◆ ice_shelf_end()

subroutine, public mom_ice_shelf::ice_shelf_end ( type(ice_shelf_cs), pointer  CS)

Deallocates all memory associated with this module.

Parameters
csA pointer to the ice shelf control structure

Definition at line 1744 of file MOM_ice_shelf.F90.

1744  type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure
1745 
1746  if (.not.associated(cs)) return
1747 
1748  call ice_shelf_state_end(cs%ISS)
1749 
1750  if (cs%active_shelf_dynamics) call ice_shelf_dyn_end(cs%dCS)
1751 
1752  deallocate(cs)
1753 

References mom_ice_shelf_dynamics::ice_shelf_dyn_end(), and mom_ice_shelf_state::ice_shelf_state_end().

Referenced by mom_main(), mom_ocean_model_nuopc::ocean_model_end(), and mom_ocean_model_mct::ocean_model_end().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ice_shelf_save_restart()

subroutine, public mom_ice_shelf::ice_shelf_save_restart ( type(ice_shelf_cs), pointer  CS,
type(time_type), intent(in)  Time,
character(len=*), intent(in), optional  directory,
logical, intent(in), optional  time_stamped,
character(len=*), intent(in), optional  filename_suffix 
)

Save the ice shelf restart file.

Parameters
csice shelf control structure
[in]timemodel time at this call
[in]directoryAn optional directory into which to write these restart files.
[in]time_stampedf true, the restart file names include a unique time stamp. The default is false.
[in]filename_suffixAn optional suffix (e.g., a time-stamp) to append to the restart file names.

Definition at line 1721 of file MOM_ice_shelf.F90.

1721  type(ice_shelf_CS), pointer :: CS !< ice shelf control structure
1722  type(time_type), intent(in) :: Time !< model time at this call
1723  character(len=*), optional, intent(in) :: directory !< An optional directory into which to write
1724  !! these restart files.
1725  logical, optional, intent(in) :: time_stamped !< f true, the restart file names include
1726  !! a unique time stamp. The default is false.
1727  character(len=*), optional, intent(in) :: filename_suffix !< An optional suffix (e.g., a
1728  !! time-stamp) to append to the restart file names.
1729  ! local variables
1730  type(ocean_grid_type), pointer :: G => null()
1731  character(len=200) :: restart_dir
1732 
1733  g => cs%grid
1734 
1735  if (present(directory)) then ; restart_dir = directory
1736  else ; restart_dir = cs%restart_output_dir ; endif
1737 
1738  call save_restart(restart_dir, time, cs%grid, cs%restart_CSp, time_stamped)
1739 

Referenced by mom_main(), mom_ocean_model_nuopc::ocean_model_restart(), mom_ocean_model_mct::ocean_model_restart(), mom_ocean_model_nuopc::ocean_model_save_restart(), mom_ocean_model_mct::ocean_model_save_restart(), and ocn_comp_mct::ocn_run_mct().

Here is the caller graph for this function:

◆ initialize_ice_shelf()

subroutine, public mom_ice_shelf::initialize_ice_shelf ( type(param_file_type), intent(in)  param_file,
type(ocean_grid_type), pointer  ocn_grid,
type(time_type), intent(inout)  Time,
type(ice_shelf_cs), pointer  CS,
type(diag_ctrl), intent(in), target  diag,
type(mech_forcing), intent(inout), optional  forces,
type(forcing), intent(inout), optional  fluxes,
type(time_type), intent(in), optional  Time_in,
logical, intent(in), optional  solo_ice_sheet_in 
)

Initializes shelf model data, parameters and diagnostics.

Parameters
[in]param_fileA structure to parse for run-time parameters
ocn_gridThe calling ocean model's horizontal grid structure
[in,out]timeThe clock that that will indicate the model time
csA pointer to the ice shelf control structure
[in]diagA structure that is used to regulate the diagnostic output.
[in,out]fluxesA structure containing pointers to any possible thermodynamic or mass-flux forcing fields.
[in,out]forcesA structure with the driving mechanical forces
[in]time_inThe time at initialization.
[in]solo_ice_sheet_inIf present, this indicates whether a solo ice-sheet driver.

Definition at line 1082 of file MOM_ice_shelf.F90.

1082  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1083  type(ocean_grid_type), pointer :: ocn_grid !< The calling ocean model's horizontal grid structure
1084  type(time_type), intent(inout) :: Time !< The clock that that will indicate the model time
1085  type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure
1086  type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate the diagnostic output.
1087  type(forcing), optional, intent(inout) :: fluxes !< A structure containing pointers to any possible
1088  !! thermodynamic or mass-flux forcing fields.
1089  type(mech_forcing), optional, intent(inout) :: forces !< A structure with the driving mechanical forces
1090  type(time_type), optional, intent(in) :: Time_in !< The time at initialization.
1091  logical, optional, intent(in) :: solo_ice_sheet_in !< If present, this indicates whether
1092  !! a solo ice-sheet driver.
1093 
1094  type(ocean_grid_type), pointer :: G => null(), og => null() ! Pointers to grids for convenience.
1095  type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
1096  ! various unit conversion factors
1097  type(ice_shelf_state), pointer :: ISS => null() !< A structure with elements that describe
1098  !! the ice-shelf state
1099  type(directories) :: dirs
1100  type(dyn_horgrid_type), pointer :: dG => null()
1101  real :: Z_rescale ! A rescaling factor for heights from the representation in
1102  ! a reastart fole to the internal representation in this run.
1103  real :: cdrag, drag_bg_vel
1104  logical :: new_sim, save_IC, var_force
1105  !This include declares and sets the variable "version".
1106 #include "version_variable.h"
1107  character(len=200) :: config
1108  character(len=200) :: IC_file,filename,inputdir
1109  character(len=40) :: mdl = "MOM_ice_shelf" ! This module's name.
1110  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq
1111  integer :: wd_halos(2)
1112  logical :: read_TideAmp, shelf_mass_is_dynamic, debug
1113  character(len=240) :: Tideamp_file
1114  real :: utide
1115  if (associated(cs)) then
1116  call mom_error(fatal, "MOM_ice_shelf.F90, initialize_ice_shelf: "// &
1117  "called with an associated control structure.")
1118  return
1119  endif
1120  allocate(cs)
1121 
1122  ! Go through all of the infrastructure initialization calls, since this is
1123  ! being treated as an independent component that just happens to use the
1124  ! MOM's grid and infrastructure.
1125  call get_mom_input(dirs=dirs)
1126 
1127  ! Determining the internal unit scaling factors for this run.
1128  call unit_scaling_init(param_file, cs%US)
1129 
1130  ! Set up the ice-shelf domain and grid
1131  wd_halos(:)=0
1132  call mom_domains_init(cs%grid%domain, param_file, min_halo=wd_halos, symmetric=grid_sym_)
1133  ! call diag_mediator_init(CS%grid,param_file,CS%diag)
1134  ! this needs to be fixed - will probably break when not using coupled driver 0
1135  call mom_grid_init(cs%grid, param_file, cs%US)
1136 
1137  call create_dyn_horgrid(dg, cs%grid%HI)
1138  call clone_mom_domain(cs%grid%Domain, dg%Domain)
1139 
1140  call set_grid_metrics(dg, param_file, cs%US)
1141  ! call set_diag_mediator_grid(CS%grid, CS%diag)
1142 
1143  ! The ocean grid possibly uses different symmetry.
1144  if (associated(ocn_grid)) then ; cs%ocn_grid => ocn_grid
1145  else ; cs%ocn_grid => cs%grid ; endif
1146 
1147  ! Convenience pointers
1148  g => cs%grid
1149  og => cs%ocn_grid
1150  us => cs%US
1151 
1152  if (is_root_pe()) then
1153  write(0,*) 'OG: ', og%isd, og%isc, og%iec, og%ied, og%jsd, og%jsc, og%jsd, og%jed
1154  write(0,*) 'IG: ', g%isd, g%isc, g%iec, g%ied, g%jsd, g%jsc, g%jsd, g%jed
1155  endif
1156 
1157  cs%Time = time ! ### This might not be in the right place?
1158  cs%diag => diag
1159 
1160  ! Are we being called from the solo ice-sheet driver? When called by the ocean
1161  ! model solo_ice_sheet_in is not preset.
1162  cs%solo_ice_sheet = .false.
1163  if (present(solo_ice_sheet_in)) cs%solo_ice_sheet = solo_ice_sheet_in
1164 
1165  if (present(time_in)) time = time_in
1166 
1167  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1168  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
1169  isdq = g%IsdB ; iedq = g%IedB ; jsdq = g%JsdB ; jedq = g%JedB
1170 
1171  cs%Lat_fusion = 3.34e5
1172  cs%override_shelf_movement = .false. ; cs%active_shelf_dynamics = .false.
1173 
1174  call log_version(param_file, mdl, version, "")
1175  call get_param(param_file, mdl, "DEBUG", debug, default=.false.)
1176  call get_param(param_file, mdl, "DEBUG_IS", cs%debug, &
1177  "If true, write verbose debugging messages for the ice shelf.", &
1178  default=debug)
1179  call get_param(param_file, mdl, "DYNAMIC_SHELF_MASS", shelf_mass_is_dynamic, &
1180  "If true, the ice sheet mass can evolve with time.", &
1181  default=.false.)
1182  if (shelf_mass_is_dynamic) then
1183  call get_param(param_file, mdl, "OVERRIDE_SHELF_MOVEMENT", cs%override_shelf_movement, &
1184  "If true, user provided code specifies the ice-shelf "//&
1185  "movement instead of the dynamic ice model.", default=.false.)
1186  cs%active_shelf_dynamics = .not.cs%override_shelf_movement
1187  call get_param(param_file, mdl, "GROUNDING_LINE_INTERPOLATE", cs%GL_regularize, &
1188  "If true, regularize the floatation condition at the "//&
1189  "grounding line as in Goldberg Holland Schoof 2009.", default=.false.)
1190  call get_param(param_file, mdl, "GROUNDING_LINE_COUPLE", cs%GL_couple, &
1191  "If true, let the floatation condition be determined by "//&
1192  "ocean column thickness. This means that update_OD_ffrac "//&
1193  "will be called. GL_REGULARIZE and GL_COUPLE are exclusive.", &
1194  default=.false., do_not_log=cs%GL_regularize)
1195  if (cs%GL_regularize) cs%GL_couple = .false.
1196  endif
1197 
1198  call get_param(param_file, mdl, "SHELF_THERMO", cs%isthermo, &
1199  "If true, use a thermodynamically interactive ice shelf.", &
1200  default=.false.)
1201  call get_param(param_file, mdl, "SHELF_THREE_EQN", cs%threeeq, &
1202  "If true, use the three equation expression of "//&
1203  "consistency to calculate the fluxes at the ice-ocean "//&
1204  "interface.", default=.true.)
1205  call get_param(param_file, mdl, "SHELF_INSULATOR", cs%insulator, &
1206  "If true, the ice shelf is a perfect insulatior "//&
1207  "(no conduction).", default=.false.)
1208  call get_param(param_file, mdl, "MELTING_CUTOFF_DEPTH", cs%cutoff_depth, &
1209  "Depth above which the melt is set to zero (it must be >= 0) "//&
1210  "Default value won't affect the solution.", default=0.0)
1211  if (cs%cutoff_depth < 0.) &
1212  call mom_error(warning,"Initialize_ice_shelf: MELTING_CUTOFF_DEPTH must be >= 0.")
1213 
1214  call get_param(param_file, mdl, "CONST_SEA_LEVEL", cs%constant_sea_level, &
1215  "If true, apply evaporative, heat and salt fluxes in "//&
1216  "the sponge region. This will avoid a large increase "//&
1217  "in sea level. This option is needed for some of the "//&
1218  "ISOMIP+ experiments (Ocean3 and Ocean4). "//&
1219  "IMPORTANT: it is not currently possible to do "//&
1220  "prefect restarts using this flag.", default=.false.)
1221 
1222  call get_param(param_file, mdl, "ISOMIP_S_SUR_SPONGE", &
1223  cs%S0, "Surface salinity in the resoring region.", &
1224  default=33.8, do_not_log=.true.)
1225 
1226  call get_param(param_file, mdl, "ISOMIP_T_SUR_SPONGE", &
1227  cs%T0, "Surface temperature in the resoring region.", &
1228  default=-1.9, do_not_log=.true.)
1229 
1230  call get_param(param_file, mdl, "SHELF_3EQ_GAMMA", cs%const_gamma, &
1231  "If true, user specifies a constant nondimensional heat-transfer coefficient "//&
1232  "(GAMMA_T_3EQ), from which the salt-transfer coefficient is then computed "//&
1233  " as GAMMA_T_3EQ/35. This is used with SHELF_THREE_EQN.", default=.false.)
1234  if (cs%const_gamma) call get_param(param_file, mdl, "SHELF_3EQ_GAMMA_T", cs%Gamma_T_3EQ, &
1235  "Nondimensional heat-transfer coefficient.",default=2.2e-2, &
1236  units="nondim.", fail_if_missing=.true.)
1237 
1238  call get_param(param_file, mdl, "ICE_SHELF_MASS_FROM_FILE", &
1239  cs%mass_from_file, "Read the mass of the "//&
1240  "ice shelf (every time step) from a file.", default=.false.)
1241 
1242  if (cs%threeeq) &
1243  call get_param(param_file, mdl, "SHELF_S_ROOT", cs%find_salt_root, &
1244  "If SHELF_S_ROOT = True, salinity at the ice/ocean interface (Sbdry) "//&
1245  "is computed from a quadratic equation. Otherwise, the previous "//&
1246  "interactive method to estimate Sbdry is used.", default=.false.)
1247  if (cs%find_salt_root) then ! read liquidus coeffs.
1248  call get_param(param_file, mdl, "TFREEZE_S0_P0",cs%lambda1, &
1249  "this is the freezing potential temperature at "//&
1250  "S=0, P=0.", units="degC", default=0.0, do_not_log=.true.)
1251  call get_param(param_file, mdl, "DTFREEZE_DS",cs%lambda1, &
1252  "this is the derivative of the freezing potential "//&
1253  "temperature with salinity.", &
1254  units="degC psu-1", default=-0.054, do_not_log=.true.)
1255  call get_param(param_file, mdl, "DTFREEZE_DP",cs%lambda3, &
1256  "this is the derivative of the freezing potential "//&
1257  "temperature with pressure.", &
1258  units="degC Pa-1", default=0.0, do_not_log=.true.)
1259 
1260  endif
1261 
1262  if (.not.cs%threeeq) &
1263  call get_param(param_file, mdl, "SHELF_2EQ_GAMMA_T", cs%gamma_t, &
1264  "If SHELF_THREE_EQN is false, this the fixed turbulent "//&
1265  "exchange velocity at the ice-ocean interface.", &
1266  units="m s-1", fail_if_missing=.true.)
1267 
1268  call get_param(param_file, mdl, "G_EARTH", cs%g_Earth, &
1269  "The gravitational acceleration of the Earth.", &
1270  units="m s-2", default = 9.80)
1271  call get_param(param_file, mdl, "C_P", cs%Cp, &
1272  "The heat capacity of sea water.", units="J kg-1 K-1", &
1273  fail_if_missing=.true.)
1274  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
1275  "The mean ocean density used with BOUSSINESQ true to "//&
1276  "calculate accelerations and the mass for conservation "//&
1277  "properties, or with BOUSSINSEQ false to convert some "//&
1278  "parameters from vertical units of m to kg m-2.", &
1279  units="kg m-3", default=1035.0) !### MAKE THIS A SEPARATE PARAMETER.
1280  call get_param(param_file, mdl, "C_P_ICE", cs%Cp_ice, &
1281  "The heat capacity of ice.", units="J kg-1 K-1", &
1282  default=2.10e3)
1283 
1284  call get_param(param_file, mdl, "ICE_SHELF_FLUX_FACTOR", cs%flux_factor, &
1285  "Non-dimensional factor applied to shelf thermodynamic "//&
1286  "fluxes.", units="none", default=1.0)
1287 
1288  call get_param(param_file, mdl, "KV_ICE", cs%kv_ice, &
1289  "The viscosity of the ice.", units="m2 s-1", default=1.0e10)
1290  call get_param(param_file, mdl, "KV_MOLECULAR", cs%kv_molec, &
1291  "The molecular kinimatic viscosity of sea water at the "//&
1292  "freezing temperature.", units="m2 s-1", default=1.95e-6)
1293  call get_param(param_file, mdl, "ICE_SHELF_SALINITY", cs%Salin_ice, &
1294  "The salinity of the ice inside the ice shelf.", units="psu", &
1295  default=0.0)
1296  call get_param(param_file, mdl, "ICE_SHELF_TEMPERATURE", cs%Temp_ice, &
1297  "The temperature at the center of the ice shelf.", &
1298  units = "degC", default=-15.0)
1299  call get_param(param_file, mdl, "KD_SALT_MOLECULAR", cs%kd_molec_salt, &
1300  "The molecular diffusivity of salt in sea water at the "//&
1301  "freezing point.", units="m2 s-1", default=8.02e-10)
1302  call get_param(param_file, mdl, "KD_TEMP_MOLECULAR", cs%kd_molec_temp, &
1303  "The molecular diffusivity of heat in sea water at the "//&
1304  "freezing point.", units="m2 s-1", default=1.41e-7)
1305  call get_param(param_file, mdl, "RHO_0", cs%density_ocean_avg, &
1306  "avg ocean density used in floatation cond", &
1307  units="kg m-3", default=1035.)
1308  call get_param(param_file, mdl, "DT_FORCING", cs%time_step, &
1309  "The time step for changing forcing, coupling with other "//&
1310  "components, or potentially writing certain diagnostics. "//&
1311  "The default value is given by DT.", units="s", default=0.0)
1312 
1313  call get_param(param_file, mdl, "COL_THICK_MELT_THRESHOLD", cs%col_thick_melt_threshold, &
1314  "The minimum ML thickness where melting is allowed.", units="m", &
1315  default=0.0)
1316 
1317  call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, &
1318  "If true, read a file (given by TIDEAMP_FILE) containing "//&
1319  "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
1320 
1321  call safe_alloc_ptr(cs%utide,isd,ied,jsd,jed) ; cs%utide(:,:) = 0.0
1322 
1323  if (read_tideamp) then
1324  call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
1325  "The path to the file containing the spatially varying "//&
1326  "tidal amplitudes.", &
1327  default="tideamp.nc")
1328  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1329  inputdir = slasher(inputdir)
1330  tideamp_file = trim(inputdir) // trim(tideamp_file)
1331  call mom_read_data(tideamp_file,'tideamp',cs%utide,g%domain,timelevel=1)
1332  else
1333  call get_param(param_file, mdl, "UTIDE", utide, &
1334  "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
1335  units="m s-1", default=0.0)
1336  cs%utide(:,:) = utide
1337  endif
1338 
1339  call eos_init(param_file, cs%eqn_of_state)
1340 
1341  !! new parameters that need to be in MOM_input
1342 
1343  if (cs%active_shelf_dynamics) then
1344 
1345  call get_param(param_file, mdl, "DENSITY_ICE", cs%density_ice, &
1346  "A typical density of ice.", units="kg m-3", default=917.0)
1347 
1348  call get_param(param_file, mdl, "INPUT_FLUX_ICE_SHELF", cs%input_flux, &
1349  "volume flux at upstream boundary", units="m2 s-1", default=0.)
1350  call get_param(param_file, mdl, "INPUT_THICK_ICE_SHELF", cs%input_thickness, &
1351  "flux thickness at upstream boundary", units="m", default=1000.)
1352  else
1353  ! This is here because of inconsistent defaults. I don't know why. RWH
1354  call get_param(param_file, mdl, "DENSITY_ICE", cs%density_ice, &
1355  "A typical density of ice.", units="kg m-3", default=900.0)
1356  endif
1357  cs%rho_ice = cs%density_ice*us%Z_to_m
1358  call get_param(param_file, mdl, "MIN_THICKNESS_SIMPLE_CALVE", &
1359  cs%min_thickness_simple_calve, &
1360  "Min thickness rule for the very simple calving law",&
1361  units="m", default=0.0, scale=us%m_to_Z)
1362 
1363  call get_param(param_file, mdl, "USTAR_SHELF_BG", cs%ustar_bg, &
1364  "The minimum value of ustar under ice sheves.", &
1365  units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
1366  call get_param(param_file, mdl, "CDRAG_SHELF", cdrag, &
1367  "CDRAG is the drag coefficient relating the magnitude of "//&
1368  "the velocity field to the surface stress.", units="nondim", &
1369  default=0.003)
1370  cs%cdrag = cdrag
1371  if (cs%ustar_bg <= 0.0) then
1372  call get_param(param_file, mdl, "DRAG_BG_VEL_SHELF", drag_bg_vel, &
1373  "DRAG_BG_VEL is either the assumed bottom velocity (with "//&
1374  "LINEAR_DRAG) or an unresolved velocity that is "//&
1375  "combined with the resolved velocity to estimate the "//&
1376  "velocity magnitude.", units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
1377  if (cs%cdrag*drag_bg_vel > 0.0) cs%ustar_bg = sqrt(cs%cdrag)*drag_bg_vel
1378  endif
1379 
1380  ! Allocate and initialize state variables to default values
1381  call ice_shelf_state_init(cs%ISS, cs%grid)
1382  iss => cs%ISS
1383 
1384  ! Allocate the arrays for passing ice-shelf data through the forcing type.
1385  if (.not. cs%solo_ice_sheet) then
1386  call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes.")
1387  ! GMM: the following assures that water/heat fluxes are just allocated
1388  ! when SHELF_THERMO = True. These fluxes are necessary if one wants to
1389  ! use either ENERGETICS_SFC_PBL (ALE mode) or BULKMIXEDLAYER (layer mode).
1390  if (present(fluxes)) &
1391  call allocate_forcing_type(cs%ocn_grid, fluxes, ustar=.true., shelf=.true., &
1392  press=.true., water=cs%isthermo, heat=cs%isthermo)
1393  if (present(forces)) &
1394  call allocate_mech_forcing(cs%ocn_grid, forces, ustar=.true., shelf=.true., press=.true.)
1395  else
1396  call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes in solo mode.")
1397  if (present(fluxes)) &
1398  call allocate_forcing_type(g, fluxes, ustar=.true., shelf=.true., press=.true.)
1399  if (present(forces)) &
1400  call allocate_mech_forcing(g, forces, ustar=.true., shelf=.true., press=.true.)
1401  endif
1402 
1403  ! Set up the bottom depth, G%D either analytically or from file
1404  call mom_initialize_topography(dg%bathyT, g%max_depth, dg, param_file)
1405  call rescale_dyn_horgrid_bathymetry(dg, us%Z_to_m)
1406 
1407  ! Set up the Coriolis parameter, G%f, usually analytically.
1408  call mom_initialize_rotation(dg%CoriolisBu, dg, param_file, us)
1409  ! This copies grid elements, including bathyT and CoriolisBu from dG to CS%grid.
1410  call copy_dyngrid_to_mom_grid(dg, cs%grid, us)
1411 
1412  call destroy_dyn_horgrid(dg)
1413 
1414  ! Set up the restarts.
1415  call restart_init(param_file, cs%restart_CSp, "Shelf.res")
1416  call register_restart_field(iss%mass_shelf, "shelf_mass", .true., cs%restart_CSp, &
1417  "Ice shelf mass", "kg m-2")
1418  call register_restart_field(iss%area_shelf_h, "shelf_area", .true., cs%restart_CSp, &
1419  "Ice shelf area in cell", "m2")
1420  call register_restart_field(iss%h_shelf, "h_shelf", .true., cs%restart_CSp, &
1421  "ice sheet/shelf thickness", "m")
1422  call register_restart_field(us%m_to_Z_restart, "m_to_Z", .false., cs%restart_CSp, &
1423  "Height unit conversion factor", "Z meter-1")
1424  if (cs%active_shelf_dynamics) then
1425  call register_restart_field(iss%hmask, "h_mask", .true., cs%restart_CSp, &
1426  "ice sheet/shelf thickness mask" ,"none")
1427  endif
1428 
1429  ! if (CS%active_shelf_dynamics) then !### Consider adding an ice shelf dynamics switch.
1430  ! Allocate CS%dCS and specify additional restarts for ice shelf dynamics
1431  call register_ice_shelf_dyn_restarts(g, param_file, cs%dCS, cs%restart_CSp)
1432  ! endif
1433 
1434  !GMM - I think we do not need to save ustar_shelf and iceshelf_melt in the restart file
1435  !if (.not. CS%solo_ice_sheet) then
1436  ! call register_restart_field(fluxes%ustar_shelf, "ustar_shelf", .false., CS%restart_CSp, &
1437  ! "Friction velocity under ice shelves", "m s-1")
1438  ! call register_restart_field(fluxes%iceshelf_melt, "iceshelf_melt", .false., CS%restart_CSp, &
1439  ! "Ice Shelf Melt Rate", "m year-1")
1440  !endif
1441 
1442  cs%restart_output_dir = dirs%restart_output_dir
1443 
1444  new_sim = .false.
1445  if ((dirs%input_filename(1:1) == 'n') .and. &
1446  (len_trim(dirs%input_filename) == 1)) new_sim = .true.
1447 
1448  if (cs%override_shelf_movement .and. cs%mass_from_file) then
1449 
1450  ! initialize the ids for reading shelf mass from a netCDF
1451  call initialize_shelf_mass(g, param_file, cs, iss)
1452 
1453  if (new_sim) then
1454  ! new simulation, initialize ice thickness as in the static case
1455  call initialize_ice_thickness(iss%h_shelf, iss%area_shelf_h, iss%hmask, g, us, param_file)
1456 
1457  ! next make sure mass is consistent with thickness
1458  do j=g%jsd,g%jed ; do i=g%isd,g%ied
1459  if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
1460  iss%mass_shelf(i,j) = iss%h_shelf(i,j)*cs%rho_ice
1461  endif
1462  enddo ; enddo
1463 
1464  if (cs%min_thickness_simple_calve > 0.0) &
1465  call ice_shelf_min_thickness_calve(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, &
1466  cs%min_thickness_simple_calve)
1467  endif
1468  endif
1469 
1470  if (cs%active_shelf_dynamics) then
1471  ! the only reason to initialize boundary conds is if the shelf is dynamic - MJH
1472 
1473  ! call initialize_ice_shelf_boundary ( CS%u_face_mask_bdry, CS%v_face_mask_bdry, &
1474  ! CS%u_flux_bdry_val, CS%v_flux_bdry_val, &
1475  ! CS%u_bdry_val, CS%v_bdry_val, CS%h_bdry_val, &
1476  ! ISS%hmask, G, param_file)
1477 
1478  endif
1479 
1480  if (new_sim .and. (.not. (cs%override_shelf_movement .and. cs%mass_from_file))) then
1481 
1482  ! This model is initialized internally or from a file.
1483  call initialize_ice_thickness(iss%h_shelf, iss%area_shelf_h, iss%hmask, g, us, param_file)
1484 
1485  ! next make sure mass is consistent with thickness
1486  do j=g%jsd,g%jed ; do i=g%isd,g%ied
1487  if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
1488  iss%mass_shelf(i,j) = iss%h_shelf(i,j)*cs%rho_ice
1489  endif
1490  enddo ; enddo
1491 
1492  ! else ! Previous block for new_sim=.T., this block restores the state.
1493  elseif (.not.new_sim) then
1494  ! This line calls a subroutine that reads the initial conditions from a restart file.
1495  call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: Restoring ice shelf from file.")
1496  call restore_state(dirs%input_filename, dirs%restart_input_dir, time, &
1497  g, cs%restart_CSp)
1498 
1499  if ((us%m_to_Z_restart /= 0.0) .and. (us%m_to_Z_restart /= us%m_to_Z)) then
1500  z_rescale = us%m_to_Z / us%m_to_Z_restart
1501  do j=g%jsc,g%jec ; do i=g%isc,g%iec
1502  iss%h_shelf(i,j) = z_rescale * iss%h_shelf(i,j)
1503  enddo ; enddo
1504  endif
1505 
1506  endif ! .not. new_sim
1507 
1508  cs%Time = time
1509 
1510  call cpu_clock_begin(id_clock_pass)
1511  call pass_var(iss%area_shelf_h, g%domain)
1512  call pass_var(iss%h_shelf, g%domain)
1513  call pass_var(iss%mass_shelf, g%domain)
1514  call pass_var(iss%hmask, g%domain)
1515  call pass_var(g%bathyT, g%domain)
1516  call cpu_clock_end(id_clock_pass)
1517 
1518  do j=jsd,jed ; do i=isd,ied
1519  if (iss%area_shelf_h(i,j) > us%L_to_m**2*g%areaT(i,j)) then
1520  call mom_error(warning,"Initialize_ice_shelf: area_shelf_h exceeds G%areaT.")
1521  iss%area_shelf_h(i,j) = us%L_to_m**2*g%areaT(i,j)
1522  endif
1523  enddo ; enddo
1524  if (present(fluxes)) then ; do j=jsd,jed ; do i=isd,ied
1525  if (g%areaT(i,j) > 0.0) fluxes%frac_shelf_h(i,j) = iss%area_shelf_h(i,j) / (us%L_to_m**2*g%areaT(i,j))
1526  enddo ; enddo ; endif
1527 
1528  if (cs%debug) then
1529  call hchksum(fluxes%frac_shelf_h, "IS init: frac_shelf_h", g%HI, haloshift=0)
1530  endif
1531 
1532  if (present(forces)) &
1533  call add_shelf_forces(g, us, cs, forces, do_shelf_area=.not.cs%solo_ice_sheet)
1534 
1535  if (present(fluxes)) call add_shelf_pressure(g, us, cs, fluxes)
1536 
1537  if (cs%active_shelf_dynamics .and. .not.cs%isthermo) then
1538  iss%water_flux(:,:) = 0.0
1539  endif
1540 
1541  if (shelf_mass_is_dynamic) &
1542  call initialize_ice_shelf_dyn(param_file, time, iss, cs%dCS, g, us, diag, new_sim, solo_ice_sheet_in)
1543 
1544  call fix_restart_unit_scaling(us)
1545 
1546  call get_param(param_file, mdl, "SAVE_INITIAL_CONDS", save_ic, &
1547  "If true, save the ice shelf initial conditions.", &
1548  default=.false.)
1549  if (save_ic) call get_param(param_file, mdl, "SHELF_IC_OUTPUT_FILE", ic_file,&
1550  "The name-root of the output file for the ice shelf "//&
1551  "initial conditions.", default="MOM_Shelf_IC")
1552 
1553  if (save_ic .and. .not.((dirs%input_filename(1:1) == 'r') .and. &
1554  (len_trim(dirs%input_filename) == 1))) then
1555  call save_restart(dirs%output_directory, cs%Time, g, &
1556  cs%restart_CSp, filename=ic_file)
1557  endif
1558 
1559 
1560  cs%id_area_shelf_h = register_diag_field('ocean_model', 'area_shelf_h', cs%diag%axesT1, cs%Time, &
1561  'Ice Shelf Area in cell', 'meter-2')
1562  cs%id_shelf_mass = register_diag_field('ocean_model', 'shelf_mass', cs%diag%axesT1, cs%Time, &
1563  'mass of shelf', 'kg/m^2')
1564  cs%id_h_shelf = register_diag_field('ocean_model', 'h_shelf', cs%diag%axesT1, cs%Time, &
1565  'ice shelf thickness', 'm', conversion=us%Z_to_m)
1566  cs%id_mass_flux = register_diag_field('ocean_model', 'mass_flux', cs%diag%axesT1,&
1567  cs%Time,'Total mass flux of freshwater across the ice-ocean interface.', 'kg/s')
1568  cs%id_melt = register_diag_field('ocean_model', 'melt', cs%diag%axesT1, cs%Time, &
1569  'Ice Shelf Melt Rate', 'm yr-1')
1570  cs%id_thermal_driving = register_diag_field('ocean_model', 'thermal_driving', cs%diag%axesT1, cs%Time, &
1571  'pot. temp. in the boundary layer minus freezing pot. temp. at the ice-ocean interface.', 'Celsius')
1572  cs%id_haline_driving = register_diag_field('ocean_model', 'haline_driving', cs%diag%axesT1, cs%Time, &
1573  'salinity in the boundary layer minus salinity at the ice-ocean interface.', 'psu')
1574  cs%id_Sbdry = register_diag_field('ocean_model', 'sbdry', cs%diag%axesT1, cs%Time, &
1575  'salinity at the ice-ocean interface.', 'psu')
1576  cs%id_u_ml = register_diag_field('ocean_model', 'u_ml', cs%diag%axesCu1, cs%Time, &
1577  'Eastward vel. in the boundary layer (used to compute ustar)', 'm s-1')
1578  cs%id_v_ml = register_diag_field('ocean_model', 'v_ml', cs%diag%axesCv1, cs%Time, &
1579  'Northward vel. in the boundary layer (used to compute ustar)', 'm s-1')
1580  cs%id_exch_vel_s = register_diag_field('ocean_model', 'exch_vel_s', cs%diag%axesT1, cs%Time, &
1581  'Sub-shelf salinity exchange velocity', 'm s-1')
1582  cs%id_exch_vel_t = register_diag_field('ocean_model', 'exch_vel_t', cs%diag%axesT1, cs%Time, &
1583  'Sub-shelf thermal exchange velocity', 'm s-1')
1584  cs%id_tfreeze = register_diag_field('ocean_model', 'tfreeze', cs%diag%axesT1, cs%Time, &
1585  'In Situ Freezing point at ice shelf interface', 'degC')
1586  cs%id_tfl_shelf = register_diag_field('ocean_model', 'tflux_shelf', cs%diag%axesT1, cs%Time, &
1587  'Heat conduction into ice shelf', 'W m-2')
1588  cs%id_ustar_shelf = register_diag_field('ocean_model', 'ustar_shelf', cs%diag%axesT1, cs%Time, &
1589  'Fric vel under shelf', 'm/s', conversion=us%Z_to_m*us%s_to_T)
1590  if (cs%active_shelf_dynamics) then
1591  cs%id_h_mask = register_diag_field('ocean_model', 'h_mask', cs%diag%axesT1, cs%Time, &
1592  'ice shelf thickness mask', 'none')
1593  endif
1594 
1595  id_clock_shelf = cpu_clock_id('Ice shelf', grain=clock_component)
1596  id_clock_pass = cpu_clock_id(' Ice shelf halo updates', grain=clock_routine)
1597 

References add_shelf_forces(), add_shelf_pressure(), mom_transcribe_grid::copy_dyngrid_to_mom_grid(), mom_eos::eos_init(), mom_unit_scaling::fix_restart_unit_scaling(), mom_get_input::get_mom_input(), mom_ice_shelf_state::ice_shelf_state_init(), id_clock_pass, id_clock_shelf, mom_ice_shelf_initialize::initialize_ice_thickness(), initialize_shelf_mass(), mom_error_handler::is_root_pe(), mom_error_handler::mom_error(), mom_grid::mom_grid_init(), mom_error_handler::mom_mesg(), mom_dyn_horgrid::rescale_dyn_horgrid_bathymetry(), mom_restart::restart_init(), mom_restart::restore_state(), mom_grid_initialize::set_grid_metrics(), and mom_unit_scaling::unit_scaling_init().

Referenced by mom_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ initialize_shelf_mass()

subroutine mom_ice_shelf::initialize_shelf_mass ( type(ocean_grid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
type(ice_shelf_cs), pointer  CS,
type(ice_shelf_state), intent(inout)  ISS,
logical, intent(in), optional  new_sim 
)
private

Initializes shelf mass based on three options (file, zero and user)

Parameters
[in]gThe ocean's grid structure.
[in]param_fileA structure to parse for run-time parameters
csA pointer to the ice shelf control structure
[in,out]issThe ice shelf state type that is being updated
[in]new_simIf present and false, this run is being restarted

Definition at line 1602 of file MOM_ice_shelf.F90.

1602 
1603  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1604  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1605  type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure
1606  type(ice_shelf_state), intent(inout) :: ISS !< The ice shelf state type that is being updated
1607  logical, optional, intent(in) :: new_sim !< If present and false, this run is being restarted
1608 
1609  integer :: i, j, is, ie, js, je
1610  logical :: read_shelf_area, new_sim_2
1611  character(len=240) :: config, inputdir, shelf_file, filename
1612  character(len=120) :: shelf_mass_var ! The name of shelf mass in the file.
1613  character(len=120) :: shelf_area_var ! The name of shelf area in the file.
1614  character(len=40) :: mdl = "MOM_ice_shelf"
1615  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1616 
1617  new_sim_2 = .true. ; if (present(new_sim)) new_sim_2 = new_sim
1618 
1619  call get_param(param_file, mdl, "ICE_SHELF_CONFIG", config, &
1620  "A string that specifies how the ice shelf is "//&
1621  "initialized. Valid options include:\n"//&
1622  " \tfile\t Read from a file.\n"//&
1623  " \tzero\t Set shelf mass to 0 everywhere.\n"//&
1624  " \tUSER\t Call USER_initialize_shelf_mass.\n", &
1625  fail_if_missing=.true.)
1626 
1627  select case ( trim(config) )
1628  case ("file")
1629 
1630  call time_interp_external_init()
1631 
1632  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1633  inputdir = slasher(inputdir)
1634 
1635  call get_param(param_file, mdl, "SHELF_FILE", shelf_file, &
1636  "If DYNAMIC_SHELF_MASS = True, OVERRIDE_SHELF_MOVEMENT = True "//&
1637  "and ICE_SHELF_MASS_FROM_FILE = True, this is the file from "//&
1638  "which to read the shelf mass and area.", &
1639  default="shelf_mass.nc")
1640  call get_param(param_file, mdl, "SHELF_MASS_VAR", shelf_mass_var, &
1641  "The variable in SHELF_FILE with the shelf mass.", &
1642  default="shelf_mass")
1643  call get_param(param_file, mdl, "READ_SHELF_AREA", read_shelf_area, &
1644  "If true, also read the area covered by ice-shelf from SHELF_FILE.", &
1645  default=.false.)
1646 
1647  filename = trim(slasher(inputdir))//trim(shelf_file)
1648  call log_param(param_file, mdl, "INPUTDIR/SHELF_FILE", filename)
1649 
1650  cs%id_read_mass = init_external_field(filename, shelf_mass_var, &
1651  domain=g%Domain%mpp_domain, verbose=cs%debug)
1652 
1653  if (read_shelf_area) then
1654  call get_param(param_file, mdl, "SHELF_AREA_VAR", shelf_area_var, &
1655  "The variable in SHELF_FILE with the shelf area.", &
1656  default="shelf_area")
1657 
1658  cs%id_read_area = init_external_field(filename,shelf_area_var, &
1659  domain=g%Domain%mpp_domain)
1660  endif
1661 
1662  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
1663  " initialize_shelf_mass: Unable to open "//trim(filename))
1664 
1665  case ("zero")
1666  do j=js,je ; do i=is,ie
1667  iss%mass_shelf(i,j) = 0.0
1668  iss%area_shelf_h(i,j) = 0.0
1669  enddo ; enddo
1670 
1671  case ("USER")
1672  call user_initialize_shelf_mass(iss%mass_shelf, iss%area_shelf_h, &
1673  iss%h_shelf, iss%hmask, g, cs%US, cs%user_CS, param_file, new_sim_2)
1674 
1675  case default ; call mom_error(fatal,"initialize_ice_shelf: "// &
1676  "Unrecognized ice shelf setup "//trim(config))
1677  end select
1678 

References mom_error_handler::mom_error().

Referenced by initialize_ice_shelf().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ shelf_calc_flux()

subroutine, public mom_ice_shelf::shelf_calc_flux ( type(surface), intent(inout)  state,
type(forcing), intent(inout)  fluxes,
type(time_type), intent(in)  Time,
real, intent(in)  time_step,
type(ice_shelf_cs), pointer  CS,
type(mech_forcing), intent(inout), optional  forces 
)

Calculates fluxes between the ocean and ice-shelf using the three-equations formulation (optional to use just two equations). See ICE_SHELF equations.

Parameters
[in,out]stateA structure containing fields that describe the surface state of the ocean. The intent is only inout to allow for halo updates.
[in,out]fluxesstructure containing pointers to any possible thermodynamic or mass-flux forcing fields.
[in]timeStart time of the fluxes.
[in]time_stepLength of time over which these fluxes will be applied [s].
csA pointer to the control structure returned by a previous call to initialize_ice_shelf.
[in,out]forcesA structure with the driving mechanical forces

Definition at line 193 of file MOM_ice_shelf.F90.

193  type(surface), intent(inout) :: state !< A structure containing fields that
194  !! describe the surface state of the ocean. The
195  !! intent is only inout to allow for halo updates.
196  type(forcing), intent(inout) :: fluxes !< structure containing pointers to any possible
197  !! thermodynamic or mass-flux forcing fields.
198  type(time_type), intent(in) :: Time !< Start time of the fluxes.
199  real, intent(in) :: time_step !< Length of time over which
200  !! these fluxes will be applied [s].
201  type(ice_shelf_CS), pointer :: CS !< A pointer to the control structure
202  !! returned by a previous call to
203  !! initialize_ice_shelf.
204  type(mech_forcing), optional, intent(inout) :: forces !< A structure with the driving mechanical forces
205 
206  type(ocean_grid_type), pointer :: G => null() ! The grid structure used by the ice shelf.
207  type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
208  ! various unit conversion factors
209  type(ice_shelf_state), pointer :: ISS => null() !< A structure with elements that describe
210  !! the ice-shelf state
211 
212  real, dimension(SZI_(CS%grid)) :: &
213  Rhoml, & !< Ocean mixed layer density [kg m-3].
214  dR0_dT, & !< Partial derivative of the mixed layer density
215  !< with temperature [kg m-3 degC-1].
216  dr0_ds, & !< Partial derivative of the mixed layer density
217  !< with salinity [kg m-3 ppt-1].
218  p_int !< The pressure at the ice-ocean interface [Pa].
219 
220  real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: &
221  exch_vel_t, & !< Sub-shelf thermal exchange velocity [m s-1]
222  exch_vel_s !< Sub-shelf salt exchange velocity [m s-1]
223 
224  real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
225  mass_flux !< total mass flux of freshwater across
226  real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
227  haline_driving !< (SSS - S_boundary) ice-ocean
228  !! interface, positive for melting and negative for freezing.
229  !! This is computed as part of the ISOMIP diagnostics.
230  real, parameter :: VK = 0.40 !< Von Karman's constant - dimensionless
231  real :: ZETA_N = 0.052 !> The fraction of the boundary layer over which the
232  !! viscosity is linearly increasing. (Was 1/8. Why?)
233  real, parameter :: RC = 0.20 ! critical flux Richardson number.
234  real :: I_ZETA_N !< The inverse of ZETA_N.
235  real :: LF, I_LF !< Latent Heat of fusion [J kg-1] and its inverse.
236  real :: I_VK !< The inverse of VK.
237  real :: PR, SC !< The Prandtl number and Schmidt number [nondim].
238 
239  ! 3 equations formulation variables
240  real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
241  Sbdry !< Salinities in the ocean at the interface with the ice shelf [ppt].
242  real :: Sbdry_it
243  real :: Sbdry1, Sbdry2, S_a, S_b, S_c ! use to find salt roots
244  real :: dS_it !< The interface salinity change during an iteration [ppt].
245  real :: hBL_neut !< The neutral boundary layer thickness [m].
246  real :: hBL_neut_h_molec !< The ratio of the neutral boundary layer thickness
247  !! to the molecular boundary layer thickness [nondim].
248  !### THESE ARE CURRENTLY POSITIVE UPWARD.
249  real :: wT_flux !< The vertical flux of heat just inside the ocean [degC m s-1].
250  real :: wB_flux !< The vertical flux of heat just inside the ocean [m2 s-3].
251  real :: dB_dS !< The derivative of buoyancy with salinity [m s-2 ppt-1].
252  real :: dB_dT !< The derivative of buoyancy with temperature [m s-2 degC-1].
253  real :: I_n_star, n_star_term, absf
254  real :: dIns_dwB !< The partial derivative of I_n_star with wB_flux, in ???.
255  real :: dT_ustar, dS_ustar
256  real :: ustar_h
257  real :: Gam_turb
258  real :: Gam_mol_t, Gam_mol_s
259  real :: RhoCp
260  real :: I_RhoLF
261  real :: ln_neut
262  real :: mass_exch
263  real :: Sb_min, Sb_max
264  real :: dS_min, dS_max
265  ! Variables used in iterating for wB_flux.
266  real :: wB_flux_new, DwB, dDwB_dwB_in
267  real :: I_Gam_T, I_Gam_S, dG_dwB, iDens
268  real :: u_at_h, v_at_h, Isqrt2
269  logical :: Sb_min_set, Sb_max_set
270  logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
271  logical :: coupled_GL ! If true, the grouding line position is determined based on
272  ! coupled ice-ocean dynamics.
273 
274  real, parameter :: c2_3 = 2.0/3.0
275  character(len=160) :: mesg ! The text of an error message
276  integer :: i, j, is, ie, js, je, ied, jed, it1, it3
277  real, parameter :: rho_fw = 1000.0 ! fresh water density
278 
279  if (.not. associated(cs)) call mom_error(fatal, "shelf_calc_flux: "// &
280  "initialize_ice_shelf must be called before shelf_calc_flux.")
281  call cpu_clock_begin(id_clock_shelf)
282 
283  g => cs%grid ; us => cs%US
284  iss => cs%ISS
285 
286  ! useful parameters
287  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; ied = g%ied ; jed = g%jed
288  i_zeta_n = 1.0 / zeta_n
289  lf = cs%Lat_fusion
290  i_rholf = 1.0/(cs%Rho0*lf)
291  i_lf = 1.0 / lf
292  sc = cs%kv_molec/cs%kd_molec_salt
293  pr = cs%kv_molec/cs%kd_molec_temp
294  i_vk = 1.0/vk
295  rhocp = cs%Rho0 * cs%Cp
296  isqrt2 = 1.0/sqrt(2.0)
297 
298  !first calculate molecular component
299  gam_mol_t = 12.5 * (pr**c2_3) - 6
300  gam_mol_s = 12.5 * (sc**c2_3) - 6
301 
302  idens = 1.0/cs%density_ocean_avg
303 
304  ! GMM, zero some fields of the ice shelf structure (ice_shelf_CS)
305  ! these fields are already set to zero during initialization
306  ! However, they seem to be changed somewhere and, for diagnostic
307  ! reasons, it is better to set them to zero again.
308  exch_vel_t(:,:) = 0.0 ; exch_vel_s(:,:) = 0.0
309  iss%tflux_shelf(:,:) = 0.0 ; iss%water_flux(:,:) = 0.0
310  iss%salt_flux(:,:) = 0.0; iss%tflux_ocn(:,:) = 0.0
311  iss%tfreeze(:,:) = 0.0
312  ! define Sbdry to avoid Run-Time Check Failure, when melt is not computed.
313  haline_driving(:,:) = 0.0
314  sbdry(:,:) = state%sss(:,:)
315 
316  !update time
317  cs%Time = time
318 
319  if (cs%override_shelf_movement) then
320  cs%time_step = time_step
321  ! update shelf mass
322  if (cs%mass_from_file) call update_shelf_mass(g, cs, iss, time)
323  endif
324 
325  if (cs%debug) then
326  call hchksum(fluxes%frac_shelf_h, "frac_shelf_h before apply melting", g%HI, haloshift=0)
327  call hchksum(state%sst, "sst before apply melting", g%HI, haloshift=0)
328  call hchksum(state%sss, "sss before apply melting", g%HI, haloshift=0)
329  call hchksum(state%u, "u_ml before apply melting", g%HI, haloshift=0)
330  call hchksum(state%v, "v_ml before apply melting", g%HI, haloshift=0)
331  call hchksum(state%ocean_mass, "ocean_mass before apply melting", g%HI, haloshift=0)
332  endif
333 
334  do j=js,je
335  ! Find the pressure at the ice-ocean interface, averaged only over the
336  ! part of the cell covered by ice shelf.
337  do i=is,ie ; p_int(i) = cs%g_Earth * iss%mass_shelf(i,j) ; enddo
338 
339  ! Calculate insitu densities and expansion coefficients
340  call calculate_density(state%sst(:,j), state%sss(:,j), p_int, &
341  rhoml(:), is, ie-is+1, cs%eqn_of_state)
342  call calculate_density_derivs(state%sst(:,j), state%sss(:,j), p_int, &
343  dr0_dt, dr0_ds, is, ie-is+1, cs%eqn_of_state)
344 
345  do i=is,ie
346  ! set ustar_shelf to zero. This is necessary if shelf_mass_is_dynamic
347  ! but it won't make a difference otherwise.
348  fluxes%ustar_shelf(i,j)= 0.0
349 
350  ! DNG - to allow this everywhere Hml>0.0 allows for melting under grounded cells
351  ! propose instead to allow where Hml > [some threshold]
352 
353  if ((idens*state%ocean_mass(i,j) > cs%col_thick_melt_threshold) .and. &
354  (iss%area_shelf_h(i,j) > 0.0) .and. &
355  (cs%isthermo) .and. (state%Hml(i,j) > 0.0) ) then
356 
357  if (cs%threeeq) then
358  ! Iteratively determine a self-consistent set of fluxes, with the ocean
359  ! salinity just below the ice-shelf as the variable that is being
360  ! iterated for.
361  ! ### SHOULD I SET USTAR_SHELF YET?
362 
363  u_at_h = state%u(i,j)
364  v_at_h = state%v(i,j)
365 
366  !### I think that CS%utide**1 should be CS%utide**2
367  ! Also I think that if taux_shelf and tauy_shelf have been calculated by the
368  ! ocean stress calculation, they should be used here or later to set ustar_shelf. - RWH
369  fluxes%ustar_shelf(i,j) = max(cs%ustar_bg, us%m_to_Z*us%T_to_s * &
370  sqrt(cs%cdrag*((u_at_h**2 + v_at_h**2) + cs%utide(i,j)**1)))
371 
372  ustar_h = us%Z_to_m*us%s_to_T*fluxes%ustar_shelf(i,j)
373 
374  ! I think that the following can be deleted without causing any problems.
375  ! if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then
376  ! ! These arrays are supposed to be stress components at C-grid points, which is
377  ! ! inconsistent with what is coded up here.
378  ! state%taux_shelf(i,j) = ustar_h*ustar_h*CS%Rho0*Isqrt2
379  ! state%tauy_shelf(i,j) = state%taux_shelf(i,j)
380  ! endif
381 
382  ! Estimate the neutral ocean boundary layer thickness as the minimum of the
383  ! reported ocean mixed layer thickness and the neutral Ekman depth.
384  absf = 0.25*us%s_to_T*((abs(us%s_to_T*g%CoriolisBu(i,j)) + abs(us%s_to_T*g%CoriolisBu(i-1,j-1))) + &
385  (abs(us%s_to_T*g%CoriolisBu(i,j-1)) + abs(us%s_to_T*g%CoriolisBu(i-1,j))))
386  if (absf*state%Hml(i,j) <= vk*ustar_h) then ; hbl_neut = state%Hml(i,j)
387  else ; hbl_neut = (vk*ustar_h) / absf ; endif
388  hbl_neut_h_molec = zeta_n * ((hbl_neut * ustar_h) / (5.0 * cs%Kv_molec))
389 
390  ! Determine the mixed layer buoyancy flux, wB_flux.
391  db_ds = (cs%g_Earth / rhoml(i)) * dr0_ds(i)
392  db_dt = (cs%g_Earth / rhoml(i)) * dr0_dt(i)
393  ln_neut = 0.0 ; if (hbl_neut_h_molec > 1.0) ln_neut = log(hbl_neut_h_molec)
394 
395  if (cs%find_salt_root) then
396  ! read liquidus parameters
397 
398  s_a = cs%lambda1 * cs%Gamma_T_3EQ * cs%Cp
399 ! S_b = -CS%Gamma_T_3EQ*(CS%lambda2-CS%lambda3*p_int(i)-state%sst(i,j)) &
400 ! -LF*CS%Gamma_T_3EQ/35.0
401 
402  s_b = cs%Gamma_T_3EQ*cs%Cp*(cs%lambda2+cs%lambda3*p_int(i)- &
403  state%sst(i,j))-lf*cs%Gamma_T_3EQ/35.0
404  s_c = lf*(cs%Gamma_T_3EQ/35.0)*state%sss(i,j)
405 
406  !### Depending on the sign of S_b, one of these will be inaccurate!
407  sbdry1 = (-s_b + sqrt(s_b*s_b-4*s_a*s_c))/(2*s_a)
408  sbdry2 = (-s_b - sqrt(s_b*s_b-4*s_a*s_c))/(2*s_a)
409  sbdry(i,j) = max(sbdry1, sbdry2)
410  ! Safety check
411  if (sbdry(i,j) < 0.) then
412  write(mesg,*) 'state%sss(i,j) = ',state%sss(i,j), 'S_a, S_b, S_c', s_a, s_b, s_c
413  call mom_error(warning, mesg, .true.)
414  write(mesg,*) 'I,J,Sbdry1,Sbdry2',i,j,sbdry1,sbdry2
415  call mom_error(warning, mesg, .true.)
416  call mom_error(fatal, "shelf_calc_flux: Negative salinity (Sbdry).")
417  endif
418  else
419  ! Guess sss as the iteration starting point for the boundary salinity.
420  sbdry(i,j) = state%sss(i,j) ; sb_max_set = .false.
421  sb_min_set = .false.
422  endif !find_salt_root
423 
424  do it1 = 1,20
425  ! Determine the potential temperature at the ice-ocean interface.
426  call calculate_tfreeze(sbdry(i,j), p_int(i), iss%tfreeze(i,j), cs%eqn_of_state)
427 
428  dt_ustar = (state%sst(i,j) - iss%tfreeze(i,j)) * ustar_h
429  ds_ustar = (state%sss(i,j) - sbdry(i,j)) * ustar_h
430 
431  ! First, determine the buoyancy flux assuming no effects of stability
432  ! on the turbulence. Following H & J '99, this limit also applies
433  ! when the buoyancy flux is destabilizing.
434 
435  if (cs%const_gamma) then ! if using a constant gamma_T
436  ! note the different form, here I_Gam_T is NOT 1/Gam_T!
437  i_gam_t = cs%Gamma_T_3EQ
438  i_gam_s = cs%Gamma_T_3EQ/35.
439  else
440  gam_turb = i_vk * (ln_neut + (0.5 * i_zeta_n - 1.0))
441  i_gam_t = 1.0 / (gam_mol_t + gam_turb)
442  i_gam_s = 1.0 / (gam_mol_s + gam_turb)
443  endif
444 
445  wt_flux = dt_ustar * i_gam_t
446  wb_flux = db_ds * (ds_ustar * i_gam_s) + db_dt * wt_flux
447 
448  if (wb_flux > 0.0) then
449  ! The buoyancy flux is stabilizing and will reduce the tubulent
450  ! fluxes, and iteration is required.
451  n_star_term = (zeta_n/rc) * (hbl_neut * vk) / ustar_h**3
452  do it3 = 1,30
453  ! n_star <= 1.0 is the ratio of working boundary layer thickness
454  ! to the neutral thickness.
455  ! hBL = n_star*hBL_neut ; hSub = 1/8*n_star*hBL
456 
457  i_n_star = sqrt(1.0 + n_star_term * wb_flux)
458  dins_dwb = 0.5 * n_star_term / i_n_star
459  if (hbl_neut_h_molec > i_n_star**2) then
460  gam_turb = i_vk * ((ln_neut - 2.0*log(i_n_star)) + &
461  (0.5*i_zeta_n*i_n_star - 1.0))
462  dg_dwb = i_vk * ( -2.0 / i_n_star + (0.5 * i_zeta_n)) * dins_dwb
463  else
464  ! The layer dominated by molecular viscosity is smaller than
465  ! the assumed boundary layer. This should be rare!
466  gam_turb = i_vk * (0.5 * i_zeta_n*i_n_star - 1.0)
467  dg_dwb = i_vk * (0.5 * i_zeta_n) * dins_dwb
468  endif
469 
470  if (cs%const_gamma) then ! if using a constant gamma_T
471  ! note the different form, here I_Gam_T is NOT 1/Gam_T!
472  i_gam_t = cs%Gamma_T_3EQ
473  i_gam_s = cs%Gamma_T_3EQ/35.
474  else
475  i_gam_t = 1.0 / (gam_mol_t + gam_turb)
476  i_gam_s = 1.0 / (gam_mol_s + gam_turb)
477  endif
478 
479  wt_flux = dt_ustar * i_gam_t
480  wb_flux_new = db_ds * (ds_ustar * i_gam_s) + db_dt * wt_flux
481 
482  ! Find the root where dwB = 0.0
483  dwb = wb_flux_new - wb_flux
484  if (abs(wb_flux_new - wb_flux) < &
485  1e-4*(abs(wb_flux_new) + abs(wb_flux))) exit
486 
487  ddwb_dwb_in = -dg_dwb * (db_ds * (ds_ustar * i_gam_s**2) + &
488  db_dt * (dt_ustar * i_gam_t**2)) - 1.0
489  ! This is Newton's method without any bounds.
490  ! ### SHOULD BOUNDS BE NEEDED?
491  wb_flux_new = wb_flux - dwb / ddwb_dwb_in
492  enddo !it3
493  endif
494 
495  iss%tflux_ocn(i,j) = rhocp * wt_flux
496  exch_vel_t(i,j) = ustar_h * i_gam_t
497  exch_vel_s(i,j) = ustar_h * i_gam_s
498 
499  !Calculate the heat flux inside the ice shelf.
500 
501  !vertical adv/diff as in H+J 1999, eqns (26) & approx from (31).
502  ! Q_ice = rho_ice * CS%CP_Ice * K_ice * dT/dz (at interface)
503  !vertical adv/diff as in H+J 199, eqs (31) & (26)...
504  ! dT/dz ~= min( (lprec/(rho_ice*K_ice))*(CS%Temp_Ice-T_freeze) , 0.0 )
505  !If this approximation is not made, iterations are required... See H+J Fig 3.
506 
507  if (iss%tflux_ocn(i,j) <= 0.0) then ! Freezing occurs, so zero ice heat flux.
508  iss%water_flux(i,j) = i_lf * iss%tflux_ocn(i,j)
509  iss%tflux_shelf(i,j) = 0.0
510  else
511  if (cs%insulator) then
512  !no conduction/perfect insulator
513  iss%tflux_shelf(i,j) = 0.0
514  iss%water_flux(i,j) = i_lf * (- iss%tflux_shelf(i,j) + iss%tflux_ocn(i,j))
515 
516  else
517  ! With melting, from H&J 1999, eqs (31) & (26)...
518  ! Q_ice ~= cp_ice * (CS%Temp_Ice-T_freeze) * lprec
519  ! RhoLF*lprec = Q_ice + ISS%tflux_ocn(i,j)
520  ! lprec = (ISS%tflux_ocn(i,j)) / (LF + cp_ice * (T_freeze-CS%Temp_Ice))
521  iss%water_flux(i,j) = iss%tflux_ocn(i,j) / &
522  (lf + cs%CP_Ice * (iss%tfreeze(i,j) - cs%Temp_Ice))
523 
524  iss%tflux_shelf(i,j) = iss%tflux_ocn(i,j) - lf*iss%water_flux(i,j)
525  endif
526 
527  endif
528  !other options: dTi/dz linear through shelf
529  ! dTi_dz = (CS%Temp_Ice - ISS%tfreeze(i,j))/G%draft(i,j)
530  ! ISS%tflux_shelf(i,j) = - Rho_Ice * CS%CP_Ice * KTI * dTi_dz
531 
532 
533  if (cs%find_salt_root) then
534  exit ! no need to do interaction, so exit loop
535  else
536 
537  mass_exch = exch_vel_s(i,j) * cs%Rho0
538  sbdry_it = (state%sss(i,j) * mass_exch + cs%Salin_ice * &
539  iss%water_flux(i,j)) / (mass_exch + iss%water_flux(i,j))
540  ds_it = sbdry_it - sbdry(i,j)
541  if (abs(ds_it) < 1e-4*(0.5*(state%sss(i,j) + sbdry(i,j) + 1.e-10))) exit
542 
543 
544  if (ds_it < 0.0) then ! Sbdry is now the upper bound.
545  if (sb_max_set .and. (sbdry(i,j) > sb_max)) &
546  call mom_error(fatal,"shelf_calc_flux: Irregular iteration for Sbdry (max).")
547  sb_max = sbdry(i,j) ; ds_max = ds_it ; sb_max_set = .true.
548  else ! Sbdry is now the lower bound.
549  if (sb_min_set .and. (sbdry(i,j) < sb_min)) &
550  call mom_error(fatal, &
551  "shelf_calc_flux: Irregular iteration for Sbdry (min).")
552  sb_min = sbdry(i,j) ; ds_min = ds_it ; sb_min_set = .true.
553  endif ! dS_it < 0.0
554 
555  if (sb_min_set .and. sb_max_set) then
556  ! Use the false position method for the next iteration.
557  sbdry(i,j) = sb_min + (sb_max-sb_min) * &
558  (ds_min / (ds_min - ds_max))
559  else
560  sbdry(i,j) = sbdry_it
561  endif ! Sb_min_set
562 
563  sbdry(i,j) = sbdry_it
564  endif ! CS%find_salt_root
565 
566  enddo !it1
567  ! Check for non-convergence and/or non-boundedness?
568 
569  else
570  ! In the 2-equation form, the mixed layer turbulent exchange velocity
571  ! is specified and large enough that the ocean salinity at the interface
572  ! is about the same as the boundary layer salinity.
573 
574  call calculate_tfreeze(state%sss(i,j), p_int(i), iss%tfreeze(i,j), cs%eqn_of_state)
575 
576  exch_vel_t(i,j) = cs%gamma_t
577  iss%tflux_ocn(i,j) = rhocp * exch_vel_t(i,j) * (state%sst(i,j) - iss%tfreeze(i,j))
578  iss%tflux_shelf(i,j) = 0.0
579  iss%water_flux(i,j) = i_lf * iss%tflux_ocn(i,j)
580  sbdry(i,j) = 0.0
581  endif
582  else !not shelf
583  iss%tflux_ocn(i,j) = 0.0
584  endif
585 
586 ! haline_driving(:,:) = state%sss(i,j) - Sbdry(i,j)
587 
588  enddo ! i-loop
589  enddo ! j-loop
590 
591  ! ISS%water_flux = net liquid water into the ocean ( kg/(m^2 s) )
592  ! We want melt in m/year
593  if (cs%const_gamma) then ! use ISOMIP+ eq. with rho_fw
594  fluxes%iceshelf_melt = iss%water_flux * (86400.0*365.0/rho_fw) * cs%flux_factor
595  else ! use original eq.
596  fluxes%iceshelf_melt = iss%water_flux * (86400.0*365.0/cs%density_ice) * cs%flux_factor
597  endif
598 
599  do j=js,je ; do i=is,ie
600  if ((idens*state%ocean_mass(i,j) > cs%col_thick_melt_threshold) .and. &
601  (iss%area_shelf_h(i,j) > 0.0) .and. &
602  (cs%isthermo) .and. (state%Hml(i,j) > 0.0) ) then
603 
604  ! Set melt to zero above a cutoff pressure
605  ! (CS%Rho0*CS%cutoff_depth*CS%g_Earth) this is needed for the isomip
606  ! test case.
607  if ((cs%g_Earth * iss%mass_shelf(i,j)) < cs%Rho0*cs%cutoff_depth* &
608  cs%g_Earth) then
609  iss%water_flux(i,j) = 0.0
610  fluxes%iceshelf_melt(i,j) = 0.0
611  endif
612  ! Compute haline driving, which is one of the diags. used in ISOMIP
613  haline_driving(i,j) = (iss%water_flux(i,j) * sbdry(i,j)) / &
614  (cs%Rho0 * exch_vel_s(i,j))
615 
616  !!!!!!!!!!!!!!!!!!!!!!!!!!!!Safety checks !!!!!!!!!!!!!!!!!!!!!!!!!
617  !1)Check if haline_driving computed above is consistent with
618  ! haline_driving = state%sss - Sbdry
619  !if (fluxes%iceshelf_melt(i,j) /= 0.0) then
620  ! if (haline_driving(i,j) /= (state%sss(i,j) - Sbdry(i,j))) then
621  ! write(mesg,*) 'at i,j=',i,j,' haline_driving, sss-Sbdry',haline_driving(i,j), &
622  ! (state%sss(i,j) - Sbdry(i,j))
623  ! call MOM_error(FATAL, &
624  ! "shelf_calc_flux: Inconsistency in melt and haline_driving"//trim(mesg))
625  ! endif
626  !endif
627 
628  ! 2) check if |melt| > 0 when ustar_shelf = 0.
629  ! this should never happen
630  if ((abs(fluxes%iceshelf_melt(i,j))>0.0) .and. (fluxes%ustar_shelf(i,j) == 0.0)) then
631  write(mesg,*) "|melt| = ",fluxes%iceshelf_melt(i,j)," > 0 and ustar_shelf = 0. at i,j", i, j
632  call mom_error(fatal, "shelf_calc_flux: "//trim(mesg))
633  endif
634  endif ! area_shelf_h
635  !!!!!!!!!!!!!!!!!!!!!!!!!!!!End of safety checks !!!!!!!!!!!!!!!!!!!
636  enddo ; enddo ! i- and j-loops
637 
638  ! mass flux [kg s-1], part of ISOMIP diags.
639  mass_flux(:,:) = 0.0
640  mass_flux(:,:) = iss%water_flux(:,:) * iss%area_shelf_h(:,:)
641 
642  if (cs%active_shelf_dynamics .or. cs%override_shelf_movement) then
643  call cpu_clock_begin(id_clock_pass)
644  call pass_var(iss%area_shelf_h, g%domain, complete=.false.)
645  call pass_var(iss%mass_shelf, g%domain)
646  call cpu_clock_end(id_clock_pass)
647  endif
648 
649  ! Melting has been computed, now is time to update thickness and mass
650  if ( cs%override_shelf_movement .and. (.not.cs%mass_from_file)) then
651  call change_thickness_using_melt(iss, g, time_step, fluxes, cs%rho_ice, cs%debug)
652 
653  if (cs%debug) then
654  call hchksum(iss%h_shelf, "h_shelf after change thickness using melt", g%HI, haloshift=0, scale=us%Z_to_m)
655  call hchksum(iss%mass_shelf, "mass_shelf after change thickness using melt", g%HI, haloshift=0)
656  endif
657  endif
658 
659  if (cs%debug) call mom_forcing_chksum("Before add shelf flux", fluxes, g, cs%US, haloshift=0)
660 
661  call add_shelf_flux(g, us, cs, state, fluxes)
662 
663  ! now the thermodynamic data is passed on... time to update the ice dynamic quantities
664 
665  if (cs%active_shelf_dynamics) then
666  update_ice_vel = .false.
667  coupled_gl = (cs%GL_couple .and. .not.cs%solo_ice_sheet)
668 
669  ! advect the ice shelf, and advance the front. Calving will be in here somewhere as well..
670  ! when we decide on how to do it
671  call update_ice_shelf(cs%dCS, iss, g, us, time_step, time, state%ocean_mass, coupled_gl)
672 
673  endif
674 
675  call enable_averaging(time_step,time,cs%diag)
676  if (cs%id_shelf_mass > 0) call post_data(cs%id_shelf_mass, iss%mass_shelf, cs%diag)
677  if (cs%id_area_shelf_h > 0) call post_data(cs%id_area_shelf_h, iss%area_shelf_h, cs%diag)
678  if (cs%id_ustar_shelf > 0) call post_data(cs%id_ustar_shelf, fluxes%ustar_shelf, cs%diag)
679  if (cs%id_melt > 0) call post_data(cs%id_melt, fluxes%iceshelf_melt, cs%diag)
680  if (cs%id_thermal_driving > 0) call post_data(cs%id_thermal_driving, (state%sst-iss%tfreeze), cs%diag)
681  if (cs%id_Sbdry > 0) call post_data(cs%id_Sbdry, sbdry, cs%diag)
682  if (cs%id_haline_driving > 0) call post_data(cs%id_haline_driving, haline_driving, cs%diag)
683  if (cs%id_mass_flux > 0) call post_data(cs%id_mass_flux, mass_flux, cs%diag)
684  if (cs%id_u_ml > 0) call post_data(cs%id_u_ml, state%u, cs%diag)
685  if (cs%id_v_ml > 0) call post_data(cs%id_v_ml, state%v, cs%diag)
686  if (cs%id_tfreeze > 0) call post_data(cs%id_tfreeze, iss%tfreeze, cs%diag)
687  if (cs%id_tfl_shelf > 0) call post_data(cs%id_tfl_shelf, iss%tflux_shelf, cs%diag)
688  if (cs%id_exch_vel_t > 0) call post_data(cs%id_exch_vel_t, exch_vel_t, cs%diag)
689  if (cs%id_exch_vel_s > 0) call post_data(cs%id_exch_vel_s, exch_vel_s, cs%diag)
690  if (cs%id_h_shelf > 0) call post_data(cs%id_h_shelf,iss%h_shelf,cs%diag)
691  if (cs%id_h_mask > 0) call post_data(cs%id_h_mask,iss%hmask,cs%diag)
692  call disable_averaging(cs%diag)
693 
694  if (present(forces)) then
695  call add_shelf_forces(g, us, cs, forces, do_shelf_area=(cs%active_shelf_dynamics .or. &
696  cs%override_shelf_movement))
697  endif
698 
699  call cpu_clock_end(id_clock_shelf)
700 
701  if (cs%debug) call mom_forcing_chksum("End of shelf calc flux", fluxes, g, cs%US, haloshift=0)
702 

References add_shelf_flux(), add_shelf_forces(), change_thickness_using_melt(), mom_diag_mediator::disable_averaging(), mom_diag_mediator::enable_averaging(), id_clock_pass, id_clock_shelf, mom_error_handler::mom_error(), and update_shelf_mass().

Referenced by mom_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ solo_time_step()

subroutine, public mom_ice_shelf::solo_time_step ( type(ice_shelf_cs), pointer  CS,
real, intent(in)  time_step,
integer, intent(inout)  nsteps,
type(time_type), intent(inout)  Time,
real, intent(in), optional  min_time_step_in 
)

This routine is for stepping a stand-alone ice shelf model without an ocean.

Parameters
csA pointer to the ice shelf control structure
[in]time_stepThe time interval for this update [s].
[in,out]nstepsThe running number of ice shelf steps.
[in,out]timeThe current model time
[in]min_time_step_inThe minimum permitted time step [s].

Definition at line 1758 of file MOM_ice_shelf.F90.

1758  type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure
1759  real, intent(in) :: time_step !< The time interval for this update [s].
1760  integer, intent(inout) :: nsteps !< The running number of ice shelf steps.
1761  type(time_type), intent(inout) :: Time !< The current model time
1762  real, optional, intent(in) :: min_time_step_in !< The minimum permitted time step [s].
1763 
1764  type(ocean_grid_type), pointer :: G => null()
1765  type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
1766  ! various unit conversion factors
1767  type(ice_shelf_state), pointer :: ISS => null() !< A structure with elements that describe
1768  !! the ice-shelf state
1769  integer :: is, iec, js, jec, i, j
1770  real :: time_step_remain
1771  real :: time_step_int, min_time_step
1772  character(len=240) :: mesg
1773  logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
1774  logical :: coupled_GL ! If true the grouding line position is determined based on
1775  ! coupled ice-ocean dynamics.
1776 
1777  g => cs%grid
1778  us => cs%US
1779  iss => cs%ISS
1780  is = g%isc ; iec = g%iec ; js = g%jsc ; jec = g%jec
1781 
1782  time_step_remain = time_step
1783  if (present (min_time_step_in)) then
1784  min_time_step = min_time_step_in
1785  else
1786  min_time_step = 1000.0 ! This is in seconds - at 1 km resolution it would imply ice is moving at ~1 meter per second
1787  endif
1788 
1789  write (mesg,*) "TIME in ice shelf call, yrs: ", time_type_to_real(time)/(365. * 86400.)
1790  call mom_mesg("solo_time_step: "//mesg)
1791 
1792  do while (time_step_remain > 0.0)
1793  nsteps = nsteps+1
1794 
1795  ! If time_step is not too long, this is unnecessary.
1796  time_step_int = min(ice_time_step_cfl(cs%dCS, iss, g), time_step)
1797 
1798  write (mesg,*) "Ice model timestep = ", time_step_int, " seconds"
1799  if (time_step_int < min_time_step) then
1800  call mom_error(fatal, "MOM_ice_shelf:solo_time_step: abnormally small timestep "//mesg)
1801  else
1802  call mom_mesg("solo_time_step: "//mesg)
1803  endif
1804 
1805  if (time_step_int >= time_step_remain) then
1806  time_step_int = time_step_remain
1807  time_step_remain = 0.0
1808  else
1809  time_step_remain = time_step_remain - time_step_int
1810  endif
1811 
1812  ! If the last mini-timestep is a day or less, we cannot expect velocities to change by much.
1813  ! Do not update the velocities if the last step is very short.
1814  update_ice_vel = ((time_step_int > min_time_step) .or. (time_step_int >= time_step))
1815  coupled_gl = .false.
1816 
1817  call update_ice_shelf(cs%dCS, iss, g, us, time_step_int, time, must_update_vel=update_ice_vel)
1818 
1819  call enable_averaging(time_step,time,cs%diag)
1820  if (cs%id_area_shelf_h > 0) call post_data(cs%id_area_shelf_h, iss%area_shelf_h, cs%diag)
1821  if (cs%id_h_shelf > 0) call post_data(cs%id_h_shelf, iss%h_shelf, cs%diag)
1822  if (cs%id_h_mask > 0) call post_data(cs%id_h_mask, iss%hmask, cs%diag)
1823  call disable_averaging(cs%diag)
1824 
1825  enddo
1826 

References mom_diag_mediator::disable_averaging(), mom_diag_mediator::enable_averaging(), mom_ice_shelf_dynamics::ice_time_step_cfl(), mom_error_handler::mom_error(), and mom_error_handler::mom_mesg().

Here is the call graph for this function:

◆ update_shelf_mass()

subroutine mom_ice_shelf::update_shelf_mass ( type(ocean_grid_type), intent(inout)  G,
type(ice_shelf_cs), intent(in)  CS,
type(ice_shelf_state), intent(inout)  ISS,
type(time_type), intent(in)  Time 
)
private

Updates the ice shelf mass using data from a file.

Parameters
[in,out]gThe ocean's grid structure.
[in]csA pointer to the ice shelf control structure
[in,out]issThe ice shelf state type that is being updated
[in]timeThe current model time

Definition at line 1683 of file MOM_ice_shelf.F90.

1683  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1684  type(ice_shelf_CS), intent(in) :: CS !< A pointer to the ice shelf control structure
1685  type(ice_shelf_state), intent(inout) :: ISS !< The ice shelf state type that is being updated
1686  type(time_type), intent(in) :: Time !< The current model time
1687 
1688  ! local variables
1689  integer :: i, j, is, ie, js, je
1690  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1691 
1692  call time_interp_external(cs%id_read_mass, time, iss%mass_shelf)
1693 
1694  do j=js,je ; do i=is,ie
1695  iss%area_shelf_h(i,j) = 0.0
1696  iss%hmask(i,j) = 0.
1697  if (iss%mass_shelf(i,j) > 0.0) then
1698  iss%area_shelf_h(i,j) = g%US%L_to_m**2*g%areaT(i,j)
1699  iss%h_shelf(i,j) = iss%mass_shelf(i,j) / cs%rho_ice
1700  iss%hmask(i,j) = 1.
1701  endif
1702  enddo ; enddo
1703 
1704  !call USER_update_shelf_mass(ISS%mass_shelf, ISS%area_shelf_h, ISS%h_shelf, &
1705  ! ISS%hmask, CS%grid, CS%user_CS, Time, .true.)
1706 
1707  if (cs%min_thickness_simple_calve > 0.0) then
1708  call ice_shelf_min_thickness_calve(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, &
1709  cs%min_thickness_simple_calve)
1710  endif
1711 
1712  call pass_var(iss%area_shelf_h, g%domain)
1713  call pass_var(iss%h_shelf, g%domain)
1714  call pass_var(iss%hmask, g%domain)
1715  call pass_var(iss%mass_shelf, g%domain)
1716 

Referenced by shelf_calc_flux().

Here is the caller graph for this function:

Variable Documentation

◆ id_clock_pass

integer mom_ice_shelf::id_clock_pass
private

CPU Clock for group pass calls.

Definition at line 185 of file MOM_ice_shelf.F90.

185 integer :: id_clock_pass !< CPU Clock for group pass calls

Referenced by initialize_ice_shelf(), and shelf_calc_flux().

◆ id_clock_shelf

integer mom_ice_shelf::id_clock_shelf

CPU Clock for the ice shelf code.

Definition at line 184 of file MOM_ice_shelf.F90.

184 integer :: id_clock_shelf !< CPU Clock for the ice shelf code

Referenced by initialize_ice_shelf(), and shelf_calc_flux().