MOM6
Functions/Subroutines | Variables
mom_state_initialization Module Reference

Detailed Description

Initialization functions for state variables, u, v, h, T and S.

Functions/Subroutines

subroutine, public mom_initialize_state (u, v, h, tv, Time, G, GV, US, PF, dirs, restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, ALE_sponge_CSp, OBC, Time_in)
 Initialize temporally evolving fields, either as initial conditions or by reading them from a restart (or saves) file. More...
 
subroutine initialize_thickness_from_file (h, G, GV, US, param_file, file_has_thickness, just_read_params)
 Reads the layer thicknesses or interface heights from a file. More...
 
subroutine adjustetatofitbathymetry (G, GV, US, eta, h)
 Adjust interface heights to fit the bathymetry and diagnose layer thickness. More...
 
subroutine initialize_thickness_uniform (h, G, GV, param_file, just_read_params)
 Initializes thickness to be uniform. More...
 
subroutine initialize_thickness_list (h, G, GV, US, param_file, just_read_params)
 Initialize thickness from a 1D list. More...
 
subroutine initialize_thickness_search
 Search density space for location of layers (not implemented!)
 
subroutine convert_thickness (h, G, GV, US, tv)
 Converts thickness from geometric to pressure units. More...
 
subroutine depress_surface (h, G, GV, US, param_file, tv, just_read_params)
 Depress the sea-surface based on an initial condition file. More...
 
subroutine trim_for_ice (PF, G, GV, US, ALE_CSp, tv, h, just_read_params)
 Adjust the layer thicknesses by cutting away the top of each model column at the depth where the hydrostatic pressure matches an imposed surface pressure read from file. More...
 
subroutine cut_off_column_top (nk, tv, GV, US, G_earth, depth, min_thickness, T, T_t, T_b, S, S_t, S_b, p_surf, h, remap_CS, z_tol, remap_answers_2018)
 Adjust the layer thicknesses by removing the top of the water column above the depth where the hydrostatic pressure matches p_surf. More...
 
subroutine initialize_velocity_from_file (u, v, G, US, param_file, just_read_params)
 Initialize horizontal velocity components from file. More...
 
subroutine initialize_velocity_zero (u, v, G, param_file, just_read_params)
 Initialize horizontal velocity components to zero. More...
 
subroutine initialize_velocity_uniform (u, v, G, US, param_file, just_read_params)
 Sets the initial velocity components to uniform. More...
 
subroutine initialize_velocity_circular (u, v, G, US, param_file, just_read_params)
 Sets the initial velocity components to be circular with no flow at edges of domain and center. More...
 
subroutine initialize_temp_salt_from_file (T, S, G, param_file, just_read_params)
 Initializes temperature and salinity from file. More...
 
subroutine initialize_temp_salt_from_profile (T, S, G, param_file, just_read_params)
 Initializes temperature and salinity from a 1D profile. More...
 
subroutine initialize_temp_salt_fit (T, S, G, GV, US, param_file, eqn_of_state, P_Ref, just_read_params)
 Initializes temperature and salinity by fitting to density. More...
 
subroutine initialize_temp_salt_linear (T, S, G, param_file, just_read_params)
 Initializes T and S with linear profiles according to reference surface layer salinity and temperature and a specified range. More...
 
subroutine initialize_sponges_file (G, GV, US, use_temperature, tv, param_file, Layer_CSp, ALE_CSp, Time)
 This subroutine sets the inverse restoration time (Idamp), and the values towards which the interface heights and an arbitrary number of tracers should be restored within each sponge. The interface height is always subject to damping, and must always be the first registered field. More...
 
subroutine set_velocity_depth_max (G)
 This subroutine sets the 4 bottom depths at velocity points to be the maximum of the adjacent depths. More...
 
subroutine compute_global_grid_integrals (G, US)
 Subroutine to pre-compute global integrals of grid quantities for later use in reporting diagnostics. More...
 
subroutine set_velocity_depth_min (G)
 This subroutine sets the 4 bottom depths at velocity points to be the minimum of the adjacent depths. More...
 
subroutine mom_temp_salt_initialize_from_z (h, tv, G, GV, US, PF, just_read_params)
 This subroutine determines the isopycnal or other coordinate interfaces and layer potential temperatures and salinities directly from a z-space file on a latitude-longitude grid. More...
 
subroutine find_interfaces (rho, zin, nk_data, Rb, depth, zi, G, US, nlevs, nkml, hml, eps_z, eps_rho)
 Find interface positions corresponding to interpolated depths in a density profile. More...
 
subroutine mom_state_init_tests (G, GV, US, tv)
 Run simple unit tests. More...
 

Variables

character(len=40) mdl = "MOM_state_initialization"
 This module's name.
 

Function/Subroutine Documentation

◆ adjustetatofitbathymetry()

subroutine mom_state_initialization::adjustetatofitbathymetry ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(inout)  eta,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  h 
)
private

Adjust interface heights to fit the bathymetry and diagnose layer thickness.

If the bottom most interface is below the topography then the bottom-most layers are contracted to GVAngstrom_m. If the bottom most interface is above the topography then the entire column is dilated (expanded) to fill the void.

Remarks
{There is a (hard-wired) "tolerance" parameter such that the criteria for adjustment must equal or exceed 10cm.}
Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in,out]etaInterface heights [Z ~> m].
[in,out]hLayer thicknesses [H ~> m or kg m-2]

Definition at line 712 of file MOM_state_initialization.F90.

713  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
714  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
715  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
716  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: eta !< Interface heights [Z ~> m].
717  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
718  ! Local variables
719  integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
720  real :: hTolerance = 0.1 !< Tolerance to exceed adjustment criteria [Z ~> m]
721  real :: hTmp, eTmp, dilate
722  character(len=100) :: mesg
723 
724  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
725  htolerance = 0.1*us%m_to_Z
726 
727  contractions = 0
728  do j=js,je ; do i=is,ie
729  if (-eta(i,j,nz+1) > g%bathyT(i,j) + htolerance) then
730  eta(i,j,nz+1) = -g%bathyT(i,j)
731  contractions = contractions + 1
732  endif
733  enddo ; enddo
734  call sum_across_pes(contractions)
735  if ((contractions > 0) .and. (is_root_pe())) then
736  write(mesg,'("Thickness initial conditions were contracted ",'// &
737  '"to fit topography in ",I8," places.")') contractions
738  call mom_error(warning, 'adjustEtaToFitBathymetry: '//mesg)
739  endif
740 
741  ! To preserve previous answers in non-Boussinesq cases, delay converting
742  ! thicknesses to units of H until the end of this routine.
743  do k=nz,1,-1 ; do j=js,je ; do i=is,ie
744  ! Collapse layers to thinnest possible if the thickness less than
745  ! the thinnest possible (or negative).
746  if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) then
747  eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
748  h(i,j,k) = gv%Angstrom_Z
749  else
750  h(i,j,k) = (eta(i,j,k) - eta(i,j,k+1))
751  endif
752  enddo ; enddo ; enddo
753 
754  dilations = 0
755  do j=js,je ; do i=is,ie
756  ! The whole column is dilated to accommodate deeper topography than
757  ! the bathymetry would indicate.
758  ! This should be... if ((G%mask2dt(i,j)*(eta(i,j,1)-eta(i,j,nz+1)) > 0.0) .and. &
759  if (-eta(i,j,nz+1) < g%bathyT(i,j) - htolerance) then
760  dilations = dilations + 1
761  if (eta(i,j,1) <= eta(i,j,nz+1)) then
762  do k=1,nz ; h(i,j,k) = (eta(i,j,1) + g%bathyT(i,j)) / real(nz) ; enddo
763  else
764  dilate = (eta(i,j,1) + g%bathyT(i,j)) / (eta(i,j,1) - eta(i,j,nz+1))
765  do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
766  endif
767  do k=nz,2,-1 ; eta(i,j,k) = eta(i,j,k+1) + h(i,j,k) ; enddo
768  endif
769  enddo ; enddo
770 
771  ! Now convert thicknesses to units of H.
772  do k=1,nz ; do j=js,je ; do i=is,ie
773  h(i,j,k) = h(i,j,k)*gv%Z_to_H
774  enddo ; enddo ; enddo
775 
776  call sum_across_pes(dilations)
777  if ((dilations > 0) .and. (is_root_pe())) then
778  write(mesg,'("Thickness initial conditions were dilated ",'// &
779  '"to fit topography in ",I8," places.")') dilations
780  call mom_error(warning, 'adjustEtaToFitBathymetry: '//mesg)
781  endif
782 

◆ compute_global_grid_integrals()

subroutine mom_state_initialization::compute_global_grid_integrals ( type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US 
)
private

Subroutine to pre-compute global integrals of grid quantities for later use in reporting diagnostics.

Parameters
[in,out]gThe ocean's grid structure
[in]usA dimensional unit scaling type

Definition at line 1932 of file MOM_state_initialization.F90.

1933  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1934  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1935  ! Local variables
1936  real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
1937  real :: area_scale
1938  integer :: i,j
1939 
1940  area_scale = us%L_to_m**2
1941  tmpforsumming(:,:) = 0.
1942  g%areaT_global = 0.0 ; g%IareaT_global = 0.0
1943  do j=g%jsc,g%jec ; do i=g%isc,g%iec
1944  tmpforsumming(i,j) = area_scale*g%areaT(i,j) * g%mask2dT(i,j)
1945  enddo ; enddo
1946  g%areaT_global = reproducing_sum(tmpforsumming)
1947  g%IareaT_global = 1. / (g%areaT_global)

◆ convert_thickness()

subroutine mom_state_initialization::convert_thickness ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  h,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(thermo_var_ptrs), intent(in)  tv 
)
private

Converts thickness from geometric to pressure units.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in,out]hInput geometric layer thicknesses being converted
[in]tvA structure pointing to various thermodynamic variables

Definition at line 922 of file MOM_state_initialization.F90.

923  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
924  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
925  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
926  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
927  intent(inout) :: h !< Input geometric layer thicknesses being converted
928  !! to layer pressure [H ~> m or kg m-2].
929  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
930  !! thermodynamic variables
931  ! Local variables
932  real, dimension(SZI_(G),SZJ_(G)) :: &
933  p_top, p_bot ! Pressure at the interfaces above and below a layer [R L2 T-2 ~> Pa]
934  real :: dz_geo(SZI_(G),SZJ_(G)) ! The change in geopotential height across a layer [L2 T-2 ~> m2 s-2]
935  real :: rho(SZI_(G)) ! The in situ density [R ~> kg m-3]
936  real :: I_gEarth ! Unit conversion factors divided by the gravitational acceleration
937  ! [H T2 R-1 L-2 ~> s2 m2 kg-1 or s2 m-1]
938  real :: HR_to_pres ! A conversion factor from the input geometric thicknesses times the layer
939  ! densities into pressure units [L2 T-2 H-1 ~> m s-2 or m4 kg-1 s-2].
940  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
941  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
942  integer :: itt, max_itt
943 
944  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
945  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
946  max_itt = 10
947 
948  if (gv%Boussinesq) then
949  call mom_error(fatal,"Not yet converting thickness with Boussinesq approx.")
950  else
951  i_gearth = gv%RZ_to_H / gv%g_Earth
952  hr_to_pres = gv%g_Earth * gv%H_to_Z
953 
954  if (associated(tv%eqn_of_state)) then
955  do j=jsq,jeq+1 ; do i=isq,ieq+1
956  p_bot(i,j) = 0.0 ; p_top(i,j) = 0.0
957  enddo ; enddo
958  eosdom(:) = eos_domain(g%HI)
959  do k=1,nz
960  do j=js,je
961  do i=is,ie ; p_top(i,j) = p_bot(i,j) ; enddo
962  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_top(:,j), rho, &
963  tv%eqn_of_state, eosdom)
964  do i=is,ie
965  p_bot(i,j) = p_top(i,j) + hr_to_pres * (h(i,j,k) * rho(i))
966  enddo
967  enddo
968 
969  do itt=1,max_itt
970  call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p_top, p_bot, 0.0, g%HI, &
971  tv%eqn_of_state, us, dz_geo)
972  if (itt < max_itt) then ; do j=js,je
973  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_bot(:,j), rho, &
974  tv%eqn_of_state, eosdom)
975  ! Use Newton's method to correct the bottom value.
976  ! The hydrostatic equation is sufficiently linear that no bounds-checking is needed.
977  do i=is,ie
978  p_bot(i,j) = p_bot(i,j) + rho(i) * (hr_to_pres*h(i,j,k) - dz_geo(i,j))
979  enddo
980  enddo ; endif
981  enddo
982 
983  do j=js,je ; do i=is,ie
984  h(i,j,k) = (p_bot(i,j) - p_top(i,j)) * i_gearth
985  enddo ; enddo
986  enddo
987  else
988  do k=1,nz ; do j=js,je ; do i=is,ie
989  h(i,j,k) = h(i,j,k) * (gv%Rlay(k) / gv%Rho0)
990  enddo ; enddo ; enddo
991  endif
992  endif
993 

◆ cut_off_column_top()

subroutine mom_state_initialization::cut_off_column_top ( integer, intent(in)  nk,
type(thermo_var_ptrs), intent(in)  tv,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, intent(in)  G_earth,
real, intent(in)  depth,
real, intent(in)  min_thickness,
real, dimension(nk), intent(inout)  T,
real, dimension(nk), intent(in)  T_t,
real, dimension(nk), intent(in)  T_b,
real, dimension(nk), intent(inout)  S,
real, dimension(nk), intent(in)  S_t,
real, dimension(nk), intent(in)  S_b,
real, intent(in)  p_surf,
real, dimension(nk), intent(inout)  h,
type(remapping_cs), pointer  remap_CS,
real, intent(in), optional  z_tol,
logical, intent(in), optional  remap_answers_2018 
)
private

Adjust the layer thicknesses by removing the top of the water column above the depth where the hydrostatic pressure matches p_surf.

Parameters
[in]nkNumber of layers
[in]tvThermodynamics structure
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]g_earthGravitational acceleration [L2 Z-1 T-2 ~> m s-2]
[in]depthDepth of ocean column [Z ~> m].
[in]min_thicknessSmallest thickness allowed [Z ~> m].
[in,out]tLayer mean temperature [degC]
[in]t_tTemperature at top of layer [degC]
[in]t_bTemperature at bottom of layer [degC]
[in,out]sLayer mean salinity [ppt]
[in]s_tSalinity at top of layer [ppt]
[in]s_bSalinity at bottom of layer [ppt]
[in]p_surfImposed pressure on ocean at surface [R L2 T-2 ~> Pa]
[in,out]hLayer thickness [H ~> m or kg m-2]
remap_csRemapping structure for remapping T and S, if associated
[in]z_tolThe tolerance with which to find the depth matching the specified pressure [Z ~> m].
[in]remap_answers_2018If true, use the order of arithmetic and expressions that recover the answers for remapping from the end of 2018. Otherwise, use more robust forms of the same expressions.

Definition at line 1171 of file MOM_state_initialization.F90.

1173  integer, intent(in) :: nk !< Number of layers
1174  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1175  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1176  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1177  real, intent(in) :: G_earth !< Gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
1178  real, intent(in) :: depth !< Depth of ocean column [Z ~> m].
1179  real, intent(in) :: min_thickness !< Smallest thickness allowed [Z ~> m].
1180  real, dimension(nk), intent(inout) :: T !< Layer mean temperature [degC]
1181  real, dimension(nk), intent(in) :: T_t !< Temperature at top of layer [degC]
1182  real, dimension(nk), intent(in) :: T_b !< Temperature at bottom of layer [degC]
1183  real, dimension(nk), intent(inout) :: S !< Layer mean salinity [ppt]
1184  real, dimension(nk), intent(in) :: S_t !< Salinity at top of layer [ppt]
1185  real, dimension(nk), intent(in) :: S_b !< Salinity at bottom of layer [ppt]
1186  real, intent(in) :: p_surf !< Imposed pressure on ocean at surface [R L2 T-2 ~> Pa]
1187  real, dimension(nk), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1188  type(remapping_CS), pointer :: remap_CS !< Remapping structure for remapping T and S,
1189  !! if associated
1190  real, optional, intent(in) :: z_tol !< The tolerance with which to find the depth
1191  !! matching the specified pressure [Z ~> m].
1192  logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic
1193  !! and expressions that recover the answers for remapping
1194  !! from the end of 2018. Otherwise, use more robust
1195  !! forms of the same expressions.
1196 
1197  ! Local variables
1198  real, dimension(nk+1) :: e ! Top and bottom edge values for reconstructions [Z ~> m]
1199  real, dimension(nk) :: h0, S0, T0, h1, S1, T1
1200  real :: P_t, P_b ! Top and bottom pressures [R L2 T-2 ~> Pa]
1201  real :: z_out, e_top
1202  logical :: answers_2018
1203  integer :: k
1204 
1205  answers_2018 = .true. ; if (present(remap_answers_2018)) answers_2018 = remap_answers_2018
1206 
1207  ! Calculate original interface positions
1208  e(nk+1) = -depth
1209  do k=nk,1,-1
1210  e(k) = e(k+1) + gv%H_to_Z*h(k)
1211  h0(k) = h(nk+1-k) ! Keep a copy to use in remapping
1212  enddo
1213 
1214  p_t = 0.
1215  e_top = e(1)
1216  do k=1,nk
1217  call find_depth_of_pressure_in_cell(t_t(k), t_b(k), s_t(k), s_b(k), e(k), e(k+1), &
1218  p_t, p_surf, gv%Rho0, g_earth, tv%eqn_of_state, &
1219  us, p_b, z_out, z_tol=z_tol)
1220  if (z_out>=e(k)) then
1221  ! Imposed pressure was less that pressure at top of cell
1222  exit
1223  elseif (z_out<=e(k+1)) then
1224  ! Imposed pressure was greater than pressure at bottom of cell
1225  e_top = e(k+1)
1226  else
1227  ! Imposed pressure was fell between pressures at top and bottom of cell
1228  e_top = z_out
1229  exit
1230  endif
1231  p_t = p_b
1232  enddo
1233  if (e_top<e(1)) then
1234  ! Clip layers from the top down, if at all
1235  do k=1,nk
1236  if (e(k) > e_top) then
1237  ! Original e(K) is too high
1238  e(k) = e_top
1239  e_top = e_top - min_thickness ! Next interface must be at least this deep
1240  endif
1241  ! This layer needs trimming
1242  h(k) = gv%Z_to_H * max( min_thickness, e(k) - e(k+1) )
1243  if (e(k) < e_top) exit ! No need to go further
1244  enddo
1245  endif
1246 
1247  ! Now we need to remap but remapping assumes the surface is at the
1248  ! same place in the two columns so we turn the column upside down.
1249  if (associated(remap_cs)) then
1250  do k=1,nk
1251  s0(k) = s(nk+1-k)
1252  t0(k) = t(nk+1-k)
1253  h1(k) = h(nk+1-k)
1254  enddo
1255  if (answers_2018) then
1256  call remapping_core_h(remap_cs, nk, h0, t0, nk, h1, t1, 1.0e-30*gv%m_to_H, 1.0e-10*gv%m_to_H)
1257  call remapping_core_h(remap_cs, nk, h0, s0, nk, h1, s1, 1.0e-30*gv%m_to_H, 1.0e-10*gv%m_to_H)
1258  else
1259  call remapping_core_h(remap_cs, nk, h0, t0, nk, h1, t1, gv%H_subroundoff, gv%H_subroundoff)
1260  call remapping_core_h(remap_cs, nk, h0, s0, nk, h1, s1, gv%H_subroundoff, gv%H_subroundoff)
1261  endif
1262  do k=1,nk
1263  s(k) = s1(nk+1-k)
1264  t(k) = t1(nk+1-k)
1265  enddo
1266  endif
1267 

◆ depress_surface()

subroutine mom_state_initialization::depress_surface ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  h,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(thermo_var_ptrs), intent(in)  tv,
logical, intent(in), optional  just_read_params 
)
private

Depress the sea-surface based on an initial condition file.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in,out]hLayer thicknesses [H ~> m or kg m-2]
[in]param_fileA structure to parse for run-time parameters
[in]tvA structure pointing to various thermodynamic variables
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 997 of file MOM_state_initialization.F90.

998  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
999  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1000  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1001  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1002  intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
1003  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1004  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
1005  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
1006  !! only read parameters without changing h.
1007  ! Local variables
1008  real, dimension(SZI_(G),SZJ_(G)) :: &
1009  eta_sfc ! The free surface height that the model should use [Z ~> m].
1010  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
1011  eta ! The free surface height that the model should use [Z ~> m].
1012  real :: dilate ! A ratio by which layers are dilated [nondim].
1013  real :: scale_factor ! A scaling factor for the eta_sfc values that are read
1014  ! in, which can be used to change units, for example.
1015  character(len=40) :: mdl = "depress_surface" ! This subroutine's name.
1016  character(len=200) :: inputdir, eta_srf_file ! Strings for file/path
1017  character(len=200) :: filename, eta_srf_var ! Strings for file/path
1018  logical :: just_read ! If true, just read parameters but set nothing.
1019  integer :: i, j, k, is, ie, js, je, nz
1020  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1021 
1022  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
1023 
1024  ! Read the surface height (or pressure) from a file.
1025 
1026  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1027  inputdir = slasher(inputdir)
1028  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_FILE", eta_srf_file,&
1029  "The initial condition file for the surface height.", &
1030  fail_if_missing=.not.just_read, do_not_log=just_read)
1031  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, &
1032  "The initial condition variable for the surface height.",&
1033  default="SSH", do_not_log=just_read)
1034  filename = trim(inputdir)//trim(eta_srf_file)
1035  if (.not.just_read) &
1036  call log_param(param_file, mdl, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
1037 
1038  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_SCALE", scale_factor, &
1039  "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into units of m", &
1040  units="variable", default=1.0, scale=us%m_to_Z, do_not_log=just_read)
1041 
1042  if (just_read) return ! All run-time parameters have been read, so return.
1043 
1044  call mom_read_data(filename, eta_srf_var, eta_sfc, g%Domain, scale=scale_factor)
1045 
1046  ! Convert thicknesses to interface heights.
1047  call find_eta(h, tv, g, gv, us, eta)
1048 
1049  do j=js,je ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
1050 ! if (eta_sfc(i,j) < eta(i,j,nz+1)) then
1051  ! Issue a warning?
1052 ! endif
1053  if (eta_sfc(i,j) > eta(i,j,1)) then
1054  ! Dilate the water column to agree, but only up to 10-fold.
1055  if (eta_sfc(i,j) - eta(i,j,nz+1) > 10.0*(eta(i,j,1) - eta(i,j,nz+1))) then
1056  dilate = 10.0
1057  call mom_error(warning, "Free surface height dilation attempted "//&
1058  "to exceed 10-fold.", all_print=.true.)
1059  else
1060  dilate = (eta_sfc(i,j) - eta(i,j,nz+1)) / (eta(i,j,1) - eta(i,j,nz+1))
1061  endif
1062  do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
1063  elseif (eta(i,j,1) > eta_sfc(i,j)) then
1064  ! Remove any mass that is above the target free surface.
1065  do k=1,nz
1066  if (eta(i,j,k) <= eta_sfc(i,j)) exit
1067  if (eta(i,j,k+1) >= eta_sfc(i,j)) then
1068  h(i,j,k) = gv%Angstrom_H
1069  else
1070  h(i,j,k) = max(gv%Angstrom_H, h(i,j,k) * &
1071  (eta_sfc(i,j) - eta(i,j,k+1)) / (eta(i,j,k) - eta(i,j,k+1)) )
1072  endif
1073  enddo
1074  endif
1075  endif ; enddo ; enddo
1076 

◆ find_interfaces()

subroutine mom_state_initialization::find_interfaces ( real, dimension( g %isd: g %ied, g %jsd: g %jed,nk_data), intent(in)  rho,
real, dimension(nk_data), intent(in)  zin,
integer, intent(in)  nk_data,
real, dimension( g %ke+1), intent(in)  Rb,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  depth,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(out)  zi,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
integer, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  nlevs,
integer, intent(in)  nkml,
real, intent(in)  hml,
real, intent(in)  eps_z,
real, intent(in)  eps_rho 
)
private

Find interface positions corresponding to interpolated depths in a density profile.

Parameters
[in]gThe ocean's grid structure
[in]nk_dataThe number of levels in the input data
[in]rhoPotential density in z-space [R ~> kg m-3]
[in]zinInput data levels [Z ~> m].
[in]rbtarget interface densities [R ~> kg m-3]
[in]depthocean depth [Z ~> m].
[out]ziThe returned interface heights [Z ~> m]
[in]usA dimensional unit scaling type
[in]nlevsnumber of valid points in each column
[in]nkmlnumber of mixed layer pieces to distribute over a depth of hml.
[in]hmlmixed layer depth [Z ~> m].
[in]eps_zA negligibly small layer thickness [Z ~> m].
[in]eps_rhoA negligibly small density difference [R ~> kg m-3].

Definition at line 2452 of file MOM_state_initialization.F90.

2454  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
2455  integer, intent(in) :: nk_data !< The number of levels in the input data
2456  real, dimension(SZI_(G),SZJ_(G),nk_data), &
2457  intent(in) :: rho !< Potential density in z-space [R ~> kg m-3]
2458  real, dimension(nk_data), intent(in) :: zin !< Input data levels [Z ~> m].
2459  real, dimension(SZK_(G)+1), intent(in) :: Rb !< target interface densities [R ~> kg m-3]
2460  real, dimension(SZI_(G),SZJ_(G)), &
2461  intent(in) :: depth !< ocean depth [Z ~> m].
2462  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
2463  intent(out) :: zi !< The returned interface heights [Z ~> m]
2464  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2465  integer, dimension(SZI_(G),SZJ_(G)), &
2466  intent(in) :: nlevs !< number of valid points in each column
2467  integer, intent(in) :: nkml !< number of mixed layer pieces to distribute over
2468  !! a depth of hml.
2469  real, intent(in) :: hml !< mixed layer depth [Z ~> m].
2470  real, intent(in) :: eps_z !< A negligibly small layer thickness [Z ~> m].
2471  real, intent(in) :: eps_rho !< A negligibly small density difference [R ~> kg m-3].
2472 
2473  ! Local variables
2474  real, dimension(nk_data) :: rho_ ! A column of densities [R ~> kg m-3]
2475  real, dimension(SZK_(G)+1) :: zi_ ! A column interface heights (negative downward) [Z ~> m].
2476  real :: slope ! The rate of change of height with density [Z R-1 ~> m4 kg-1]
2477  real :: drhodz ! A local vertical density gradient [R Z-1 ~> kg m-4]
2478  real, parameter :: zoff=0.999
2479  logical :: unstable ! True if the column is statically unstable anywhere.
2480  integer :: nlevs_data ! The number of data values in a column.
2481  logical :: work_down ! This indicates whether this pass goes up or down the water column.
2482  integer :: k_int, lo_int, hi_int, mid
2483  integer :: i, j, k, is, ie, js, je, nz
2484 
2485  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2486 
2487  zi(:,:,:) = 0.0
2488 
2489  do j=js,je ; do i=is,ie
2490  nlevs_data = nlevs(i,j)
2491  do k=1,nlevs_data ; rho_(k) = rho(i,j,k) ; enddo
2492 
2493  unstable=.true.
2494  work_down = .true.
2495  do while (unstable)
2496  ! Modifiy the input profile until it no longer has densities that decrease with depth.
2497  unstable=.false.
2498  if (work_down) then
2499  do k=2,nlevs_data-1 ; if (rho_(k) - rho_(k-1) < 0.0 ) then
2500  if (k == 2) then
2501  rho_(k-1) = rho_(k) - eps_rho
2502  else
2503  drhodz = (rho_(k+1)-rho_(k-1)) / (zin(k+1)-zin(k-1))
2504  if (drhodz < 0.0) unstable=.true.
2505  rho_(k) = rho_(k-1) + drhodz*zoff*(zin(k)-zin(k-1))
2506  endif
2507  endif ; enddo
2508  work_down = .false.
2509  else
2510  do k=nlevs_data-1,2,-1 ; if (rho_(k+1) - rho_(k) < 0.0) then
2511  if (k == nlevs_data-1) then
2512  rho_(k+1) = rho_(k-1) + eps_rho !### This should be rho_(k) + eps_rho
2513  else
2514  drhodz = (rho_(k+1)-rho_(k-1)) / (zin(k+1)-zin(k-1))
2515  if (drhodz < 0.0) unstable=.true.
2516  rho_(k) = rho_(k+1) - drhodz*(zin(k+1)-zin(k))
2517  endif
2518  endif ; enddo
2519  work_down = .true.
2520  endif
2521  enddo
2522 
2523  ! Find and store the interface depths.
2524  zi_(1) = 0.0
2525  do k=2,nz
2526  ! Find the value of k_int in the list of rho_ where rho_(k_int) <= Rb(K) < rho_(k_int+1).
2527  ! This might be made a little faster by exploiting the fact that Rb is
2528  ! monotonically increasing and not resetting lo_int back to 1 inside the K loop.
2529  lo_int = 1 ; hi_int = nlevs_data
2530  do while (lo_int < hi_int)
2531  mid = (lo_int+hi_int) / 2
2532  if (rb(k) < rho_(mid)) then ; hi_int = mid
2533  else ; lo_int = mid+1 ; endif
2534  enddo
2535  k_int = max(1, lo_int-1)
2536 
2537  ! Linearly interpolate to find the depth, zi_, where Rb would be found.
2538  slope = (zin(k_int+1) - zin(k_int)) / max(rho_(k_int+1) - rho_(k_int), eps_rho)
2539  zi_(k) = -1.0*(zin(k_int) + slope*(rb(k)-rho_(k_int)))
2540  zi_(k) = min(max(zi_(k), -depth(i,j)), -1.0*hml)
2541  enddo
2542  zi_(nz+1) = -depth(i,j)
2543  if (nkml > 0) then ; do k=2,nkml+1
2544  zi_(k) = max(hml*((1.0-real(k))/real(nkml)), -depth(i,j))
2545  enddo ; endif
2546  do k=nz,max(nkml+2,2),-1
2547  if (zi_(k) < zi_(k+1) + eps_z) zi_(k) = zi_(k+1) + eps_z
2548  if (zi_(k) > -1.0*hml) zi_(k) = max(-1.0*hml, -depth(i,j))
2549  enddo
2550 
2551  do k=1,nz+1
2552  zi(i,j,k) = zi_(k)
2553  enddo
2554  enddo ; enddo ! i- and j- loops
2555 

◆ initialize_sponges_file()

subroutine mom_state_initialization::initialize_sponges_file ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
logical, intent(in)  use_temperature,
type(thermo_var_ptrs), intent(in)  tv,
type(param_file_type), intent(in)  param_file,
type(sponge_cs), pointer  Layer_CSp,
type(ale_sponge_cs), pointer  ALE_CSp,
type(time_type), intent(in)  Time 
)
private

This subroutine sets the inverse restoration time (Idamp), and the values towards which the interface heights and an arbitrary number of tracers should be restored within each sponge. The interface height is always subject to damping, and must always be the first registered field.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]use_temperatureIf true, T & S are state variables.
[in]tvA structure pointing to various thermodynamic variables.
[in]param_fileA structure to parse for run-time parameters.
layer_cspA pointer that is set to point to the control structure for this module (in layered mode).
ale_cspA pointer that is set to point to the control structure for this module (in ALE mode).
[in]timeTime at the start of the run segment. Time_in overrides any value set for Time.

Definition at line 1712 of file MOM_state_initialization.F90.

1713  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1714  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1715  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1716  logical, intent(in) :: use_temperature !< If true, T & S are state variables.
1717  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic
1718  !! variables.
1719  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
1720  type(sponge_CS), pointer :: Layer_CSp !< A pointer that is set to point to the control
1721  !! structure for this module (in layered mode).
1722  type(ALE_sponge_CS), pointer :: ALE_CSp !< A pointer that is set to point to the control
1723  !! structure for this module (in ALE mode).
1724  type(time_type), intent(in) :: Time !< Time at the start of the run segment. Time_in
1725  !! overrides any value set for Time.
1726  ! Local variables
1727  real, allocatable, dimension(:,:,:) :: eta ! The target interface heights [Z ~> m].
1728  real, allocatable, dimension(:,:,:) :: h ! The target interface thicknesses [H ~> m or kg m-2].
1729 
1730  real, dimension (SZI_(G),SZJ_(G),SZK_(G)) :: &
1731  tmp, tmp2 ! A temporary array for tracers.
1732  real, dimension (SZI_(G),SZJ_(G)) :: &
1733  tmp_2d ! A temporary array for tracers.
1734  real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading sponge fields
1735 
1736  real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate [T-1 ~> s-1].
1737  real :: pres(SZI_(G)) ! An array of the reference pressure [R L2 T-2 ~> Pa].
1738 
1739  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1740  integer :: i, j, k, is, ie, js, je, nz
1741  integer :: isd, ied, jsd, jed
1742  integer, dimension(4) :: siz
1743  integer :: nz_data ! The size of the sponge source grid
1744  character(len=40) :: potemp_var, salin_var, Idamp_var, eta_var
1745  character(len=40) :: mdl = "initialize_sponges_file"
1746  character(len=200) :: damping_file, state_file ! Strings for filenames
1747  character(len=200) :: filename, inputdir ! Strings for file/path and path.
1748 
1749  logical :: use_ALE ! True if ALE is being used, False if in layered mode
1750  logical :: time_space_interp_sponge ! True if using sponge data which
1751  ! need to be interpolated from in both the horizontal dimension and in
1752  ! time prior to vertical remapping.
1753 
1754 
1755  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1756  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1757 
1758  pres(:) = 0.0 ; tmp(:,:,:) = 0.0 ; idamp(:,:) = 0.0
1759 
1760  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1761  inputdir = slasher(inputdir)
1762  call get_param(param_file, mdl, "SPONGE_DAMPING_FILE", damping_file, &
1763  "The name of the file with the sponge damping rates.", &
1764  fail_if_missing=.true.)
1765  call get_param(param_file, mdl, "SPONGE_STATE_FILE", state_file, &
1766  "The name of the file with the state to damp toward.", &
1767  default=damping_file)
1768  call get_param(param_file, mdl, "SPONGE_PTEMP_VAR", potemp_var, &
1769  "The name of the potential temperature variable in "//&
1770  "SPONGE_STATE_FILE.", default="PTEMP")
1771  call get_param(param_file, mdl, "SPONGE_SALT_VAR", salin_var, &
1772  "The name of the salinity variable in "//&
1773  "SPONGE_STATE_FILE.", default="SALT")
1774  call get_param(param_file, mdl, "SPONGE_ETA_VAR", eta_var, &
1775  "The name of the interface height variable in "//&
1776  "SPONGE_STATE_FILE.", default="ETA")
1777  call get_param(param_file, mdl, "SPONGE_IDAMP_VAR", idamp_var, &
1778  "The name of the inverse damping rate variable in "//&
1779  "SPONGE_DAMPING_FILE.", default="Idamp")
1780  call get_param(param_file, mdl, "USE_REGRIDDING", use_ale, do_not_log = .true.)
1781  time_space_interp_sponge = .false.
1782  call get_param(param_file, mdl, "NEW_SPONGES", time_space_interp_sponge, &
1783  "Set True if using the newer sponging code which "//&
1784  "performs on-the-fly regridding in lat-lon-time.",&
1785  "of sponge restoring data.", default=.false.)
1786  if (time_space_interp_sponge) then
1787  call mom_error(warning, " initialize_sponges: NEW_SPONGES has been deprecated. "//&
1788  "Please use INTERPOLATE_SPONGE_TIME_SPACE instead. Setting "//&
1789  "INTERPOLATE_SPONGE_TIME_SPACE = True.")
1790  endif
1791  call get_param(param_file, mdl, "INTERPOLATE_SPONGE_TIME_SPACE", time_space_interp_sponge, &
1792  "Set True if using the newer sponging code which "//&
1793  "performs on-the-fly regridding in lat-lon-time.",&
1794  "of sponge restoring data.", default=time_space_interp_sponge)
1795 
1796  filename = trim(inputdir)//trim(damping_file)
1797  call log_param(param_file, mdl, "INPUTDIR/SPONGE_DAMPING_FILE", filename)
1798  if (.not.file_exists(filename, g%Domain)) &
1799  call mom_error(fatal, " initialize_sponges: Unable to open "//trim(filename))
1800 
1801  if (time_space_interp_sponge .and. .not. use_ale) &
1802  call mom_error(fatal, " initialize_sponges: Time-varying sponges are currently unavailable in layered mode ")
1803 
1804  call mom_read_data(filename, idamp_var, idamp(:,:), g%Domain, scale=us%T_to_s)
1805 
1806  ! Now register all of the fields which are damped in the sponge.
1807  ! By default, momentum is advected vertically within the sponge, but
1808  ! momentum is typically not damped within the sponge.
1809 
1810  filename = trim(inputdir)//trim(state_file)
1811  call log_param(param_file, mdl, "INPUTDIR/SPONGE_STATE_FILE", filename)
1812  if (.not.file_exists(filename, g%Domain)) &
1813  call mom_error(fatal, " initialize_sponges: Unable to open "//trim(filename))
1814 
1815 
1816  if (.not. use_ale) then
1817  ! The first call to set_up_sponge_field is for the interface heights if in layered mode.
1818  allocate(eta(isd:ied,jsd:jed,nz+1)); eta(:,:,:) = 0.0
1819  call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
1820 
1821  do j=js,je ; do i=is,ie
1822  eta(i,j,nz+1) = -g%bathyT(i,j)
1823  enddo ; enddo
1824  do k=nz,1,-1 ; do j=js,je ; do i=is,ie
1825  if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) &
1826  eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
1827  enddo ; enddo ; enddo
1828  ! Set the inverse damping rates so that the model will know where to
1829  ! apply the sponges, along with the interface heights.
1830  call initialize_sponge(idamp, eta, g, param_file, layer_csp, gv)
1831  deallocate(eta)
1832 
1833  if ( gv%nkml>0) then
1834  ! This call to set_up_sponge_ML_density registers the target values of the
1835  ! mixed layer density, which is used in determining which layers can be
1836  ! inflated without causing static instabilities.
1837  do i=is-1,ie ; pres(i) = tv%P_Ref ; enddo
1838  eosdom(:) = eos_domain(g%HI)
1839 
1840  call mom_read_data(filename, potemp_var, tmp(:,:,:), g%Domain)
1841  call mom_read_data(filename, salin_var, tmp2(:,:,:), g%Domain)
1842 
1843  do j=js,je
1844  call calculate_density(tmp(:,j,1), tmp2(:,j,1), pres, tmp_2d(:,j), tv%eqn_of_state, eosdom)
1845  enddo
1846 
1847  call set_up_sponge_ml_density(tmp_2d, g, layer_csp)
1848  endif
1849 
1850  ! Now register all of the tracer fields which are damped in the
1851  ! sponge. By default, momentum is advected vertically within the
1852  ! sponge, but momentum is typically not damped within the sponge.
1853 
1854 
1855  ! The remaining calls to set_up_sponge_field can be in any order.
1856  if ( use_temperature) then
1857  call mom_read_data(filename, potemp_var, tmp(:,:,:), g%Domain)
1858  call set_up_sponge_field(tmp, tv%T, g, nz, layer_csp)
1859  call mom_read_data(filename, salin_var, tmp(:,:,:), g%Domain)
1860  call set_up_sponge_field(tmp, tv%S, g, nz, layer_csp)
1861  endif
1862 
1863  endif
1864 
1865 
1866  if (use_ale) then
1867  if (.not. time_space_interp_sponge) then ! ALE mode
1868  call field_size(filename,eta_var,siz,no_domain=.true.)
1869  if (siz(1) /= g%ieg-g%isg+1 .or. siz(2) /= g%jeg-g%jsg+1) &
1870  call mom_error(fatal,"initialize_sponge_file: Array size mismatch for sponge data.")
1871  nz_data = siz(3)-1
1872  allocate(eta(isd:ied,jsd:jed,nz_data+1))
1873  allocate(h(isd:ied,jsd:jed,nz_data))
1874  call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
1875  do j=js,je ; do i=is,ie
1876  eta(i,j,nz+1) = -g%bathyT(i,j)
1877  enddo ; enddo
1878  do k=nz,1,-1 ; do j=js,je ; do i=is,ie
1879  if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) &
1880  eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
1881  enddo ; enddo ; enddo
1882  do k=1,nz; do j=js,je ; do i=is,ie
1883  h(i,j,k) = gv%Z_to_H*(eta(i,j,k)-eta(i,j,k+1))
1884  enddo ; enddo ; enddo
1885  call initialize_ale_sponge(idamp, g, param_file, ale_csp, h, nz_data)
1886  deallocate(eta)
1887  deallocate(h)
1888  if (use_temperature) then
1889  allocate(tmp_tr(isd:ied,jsd:jed,nz_data))
1890  call mom_read_data(filename, potemp_var, tmp_tr(:,:,:), g%Domain)
1891  call set_up_ale_sponge_field(tmp_tr, g, tv%T, ale_csp)
1892  call mom_read_data(filename, salin_var, tmp_tr(:,:,:), g%Domain)
1893  call set_up_ale_sponge_field(tmp_tr, g, tv%S, ale_csp)
1894  deallocate(tmp_tr)
1895  endif
1896  else
1897  ! Initialize sponges without supplying sponge grid
1898  call initialize_ale_sponge(idamp, g, param_file, ale_csp)
1899  ! The remaining calls to set_up_sponge_field can be in any order.
1900  if ( use_temperature) then
1901  call set_up_ale_sponge_field(filename, potemp_var, time, g, gv, us, tv%T, ale_csp)
1902  call set_up_ale_sponge_field(filename, salin_var, time, g, gv, us, tv%S, ale_csp)
1903  endif
1904  endif
1905 
1906 
1907 
1908  endif
1909 
1910 

◆ initialize_temp_salt_fit()

subroutine mom_state_initialization::initialize_temp_salt_fit ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  T,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  S,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(eos_type), pointer  eqn_of_state,
real, intent(in)  P_Ref,
logical, intent(in), optional  just_read_params 
)
private

Initializes temperature and salinity by fitting to density.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[out]tThe potential temperature that is being initialized [degC].
[out]sThe salinity that is being initialized [ppt].
[in]usA dimensional unit scaling type
[in]param_fileA structure to parse for run-time parameters.
eqn_of_stateEquation of state structure
[in]p_refThe coordinate-density reference pressure [R L2 T-2 ~> Pa].
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 1551 of file MOM_state_initialization.F90.

1552  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1553  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1554  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: T !< The potential temperature that is
1555  !! being initialized [degC].
1556  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: S !< The salinity that is being
1557  !! initialized [ppt].
1558  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1559  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1560  !! parameters.
1561  type(EOS_type), pointer :: eqn_of_state !< Equation of state structure
1562  real, intent(in) :: P_Ref !< The coordinate-density reference pressure
1563  !! [R L2 T-2 ~> Pa].
1564  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
1565  !! only read parameters without changing h.
1566  ! Local variables
1567  real :: T0(SZK_(G)) ! Layer potential temperatures [degC]
1568  real :: S0(SZK_(G)) ! Layer salinities [degC]
1569  real :: T_Ref ! Reference Temperature [degC]
1570  real :: S_Ref ! Reference Salinity [ppt]
1571  real :: pres(SZK_(G)) ! An array of the reference pressure [R L2 T-2 ~> Pa].
1572  real :: drho_dT(SZK_(G)) ! Derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
1573  real :: drho_dS(SZK_(G)) ! Derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
1574  real :: rho_guess(SZK_(G)) ! Potential density at T0 & S0 [R ~> kg m-3].
1575  logical :: fit_salin ! If true, accept the prescribed temperature and fit the salinity.
1576  logical :: just_read ! If true, just read parameters but set nothing.
1577  character(len=40) :: mdl = "initialize_temp_salt_fit" ! This subroutine's name.
1578  integer :: i, j, k, itt, nz
1579  nz = g%ke
1580 
1581  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
1582 
1583  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1584 
1585  call get_param(param_file, mdl, "T_REF", t_ref, &
1586  "A reference temperature used in initialization.", &
1587  units="degC", fail_if_missing=.not.just_read, do_not_log=just_read)
1588  call get_param(param_file, mdl, "S_REF", s_ref, &
1589  "A reference salinity used in initialization.", units="PSU", &
1590  default=35.0, do_not_log=just_read)
1591  call get_param(param_file, mdl, "FIT_SALINITY", fit_salin, &
1592  "If true, accept the prescribed temperature and fit the "//&
1593  "salinity; otherwise take salinity and fit temperature.", &
1594  default=.false., do_not_log=just_read)
1595 
1596  if (just_read) return ! All run-time parameters have been read, so return.
1597 
1598  do k=1,nz
1599  pres(k) = p_ref ; s0(k) = s_ref
1600  t0(k) = t_ref
1601  enddo
1602 
1603  call calculate_density(t0(1), s0(1), pres(1), rho_guess(1), eqn_of_state)
1604  call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state, (/1,1/) )
1605 
1606  if (fit_salin) then
1607  ! A first guess of the layers' temperatures.
1608  do k=nz,1,-1
1609  s0(k) = max(0.0, s0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_ds(1))
1610  enddo
1611  ! Refine the guesses for each layer.
1612  do itt=1,6
1613  call calculate_density(t0, s0, pres, rho_guess, eqn_of_state)
1614  call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state)
1615  do k=1,nz
1616  s0(k) = max(0.0, s0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_ds(k))
1617  enddo
1618  enddo
1619  else
1620  ! A first guess of the layers' temperatures.
1621  do k=nz,1,-1
1622  t0(k) = t0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_dt(1)
1623  enddo
1624  do itt=1,6
1625  call calculate_density(t0, s0, pres, rho_guess, eqn_of_state)
1626  call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state)
1627  do k=1,nz
1628  t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
1629  enddo
1630  enddo
1631  endif
1632 
1633  do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
1634  t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
1635  enddo ; enddo ; enddo
1636 
1637  call calltree_leave(trim(mdl)//'()')

◆ initialize_temp_salt_from_file()

subroutine mom_state_initialization::initialize_temp_salt_from_file ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(out)  T,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(out)  S,
type(ocean_grid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
logical, intent(in), optional  just_read_params 
)
private

Initializes temperature and salinity from file.

Parameters
[in]gThe ocean's grid structure
[out]tThe potential temperature that is being initialized [degC]
[out]sThe salinity that is being initialized [ppt]
[in]param_fileA structure to parse for run-time parameters
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 1450 of file MOM_state_initialization.F90.

1451  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1452  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: T !< The potential temperature that is
1453  !! being initialized [degC]
1454  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: S !< The salinity that is
1455  !! being initialized [ppt]
1456  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1457  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
1458  !! only read parameters without changing h.
1459  ! Local variables
1460  logical :: just_read ! If true, just read parameters but set nothing.
1461  character(len=200) :: filename, salt_filename ! Full paths to input files
1462  character(len=200) :: ts_file, salt_file, inputdir ! Strings for file/path
1463  character(len=40) :: mdl = "initialize_temp_salt_from_file"
1464  character(len=64) :: temp_var, salt_var ! Temperature and salinity names in files
1465 
1466  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
1467 
1468  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1469 
1470  call get_param(param_file, mdl, "TS_FILE", ts_file, &
1471  "The initial condition file for temperature.", &
1472  fail_if_missing=.not.just_read, do_not_log=just_read)
1473  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1474  inputdir = slasher(inputdir)
1475 
1476  filename = trim(inputdir)//trim(ts_file)
1477  if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/TS_FILE", filename)
1478  call get_param(param_file, mdl, "TEMP_IC_VAR", temp_var, &
1479  "The initial condition variable for potential temperature.", &
1480  default="PTEMP", do_not_log=just_read)
1481  call get_param(param_file, mdl, "SALT_IC_VAR", salt_var, &
1482  "The initial condition variable for salinity.", &
1483  default="SALT", do_not_log=just_read)
1484  call get_param(param_file, mdl, "SALT_FILE", salt_file, &
1485  "The initial condition file for salinity.", &
1486  default=trim(ts_file), do_not_log=just_read)
1487 
1488  if (just_read) return ! All run-time parameters have been read, so return.
1489 
1490  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
1491  " initialize_temp_salt_from_file: Unable to open "//trim(filename))
1492 
1493  ! Read the temperatures and salinities from netcdf files.
1494  call mom_read_data(filename, temp_var, t(:,:,:), g%Domain)
1495 
1496  salt_filename = trim(inputdir)//trim(salt_file)
1497  if (.not.file_exists(salt_filename, g%Domain)) call mom_error(fatal, &
1498  " initialize_temp_salt_from_file: Unable to open "//trim(salt_filename))
1499 
1500  call mom_read_data(salt_filename, salt_var, s(:,:,:), g%Domain)
1501 
1502  call calltree_leave(trim(mdl)//'()')

◆ initialize_temp_salt_from_profile()

subroutine mom_state_initialization::initialize_temp_salt_from_profile ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  T,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  S,
type(ocean_grid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
logical, intent(in), optional  just_read_params 
)
private

Initializes temperature and salinity from a 1D profile.

Parameters
[in]gThe ocean's grid structure
[out]tThe potential temperature that is being initialized [degC]
[out]sThe salinity that is being initialized [ppt]
[in]param_fileA structure to parse for run-time parameters
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 1506 of file MOM_state_initialization.F90.

1507  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1508  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: T !< The potential temperature that is
1509  !! being initialized [degC]
1510  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: S !< The salinity that is
1511  !! being initialized [ppt]
1512  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1513  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
1514  !! only read parameters without changing h.
1515  ! Local variables
1516  real, dimension(SZK_(G)) :: T0, S0
1517  integer :: i, j, k
1518  logical :: just_read ! If true, just read parameters but set nothing.
1519  character(len=200) :: filename, ts_file, inputdir ! Strings for file/path
1520  character(len=40) :: mdl = "initialize_temp_salt_from_profile"
1521 
1522  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
1523 
1524  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1525 
1526  call get_param(param_file, mdl, "TS_FILE", ts_file, &
1527  "The file with the reference profiles for temperature "//&
1528  "and salinity.", fail_if_missing=.not.just_read, do_not_log=just_read)
1529 
1530  if (just_read) return ! All run-time parameters have been read, so return.
1531 
1532  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1533  inputdir = slasher(inputdir)
1534  filename = trim(inputdir)//trim(ts_file)
1535  call log_param(param_file, mdl, "INPUTDIR/TS_FILE", filename)
1536  if (.not.file_exists(filename)) call mom_error(fatal, &
1537  " initialize_temp_salt_from_profile: Unable to open "//trim(filename))
1538 
1539  ! Read the temperatures and salinities from a netcdf file.
1540  call mom_read_data(filename, "PTEMP", t0(:))
1541  call mom_read_data(filename, "SALT", s0(:))
1542 
1543  do k=1,g%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
1544  t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
1545  enddo ; enddo ; enddo
1546 
1547  call calltree_leave(trim(mdl)//'()')

◆ initialize_temp_salt_linear()

subroutine mom_state_initialization::initialize_temp_salt_linear ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  T,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  S,
type(ocean_grid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
logical, intent(in), optional  just_read_params 
)
private

Initializes T and S with linear profiles according to reference surface layer salinity and temperature and a specified range.

Remarks
Note that the linear distribution is set up with respect to the layer number, not the physical position).
Parameters
[in]gThe ocean's grid structure
[out]tThe potential temperature that is being initialized [degC]
[out]sThe salinity that is being initialized [ppt]
[in]param_fileA structure to parse for run-time parameters
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 1645 of file MOM_state_initialization.F90.

1646  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1647  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: T !< The potential temperature that is
1648  !! being initialized [degC]
1649  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: S !< The salinity that is
1650  !! being initialized [ppt]
1651  type(param_file_type), intent(in) :: param_file !< A structure to parse for
1652  !! run-time parameters
1653  logical, optional, intent(in) :: just_read_params !< If present and true,
1654  !! this call will only read
1655  !! parameters without
1656  !! changing h.
1657 
1658  integer :: k
1659  real :: delta_S, delta_T
1660  real :: S_top, T_top ! Reference salinity and temerature within surface layer
1661  real :: S_range, T_range ! Range of salinities and temperatures over the vertical
1662  real :: delta
1663  logical :: just_read ! If true, just read parameters but set nothing.
1664  character(len=40) :: mdl = "initialize_temp_salt_linear" ! This subroutine's name.
1665 
1666  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
1667 
1668  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1669  call get_param(param_file, mdl, "T_TOP", t_top, &
1670  "Initial temperature of the top surface.", &
1671  units="degC", fail_if_missing=.not.just_read, do_not_log=just_read)
1672  call get_param(param_file, mdl, "T_RANGE", t_range, &
1673  "Initial temperature difference (top-bottom).", &
1674  units="degC", fail_if_missing=.not.just_read, do_not_log=just_read)
1675  call get_param(param_file, mdl, "S_TOP", s_top, &
1676  "Initial salinity of the top surface.", &
1677  units="PSU", fail_if_missing=.not.just_read, do_not_log=just_read)
1678  call get_param(param_file, mdl, "S_RANGE", s_range, &
1679  "Initial salinity difference (top-bottom).", &
1680  units="PSU", fail_if_missing=.not.just_read, do_not_log=just_read)
1681 
1682  if (just_read) return ! All run-time parameters have been read, so return.
1683 
1684  ! Prescribe salinity
1685 ! delta_S = S_range / ( G%ke - 1.0 )
1686 ! S(:,:,1) = S_top
1687 ! do k = 2,G%ke
1688 ! S(:,:,k) = S(:,:,k-1) + delta_S
1689 ! enddo
1690  do k = 1,g%ke
1691  s(:,:,k) = s_top - s_range*((real(k)-0.5)/real(g%ke))
1692  t(:,:,k) = t_top - t_range*((real(k)-0.5)/real(g%ke))
1693  enddo
1694 
1695  ! Prescribe temperature
1696 ! delta_T = T_range / ( G%ke - 1.0 )
1697 ! T(:,:,1) = T_top
1698 ! do k = 2,G%ke
1699 ! T(:,:,k) = T(:,:,k-1) + delta_T
1700 ! enddo
1701 ! delta = 1
1702 ! T(:,:,G%ke/2 - (delta-1):G%ke/2 + delta) = 1.0
1703 
1704  call calltree_leave(trim(mdl)//'()')

◆ initialize_thickness_from_file()

subroutine mom_state_initialization::initialize_thickness_from_file ( real, dimension(szi_(g),szj_(g),szk_(gv)), intent(out)  h,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
logical, intent(in)  file_has_thickness,
logical, intent(in), optional  just_read_params 
)
private

Reads the layer thicknesses or interface heights from a file.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[out]hThe thickness that is being initialized [H ~> m or kg m-2].
[in]param_fileA structure indicating the open file to parse for model parameter values.
[in]file_has_thicknessIf true, this file contains layer thicknesses; otherwise it contains interface heights.
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 619 of file MOM_state_initialization.F90.

621  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
622  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
623  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
624  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
625  intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2].
626  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
627  !! to parse for model parameter values.
628  logical, intent(in) :: file_has_thickness !< If true, this file contains layer
629  !! thicknesses; otherwise it contains
630  !! interface heights.
631  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
632  !! only read parameters without changing h.
633 
634  ! Local variables
635  real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1) ! Interface heights, in depth units.
636  integer :: inconsistent = 0
637  logical :: correct_thickness
638  logical :: just_read ! If true, just read parameters but set nothing.
639  character(len=40) :: mdl = "initialize_thickness_from_file" ! This subroutine's name.
640  character(len=200) :: filename, thickness_file, inputdir, mesg ! Strings for file/path
641  integer :: i, j, k, is, ie, js, je, nz
642 
643  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
644 
645  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
646 
647  if (.not.just_read) &
648  call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
649 
650  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".", do_not_log=just_read)
651  inputdir = slasher(inputdir)
652  call get_param(param_file, mdl, "THICKNESS_FILE", thickness_file, &
653  "The name of the thickness file.", &
654  fail_if_missing=.not.just_read, do_not_log=just_read)
655 
656  filename = trim(inputdir)//trim(thickness_file)
657  if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/THICKNESS_FILE", filename)
658 
659  if ((.not.just_read) .and. (.not.file_exists(filename, g%Domain))) call mom_error(fatal, &
660  " initialize_thickness_from_file: Unable to open "//trim(filename))
661 
662  if (file_has_thickness) then
663  !### Consider adding a parameter to use to rescale h.
664  if (just_read) return ! All run-time parameters have been read, so return.
665  call mom_read_data(filename, "h", h(:,:,:), g%Domain, scale=gv%m_to_H)
666  else
667  call get_param(param_file, mdl, "ADJUST_THICKNESS", correct_thickness, &
668  "If true, all mass below the bottom removed if the "//&
669  "topography is shallower than the thickness input file "//&
670  "would indicate.", default=.false., do_not_log=just_read)
671  if (just_read) return ! All run-time parameters have been read, so return.
672 
673  call mom_read_data(filename, "eta", eta(:,:,:), g%Domain, scale=us%m_to_Z)
674 
675  if (correct_thickness) then
676  call adjustetatofitbathymetry(g, gv, us, eta, h)
677  else
678  do k=nz,1,-1 ; do j=js,je ; do i=is,ie
679  if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) then
680  eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
681  h(i,j,k) = gv%Angstrom_H
682  else
683  h(i,j,k) = gv%Z_to_H * (eta(i,j,k) - eta(i,j,k+1))
684  endif
685  enddo ; enddo ; enddo
686 
687  do j=js,je ; do i=is,ie
688  if (abs(eta(i,j,nz+1) + g%bathyT(i,j)) > 1.0*us%m_to_Z) &
689  inconsistent = inconsistent + 1
690  enddo ; enddo
691  call sum_across_pes(inconsistent)
692 
693  if ((inconsistent > 0) .and. (is_root_pe())) then
694  write(mesg,'("Thickness initial conditions are inconsistent ",'// &
695  '"with topography in ",I8," places.")') inconsistent
696  call mom_error(warning, mesg)
697  endif
698  endif
699 
700  endif
701  call calltree_leave(trim(mdl)//'()')

◆ initialize_thickness_list()

subroutine mom_state_initialization::initialize_thickness_list ( real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(out)  h,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
logical, intent(in), optional  just_read_params 
)
private

Initialize thickness from a 1D list.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[out]hThe thickness that is being initialized [H ~> m or kg m-2].
[in]param_fileA structure indicating the open file to parse for model parameter values.
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 841 of file MOM_state_initialization.F90.

842  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
843  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
844  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
845  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
846  intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2].
847  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
848  !! to parse for model parameter values.
849  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
850  !! only read parameters without changing h.
851  ! Local variables
852  character(len=40) :: mdl = "initialize_thickness_list" ! This subroutine's name.
853  real :: e0(SZK_(G)+1) ! The resting interface heights, in depth units [Z ~> m],
854  ! usually negative because it is positive upward.
855  real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface
856  ! positive upward, in depth units [Z ~> m].
857  logical :: just_read ! If true, just read parameters but set nothing.
858  character(len=200) :: filename, eta_file, inputdir ! Strings for file/path
859  character(len=72) :: eta_var
860  integer :: i, j, k, is, ie, js, je, nz
861 
862  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
863 
864  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
865 
866  call get_param(param_file, mdl, "INTERFACE_IC_FILE", eta_file, &
867  "The file from which horizontal mean initial conditions "//&
868  "for interface depths can be read.", fail_if_missing=.true.)
869  call get_param(param_file, mdl, "INTERFACE_IC_VAR", eta_var, &
870  "The variable name for horizontal mean initial conditions "//&
871  "for interface depths relative to mean sea level.", &
872  default="eta")
873 
874  if (just_read) return
875 
876  call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
877 
878  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
879  filename = trim(slasher(inputdir))//trim(eta_file)
880  call log_param(param_file, mdl, "INPUTDIR/INTERFACE_IC_FILE", filename)
881 
882  e0(:) = 0.0
883  call mom_read_data(filename, eta_var, e0(:), scale=us%m_to_Z)
884 
885  if ((abs(e0(1)) - 0.0) > 0.001) then
886  ! This list probably starts with the interior interface, so shift it up.
887  do k=nz+1,2,-1 ; e0(k) = e0(k-1) ; enddo
888  e0(1) = 0.0
889  endif
890 
891  if (e0(2) > e0(1)) then ! Switch to the convention for interface heights increasing upward.
892  do k=1,nz ; e0(k) = -e0(k) ; enddo
893  endif
894 
895  do j=js,je ; do i=is,ie
896  ! This sets the initial thickness (in m) of the layers. The
897  ! thicknesses are set to insure that: 1. each layer is at least an
898  ! Angstrom thick, and 2. the interfaces are where they should be
899  ! based on the resting depths and interface height perturbations,
900  ! as long at this doesn't interfere with 1.
901  eta1d(nz+1) = -g%bathyT(i,j)
902  do k=nz,1,-1
903  eta1d(k) = e0(k)
904  if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
905  eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
906  h(i,j,k) = gv%Angstrom_H
907  else
908  h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
909  endif
910  enddo
911  enddo ; enddo
912 
913  call calltree_leave(trim(mdl)//'()')

◆ initialize_thickness_uniform()

subroutine mom_state_initialization::initialize_thickness_uniform ( real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(out)  h,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(param_file_type), intent(in)  param_file,
logical, intent(in), optional  just_read_params 
)
private

Initializes thickness to be uniform.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[out]hThe thickness that is being initialized [H ~> m or kg m-2].
[in]param_fileA structure indicating the open file to parse for model parameter values.
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 786 of file MOM_state_initialization.F90.

787  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
788  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
789  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
790  intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2].
791  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
792  !! to parse for model parameter values.
793  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
794  !! only read parameters without changing h.
795  ! Local variables
796  character(len=40) :: mdl = "initialize_thickness_uniform" ! This subroutine's name.
797  real :: e0(SZK_(G)+1) ! The resting interface heights, in depth units, usually
798  ! negative because it is positive upward.
799  real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface
800  ! positive upward, in depth units.
801  logical :: just_read ! If true, just read parameters but set nothing.
802  integer :: i, j, k, is, ie, js, je, nz
803 
804  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
805 
806  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
807 
808  if (just_read) return ! This subroutine has no run-time parameters.
809 
810  call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
811 
812  if (g%max_depth<=0.) call mom_error(fatal,"initialize_thickness_uniform: "// &
813  "MAXIMUM_DEPTH has a non-sensical value! Was it set?")
814 
815  do k=1,nz
816  e0(k) = -g%max_depth * real(k-1) / real(nz)
817  enddo
818 
819  do j=js,je ; do i=is,ie
820  ! This sets the initial thickness (in m) of the layers. The
821  ! thicknesses are set to insure that: 1. each layer is at least an
822  ! Angstrom thick, and 2. the interfaces are where they should be
823  ! based on the resting depths and interface height perturbations,
824  ! as long at this doesn't interfere with 1.
825  eta1d(nz+1) = -g%bathyT(i,j)
826  do k=nz,1,-1
827  eta1d(k) = e0(k)
828  if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
829  eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
830  h(i,j,k) = gv%Angstrom_H
831  else
832  h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
833  endif
834  enddo
835  enddo ; enddo
836 
837  call calltree_leave(trim(mdl)//'()')

◆ initialize_velocity_circular()

subroutine mom_state_initialization::initialize_velocity_circular ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out)  v,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
logical, intent(in), optional  just_read_params 
)
private

Sets the initial velocity components to be circular with no flow at edges of domain and center.

Parameters
[in]gThe ocean's grid structure
[out]uThe zonal velocity that is being initialized [L T-1 ~> m s-1]
[out]vThe meridional velocity that is being initialized [L T-1 ~> m s-1]
[in]usA dimensional unit scaling type
[in]param_fileA structure indicating the open file to parse for model parameter values.
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 1387 of file MOM_state_initialization.F90.

1388  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1389  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1390  intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1391  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1392  intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1393  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1394  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1395  !! parse for model parameter values.
1396  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
1397  !! only read parameters without changing h.
1398  ! Local variables
1399  character(len=200) :: mdl = "initialize_velocity_circular"
1400  real :: circular_max_u ! The amplitude of the zonal flow [L T-1 ~> m s-1]
1401  real :: dpi ! A local variable storing pi = 3.14159265358979...
1402  real :: psi1, psi2 ! Values of the streamfunction at two points [L2 T-1 ~> m2 s-1]
1403  logical :: just_read ! If true, just read parameters but set nothing.
1404  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1405  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1406  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1407 
1408  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
1409 
1410  call get_param(param_file, mdl, "CIRCULAR_MAX_U", circular_max_u, &
1411  "The amplitude of zonal flow from which to scale the "// &
1412  "circular stream function [m s-1].", &
1413  units="m s-1", default=0., scale=us%m_s_to_L_T, do_not_log=just_read)
1414 
1415  if (just_read) return ! All run-time parameters have been read, so return.
1416 
1417  dpi=acos(0.0)*2.0 ! pi
1418 
1419  do k=1,nz ; do j=js,je ; do i=isq,ieq
1420  psi1 = my_psi(i,j)
1421  psi2 = my_psi(i,j-1)
1422  u(i,j,k) = (psi1 - psi2) / g%dy_Cu(i,j) ! *(circular_max_u*G%len_lon/(2.0*dpi))
1423  enddo ; enddo ; enddo
1424  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1425  psi1 = my_psi(i,j)
1426  psi2 = my_psi(i-1,j)
1427  v(i,j,k) = (psi2 - psi1) / g%dx_Cv(i,j) ! *(circular_max_u*G%len_lon/(2.0*dpi))
1428  enddo ; enddo ; enddo
1429 
1430  contains
1431 
1432  !> Returns the value of a circular stream function at (ig,jg) in [L2 T-1 ~> m2 s-1]
1433  real function my_psi(ig,jg)
1434  integer :: ig !< Global i-index
1435  integer :: jg !< Global j-index
1436  ! Local variables
1437  real :: x, y, r ! [nondim]
1438 
1439  x = 2.0*(g%geoLonBu(ig,jg)-g%west_lon) / g%len_lon - 1.0 ! -1<x<1
1440  y = 2.0*(g%geoLatBu(ig,jg)-g%south_lat) / g%len_lat - 1.0 ! -1<y<1
1441  r = sqrt( x**2 + y**2 ) ! Circular stream function is a function of radius only
1442  r = min(1.0, r) ! Flatten stream function in corners of box
1443  my_psi = 0.5*(1.0 - cos(dpi*r))
1444  my_psi = my_psi * (circular_max_u * g%US%m_to_L*g%len_lon*1e3 / dpi) ! len_lon is in km
1445  end function my_psi
1446 

◆ initialize_velocity_from_file()

subroutine mom_state_initialization::initialize_velocity_from_file ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out)  v,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
logical, intent(in), optional  just_read_params 
)
private

Initialize horizontal velocity components from file.

Parameters
[in]gThe ocean's grid structure
[out]uThe zonal velocity that is being initialized [L T-1 ~> m s-1]
[out]vThe meridional velocity that is being initialized [L T-1 ~> m s-1]
[in]usA dimensional unit scaling type
[in]param_fileA structure indicating the open file to parse for modelparameter values.
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 1271 of file MOM_state_initialization.F90.

1272  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1273  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1274  intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1275  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1276  intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1277  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1278  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1279  !! parse for modelparameter values.
1280  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
1281  !! only read parameters without changing h.
1282  ! Local variables
1283  character(len=40) :: mdl = "initialize_velocity_from_file" ! This subroutine's name.
1284  character(len=200) :: filename,velocity_file,inputdir ! Strings for file/path
1285  logical :: just_read ! If true, just read parameters but set nothing.
1286 
1287  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
1288 
1289  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1290 
1291  call get_param(param_file, mdl, "VELOCITY_FILE", velocity_file, &
1292  "The name of the velocity initial condition file.", &
1293  fail_if_missing=.not.just_read, do_not_log=just_read)
1294  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1295  inputdir = slasher(inputdir)
1296 
1297  if (just_read) return ! All run-time parameters have been read, so return.
1298 
1299  filename = trim(inputdir)//trim(velocity_file)
1300  call log_param(param_file, mdl, "INPUTDIR/VELOCITY_FILE", filename)
1301 
1302  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
1303  " initialize_velocity_from_file: Unable to open "//trim(filename))
1304 
1305  ! Read the velocities from a netcdf file.
1306  call mom_read_vector(filename, "u", "v", u(:,:,:), v(:,:,:), g%Domain, scale=us%m_s_to_L_T)
1307 
1308  call calltree_leave(trim(mdl)//'()')

◆ initialize_velocity_uniform()

subroutine mom_state_initialization::initialize_velocity_uniform ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out)  v,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
logical, intent(in), optional  just_read_params 
)
private

Sets the initial velocity components to uniform.

Parameters
[in]gThe ocean's grid structure
[out]uThe zonal velocity that is being initialized [L T-1 ~> m s-1]
[out]vThe meridional velocity that is being initialized [L T-1 ~> m s-1]
[in]usA dimensional unit scaling type
[in]param_fileA structure indicating the open file to parse for modelparameter values.
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 1346 of file MOM_state_initialization.F90.

1347  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1348  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1349  intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1350  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1351  intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1352  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1353  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1354  !! parse for modelparameter values.
1355  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
1356  !! only read parameters without changing h.
1357  ! Local variables
1358  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1359  real :: initial_u_const, initial_v_const
1360  logical :: just_read ! If true, just read parameters but set nothing.
1361  character(len=200) :: mdl = "initialize_velocity_uniform" ! This subroutine's name.
1362  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1363  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1364 
1365  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
1366 
1367  call get_param(param_file, mdl, "INITIAL_U_CONST", initial_u_const, &
1368  "A initial uniform value for the zonal flow.", &
1369  default=0.0, units="m s-1", scale=us%m_s_to_L_T, do_not_log=just_read)
1370  call get_param(param_file, mdl, "INITIAL_V_CONST", initial_v_const, &
1371  "A initial uniform value for the meridional flow.", &
1372  default=0.0, units="m s-1", scale=us%m_s_to_L_T, do_not_log=just_read)
1373 
1374  if (just_read) return ! All run-time parameters have been read, so return.
1375 
1376  do k=1,nz ; do j=js,je ; do i=isq,ieq
1377  u(i,j,k) = initial_u_const
1378  enddo ; enddo ; enddo
1379  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1380  v(i,j,k) = initial_v_const
1381  enddo ; enddo ; enddo
1382 

◆ initialize_velocity_zero()

subroutine mom_state_initialization::initialize_velocity_zero ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out)  v,
type(ocean_grid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
logical, intent(in), optional  just_read_params 
)
private

Initialize horizontal velocity components to zero.

Parameters
[in]gThe ocean's grid structure
[out]uThe zonal velocity that is being initialized [L T-1 ~> m s-1]
[out]vThe meridional velocity that is being initialized [L T-1 ~> m s-1]
[in]param_fileA structure indicating the open file to parse for modelparameter values.
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 1312 of file MOM_state_initialization.F90.

1313  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1314  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1315  intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1316  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1317  intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1318  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1319  !! parse for modelparameter values.
1320  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
1321  !! only read parameters without changing h.
1322  ! Local variables
1323  character(len=200) :: mdl = "initialize_velocity_zero" ! This subroutine's name.
1324  logical :: just_read ! If true, just read parameters but set nothing.
1325  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1326  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1327  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1328 
1329  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
1330 
1331  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1332 
1333  if (just_read) return ! All run-time parameters have been read, so return.
1334 
1335  do k=1,nz ; do j=js,je ; do i=isq,ieq
1336  u(i,j,k) = 0.0
1337  enddo ; enddo ; enddo
1338  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1339  v(i,j,k) = 0.0
1340  enddo ; enddo ; enddo
1341 
1342  call calltree_leave(trim(mdl)//'()')

◆ mom_initialize_state()

subroutine, public mom_state_initialization::mom_initialize_state ( real, dimension(szib_(g),szj_(g),szk_(g)), intent(out)  u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(out)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(out)  h,
type(thermo_var_ptrs), intent(inout)  tv,
type(time_type), intent(inout)  Time,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  PF,
type(directories), intent(in)  dirs,
type(mom_restart_cs), pointer  restart_CS,
type(ale_cs), pointer  ALE_CSp,
type(tracer_registry_type), pointer  tracer_Reg,
type(sponge_cs), pointer  sponge_CSp,
type(ale_sponge_cs), pointer  ALE_sponge_CSp,
type(ocean_obc_type), pointer  OBC,
type(time_type), intent(in), optional  Time_in 
)

Initialize temporally evolving fields, either as initial conditions or by reading them from a restart (or saves) file.

Parameters
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[out]uThe zonal velocity that is being
[out]vThe meridional velocity that is being
[out]hLayer thicknesses [H ~> m or kg m-2]
[in,out]tvA structure pointing to various thermodynamic variables
[in,out]timeTime at the start of the run segment.
[in]pfA structure indicating the open file to parse for model parameter values.
[in]dirsA structure containing several relevant directory paths.
restart_csA pointer to the restart control structure.
ale_cspThe ALE control structure for remapping
tracer_regA pointer to the tracer registry
sponge_cspThe layerwise sponge control structure.
ale_sponge_cspThe ALE sponge control structure.
obcThe open boundary condition control structure.
[in]time_inTime at the start of the run segment. Time_in overrides any value set for Time.

Definition at line 114 of file MOM_state_initialization.F90.

117  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
118  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
119  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
120  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
121  intent(out) :: u !< The zonal velocity that is being
122  !! initialized [L T-1 ~> m s-1]
123  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
124  intent(out) :: v !< The meridional velocity that is being
125  !! initialized [L T-1 ~> m s-1]
126  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
127  intent(out) :: h !< Layer thicknesses [H ~> m or kg m-2]
128  type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic
129  !! variables
130  type(time_type), intent(inout) :: Time !< Time at the start of the run segment.
131  type(param_file_type), intent(in) :: PF !< A structure indicating the open file to parse
132  !! for model parameter values.
133  type(directories), intent(in) :: dirs !< A structure containing several relevant
134  !! directory paths.
135  type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control
136  !! structure.
137  type(ALE_CS), pointer :: ALE_CSp !< The ALE control structure for remapping
138  type(tracer_registry_type), pointer :: tracer_Reg !< A pointer to the tracer registry
139  type(sponge_CS), pointer :: sponge_CSp !< The layerwise sponge control structure.
140  type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< The ALE sponge control structure.
141  type(ocean_OBC_type), pointer :: OBC !< The open boundary condition control structure.
142  type(time_type), optional, intent(in) :: Time_in !< Time at the start of the run segment.
143  !! Time_in overrides any value set for Time.
144  ! Local variables
145  character(len=200) :: filename ! The name of an input file.
146  character(len=200) :: filename2 ! The name of an input files.
147  character(len=200) :: inputdir ! The directory where NetCDF input files are.
148  character(len=200) :: config
149  real :: H_rescale ! A rescaling factor for thicknesses from the representation in
150  ! a restart file to the internal representation in this run.
151  real :: vel_rescale ! A rescaling factor for velocities from the representation in
152  ! a restart file to the internal representation in this run.
153  real :: dt ! The baroclinic dynamics timestep for this run [T ~> s].
154  logical :: from_Z_file, useALE
155  logical :: new_sim
156  integer :: write_geom
157  logical :: use_temperature, use_sponge, use_OBC
158  logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
159  logical :: depress_sfc ! If true, remove the mass that would be displaced
160  ! by a large surface pressure by squeezing the column.
161  logical :: trim_ic_for_p_surf ! If true, remove the mass that would be displaced
162  ! by a large surface pressure, such as with an ice sheet.
163  logical :: regrid_accelerate
164  integer :: regrid_iterations
165 ! logical :: Analytic_FV_PGF, obsol_test
166  logical :: convert
167  logical :: just_read ! If true, only read the parameters because this
168  ! is a run from a restart file; this option
169  ! allows the use of Fatal unused parameters.
170  type(EOS_type), pointer :: eos => null()
171  logical :: debug ! If true, write debugging output.
172  logical :: debug_obc ! If true, do debugging calls related to OBCs.
173  logical :: debug_layers = .false.
174  character(len=80) :: mesg
175 ! This include declares and sets the variable "version".
176 #include "version_variable.h"
177  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
178  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
179 
180  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
181  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
182  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
183  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
184 
185  call calltree_enter("MOM_initialize_state(), MOM_state_initialization.F90")
186  call log_version(pf, mdl, version, "")
187  call get_param(pf, mdl, "DEBUG", debug, default=.false.)
188  call get_param(pf, mdl, "DEBUG_OBC", debug_obc, default=.false.)
189 
190  new_sim = determine_is_new_run(dirs%input_filename, dirs%restart_input_dir, &
191  g, restart_cs)
192  just_read = .not.new_sim
193 
194  call get_param(pf, mdl, "INPUTDIR", inputdir, &
195  "The directory in which input files are found.", default=".")
196  inputdir = slasher(inputdir)
197 
198  use_temperature = associated(tv%T)
199  useale = associated(ale_csp)
200  use_eos = associated(tv%eqn_of_state)
201  use_obc = associated(obc)
202  if (use_eos) eos => tv%eqn_of_state
203 
204  !====================================================================
205  ! Initialize temporally evolving fields, either as initial
206  ! conditions or by reading them from a restart (or saves) file.
207  !====================================================================
208 
209  if (new_sim) then
210  call mom_mesg("Run initialized internally.", 3)
211 
212  if (present(time_in)) time = time_in
213  ! Otherwise leave Time at its input value.
214 
215  ! This initialization should not be needed. Certainly restricting it
216  ! to the computational domain helps detect possible uninitialized
217  ! data in halos which should be covered by the pass_var(h) later.
218  !do k = 1, nz; do j = js, je; do i = is, ie
219  ! h(i,j,k) = 0.
220  !enddo
221  endif
222 
223  ! The remaining initialization calls are done, regardless of whether the
224  ! fields are actually initialized here (if just_read=.false.) or whether it
225  ! is just to make sure that all valid parameters are read to enable the
226  ! detection of unused parameters.
227  call get_param(pf, mdl, "INIT_LAYERS_FROM_Z_FILE", from_z_file, &
228  "If true, initialize the layer thicknesses, temperatures, "//&
229  "and salinities from a Z-space file on a latitude-longitude "//&
230  "grid.", default=.false., do_not_log=just_read)
231 
232  if (from_z_file) then
233  ! Initialize thickness and T/S from z-coordinate data in a file.
234  if (.NOT.use_temperature) call mom_error(fatal,"MOM_initialize_state : "//&
235  "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true")
236 
237  call mom_temp_salt_initialize_from_z(h, tv, g, gv, us, pf, just_read_params=just_read)
238 
239  else
240  ! Initialize thickness, h.
241  call get_param(pf, mdl, "THICKNESS_CONFIG", config, &
242  "A string that determines how the initial layer "//&
243  "thicknesses are specified for a new run: \n"//&
244  " \t file - read interface heights from the file specified \n"//&
245  " \t thickness_file - read thicknesses from the file specified \n"//&
246  " \t\t by (THICKNESS_FILE).\n"//&
247  " \t coord - determined by ALE coordinate.\n"//&
248  " \t uniform - uniform thickness layers evenly distributed \n"//&
249  " \t\t between the surface and MAXIMUM_DEPTH. \n"//&
250  " \t list - read a list of positive interface depths. \n"//&
251  " \t DOME - use a slope and channel configuration for the \n"//&
252  " \t\t DOME sill-overflow test case. \n"//&
253  " \t ISOMIP - use a configuration for the \n"//&
254  " \t\t ISOMIP test case. \n"//&
255  " \t benchmark - use the benchmark test case thicknesses. \n"//&
256  " \t Neverworld - use the Neverworld test case thicknesses. \n"//&
257  " \t search - search a density profile for the interface \n"//&
258  " \t\t densities. This is not yet implemented. \n"//&
259  " \t circle_obcs - the circle_obcs test case is used. \n"//&
260  " \t DOME2D - 2D version of DOME initialization. \n"//&
261  " \t adjustment2d - 2D lock exchange thickness ICs. \n"//&
262  " \t sloshing - sloshing gravity thickness ICs. \n"//&
263  " \t seamount - no motion test with seamount ICs. \n"//&
264  " \t dumbbell - sloshing channel ICs. \n"//&
265  " \t soliton - Equatorial Rossby soliton. \n"//&
266  " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
267  " \t USER - call a user modified routine.", &
268  default="uniform", do_not_log=just_read)
269  select case (trim(config))
270  case ("file")
271  call initialize_thickness_from_file(h, g, gv, us, pf, .false., just_read_params=just_read)
272  case ("thickness_file")
273  call initialize_thickness_from_file(h, g, gv, us, pf, .true., just_read_params=just_read)
274  case ("coord")
275  if (new_sim .and. useale) then
276  call ale_initthicknesstocoord( ale_csp, g, gv, h )
277  elseif (new_sim) then
278  call mom_error(fatal, "MOM_initialize_state: USE_REGRIDDING must be True "//&
279  "for THICKNESS_CONFIG of 'coord'")
280  endif
281  case ("uniform"); call initialize_thickness_uniform(h, g, gv, pf, &
282  just_read_params=just_read)
283  case ("list"); call initialize_thickness_list(h, g, gv, us, pf, &
284  just_read_params=just_read)
285  case ("DOME"); call dome_initialize_thickness(h, g, gv, pf, &
286  just_read_params=just_read)
287  case ("ISOMIP"); call isomip_initialize_thickness(h, g, gv, us, pf, tv, &
288  just_read_params=just_read)
289  case ("benchmark"); call benchmark_initialize_thickness(h, g, gv, us, pf, &
290  tv%eqn_of_state, tv%P_Ref, just_read_params=just_read)
291  case ("Neverwoorld","Neverland"); call neverworld_initialize_thickness(h, g, gv, us, pf, &
292  tv%eqn_of_state, tv%P_Ref)
293  case ("search"); call initialize_thickness_search
294  case ("circle_obcs"); call circle_obcs_initialize_thickness(h, g, gv, pf, &
295  just_read_params=just_read)
296  case ("lock_exchange"); call lock_exchange_initialize_thickness(h, g, gv, us, &
297  pf, just_read_params=just_read)
298  case ("external_gwave"); call external_gwave_initialize_thickness(h, g, gv, us, &
299  pf, just_read_params=just_read)
300  case ("DOME2D"); call dome2d_initialize_thickness(h, g, gv, us, pf, &
301  just_read_params=just_read)
302  case ("adjustment2d"); call adjustment_initialize_thickness(h, g, gv, us, &
303  pf, just_read_params=just_read)
304  case ("sloshing"); call sloshing_initialize_thickness(h, g, gv, us, pf, &
305  just_read_params=just_read)
306  case ("seamount"); call seamount_initialize_thickness(h, g, gv, us, pf, &
307  just_read_params=just_read)
308  case ("dumbbell"); call dumbbell_initialize_thickness(h, g, gv, us, pf, &
309  just_read_params=just_read)
310  case ("soliton"); call soliton_initialize_thickness(h, g, gv, us)
311  case ("phillips"); call phillips_initialize_thickness(h, g, gv, us, pf, &
312  just_read_params=just_read)
313  case ("rossby_front"); call rossby_front_initialize_thickness(h, g, gv, us, &
314  pf, just_read_params=just_read)
315  case ("USER"); call user_initialize_thickness(h, g, gv, pf, &
316  just_read_params=just_read)
317  case default ; call mom_error(fatal, "MOM_initialize_state: "//&
318  "Unrecognized layer thickness configuration "//trim(config))
319  end select
320 
321  ! Initialize temperature and salinity (T and S).
322  if ( use_temperature ) then
323  call get_param(pf, mdl, "TS_CONFIG", config, &
324  "A string that determines how the initial tempertures "//&
325  "and salinities are specified for a new run: \n"//&
326  " \t file - read velocities from the file specified \n"//&
327  " \t\t by (TS_FILE). \n"//&
328  " \t fit - find the temperatures that are consistent with \n"//&
329  " \t\t the layer densities and salinity S_REF. \n"//&
330  " \t TS_profile - use temperature and salinity profiles \n"//&
331  " \t\t (read from TS_FILE) to set layer densities. \n"//&
332  " \t benchmark - use the benchmark test case T & S. \n"//&
333  " \t linear - linear in logical layer space. \n"//&
334  " \t DOME2D - 2D DOME initialization. \n"//&
335  " \t ISOMIP - ISOMIP initialization. \n"//&
336  " \t adjustment2d - 2d lock exchange T/S ICs. \n"//&
337  " \t sloshing - sloshing mode T/S ICs. \n"//&
338  " \t seamount - no motion test with seamount ICs. \n"//&
339  " \t dumbbell - sloshing channel ICs. \n"//&
340  " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
341  " \t SCM_CVMix_tests - used in the SCM CVMix tests.\n"//&
342  " \t USER - call a user modified routine.", &
343  fail_if_missing=new_sim, do_not_log=just_read)
344 ! " \t baroclinic_zone - an analytic baroclinic zone. \n"//&
345  select case (trim(config))
346  case ("fit"); call initialize_temp_salt_fit(tv%T, tv%S, g, gv, us, pf, &
347  eos, tv%P_Ref, just_read_params=just_read)
348  case ("file"); call initialize_temp_salt_from_file(tv%T, tv%S, g, &
349  pf, just_read_params=just_read)
350  case ("benchmark"); call benchmark_init_temperature_salinity(tv%T, tv%S, &
351  g, gv, us, pf, eos, tv%P_Ref, just_read_params=just_read)
352  case ("TS_profile") ; call initialize_temp_salt_from_profile(tv%T, tv%S, &
353  g, pf, just_read_params=just_read)
354  case ("linear"); call initialize_temp_salt_linear(tv%T, tv%S, g, pf, &
355  just_read_params=just_read)
356  case ("DOME2D"); call dome2d_initialize_temperature_salinity ( tv%T, &
357  tv%S, h, g, gv, pf, eos, just_read_params=just_read)
358  case ("ISOMIP"); call isomip_initialize_temperature_salinity ( tv%T, &
359  tv%S, h, g, gv, us, pf, eos, just_read_params=just_read)
360  case ("adjustment2d"); call adjustment_initialize_temperature_salinity ( tv%T, &
361  tv%S, h, g, gv, pf, eos, just_read_params=just_read)
362  case ("baroclinic_zone"); call baroclinic_zone_init_temperature_salinity( tv%T, &
363  tv%S, h, g, gv, us, pf, just_read_params=just_read)
364  case ("sloshing"); call sloshing_initialize_temperature_salinity(tv%T, &
365  tv%S, h, g, gv, pf, eos, just_read_params=just_read)
366  case ("seamount"); call seamount_initialize_temperature_salinity(tv%T, &
367  tv%S, h, g, gv, pf, eos, just_read_params=just_read)
368  case ("dumbbell"); call dumbbell_initialize_temperature_salinity(tv%T, &
369  tv%S, h, g, gv, pf, eos, just_read_params=just_read)
370  case ("rossby_front"); call rossby_front_initialize_temperature_salinity ( tv%T, &
371  tv%S, h, g, gv, pf, eos, just_read_params=just_read)
372  case ("SCM_CVMix_tests"); call scm_cvmix_tests_ts_init(tv%T, tv%S, h, &
373  g, gv, us, pf, just_read_params=just_read)
374  case ("dense"); call dense_water_initialize_ts(g, gv, pf, eos, tv%T, tv%S, &
375  h, just_read_params=just_read)
376  case ("USER"); call user_init_temperature_salinity(tv%T, tv%S, g, pf, eos, &
377  just_read_params=just_read)
378  case default ; call mom_error(fatal, "MOM_initialize_state: "//&
379  "Unrecognized Temp & salt configuration "//trim(config))
380  end select
381  endif
382  endif ! not from_Z_file.
383  if (use_temperature .and. use_obc) &
384  call fill_temp_salt_segments(g, obc, tv)
385 
386  ! The thicknesses in halo points might be needed to initialize the velocities.
387  if (new_sim) call pass_var(h, g%Domain)
388 
389  ! Initialize velocity components, u and v
390  call get_param(pf, mdl, "VELOCITY_CONFIG", config, &
391  "A string that determines how the initial velocities "//&
392  "are specified for a new run: \n"//&
393  " \t file - read velocities from the file specified \n"//&
394  " \t\t by (VELOCITY_FILE). \n"//&
395  " \t zero - the fluid is initially at rest. \n"//&
396  " \t uniform - the flow is uniform (determined by\n"//&
397  " \t\t parameters INITIAL_U_CONST and INITIAL_V_CONST).\n"//&
398  " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
399  " \t soliton - Equatorial Rossby soliton.\n"//&
400  " \t USER - call a user modified routine.", default="zero", &
401  do_not_log=just_read)
402  select case (trim(config))
403  case ("file"); call initialize_velocity_from_file(u, v, g, us, pf, &
404  just_read_params=just_read)
405  case ("zero"); call initialize_velocity_zero(u, v, g, pf, &
406  just_read_params=just_read)
407  case ("uniform"); call initialize_velocity_uniform(u, v, g, us, pf, &
408  just_read_params=just_read)
409  case ("circular"); call initialize_velocity_circular(u, v, g, us, pf, &
410  just_read_params=just_read)
411  case ("phillips"); call phillips_initialize_velocity(u, v, g, gv, us, pf, &
412  just_read_params=just_read)
413  case ("rossby_front"); call rossby_front_initialize_velocity(u, v, h, &
414  g, gv, us, pf, just_read_params=just_read)
415  case ("soliton"); call soliton_initialize_velocity(u, v, h, g, us)
416  case ("USER"); call user_initialize_velocity(u, v, g, us, pf, &
417  just_read_params=just_read)
418  case default ; call mom_error(fatal, "MOM_initialize_state: "//&
419  "Unrecognized velocity configuration "//trim(config))
420  end select
421 
422  if (new_sim) call pass_vector(u, v, g%Domain)
423  if (debug .and. new_sim) then
424  call uvchksum("MOM_initialize_state [uv]", u, v, g%HI, haloshift=1, scale=us%m_s_to_L_T)
425  endif
426 
427  ! Optionally convert the thicknesses from m to kg m-2. This is particularly
428  ! useful in a non-Boussinesq model.
429  call get_param(pf, mdl, "CONVERT_THICKNESS_UNITS", convert, &
430  "If true, convert the thickness initial conditions from "//&
431  "units of m to kg m-2 or vice versa, depending on whether "//&
432  "BOUSSINESQ is defined. This does not apply if a restart "//&
433  "file is read.", default=.not.gv%Boussinesq, do_not_log=just_read)
434 
435  if (new_sim .and. convert .and. .not.gv%Boussinesq) &
436  ! Convert thicknesses from geomtric distances to mass-per-unit-area.
437  call convert_thickness(h, g, gv, us, tv)
438 
439  ! Remove the mass that would be displaced by an ice shelf or inverse barometer.
440  call get_param(pf, mdl, "DEPRESS_INITIAL_SURFACE", depress_sfc, &
441  "If true, depress the initial surface to avoid huge "//&
442  "tsunamis when a large surface pressure is applied.", &
443  default=.false., do_not_log=just_read)
444  call get_param(pf, mdl, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, &
445  "If true, cuts way the top of the column for initial conditions "//&
446  "at the depth where the hydrostatic pressure matches the imposed "//&
447  "surface pressure which is read from file.", default=.false., &
448  do_not_log=just_read)
449  if (depress_sfc .and. trim_ic_for_p_surf) call mom_error(fatal, "MOM_initialize_state: "//&
450  "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True")
451  if (new_sim .and. debug .and. (depress_sfc .or. trim_ic_for_p_surf)) &
452  call hchksum(h, "Pre-depress: h ", g%HI, haloshift=1, scale=gv%H_to_m)
453  if (depress_sfc) call depress_surface(h, g, gv, us, pf, tv, just_read_params=just_read)
454  if (trim_ic_for_p_surf) call trim_for_ice(pf, g, gv, us, ale_csp, tv, h, just_read_params=just_read)
455 
456  ! Perhaps we want to run the regridding coordinate generator for multiple
457  ! iterations here so the initial grid is consistent with the coordinate
458  if (useale) then
459  call get_param(pf, mdl, "REGRID_ACCELERATE_INIT", regrid_accelerate, &
460  "If true, runs REGRID_ACCELERATE_ITERATIONS iterations of the regridding "//&
461  "algorithm to push the initial grid to be consistent with the initial "//&
462  "condition. Useful only for state-based and iterative coordinates.", &
463  default=.false., do_not_log=just_read)
464  if (regrid_accelerate) then
465  call get_param(pf, mdl, "REGRID_ACCELERATE_ITERATIONS", regrid_iterations, &
466  "The number of regridding iterations to perform to generate "//&
467  "an initial grid that is consistent with the initial conditions.", &
468  default=1, do_not_log=just_read)
469 
470  call get_param(pf, mdl, "DT", dt, "Timestep", fail_if_missing=.true., scale=us%s_to_T)
471 
472  if (new_sim .and. debug) &
473  call hchksum(h, "Pre-ALE_regrid: h ", g%HI, haloshift=1, scale=gv%H_to_m)
474  call ale_regrid_accelerated(ale_csp, g, gv, h, tv, regrid_iterations, u, v, obc, tracer_reg, &
475  dt=dt, initial=.true.)
476  endif
477  endif
478  ! This is the end of the block of code that might have initialized fields
479  ! internally at the start of a new run.
480 
481  if (.not.new_sim) then ! This block restores the state from a restart file.
482  ! This line calls a subroutine that reads the initial conditions
483  ! from a previously generated file.
484  call restore_state(dirs%input_filename, dirs%restart_input_dir, time, &
485  g, restart_cs)
486  if (present(time_in)) time = time_in
487  if ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
488  h_rescale = gv%m_to_H / gv%m_to_H_restart
489  do k=1,nz ; do j=js,je ; do i=is,ie ; h(i,j,k) = h_rescale * h(i,j,k) ; enddo ; enddo ; enddo
490  endif
491  if ( (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
492  ((us%m_to_L * us%s_to_T_restart) /= (us%m_to_L_restart * us%s_to_T)) ) then
493  vel_rescale = (us%m_to_L * us%s_to_T_restart) / (us%m_to_L_restart * us%s_to_T)
494  do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb ; u(i,j,k) = vel_rescale * u(i,j,k) ; enddo ; enddo ; enddo
495  do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied ; v(i,j,k) = vel_rescale * v(i,j,k) ; enddo ; enddo ; enddo
496  endif
497  endif
498 
499  if ( use_temperature ) then
500  call pass_var(tv%T, g%Domain, complete=.false.)
501  call pass_var(tv%S, g%Domain, complete=.false.)
502  endif
503  call pass_var(h, g%Domain)
504 
505  if (debug) then
506  call hchksum(h, "MOM_initialize_state: h ", g%HI, haloshift=1, scale=gv%H_to_m)
507  if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", g%HI, haloshift=1)
508  if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", g%HI, haloshift=1)
509  if ( use_temperature .and. debug_layers) then ; do k=1,nz
510  write(mesg,'("MOM_IS: T[",I2,"]")') k
511  call hchksum(tv%T(:,:,k), mesg, g%HI, haloshift=1)
512  write(mesg,'("MOM_IS: S[",I2,"]")') k
513  call hchksum(tv%S(:,:,k), mesg, g%HI, haloshift=1)
514  enddo ; endif
515  endif
516 
517  call get_param(pf, mdl, "SPONGE", use_sponge, &
518  "If true, sponges may be applied anywhere in the domain. "//&
519  "The exact location and properties of those sponges are "//&
520  "specified via SPONGE_CONFIG.", default=.false.)
521  if ( use_sponge ) then
522  call get_param(pf, mdl, "SPONGE_CONFIG", config, &
523  "A string that sets how the sponges are configured: \n"//&
524  " \t file - read sponge properties from the file \n"//&
525  " \t\t specified by (SPONGE_FILE).\n"//&
526  " \t ISOMIP - apply ale sponge in the ISOMIP case \n"//&
527  " \t RGC - apply sponge in the rotating_gravity_current case \n"//&
528  " \t DOME - use a slope and channel configuration for the \n"//&
529  " \t\t DOME sill-overflow test case. \n"//&
530  " \t BFB - Sponge at the southern boundary of the domain\n"//&
531  " \t\t for buoyancy-forced basin case.\n"//&
532  " \t USER - call a user modified routine.", default="file")
533  select case (trim(config))
534  case ("DOME"); call dome_initialize_sponges(g, gv, us, tv, pf, sponge_csp)
535  case ("DOME2D"); call dome2d_initialize_sponges(g, gv, us, tv, pf, useale, &
536  sponge_csp, ale_sponge_csp)
537  case ("ISOMIP"); call isomip_initialize_sponges(g, gv, us, tv, pf, useale, &
538  sponge_csp, ale_sponge_csp)
539  case("RGC"); call rgc_initialize_sponges(g, gv, us, tv, u, v, pf, useale, &
540  sponge_csp, ale_sponge_csp)
541  case ("USER"); call user_initialize_sponges(g, gv, use_temperature, tv, pf, sponge_csp, h)
542  case ("BFB"); call bfb_initialize_sponges_southonly(g, gv, us, use_temperature, tv, pf, &
543  sponge_csp, h)
544  case ("DUMBBELL"); call dumbbell_initialize_sponges(g, gv, us, tv, pf, useale, &
545  sponge_csp, ale_sponge_csp)
546  case ("phillips"); call phillips_initialize_sponges(g, gv, us, tv, pf, sponge_csp, h)
547  case ("dense"); call dense_water_initialize_sponges(g, gv, us, tv, pf, useale, &
548  sponge_csp, ale_sponge_csp)
549  case ("file"); call initialize_sponges_file(g, gv, us, use_temperature, tv, pf, &
550  sponge_csp, ale_sponge_csp, time)
551  case default ; call mom_error(fatal, "MOM_initialize_state: "//&
552  "Unrecognized sponge configuration "//trim(config))
553  end select
554  endif
555 
556  ! Reads OBC parameters not pertaining to the location of the boundaries
557  call open_boundary_init(g, gv, us, pf, obc, restart_cs)
558 
559  ! This controls user code for setting open boundary data
560  if (associated(obc)) then
561  call initialize_segment_data(g, obc, pf) ! call initialize_segment_data(G, OBC, param_file)
562 ! call open_boundary_config(G, US, PF, OBC)
563  ! Call this once to fill boundary arrays from fixed values
564  if (.not. obc%needs_IO_for_data) &
565  call update_obc_segment_data(g, gv, us, obc, tv, h, time)
566 
567  call get_param(pf, mdl, "OBC_USER_CONFIG", config, &
568  "A string that sets how the user code is invoked to set open boundary data: \n"//&
569  " DOME - specified inflow on northern boundary\n"//&
570  " dyed_channel - supercritical with dye on the inflow boundary\n"//&
571  " dyed_obcs - circle_obcs with dyes on the open boundaries\n"//&
572  " Kelvin - barotropic Kelvin wave forcing on the western boundary\n"//&
573  " shelfwave - Flather with shelf wave forcing on western boundary\n"//&
574  " supercritical - now only needed here for the allocations\n"//&
575  " tidal_bay - Flather with tidal forcing on eastern boundary\n"//&
576  " USER - user specified", default="none")
577  if (trim(config) == "DOME") then
578  call dome_set_obc_data(obc, tv, g, gv, us, pf, tracer_reg)
579  elseif (trim(config) == "dyed_channel") then
580  call dyed_channel_set_obc_tracer_data(obc, g, gv, pf, tracer_reg)
581  obc%update_OBC = .true.
582  elseif (trim(config) == "dyed_obcs") then
583  call dyed_obcs_set_obc_data(obc, g, gv, pf, tracer_reg)
584  elseif (trim(config) == "Kelvin") then
585  obc%update_OBC = .true.
586  elseif (trim(config) == "shelfwave") then
587  obc%update_OBC = .true.
588  elseif (lowercase(trim(config)) == "supercritical") then
589  call supercritical_set_obc_data(obc, g, pf)
590  elseif (trim(config) == "tidal_bay") then
591  obc%update_OBC = .true.
592  elseif (trim(config) == "USER") then
593  call user_set_obc_data(obc, tv, g, pf, tracer_reg)
594  elseif (.not. trim(config) == "none") then
595  call mom_error(fatal, "The open boundary conditions specified by "//&
596  "OBC_USER_CONFIG = "//trim(config)//" have not been fully implemented.")
597  endif
598  if (open_boundary_query(obc, apply_open_obc=.true.)) then
599  call set_tracer_data(obc, tv, h, g, pf, tracer_reg)
600  endif
601  endif
602 ! if (open_boundary_query(OBC, apply_nudged_OBC=.true.)) then
603 ! call set_3D_OBC_data(OBC, tv, h, G, PF, tracer_Reg)
604 ! endif
605  ! Still need a way to specify the boundary values
606  if (debug.and.associated(obc)) then
607  call hchksum(g%mask2dT, 'MOM_initialize_state: mask2dT ', g%HI)
608  call uvchksum('MOM_initialize_state: mask2dC[uv]', g%mask2dCu, &
609  g%mask2dCv, g%HI)
610  call qchksum(g%mask2dBu, 'MOM_initialize_state: mask2dBu ', g%HI)
611  endif
612 
613  if (debug_obc) call open_boundary_test_extern_h(g, gv, obc, h)
614  call calltree_leave('MOM_initialize_state()')
615 

◆ mom_state_init_tests()

subroutine mom_state_initialization::mom_state_init_tests ( type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(thermo_var_ptrs), intent(in)  tv 
)
private

Run simple unit tests.

Parameters
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]tvThermodynamics structure.

Definition at line 2559 of file MOM_state_initialization.F90.

2560  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2561  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2562  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2563  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
2564 
2565  ! Local variables
2566  integer, parameter :: nk=5
2567  real, dimension(nk) :: T, T_t, T_b ! Temperatures [degC]
2568  real, dimension(nk) :: S, S_t, S_b ! Salinities [ppt]
2569  real, dimension(nk) :: rho ! Layer density [R ~> kg m-3]
2570  real, dimension(nk) :: h ! Layer thicknesses [H ~> m or kg m-2]
2571  real, dimension(nk) :: z ! Height of layer center [Z ~> m]
2572  real, dimension(nk+1) :: e ! Interface heights [Z ~> m]
2573  integer :: k
2574  real :: P_tot, P_t, P_b ! Pressures [R L2 T-2 ~> Pa]
2575  real :: z_out ! Output height [Z ~> m]
2576  real :: I_z_scale ! The inverse of the height scale for prescribed gradients [Z-1 ~> m-1]
2577  type(remapping_CS), pointer :: remap_CS => null()
2578 
2579  i_z_scale = 1.0 / (500.0*us%m_to_Z)
2580  do k = 1, nk
2581  h(k) = 100.0*gv%m_to_H
2582  enddo
2583  e(1) = 0.
2584  do k = 1, nk
2585  e(k+1) = e(k) - gv%H_to_Z * h(k)
2586  enddo
2587  p_tot = 0.
2588  do k = 1, nk
2589  z(k) = 0.5 * ( e(k) + e(k+1) )
2590  t_t(k) = 20. + (0. * i_z_scale) * e(k)
2591  t(k) = 20. + (0. * i_z_scale)*z(k)
2592  t_b(k) = 20. + (0. * i_z_scale)*e(k+1)
2593  s_t(k) = 35. - (0. * i_z_scale)*e(k)
2594  s(k) = 35. + (0. * i_z_scale)*z(k)
2595  s_b(k) = 35. - (0. * i_z_scale)*e(k+1)
2596  call calculate_density(0.5*(t_t(k)+t_b(k)), 0.5*(s_t(k)+s_b(k)), -gv%Rho0*gv%g_Earth*us%m_to_Z*z(k), &
2597  rho(k), tv%eqn_of_state)
2598  p_tot = p_tot + gv%g_Earth * rho(k) * gv%H_to_Z*h(k)
2599  enddo
2600 
2601  p_t = 0.
2602  do k = 1, nk
2603  call find_depth_of_pressure_in_cell(t_t(k), t_b(k), s_t(k), s_b(k), e(k), e(k+1), p_t, 0.5*p_tot, &
2604  gv%Rho0, gv%g_Earth, tv%eqn_of_state, us, p_b, z_out)
2605  write(0,*) k, us%RL2_T2_to_Pa*p_t, us%RL2_T2_to_Pa*p_b, 0.5*us%RL2_T2_to_Pa*p_tot, &
2606  us%Z_to_m*e(k), us%Z_to_m*e(k+1), us%Z_to_m*z_out
2607  p_t = p_b
2608  enddo
2609  write(0,*) us%RL2_T2_to_Pa*p_b, us%RL2_T2_to_Pa*p_tot
2610 
2611  write(0,*) ''
2612  write(0,*) ' ==================================================================== '
2613  write(0,*) ''
2614  write(0,*) gv%H_to_m*h
2615  call cut_off_column_top(nk, tv, gv, us, gv%g_Earth, -e(nk+1), gv%Angstrom_Z, &
2616  t, t_t, t_b, s, s_t, s_b, 0.5*p_tot, h, remap_cs)
2617  write(0,*) gv%H_to_m*h
2618 

◆ mom_temp_salt_initialize_from_z()

subroutine mom_state_initialization::mom_temp_salt_initialize_from_z ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(out)  h,
type(thermo_var_ptrs), intent(inout)  tv,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  PF,
logical, intent(in), optional  just_read_params 
)
private

This subroutine determines the isopycnal or other coordinate interfaces and layer potential temperatures and salinities directly from a z-space file on a latitude-longitude grid.

Parameters
[in,out]gThe ocean's grid structure
[out]hLayer thicknesses being initialized [H ~> m or kg m-2]
[in,out]tvA structure pointing to various thermodynamic variables including temperature and salinity
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in]pfA structure indicating the open file to parse for model parameter values.
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 1970 of file MOM_state_initialization.F90.

1971  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1972  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1973  intent(out) :: h !< Layer thicknesses being initialized [H ~> m or kg m-2]
1974  type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic
1975  !! variables including temperature and salinity
1976  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1977  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1978  type(param_file_type), intent(in) :: PF !< A structure indicating the open file
1979  !! to parse for model parameter values.
1980  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
1981  !! only read parameters without changing h.
1982 
1983  ! Local variables
1984  character(len=200) :: filename !< The name of an input file containing temperature
1985  !! and salinity in z-space; also used for ice shelf area.
1986  character(len=200) :: tfilename !< The name of an input file containing only temperature
1987  !! in z-space.
1988  character(len=200) :: sfilename !< The name of an input file containing only salinity
1989  !! in z-space.
1990  character(len=200) :: shelf_file !< The name of an input file used for ice shelf area.
1991  character(len=200) :: inputdir !! The directory where NetCDF input filesare.
1992  character(len=200) :: mesg, area_varname, ice_shelf_file
1993 
1994  type(EOS_type), pointer :: eos => null()
1995  type(thermo_var_ptrs) :: tv_loc ! A temporary thermo_var container
1996  type(verticalGrid_type) :: GV_loc ! A temporary vertical grid structure
1997  ! This include declares and sets the variable "version".
1998 # include "version_variable.h"
1999  character(len=40) :: mdl = "MOM_initialize_layers_from_Z" ! This module's name.
2000 
2001  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
2002  integer :: is, ie, js, je, nz ! compute domain indices
2003  integer :: isc,iec,jsc,jec ! global compute domain indices
2004  integer :: isg, ieg, jsg, jeg ! global extent
2005  integer :: isd, ied, jsd, jed ! data domain indices
2006 
2007  integer :: i, j, k, ks, np, ni, nj
2008  integer :: nkml ! The number of layers in the mixed layer.
2009 
2010  integer :: kd, inconsistent
2011  integer :: nkd ! number of levels to use for regridding input arrays
2012  real :: eps_Z ! A negligibly thin layer thickness [Z ~> m].
2013  real :: eps_rho ! A negligibly small density difference [R ~> kg m-3].
2014  real :: PI_180 ! for conversion from degrees to radians
2015 
2016  real, dimension(:,:), pointer :: shelf_area => null()
2017  real :: Hmix_default ! The default initial mixed layer depth [m].
2018  real :: Hmix_depth ! The mixed layer depth in the initial condition [Z ~> m].
2019  real :: dilate ! A dilation factor to match topography [nondim]
2020  real :: missing_value_temp, missing_value_salt
2021  logical :: correct_thickness
2022  character(len=40) :: potemp_var, salin_var
2023  character(len=8) :: laynum
2024 
2025  integer, parameter :: niter=10 ! number of iterations for t/s adjustment to layer density
2026  logical :: just_read ! If true, just read parameters but set nothing.
2027  logical :: adjust_temperature = .true. ! fit t/s to target densities
2028  real, parameter :: missing_value = -1.e20
2029  real, parameter :: temp_land_fill = 0.0, salt_land_fill = 35.0
2030  logical :: reentrant_x, tripolar_n,dbg
2031  logical :: debug = .false. ! manually set this to true for verbose output
2032 
2033  ! data arrays
2034  real, dimension(:), allocatable :: z_edges_in, z_in ! Interface heights [Z ~> m]
2035  real, dimension(:), allocatable :: Rb ! Interface densities [R ~> kg m-3]
2036  real, dimension(:,:,:), allocatable, target :: temp_z, salt_z, mask_z
2037  real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3]
2038  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: zi ! Interface heights [Z ~> m].
2039  integer, dimension(SZI_(G),SZJ_(G)) :: nlevs
2040  real, dimension(SZI_(G)) :: press ! Pressures [R L2 T-2 ~> Pa].
2041 
2042  ! Local variables for ALE remapping
2043  real, dimension(:), allocatable :: hTarget ! Target thicknesses [Z ~> m].
2044  real, dimension(:,:), allocatable :: area_shelf_h ! Shelf-covered area per grid cell [L2 ~> m2]
2045  real, dimension(:,:), allocatable, target :: frac_shelf_h ! Fractional shelf area per grid cell [nondim]
2046  real, dimension(:,:,:), allocatable, target :: tmpT1dIn, tmpS1dIn
2047  real, dimension(:,:,:), allocatable :: tmp_mask_in
2048  real, dimension(:,:,:), allocatable :: h1 ! Thicknesses [H ~> m or kg m-2].
2049  real, dimension(:,:,:), allocatable :: dz_interface ! Change in position of interface due to regridding
2050  real :: zTopOfCell, zBottomOfCell ! Heights in Z units [Z ~> m].
2051  type(regridding_CS) :: regridCS ! Regridding parameters and work arrays
2052  type(remapping_CS) :: remapCS ! Remapping parameters and work arrays
2053 
2054  logical :: homogenize, useALEremapping, remap_full_column, remap_general, remap_old_alg
2055  logical :: answers_2018, default_2018_answers, hor_regrid_answers_2018
2056  logical :: use_ice_shelf
2057  logical :: pre_gridded
2058  logical :: separate_mixed_layer ! If true, handle the mixed layers differently.
2059  character(len=10) :: remappingScheme
2060  real :: tempAvg, saltAvg
2061  integer :: nPoints, ans
2062  integer :: id_clock_routine, id_clock_read, id_clock_interp, id_clock_fill, id_clock_ALE
2063 
2064  id_clock_routine = cpu_clock_id('(Initialize from Z)', grain=clock_routine)
2065  id_clock_ale = cpu_clock_id('(Initialize from Z) ALE', grain=clock_loop)
2066 
2067  call cpu_clock_begin(id_clock_routine)
2068 
2069  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2070  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2071  isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
2072 
2073  pi_180=atan(1.0)/45.
2074 
2075  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
2076 
2077  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
2078  if (.not.just_read) call log_version(pf, mdl, version, "")
2079 
2080  inputdir = "." ; call get_param(pf, mdl, "INPUTDIR", inputdir)
2081  inputdir = slasher(inputdir)
2082 
2083  eos => tv%eqn_of_state
2084 
2085  reentrant_x = .false. ; call get_param(pf, mdl, "REENTRANT_X", reentrant_x, default=.true.)
2086  tripolar_n = .false. ; call get_param(pf, mdl, "TRIPOLAR_N", tripolar_n, default=.false.)
2087 
2088  call get_param(pf, mdl, "TEMP_SALT_Z_INIT_FILE",filename, &
2089  "The name of the z-space input file used to initialize "//&
2090  "temperatures (T) and salinities (S). If T and S are not "//&
2091  "in the same file, TEMP_Z_INIT_FILE and SALT_Z_INIT_FILE "//&
2092  "must be set.",default="temp_salt_z.nc",do_not_log=just_read)
2093  call get_param(pf, mdl, "TEMP_Z_INIT_FILE",tfilename, &
2094  "The name of the z-space input file used to initialize "//&
2095  "temperatures, only.", default=trim(filename),do_not_log=just_read)
2096  call get_param(pf, mdl, "SALT_Z_INIT_FILE",sfilename, &
2097  "The name of the z-space input file used to initialize "//&
2098  "temperatures, only.", default=trim(filename),do_not_log=just_read)
2099  filename = trim(inputdir)//trim(filename)
2100  tfilename = trim(inputdir)//trim(tfilename)
2101  sfilename = trim(inputdir)//trim(sfilename)
2102  call get_param(pf, mdl, "Z_INIT_FILE_PTEMP_VAR", potemp_var, &
2103  "The name of the potential temperature variable in "//&
2104  "TEMP_Z_INIT_FILE.", default="ptemp",do_not_log=just_read)
2105  call get_param(pf, mdl, "Z_INIT_FILE_SALT_VAR", salin_var, &
2106  "The name of the salinity variable in "//&
2107  "SALT_Z_INIT_FILE.", default="salt",do_not_log=just_read)
2108  call get_param(pf, mdl, "Z_INIT_HOMOGENIZE", homogenize, &
2109  "If True, then horizontally homogenize the interpolated "//&
2110  "initial conditions.", default=.false., do_not_log=just_read)
2111  call get_param(pf, mdl, "Z_INIT_ALE_REMAPPING", usealeremapping, &
2112  "If True, then remap straight to model coordinate from file.", &
2113  default=.false., do_not_log=just_read)
2114  call get_param(pf, mdl, "Z_INIT_REMAPPING_SCHEME", remappingscheme, &
2115  "The remapping scheme to use if using Z_INIT_ALE_REMAPPING "//&
2116  "is True.", default="PPM_IH4", do_not_log=just_read)
2117  call get_param(pf, mdl, "Z_INIT_REMAP_GENERAL", remap_general, &
2118  "If false, only initializes to z* coordinates. "//&
2119  "If true, allows initialization directly to general coordinates.",&
2120  default=.false., do_not_log=just_read)
2121  call get_param(pf, mdl, "Z_INIT_REMAP_FULL_COLUMN", remap_full_column, &
2122  "If false, only reconstructs profiles for valid data points. "//&
2123  "If true, inserts vanished layers below the valid data.", &
2124  default=remap_general, do_not_log=just_read)
2125  call get_param(pf, mdl, "Z_INIT_REMAP_OLD_ALG", remap_old_alg, &
2126  "If false, uses the preferred remapping algorithm for initialization. "//&
2127  "If true, use an older, less robust algorithm for remapping.", &
2128  default=.false., do_not_log=just_read)
2129  call get_param(pf, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
2130  "This sets the default value for the various _2018_ANSWERS parameters.", &
2131  default=.false.)
2132  call get_param(pf, mdl, "TEMP_SALT_INIT_VERTICAL_REMAP_ONLY", pre_gridded, &
2133  "If true, initial conditions are on the model horizontal grid. " //&
2134  "Extrapolation over missing ocean values is done using an ICE-9 "//&
2135  "procedure with vertical ALE remapping .", &
2136  default=.false.)
2137  if (usealeremapping) then
2138  call get_param(pf, mdl, "REMAPPING_2018_ANSWERS", answers_2018, &
2139  "If true, use the order of arithmetic and expressions that recover the "//&
2140  "answers from the end of 2018. Otherwise, use updated and more robust "//&
2141  "forms of the same expressions.", default=default_2018_answers)
2142  endif
2143  call get_param(pf, mdl, "HOR_REGRID_2018_ANSWERS", hor_regrid_answers_2018, &
2144  "If true, use the order of arithmetic for horizonal regridding that recovers "//&
2145  "the answers from the end of 2018. Otherwise, use rotationally symmetric "//&
2146  "forms of the same expressions.", default=default_2018_answers)
2147  call get_param(pf, mdl, "ICE_SHELF", use_ice_shelf, default=.false.)
2148  if (use_ice_shelf) then
2149  call get_param(pf, mdl, "ICE_THICKNESS_FILE", ice_shelf_file, &
2150  "The file from which the ice bathymetry and area are read.", &
2151  fail_if_missing=.not.just_read, do_not_log=just_read)
2152  shelf_file = trim(inputdir)//trim(ice_shelf_file)
2153  if (.not.just_read) call log_param(pf, mdl, "INPUTDIR/THICKNESS_FILE", shelf_file)
2154  call get_param(pf, mdl, "ICE_AREA_VARNAME", area_varname, &
2155  "The name of the area variable in ICE_THICKNESS_FILE.", &
2156  fail_if_missing=.not.just_read, do_not_log=just_read)
2157  endif
2158  if (.not.usealeremapping) then
2159  call get_param(pf, mdl, "ADJUST_THICKNESS", correct_thickness, &
2160  "If true, all mass below the bottom removed if the "//&
2161  "topography is shallower than the thickness input file "//&
2162  "would indicate.", default=.false., do_not_log=just_read)
2163 
2164  call get_param(pf, mdl, "FIT_TO_TARGET_DENSITY_IC", adjust_temperature, &
2165  "If true, all the interior layers are adjusted to "//&
2166  "their target densities using mostly temperature "//&
2167  "This approach can be problematic, particularly in the "//&
2168  "high latitudes.", default=.true., do_not_log=just_read)
2169  call get_param(pf, mdl, "Z_INIT_SEPARATE_MIXED_LAYER", separate_mixed_layer, &
2170  "If true, distribute the topmost Z_INIT_HMIX_DEPTH of water over NKML layers, "//&
2171  "and do not correct the density of the topmost NKML+NKBL layers. Otherwise "//&
2172  "all layers are initialized based on the depths of their target densities.", &
2173  default=.false., do_not_log=just_read.or.(gv%nkml==0))
2174  if (gv%nkml == 0) separate_mixed_layer = .false.
2175  call get_param(pf, mdl, "MINIMUM_DEPTH", hmix_default, default=0.0)
2176  call get_param(pf, mdl, "Z_INIT_HMIX_DEPTH", hmix_depth, &
2177  "The mixed layer depth in the initial conditions when Z_INIT_SEPARATE_MIXED_LAYER "//&
2178  "is set to true.", default=hmix_default, units="m", scale=us%m_to_Z, &
2179  do_not_log=(just_read .or. .not.separate_mixed_layer))
2180  ! Reusing MINIMUM_DEPTH for the default mixed layer depth may be a strange choice, but
2181  ! it reproduces previous answers.
2182  endif
2183  if (just_read) then
2184  call cpu_clock_end(id_clock_routine)
2185  return ! All run-time parameters have been read, so return.
2186  endif
2187 
2188  eps_z = gv%Angstrom_Z
2189  eps_rho = 1.0e-10*us%kg_m3_to_R
2190 
2191  ! Read input grid coordinates for temperature and salinity field
2192  ! in z-coordinate dataset. The file is REQUIRED to contain the
2193  ! following:
2194  !
2195  ! dimension variables:
2196  ! lon (degrees_E), lat (degrees_N), depth(meters)
2197  ! variables:
2198  ! ptemp(lon,lat,depth) : degC, potential temperature
2199  ! salt (lon,lat,depth) : ppt, salinity
2200  !
2201  ! The first record will be read if there are multiple time levels.
2202  ! The observation grid MUST tile the model grid. If the model grid extends
2203  ! to the North/South Pole past the limits of the input data, they are extrapolated using the average
2204  ! value at the northernmost/southernmost latitude.
2205 
2206  call horiz_interp_and_extrap_tracer(tfilename, potemp_var, 1.0, 1, &
2207  g, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, &
2208  tripolar_n, homogenize, m_to_z=us%m_to_Z, answers_2018=hor_regrid_answers_2018, ongrid=pre_gridded)
2209 
2210  call horiz_interp_and_extrap_tracer(sfilename, salin_var, 1.0, 1, &
2211  g, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, &
2212  tripolar_n, homogenize, m_to_z=us%m_to_Z, answers_2018=hor_regrid_answers_2018, ongrid=pre_gridded)
2213 
2214  kd = size(z_in,1)
2215 
2216  ! Convert the sign convention of Z_edges_in.
2217  do k=1,size(z_edges_in,1) ; z_edges_in(k) = -z_edges_in(k) ; enddo
2218 
2219  allocate(rho_z(isd:ied,jsd:jed,kd))
2220  allocate(area_shelf_h(isd:ied,jsd:jed))
2221  allocate(frac_shelf_h(isd:ied,jsd:jed))
2222 
2223  ! Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 or NEMO
2224  call convert_temp_salt_for_teos10(temp_z, salt_z, g%HI, kd, mask_z, eos)
2225 
2226  press(:) = tv%P_Ref
2227  eosdom(:) = eos_domain(g%HI)
2228  do k=1,kd ; do j=js,je
2229  call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), eos, eosdom)
2230  enddo ; enddo
2231 
2232  call pass_var(temp_z,g%Domain)
2233  call pass_var(salt_z,g%Domain)
2234  call pass_var(mask_z,g%Domain)
2235  call pass_var(rho_z,g%Domain)
2236 
2237  ! This is needed for building an ALE grid under ice shelves
2238  if (use_ice_shelf) then
2239  if (.not.file_exists(shelf_file, g%Domain)) call mom_error(fatal, &
2240  "MOM_temp_salt_initialize_from_Z: Unable to open shelf file "//trim(shelf_file))
2241 
2242  call mom_read_data(shelf_file, trim(area_varname), area_shelf_h, g%Domain, scale=us%m_to_L**2)
2243 
2244  ! Initialize frac_shelf_h with zeros (open water everywhere)
2245  frac_shelf_h(:,:) = 0.0
2246  ! Compute fractional ice shelf coverage of h
2247  do j=jsd,jed ; do i=isd,ied
2248  if (g%areaT(i,j) > 0.0) &
2249  frac_shelf_h(i,j) = area_shelf_h(i,j) / g%areaT(i,j)
2250  enddo ; enddo
2251  ! Pass to the pointer for use as an argument to regridding_main
2252  shelf_area => frac_shelf_h
2253 
2254  endif
2255 
2256  ! Done with horizontal interpolation.
2257  ! Now remap to model coordinates
2258  if (usealeremapping) then
2259  call cpu_clock_begin(id_clock_ale)
2260  nkd = max(gv%ke, kd)
2261 
2262  ! Build the source grid and copy data onto model-shaped arrays with vanished layers
2263  allocate( tmp_mask_in(isd:ied,jsd:jed,nkd) ) ; tmp_mask_in(:,:,:) = 0.
2264  allocate( h1(isd:ied,jsd:jed,nkd) ) ; h1(:,:,:) = 0.
2265  allocate( tmpt1din(isd:ied,jsd:jed,nkd) ) ; tmpt1din(:,:,:) = 0.
2266  allocate( tmps1din(isd:ied,jsd:jed,nkd) ) ; tmps1din(:,:,:) = 0.
2267  do j = js, je ; do i = is, ie
2268  if (g%mask2dT(i,j)>0.) then
2269  ztopofcell = 0. ; zbottomofcell = 0.
2270  tmp_mask_in(i,j,1:kd) = mask_z(i,j,:)
2271  do k = 1, nkd
2272  if (tmp_mask_in(i,j,k)>0. .and. k<=kd) then
2273  zbottomofcell = max( z_edges_in(k+1), -g%bathyT(i,j) )
2274  tmpt1din(i,j,k) = temp_z(i,j,k)
2275  tmps1din(i,j,k) = salt_z(i,j,k)
2276  elseif (k>1) then
2277  zbottomofcell = -g%bathyT(i,j)
2278  tmpt1din(i,j,k) = tmpt1din(i,j,k-1)
2279  tmps1din(i,j,k) = tmps1din(i,j,k-1)
2280  else ! This next block should only ever be reached over land
2281  tmpt1din(i,j,k) = -99.9
2282  tmps1din(i,j,k) = -99.9
2283  endif
2284  h1(i,j,k) = gv%Z_to_H * (ztopofcell - zbottomofcell)
2285  ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
2286  enddo
2287  h1(i,j,kd) = h1(i,j,kd) + gv%Z_to_H * max(0., ztopofcell + g%bathyT(i,j) )
2288  ! The max here is in case the data data is shallower than model
2289  endif ! mask2dT
2290  enddo ; enddo
2291  deallocate( tmp_mask_in )
2292  call pass_var(h1, g%Domain)
2293  call pass_var(tmpt1din, g%Domain)
2294  call pass_var(tmps1din, g%Domain)
2295 
2296  ! Build the target grid (and set the model thickness to it)
2297  ! This call can be more general but is hard-coded for z* coordinates... ????
2298  call ale_initregridding( gv, us, g%max_depth, pf, mdl, regridcs ) ! sets regridCS
2299 
2300  if (.not. remap_general) then
2301  ! This is the old way of initializing to z* coordinates only
2302  allocate( htarget(nz) )
2303  htarget = getcoordinateresolution( regridcs )
2304  do j = js, je ; do i = is, ie
2305  h(i,j,:) = 0.
2306  if (g%mask2dT(i,j)>0.) then
2307  ! Build the target grid combining hTarget and topography
2308  ztopofcell = 0. ; zbottomofcell = 0.
2309  do k = 1, nz
2310  zbottomofcell = max( ztopofcell - htarget(k), -g%bathyT(i,j) )
2311  h(i,j,k) = gv%Z_to_H * (ztopofcell - zbottomofcell)
2312  ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
2313  enddo
2314  else
2315  h(i,j,:) = 0.
2316  endif ! mask2dT
2317  enddo ; enddo
2318  call pass_var(h, g%Domain)
2319  deallocate( htarget )
2320  endif
2321 
2322  ! Now remap from source grid to target grid, first setting reconstruction parameters
2323  call initialize_remapping( remapcs, remappingscheme, boundary_extrapolation=.false., answers_2018=answers_2018 )
2324  if (remap_general) then
2325  call set_regrid_params( regridcs, min_thickness=0. )
2326  tv_loc = tv
2327  tv_loc%T => tmpt1din
2328  tv_loc%S => tmps1din
2329  gv_loc = gv
2330  gv_loc%ke = nkd
2331  allocate( dz_interface(isd:ied,jsd:jed,nkd+1) ) ! Need for argument to regridding_main() but is not used
2332  if (use_ice_shelf) then
2333  call regridding_main( remapcs, regridcs, g, gv_loc, h1, tv_loc, h, dz_interface, shelf_area )
2334  else
2335  call regridding_main( remapcs, regridcs, g, gv_loc, h1, tv_loc, h, dz_interface )
2336  endif
2337  deallocate( dz_interface )
2338  endif
2339  call ale_remap_scalar(remapcs, g, gv, nkd, h1, tmpt1din, h, tv%T, all_cells=remap_full_column, &
2340  old_remap=remap_old_alg, answers_2018=answers_2018 )
2341  call ale_remap_scalar(remapcs, g, gv, nkd, h1, tmps1din, h, tv%S, all_cells=remap_full_column, &
2342  old_remap=remap_old_alg, answers_2018=answers_2018 )
2343  deallocate( h1 )
2344  deallocate( tmpt1din )
2345  deallocate( tmps1din )
2346 
2347  call cpu_clock_end(id_clock_ale)
2348 
2349  else ! remap to isopycnal layer space
2350 
2351  ! Next find interface positions using local arrays
2352  ! nlevs contains the number of valid data points in each column
2353  nlevs = int(sum(mask_z,dim=3))
2354 
2355  ! Rb contains the layer interface densities
2356  allocate(rb(nz+1))
2357  do k=2,nz ; rb(k) = 0.5*(gv%Rlay(k-1)+gv%Rlay(k)) ; enddo
2358  rb(1) = 0.0
2359  if (nz>1) then
2360  rb(nz+1) = 2.0*gv%Rlay(nz) - gv%Rlay(nz-1)
2361  else
2362  rb(nz+1) = 2.0 * gv%Rlay(1)
2363  endif
2364 
2365  nkml = 0 ; if (separate_mixed_layer) nkml = gv%nkml
2366 
2367  call find_interfaces(rho_z, z_in, kd, rb, g%bathyT, zi, g, us, &
2368  nlevs, nkml, hml=hmix_depth, eps_z=eps_z, eps_rho=eps_rho)
2369 
2370  if (correct_thickness) then
2371  call adjustetatofitbathymetry(g, gv, us, zi, h)
2372  else
2373  do k=nz,1,-1 ; do j=js,je ; do i=is,ie
2374  if (zi(i,j,k) < (zi(i,j,k+1) + gv%Angstrom_Z)) then
2375  zi(i,j,k) = zi(i,j,k+1) + gv%Angstrom_Z
2376  h(i,j,k) = gv%Angstrom_H
2377  else
2378  h(i,j,k) = gv%Z_to_H * (zi(i,j,k) - zi(i,j,k+1))
2379  endif
2380  enddo ; enddo ; enddo
2381  inconsistent=0
2382  do j=js,je ; do i=is,ie
2383  if (abs(zi(i,j,nz+1) + g%bathyT(i,j)) > 1.0*us%m_to_Z) &
2384  inconsistent = inconsistent + 1
2385  enddo ; enddo
2386  call sum_across_pes(inconsistent)
2387 
2388  if ((inconsistent > 0) .and. (is_root_pe())) then
2389  write(mesg, '("Thickness initial conditions are inconsistent ",'// &
2390  '"with topography in ",I5," places.")') inconsistent
2391  call mom_error(warning, mesg)
2392  endif
2393  endif
2394 
2395  call tracer_z_init_array(temp_z, z_edges_in, kd, zi, missing_value, g, nz, nlevs, eps_z, tv%T)
2396  call tracer_z_init_array(salt_z, z_edges_in, kd, zi, missing_value, g, nz, nlevs, eps_z, tv%S)
2397 
2398  do k=1,nz
2399  npoints = 0 ; tempavg = 0. ; saltavg = 0.
2400  do j=js,je ; do i=is,ie ; if (g%mask2dT(i,j) >= 1.0) then
2401  npoints = npoints + 1
2402  tempavg = tempavg + tv%T(i,j,k)
2403  saltavg = saltavg + tv%S(i,j,k)
2404  endif ; enddo ; enddo
2405 
2406  ! Horizontally homogenize data to produce perfectly "flat" initial conditions
2407  if (homogenize) then
2408  call sum_across_pes(npoints)
2409  call sum_across_pes(tempavg)
2410  call sum_across_pes(saltavg)
2411  if (npoints>0) then
2412  tempavg = tempavg / real(npoints)
2413  saltavg = saltavg / real(npoints)
2414  endif
2415  tv%T(:,:,k) = tempavg
2416  tv%S(:,:,k) = saltavg
2417  endif
2418  enddo
2419 
2420  endif ! useALEremapping
2421 
2422  ! Fill land values
2423  do k=1,nz ; do j=js,je ; do i=is,ie
2424  if (tv%T(i,j,k) == missing_value) then
2425  tv%T(i,j,k) = temp_land_fill
2426  tv%S(i,j,k) = salt_land_fill
2427  endif
2428  enddo ; enddo ; enddo
2429 
2430 
2431  if (adjust_temperature .and. .not. usealeremapping) then
2432  ! Finally adjust to target density
2433  ks = 1 ; if (separate_mixed_layer) ks = gv%nk_rho_varies + 1
2434  call determine_temperature(tv%T, tv%S, gv%Rlay(1:nz), tv%P_Ref, niter, &
2435  missing_value, h, ks, g, us, eos)
2436  endif
2437 
2438  deallocate(z_in, z_edges_in, temp_z, salt_z, mask_z)
2439  deallocate(rho_z, area_shelf_h, frac_shelf_h)
2440 
2441  call pass_var(h, g%Domain)
2442  call pass_var(tv%T, g%Domain)
2443  call pass_var(tv%S, g%Domain)
2444 
2445  call calltree_leave(trim(mdl)//'()')
2446  call cpu_clock_end(id_clock_routine)
2447 

◆ set_velocity_depth_max()

subroutine mom_state_initialization::set_velocity_depth_max ( type(ocean_grid_type), intent(inout)  G)
private

This subroutine sets the 4 bottom depths at velocity points to be the maximum of the adjacent depths.

Parameters
[in,out]gThe ocean's grid structure

Definition at line 1915 of file MOM_state_initialization.F90.

1916  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1917  ! Local variables
1918  integer :: i, j
1919 
1920  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
1921  g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
1922  g%Dopen_u(i,j) = g%Dblock_u(i,j)
1923  enddo ; enddo
1924  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
1925  g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
1926  g%Dopen_v(i,j) = g%Dblock_v(i,j)
1927  enddo ; enddo

◆ set_velocity_depth_min()

subroutine mom_state_initialization::set_velocity_depth_min ( type(ocean_grid_type), intent(inout)  G)
private

This subroutine sets the 4 bottom depths at velocity points to be the minimum of the adjacent depths.

Parameters
[in,out]gThe ocean's grid structure

Definition at line 1952 of file MOM_state_initialization.F90.

1953  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1954  ! Local variables
1955  integer :: i, j
1956 
1957  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
1958  g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
1959  g%Dopen_u(i,j) = g%Dblock_u(i,j)
1960  enddo ; enddo
1961  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
1962  g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
1963  g%Dopen_v(i,j) = g%Dblock_v(i,j)
1964  enddo ; enddo

◆ trim_for_ice()

subroutine mom_state_initialization::trim_for_ice ( type(param_file_type), intent(in)  PF,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(ale_cs), pointer  ALE_CSp,
type(thermo_var_ptrs), intent(inout)  tv,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  h,
logical, intent(in), optional  just_read_params 
)
private

Adjust the layer thicknesses by cutting away the top of each model column at the depth where the hydrostatic pressure matches an imposed surface pressure read from file.

Parameters
[in]pfParameter file structure
[in]gOcean grid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
ale_cspALE control structure
[in,out]tvThermodynamics structure
[in,out]hLayer thickness [H ~> m or kg m-2]
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 1081 of file MOM_state_initialization.F90.

1082  type(param_file_type), intent(in) :: PF !< Parameter file structure
1083  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1084  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1085  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1086  type(ALE_CS), pointer :: ALE_CSp !< ALE control structure
1087  type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure
1088  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1089  intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1090  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
1091  !! only read parameters without changing h.
1092  ! Local variables
1093  character(len=200) :: mdl = "trim_for_ice"
1094  real, dimension(SZI_(G),SZJ_(G)) :: p_surf ! Imposed pressure on ocean at surface [R L2 T-2 ~> Pa]
1095  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: S_t, S_b ! Top and bottom edge values for reconstructions
1096  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: T_t, T_b ! of salinity [ppt] and temperature [degC] within each layer.
1097  character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path
1098  real :: scale_factor ! A file-dependent scaling factor for the input pressure.
1099  real :: min_thickness ! The minimum layer thickness, recast into Z units [Z ~> m].
1100  integer :: i, j, k
1101  logical :: default_2018_answers, remap_answers_2018
1102  logical :: just_read ! If true, just read parameters but set nothing.
1103  logical :: use_remapping ! If true, remap the initial conditions.
1104  type(remapping_CS), pointer :: remap_CS => null()
1105 
1106  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
1107 
1108  call get_param(pf, mdl, "SURFACE_PRESSURE_FILE", p_surf_file, &
1109  "The initial condition file for the surface pressure exerted by ice.", &
1110  fail_if_missing=.not.just_read, do_not_log=just_read)
1111  call get_param(pf, mdl, "SURFACE_PRESSURE_VAR", p_surf_var, &
1112  "The initial condition variable for the surface pressure exerted by ice.", &
1113  units="Pa", default="", do_not_log=just_read)
1114  call get_param(pf, mdl, "INPUTDIR", inputdir, default=".", do_not_log=.true.)
1115  filename = trim(slasher(inputdir))//trim(p_surf_file)
1116  if (.not.just_read) call log_param(pf, mdl, "!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
1117 
1118  call get_param(pf, mdl, "SURFACE_PRESSURE_SCALE", scale_factor, &
1119  "A scaling factor to convert SURFACE_PRESSURE_VAR from "//&
1120  "file SURFACE_PRESSURE_FILE into a surface pressure.", &
1121  units="file dependent", default=1., do_not_log=just_read)
1122  call get_param(pf, mdl, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', &
1123  units='m', default=1.e-3, do_not_log=just_read, scale=us%m_to_Z)
1124  call get_param(pf, mdl, "TRIMMING_USES_REMAPPING", use_remapping, &
1125  'When trimming the column, also remap T and S.', &
1126  default=.false., do_not_log=just_read)
1127  remap_answers_2018 = .true.
1128  if (use_remapping) then
1129  call get_param(pf, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
1130  "This sets the default value for the various _2018_ANSWERS parameters.", &
1131  default=.false.)
1132  call get_param(pf, mdl, "REMAPPING_2018_ANSWERS", remap_answers_2018, &
1133  "If true, use the order of arithmetic and expressions that recover the "//&
1134  "answers from the end of 2018. Otherwise, use updated and more robust "//&
1135  "forms of the same expressions.", default=default_2018_answers)
1136  endif
1137 
1138  if (just_read) return ! All run-time parameters have been read, so return.
1139 
1140  call mom_read_data(filename, p_surf_var, p_surf, g%Domain, &
1141  scale=scale_factor*us%kg_m3_to_R*us%m_s_to_L_T**2)
1142 
1143  if (use_remapping) then
1144  allocate(remap_cs)
1145  call initialize_remapping(remap_cs, 'PLM', boundary_extrapolation=.true.)
1146  endif
1147 
1148  ! Find edge values of T and S used in reconstructions
1149  if ( associated(ale_csp) ) then ! This should only be associated if we are in ALE mode
1150  call ts_plm_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, .true.)
1151  else
1152 ! call MOM_error(FATAL, "trim_for_ice: Does not work without ALE mode")
1153  do k=1,g%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
1154  t_t(i,j,k) = tv%T(i,j,k) ; t_b(i,j,k) = tv%T(i,j,k)
1155  s_t(i,j,k) = tv%S(i,j,k) ; s_b(i,j,k) = tv%S(i,j,k)
1156  enddo ; enddo ; enddo
1157  endif
1158 
1159  do j=g%jsc,g%jec ; do i=g%isc,g%iec
1160  call cut_off_column_top(gv%ke, tv, gv, us, gv%g_Earth, g%bathyT(i,j), &
1161  min_thickness, tv%T(i,j,:), t_t(i,j,:), t_b(i,j,:), &
1162  tv%S(i,j,:), s_t(i,j,:), s_b(i,j,:), p_surf(i,j), h(i,j,:), remap_cs, &
1163  z_tol=1.0e-5*us%m_to_Z, remap_answers_2018=remap_answers_2018)
1164  enddo ; enddo
1165 
mom_coms::reproducing_sum
Find an accurate and order-invariant sum of a distributed 2d or 3d field.
Definition: MOM_coms.F90:53