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!) More...
 
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)
 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, 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 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. More...
 

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 719 of file MOM_state_initialization.F90.

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

Referenced by initialize_thickness_from_file(), and mom_temp_salt_initialize_from_z().

Here is the caller graph for this function:

◆ 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 1893 of file MOM_state_initialization.F90.

1893  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1894  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1895  ! Local variables
1896  real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
1897  real :: area_scale
1898  integer :: i,j
1899 
1900  area_scale = us%L_to_m**2
1901  tmpforsumming(:,:) = 0.
1902  g%areaT_global = 0.0 ; g%IareaT_global = 0.0
1903  do j=g%jsc,g%jec ; do i=g%isc,g%iec
1904  tmpforsumming(i,j) = area_scale*g%areaT(i,j) * g%mask2dT(i,j)
1905  enddo ; enddo
1906  g%areaT_global = reproducing_sum(tmpforsumming)
1907  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 929 of file MOM_state_initialization.F90.

929  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
930  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
931  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
932  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
933  intent(inout) :: h !< Input geometric layer thicknesses being converted
934  !! to layer pressure [H ~> m or kg m-2].
935  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
936  !! thermodynamic variables
937  ! Local variables
938  real, dimension(SZI_(G),SZJ_(G)) :: &
939  p_top, p_bot
940  real :: dz_geo(SZI_(G),SZJ_(G)) ! The change in geopotential height
941  ! across a layer [m2 s-2].
942  real :: rho(SZI_(G))
943  real :: I_gEarth
944  real :: Hm_rho_to_Pa ! A conversion factor from the input geometric thicknesses times the
945  ! layer densities into Pa [Pa m3 H-1 kg-1 ~> s-2 m2 or s-2 m5 kg-1].
946  logical :: Boussinesq
947  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
948  integer :: itt, max_itt
949 
950  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
951  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
952  max_itt = 10
953  boussinesq = gv%Boussinesq
954  i_gearth = 1.0 / (gv%mks_g_Earth)
955  hm_rho_to_pa = gv%mks_g_Earth * gv%H_to_m ! = GV%H_to_Pa / (US%R_to_kg_m3*GV%Rho0)
956 
957  if (boussinesq) then
958  call mom_error(fatal,"Not yet converting thickness with Boussinesq approx.")
959  else
960  if (associated(tv%eqn_of_state)) then
961  do j=jsq,jeq+1 ; do i=isq,ieq+1
962  p_bot(i,j) = 0.0 ; p_top(i,j) = 0.0
963  enddo ; enddo
964  do k=1,nz
965  do j=js,je
966  do i=is,ie ; p_top(i,j) = p_bot(i,j) ; enddo
967  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_top(:,j), rho, &
968  is, ie-is+1, tv%eqn_of_state)
969  do i=is,ie
970  p_bot(i,j) = p_top(i,j) + hm_rho_to_pa * (h(i,j,k) * rho(i))
971  enddo
972  enddo
973 
974  do itt=1,max_itt
975  call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p_top, p_bot, &
976  0.0, g%HI, tv%eqn_of_state, dz_geo)
977  if (itt < max_itt) then ; do j=js,je
978  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_bot(:,j), rho, &
979  is, ie-is+1, tv%eqn_of_state)
980  ! Use Newton's method to correct the bottom value.
981  ! The hydrostatic equation is linear to such a
982  ! high degree that no bounds-checking is needed.
983  do i=is,ie
984  p_bot(i,j) = p_bot(i,j) + rho(i) * &
985  (hm_rho_to_pa*h(i,j,k) - dz_geo(i,j))
986  enddo
987  enddo ; endif
988  enddo
989 
990  do j=js,je ; do i=is,ie
991  h(i,j,k) = (p_bot(i,j) - p_top(i,j)) * gv%kg_m2_to_H * i_gearth
992  enddo ; enddo
993  enddo
994  else
995  do k=1,nz ; do j=js,je ; do i=is,ie
996  h(i,j,k) = h(i,j,k) * (gv%Rlay(k) / gv%Rho0)
997  enddo ; enddo ; enddo
998  endif
999  endif
1000 

References mom_eos::int_specific_vol_dp().

Referenced by mom_initialize_state().

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

◆ 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 
)
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 [m2 Z-1 s-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 [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].

Definition at line 1168 of file MOM_state_initialization.F90.

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

References mom_remapping::remapping_core_h().

Referenced by mom_state_init_tests(), and trim_for_ice().

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

◆ 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 1005 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

Here is the caller graph for this function:

◆ 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  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.
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 1694 of file MOM_state_initialization.F90.

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

References mom_error_handler::mom_error().

Referenced by mom_initialize_state().

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

◆ 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_stateInteger that selects the equatio of state.
[in]p_refThe coordinate-density reference pressure [Pa].
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 1533 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

Here is the caller graph for this function:

◆ 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 1432 of file MOM_state_initialization.F90.

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

References mom_error_handler::mom_error().

Referenced by mom_initialize_state().

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

◆ 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 1488 of file MOM_state_initialization.F90.

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

References mom_error_handler::mom_error().

Referenced by mom_initialize_state().

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

◆ 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 1627 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

Here is the caller graph for this function:

◆ 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 627 of file MOM_state_initialization.F90.

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

References adjustetatofitbathymetry(), mom_error_handler::calltree_enter(), and mom_error_handler::calltree_leave().

Referenced by mom_initialize_state().

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

◆ 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 848 of file MOM_state_initialization.F90.

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

References mom_error_handler::calltree_enter(), and mom_error_handler::calltree_leave().

Referenced by mom_initialize_state().

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

◆ initialize_thickness_search()

subroutine mom_state_initialization::initialize_thickness_search ( )
private

Search density space for location of layers (not implemented!)

Definition at line 924 of file MOM_state_initialization.F90.

924  call mom_error(fatal," MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED")

Referenced by mom_initialize_state().

Here is the caller graph for this function:

◆ 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 793 of file MOM_state_initialization.F90.

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

References mom_error_handler::calltree_enter(), and mom_error_handler::calltree_leave().

Referenced by mom_initialize_state().

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

◆ 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 1370 of file MOM_state_initialization.F90.

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

References my_psi().

Referenced by mom_initialize_state().

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

◆ 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 1254 of file MOM_state_initialization.F90.

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

References mom_error_handler::calltree_enter(), and mom_error_handler::calltree_leave().

Referenced by mom_initialize_state().

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

◆ 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 1329 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

Here is the caller graph for this function:

◆ 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 1295 of file MOM_state_initialization.F90.

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

References mom_error_handler::calltree_enter(), and mom_error_handler::calltree_leave().

Referenced by mom_initialize_state().

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

◆ 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 125 of file MOM_state_initialization.F90.

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

References adjustment_initialization::adjustment_initialize_temperature_salinity(), mom_ale::ale_regrid_accelerated(), baroclinic_zone_initialization::baroclinic_zone_init_temperature_salinity(), benchmark_initialization::benchmark_init_temperature_salinity(), bfb_initialization::bfb_initialize_sponges_southonly(), mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), circle_obcs_initialization::circle_obcs_initialize_thickness(), convert_thickness(), dense_water_initialization::dense_water_initialize_sponges(), depress_surface(), mom_restart::determine_is_new_run(), dome2d_initialization::dome2d_initialize_sponges(), dome_initialization::dome_initialize_sponges(), dumbbell_initialization::dumbbell_initialize_sponges(), dyed_channel_initialization::dyed_channel_set_obc_tracer_data(), dyed_obcs_initialization::dyed_obcs_set_obc_data(), external_gwave_initialization::external_gwave_initialize_thickness(), initialize_sponges_file(), initialize_temp_salt_fit(), initialize_temp_salt_from_file(), initialize_temp_salt_from_profile(), initialize_temp_salt_linear(), initialize_thickness_from_file(), initialize_thickness_list(), initialize_thickness_search(), initialize_thickness_uniform(), initialize_velocity_circular(), initialize_velocity_from_file(), initialize_velocity_uniform(), initialize_velocity_zero(), isomip_initialization::isomip_initialize_temperature_salinity(), lock_exchange_initialization::lock_exchange_initialize_thickness(), mom_string_functions::lowercase(), mdl, mom_temp_salt_initialize_from_z(), neverland_initialization::neverland_initialize_thickness(), phillips_initialization::phillips_initialize_sponges(), mom_restart::restore_state(), rgc_initialization::rgc_initialize_sponges(), rossby_front_2d_initialization::rossby_front_initialize_velocity(), scm_cvmix_tests::scm_cvmix_tests_ts_init(), seamount_initialization::seamount_initialize_temperature_salinity(), sloshing_initialization::sloshing_initialize_temperature_salinity(), soliton_initialization::soliton_initialize_thickness(), supercritical_initialization::supercritical_set_obc_data(), trim_for_ice(), mom_open_boundary::update_obc_segment_data(), and user_initialization::user_initialize_sponges().

Referenced by mom::initialize_mom().

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

◆ 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 2388 of file MOM_state_initialization.F90.

2388  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2389  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2390  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2391  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
2392  ! Local variables
2393  integer, parameter :: nk=5
2394  real, dimension(nk) :: T, T_t, T_b, S, S_t, S_b, rho, h, z
2395  real, dimension(nk+1) :: e
2396  integer :: k
2397  real :: P_tot, P_t, P_b, z_out
2398  type(remapping_CS), pointer :: remap_CS => null()
2399 
2400  do k = 1, nk
2401  h(k) = 100.
2402  enddo
2403  e(1) = 0.
2404  do k = 1, nk
2405  e(k+1) = e(k) - h(k)
2406  enddo
2407  p_tot = 0.
2408  do k = 1, nk
2409  z(k) = 0.5 * ( e(k) + e(k+1) )
2410  t_t(k) = 20.+(0./500.)*e(k)
2411  t(k) = 20.+(0./500.)*z(k)
2412  t_b(k) = 20.+(0./500.)*e(k+1)
2413  s_t(k) = 35.-(0./500.)*e(k)
2414  s(k) = 35.+(0./500.)*z(k)
2415  s_b(k) = 35.-(0./500.)*e(k+1)
2416  call calculate_density(0.5*(t_t(k)+t_b(k)), 0.5*(s_t(k)+s_b(k)), -us%R_to_kg_m3*gv%Rho0*gv%mks_g_Earth*z(k), &
2417  rho(k), tv%eqn_of_state)
2418  p_tot = p_tot + gv%mks_g_Earth * rho(k) * h(k)
2419  enddo
2420 
2421  p_t = 0.
2422  do k = 1, nk
2423  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, &
2424  us%R_to_kg_m3*gv%Rho0, gv%mks_g_Earth, tv%eqn_of_state, p_b, z_out)
2425  write(0,*) k,p_t,p_b,0.5*p_tot,e(k),e(k+1),z_out
2426  p_t = p_b
2427  enddo
2428  write(0,*) p_b,p_tot
2429 
2430  write(0,*) ''
2431  write(0,*) ' ==================================================================== '
2432  write(0,*) ''
2433  write(0,*) h
2434  call cut_off_column_top(nk, tv, gv, us, gv%mks_g_Earth, -e(nk+1), gv%Angstrom_H, &
2435  t, t_t, t_b, s, s_t, s_b, 0.5*p_tot, h, remap_cs)
2436  write(0,*) h
2437 

References cut_off_column_top().

Here is the call graph for this function:

◆ 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 1931 of file MOM_state_initialization.F90.

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

References adjustetatofitbathymetry(), and mom_error_handler::mom_error().

Referenced by mom_initialize_state().

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

◆ 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 1876 of file MOM_state_initialization.F90.

1876  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1877  ! Local variables
1878  integer :: i, j
1879 
1880  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
1881  g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
1882  g%Dopen_u(i,j) = g%Dblock_u(i,j)
1883  enddo ; enddo
1884  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
1885  g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
1886  g%Dopen_v(i,j) = g%Dblock_v(i,j)
1887  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 1913 of file MOM_state_initialization.F90.

1913  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1914  ! Local variables
1915  integer :: i, j
1916 
1917  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
1918  g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
1919  g%Dopen_u(i,j) = g%Dblock_u(i,j)
1920  enddo ; enddo
1921  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
1922  g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
1923  g%Dopen_v(i,j) = g%Dblock_v(i,j)
1924  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 1089 of file MOM_state_initialization.F90.

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

References cut_off_column_top().

Referenced by mom_initialize_state().

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

Variable Documentation

◆ mdl

character(len=40) mom_state_initialization::mdl = "MOM_state_initialization"
private

This module's name.

Definition at line 116 of file MOM_state_initialization.F90.

116 character(len=40) :: mdl = "MOM_state_initialization" !< This module's name.

Referenced by mom_initialize_state().

my_psi
real function my_psi(ig, jg)
Returns the value of a circular stream function at (ig,jg)
Definition: MOM_state_initialization.F90:1415
mom_coms::reproducing_sum
Find an accurate and order-invariant sum of distributed 2d or 3d fields.
Definition: MOM_coms.F90:54