mom_state_initialization Module Reference

Detailed Description

Initialize state variables, u, v, h, T and S.

Functions/Subroutines
subroutine, public	mom_initialize_state (u, v, h, tv, Time, G, GV, PF, dirs, restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, ALE_sponge_CSp, OBC, Time_in)
subroutine	initialize_thickness_from_file (h, G, GV, param_file, file_has_thickness, just_read_params)
	This subroutine reads the layer thicknesses or interface heights from a file. More...
subroutine	adjustetatofitbathymetry (G, GV, eta, h)
	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_z. If the bottom most interface is above the topography then the entire column is dilated (expanded) to fill the void. More...
subroutine	initialize_thickness_uniform (h, G, GV, param_file, just_read_params)
subroutine	initialize_thickness_search
subroutine	convert_thickness (h, G, GV, tv)
subroutine	depress_surface (h, G, GV, param_file, tv, just_read_params)
subroutine	trim_for_ice (PF, G, GV, 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, Rho0, G_earth, depth, min_thickness, T, T_t, T_b, S, S_t, S_b, p_surf, h, remap_CS)
	Adjust the layer thicknesses by cutting away the top at the depth where the hydrostatic pressure matches p_surf. More...
subroutine	initialize_velocity_from_file (u, v, G, param_file, just_read_params)
subroutine	initialize_velocity_zero (u, v, G, param_file, just_read_params)
subroutine	initialize_velocity_uniform (u, v, G, param_file, just_read_params)
subroutine	initialize_velocity_circular (u, v, G, param_file, just_read_params)
subroutine	initialize_temp_salt_from_file (T, S, G, param_file, just_read_params)
subroutine	initialize_temp_salt_from_profile (T, S, G, param_file, just_read_params)
subroutine	initialize_temp_salt_fit (T, S, G, GV, param_file, eqn_of_state, P_Ref, just_read_params)
subroutine	initialize_temp_salt_linear (T, S, G, param_file, just_read_params)
subroutine	initialize_sponges_file (G, GV, use_temperature, tv, param_file, CSp)
	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)
subroutine	compute_global_grid_integrals (G)
subroutine	set_velocity_depth_min (G)
subroutine	mom_temp_salt_initialize_from_z (h, tv, G, GV, 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, tv)
	Run simple unit tests. More...

Variables
character(len=40)	mdl = "MOM_state_initialization"

Function/Subroutine Documentation

adjustetatofitbathymetry()

subroutine mom_state_initialization::adjustetatofitbathymetry ( type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, 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_z. 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]	G	Grid type
[in,out]	eta	Interface heights
[out]	h	Layer thicknesses
[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in,out]	h	Layer thicknesses, in m

Definition at line 684 of file MOM_state_initialization.F90.

Referenced by initialize_thickness_from_file(), and mom_temp_salt_initialize_from_z().

  type(ocean_grid_type),                          intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type),                        intent(in)    :: gv   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)+1),    intent(inout) :: eta
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)),      intent(inout) :: h    !< Layer thicknesses, in m
  ! Local variables
  integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
  real, parameter :: htolerance = 0.1 !<  Tolerance to exceed adjustment criteria (m)
  real :: htmp, etmp, dilate
  character(len=100) :: mesg

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  contractions = 0
  do j=js,je ; do i=is,ie
    if (-eta(i,j,nz+1) > g%bathyT(i,j) + htolerance) then
      eta(i,j,nz+1) = -g%bathyT(i,j)
      contractions = contractions + 1
    endif
  enddo ; enddo
  call sum_across_pes(contractions)
  if ((contractions > 0) .and. (is_root_pe())) then
    write(mesg,'("Thickness initial conditions were contracted ",'// &
               '"to fit topography in ",I8," places.")') contractions
    call mom_error(warning, 'adjustEtaToFitBathymetry: '//mesg)
  endif

  do k=nz,1,-1 ; do j=js,je ; do i=is,ie
    ! Collapse layers to thinnest possible if the thickness less than
    ! the thinnest possible (or negative).
    if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_z)) then
      eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_z
      h(i,j,k) = gv%Angstrom_z
    else
      h(i,j,k) = eta(i,j,k) - eta(i,j,k+1)
    endif
  enddo ; enddo ; enddo

  dilations = 0
  do j=js,je ; do i=is,ie
    !   The whole column is dilated to accommodate deeper topography than
    ! the bathymetry would indicate.
    ! This should be...  if ((G%mask2dt(i,j)*(eta(i,j,1)-eta(i,j,nz+1)) > 0.0) .and. &
    if (-eta(i,j,nz+1) < g%bathyT(i,j) - htolerance) then
      dilations = dilations + 1
      if (eta(i,j,1) <= eta(i,j,nz+1)) then
        do k=1,nz ; h(i,j,k) = (eta(i,j,1)+g%bathyT(i,j)) / real(nz) ; enddo
      else
        dilate = (eta(i,j,1)+g%bathyT(i,j)) / (eta(i,j,1)-eta(i,j,nz+1))
        do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
      endif
      do k=nz, 2, -1; eta(i,j,k) = eta(i,j,k+1) + h(i,j,k); enddo
    endif
  enddo ; enddo
  call sum_across_pes(dilations)
  if ((dilations > 0) .and. (is_root_pe())) then
    write(mesg,'("Thickness initial conditions were dilated ",'// &
               '"to fit topography in ",I8," places.")') dilations
    call mom_error(warning, 'adjustEtaToFitBathymetry: '//mesg)
  endif

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 )

private

Parameters

[in,out]

g

The ocean's grid structure

Definition at line 1777 of file MOM_state_initialization.F90.

  type(ocean_grid_type), intent(inout) :: g    !< The ocean's grid structure
  ! Subroutine to pre-compute global integrals of grid quantities for
  ! later use in reporting diagnostics
  real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpforsumming
  integer :: i,j

  tmpforsumming(:,:) = 0.
  g%areaT_global = 0.0 ; g%IareaT_global = 0.0
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    tmpforsumming(i,j) = g%areaT(i,j) * g%mask2dT(i,j)
  enddo ; enddo
  g%areaT_global = reproducing_sum(tmpforsumming)
  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(thermo_var_ptrs), intent(in)  tv  )

private

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in,out]	h	Layer thicknesses, being
[in]	tv	A structure pointing to various thermodynamic variables

Definition at line 818 of file MOM_state_initialization.F90.

References mom_eos::int_specific_vol_dp().

Referenced by mom_initialize_state().

  type(ocean_grid_type),                  intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type),                intent(in)    :: gv   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                                          intent(inout) :: h    !< Layer thicknesses, being
                                                                !! converted from m to H (m or kg
                                                                !! m-2)
  type(thermo_var_ptrs),                  intent(in)    :: tv   !< A structure pointing to various
                                                                !! thermodynamic variables
! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    p_top, p_bot
  real :: dz_geo(szi_(g),szj_(g))      ! The change in geopotential height
                                       ! across a layer, in m2 s-2.
  real :: rho(szi_(g))
  real :: i_gearth
  logical :: boussinesq
  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
  integer :: itt, max_itt

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  max_itt = 10
  boussinesq = gv%Boussinesq
  i_gearth = 1.0 / gv%g_Earth

  if (boussinesq) then
    call mom_error(fatal,"Not yet converting thickness with Boussinesq approx.")
  else
    if (associated(tv%eqn_of_state)) then
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        p_bot(i,j) = 0.0 ; p_top(i,j) = 0.0
      enddo ; enddo
      do k=1,nz
        do j=js,je
          do i=is,ie ; p_top(i,j) = p_bot(i,j) ; enddo
          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_top(:,j), rho, &
                                 is, ie-is+1, tv%eqn_of_state)
          do i=is,ie
            p_bot(i,j) = p_top(i,j) + gv%g_Earth * h(i,j,k) * rho(i)
          enddo
        enddo

        do itt=1,max_itt
          call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p_top, p_bot, &
                                   0.0, g%HI, tv%eqn_of_state, dz_geo)
          if (itt < max_itt) then ; do j=js,je
            call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_bot(:,j), rho, &
                                   is, ie-is+1, tv%eqn_of_state)
            ! Use Newton's method to correct the bottom value.
            !   The hydrostatic equation is linear to such a
            ! high degree that no bounds-checking is needed.
            do i=is,ie
              p_bot(i,j) = p_bot(i,j) + rho(i) * (gv%g_Earth*h(i,j,k) - dz_geo(i,j))
            enddo
          enddo ; endif
        enddo

        do j=js,je ; do i=is,ie
          h(i,j,k) = (p_bot(i,j) - p_top(i,j)) * gv%kg_m2_to_H * i_gearth
        enddo ; enddo
      enddo
    else
      do k=1,nz ; do j=js,je ; do i=is,ie
        h(i,j,k) = h(i,j,k) * gv%Rlay(k) * gv%kg_m2_to_H
      enddo ; enddo ; enddo
    endif
  endif

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, real, intent(in)  Rho0, 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  )

private

Adjust the layer thicknesses by cutting away the top at the depth where the hydrostatic pressure matches p_surf.

Parameters

[in]	nk	Number of layers
[in]	tv	Thermodynamics structure
[in]	rho0	Reference density (kg/m3)
[in]	g_earth	Gravitational acceleration (m/s2)
[in]	depth	Depth of ocean column (m)
[in]	min_thickness	Smallest thickness allowed (m)
[in,out]	t	Layer mean temperature
[in]	t_t	Temperature at top of layer
[in]	t_b	Temperature at bottom of layer
[in,out]	s	Layer mean salinity
[in]	s_t	Salinity at top of layer
[in]	s_b	Salinity at bottom of layer
[in]	p_surf	Imposed pressure on ocean at surface (Pa)
[in,out]	h	Layer thickness (H units, m or Pa)

Definition at line 1064 of file MOM_state_initialization.F90.

References mom_remapping::remapping_core_h().

Referenced by mom_state_init_tests(), and trim_for_ice().

  integer,               intent(in)    :: nk !< Number of layers
  type(thermo_var_ptrs), intent(in)    :: tv !< Thermodynamics structure
  real,                  intent(in)    :: rho0 !< Reference density (kg/m3)
  real,                  intent(in)    :: g_earth !< Gravitational acceleration (m/s2)
  real,                  intent(in)    :: depth !< Depth of ocean column (m)
  real,                  intent(in)    :: min_thickness !< Smallest thickness allowed (m)
  real, dimension(nk),   intent(inout) :: t !< Layer mean temperature
  real, dimension(nk),   intent(in)    :: t_t !< Temperature at top of layer
  real, dimension(nk),   intent(in)    :: t_b !< Temperature at bottom of layer
  real, dimension(nk),   intent(inout) :: s !< Layer mean salinity
  real, dimension(nk),   intent(in)    :: s_t !< Salinity at top of layer
  real, dimension(nk),   intent(in)    :: s_b !< Salinity at bottom of layer
  real,                  intent(in)    :: p_surf !< Imposed pressure on ocean at surface (Pa)
  real, dimension(nk),   intent(inout) :: h !< Layer thickness (H units, m or Pa)
  type(remapping_cs),    pointer       :: remap_cs ! Remapping structure for remapping T and S, if associated
  ! Local variables
  real, dimension(nk+1) :: e ! Top and bottom edge values for reconstructions
  real, dimension(nk) :: h0, s0, t0, h1, s1, t1
  real :: p_t, p_b, z_out, e_top
  integer :: k

  ! Calculate original interface positions
  e(nk+1) = -depth
  do k=nk,1,-1
    e(k) = e(k+1) + h(k)
    h0(k) = h(nk+1-k) ! Keep a copy to use in remapping
  enddo

  p_t = 0.
  e_top = e(1)
  do k=1,nk
    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, p_surf, rho0, g_earth, tv%eqn_of_state, p_b, z_out)
    if (z_out>=e(k)) then
      ! Imposed pressure was less that pressure at top of cell
      exit
    elseif (z_out<=e(k+1)) then
      ! Imposed pressure was greater than pressure at bottom of cell
      e_top = e(k+1)
    else
      ! Imposed pressure was fell between pressures at top and bottom of cell
      e_top = z_out
      exit
    endif
    p_t = p_b
  enddo
  if (e_top<e(1)) then
    ! Clip layers from the top down, if at all
    do k=1,nk
      if (e(k)>e_top) then
        ! Original e(K) is too high
        e(k) = e_top
        e_top = e_top - min_thickness ! Next interface must be at least this deep
      endif
      ! This layer needs trimming
      h(k) = max( min_thickness, e(k) - e(k+1) )
      if (e(k)<e_top) exit ! No need to go further
    enddo
  endif

  ! Now we need to remap but remapping assumes the surface is at the
  ! same place in the two columns so we turn the column upside down.
  if (associated(remap_cs)) then
    do k=1,nk
      s0(k) = s(nk+1-k)
      t0(k) = t(nk+1-k)
      h1(k) = h(nk+1-k)
    enddo
    call remapping_core_h(remap_cs, nk, h0, t0, nk, h1, t1)
    call remapping_core_h(remap_cs, nk, h0, s0, nk, h1, s1)
    do k=1,nk
      s(k) = s1(nk+1-k)
      t(k) = t1(nk+1-k)
    enddo
  endif

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(param_file_type), intent(in)  param_file, type(thermo_var_ptrs), intent(in)  tv, logical, intent(in), optional  just_read_params  )

private

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in,out]	h	Layer thicknesses, in H (usually m or kg m-2)
[in]	param_file	A structure to parse for run-time parameters
[in]	tv	A structure pointing to various thermodynamic variables
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 891 of file MOM_state_initialization.F90.

Referenced by mom_initialize_state().

  type(ocean_grid_type),   intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type), intent(in)    :: gv   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                           intent(inout) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure pointing to various thermodynamic variables
  logical,       optional, intent(in)    :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

  real, dimension(SZI_(G),SZJ_(G)) :: &
    eta_sfc  ! The free surface height that the model should use, in m.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
    eta  ! The free surface height that the model should use, in m.
  real :: dilate  ! A ratio by which layers are dilated, nondim.
  real :: scale_factor ! A scaling factor for the eta_sfc values that are read
                       ! in, which can be used to change units, for example.
  character(len=40)  :: mdl = "depress_surface" ! This subroutine's name.
  character(len=200) :: inputdir, eta_srf_file ! Strings for file/path
  character(len=200) :: filename, eta_srf_var  ! Strings for file/path
  logical :: just_read    ! If true, just read parameters but set nothing.
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  ! Read the surface height (or pressure) from a file.

  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_FILE", eta_srf_file,&
                 "The initial condition file for the surface height.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, &
                 "The initial condition variable for the surface height.",&
                 default="SSH", do_not_log=just_read)
  filename = trim(inputdir)//trim(eta_srf_file)
  if (.not.just_read) &
    call log_param(param_file,  mdl, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)

  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_SCALE", scale_factor, &
                 "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into \n"//&
                 "units of m", units="variable", default=1.0, do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  call read_data(filename,eta_srf_var,eta_sfc,domain=g%Domain%mpp_domain)

  if (scale_factor /= 1.0) then ; do j=js,je ; do i=is,ie
    eta_sfc(i,j) = eta_sfc(i,j) * scale_factor
  enddo ; enddo ; endif

  ! Convert thicknesses to interface heights.
  call find_eta(h, tv, gv%g_Earth, g, gv, eta)

  do j=js,je ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
!    if (eta_sfc(i,j) < eta(i,j,nz+1)) then
      ! Issue a warning?
!    endif
    if (eta_sfc(i,j) > eta(i,j,1)) then
      ! Dilate the water column to agree, but only up to 10-fold.
      if (eta_sfc(i,j) - eta(i,j,nz+1) > 10.0*(eta(i,j,1) - eta(i,j,nz+1))) then
        dilate = 10.0
        call mom_error(warning, "Free surface height dilation attempted "//&
               "to exceed 10-fold.", all_print=.true.)
      else
        dilate = (eta_sfc(i,j) - eta(i,j,nz+1)) / (eta(i,j,1) - eta(i,j,nz+1))
      endif
      do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
    elseif (eta(i,j,1) > eta_sfc(i,j)) then
      ! Remove any mass that is above the target free surface.
      do k=1,nz
        if (eta(i,j,k) <= eta_sfc(i,j)) exit
        if (eta(i,j,k+1) >= eta_sfc(i,j)) then
          h(i,j,k) = gv%Angstrom
        else
          h(i,j,k) = max(gv%Angstrom, h(i,j,k) * &
              (eta_sfc(i,j) - eta(i,j,k+1)) / (eta(i,j,k) - eta(i,j,k+1)) )
        endif
      enddo
    endif
  endif ; enddo ; enddo

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, 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  )

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]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	tv	A structure pointing to various thermodynamic variables.
[in]	param_file	A structure to parse for run-time parameters.

Definition at line 1621 of file MOM_state_initialization.F90.

References mom_sponge::initialize_sponge().

Referenced by mom_initialize_state().

  type(ocean_grid_type),   intent(in) :: g    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in) :: gv   !< The ocean's vertical grid structure.
  logical,                 intent(in) :: use_temperature
  type(thermo_var_ptrs),   intent(in) :: tv   !< A structure pointing to various thermodynamic
                                              !! variables.
  type(param_file_type),   intent(in) :: param_file !< A structure to parse for run-time parameters.
  type(sponge_cs),         pointer    :: csp  !! A pointer that is set to point to the control
                                              !! structure for this module.
!   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.                                        !

! Arguments: from_file - .true. if the variables that are used here are to
!                        be read from a file; .false. to be set internally.
!  (in)      filename - The name of the file to read for all fields
!                       except the inverse damping rate.
!  (in)      damp_file - The name of the file from which to read the
!                        inverse damping rate.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      use_temperature - If true, T & S are state variables.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in/out)  CSp - A pointer that is set to point to the control structure
!                  for this module

  real :: eta(szi_(g),szj_(g),szk_(g)+1) ! The target interface heights, in m.
  real, dimension (SZI_(G),SZJ_(G),SZK_(G)) :: &
    tmp, tmp2 ! A temporary array for tracers.
  real, dimension (SZI_(G),SZJ_(G)) :: &
    tmp_2d ! A temporary array for tracers.

  real :: idamp(szi_(g),szj_(g))    ! The inverse damping rate, in s-1.
  real :: pres(szi_(g))     ! An array of the reference pressure, in Pa.

  integer :: i, j, k, is, ie, js, je, nz
  character(len=40) :: potemp_var, salin_var, idamp_var, eta_var
  character(len=40) :: mdl = "initialize_sponges_file"
  character(len=200) :: damping_file, state_file  ! Strings for filenames
  character(len=200) :: filename, inputdir ! Strings for file/path and path.
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  pres(:) = 0.0 ; eta(:,:,:) = 0.0 ; tmp(:,:,:) = 0.0 ; idamp(:,:) = 0.0

  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  call get_param(param_file, mdl, "SPONGE_DAMPING_FILE", damping_file, &
                 "The name of the file with the sponge damping rates.", &
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "SPONGE_STATE_FILE", state_file, &
                 "The name of the file with the state to damp toward.", &
                 default=damping_file)

  call get_param(param_file, mdl, "SPONGE_PTEMP_VAR", potemp_var, &
                 "The name of the potential temperature variable in \n"//&
                 "SPONGE_STATE_FILE.", default="PTEMP")
  call get_param(param_file, mdl, "SPONGE_SALT_VAR", salin_var, &
                 "The name of the salinity variable in \n"//&
                 "SPONGE_STATE_FILE.", default="SALT")
  call get_param(param_file, mdl, "SPONGE_ETA_VAR", eta_var, &
                 "The name of the interface height variable in \n"//&
                 "SPONGE_STATE_FILE.", default="ETA")
  call get_param(param_file, mdl, "SPONGE_IDAMP_VAR", idamp_var, &
                 "The name of the inverse damping rate variable in \n"//&
                 "SPONGE_DAMPING_FILE.", default="IDAMP")

  filename = trim(inputdir)//trim(damping_file)
  call log_param(param_file, mdl, "INPUTDIR/SPONGE_DAMPING_FILE", filename)
  if (.not.file_exists(filename, g%Domain)) &
    call mom_error(fatal, " initialize_sponges: Unable to open "//trim(filename))


  call read_data(filename,"Idamp",idamp(:,:), domain=g%Domain%mpp_domain)

! Now register all of the fields which are damped in the sponge.     !
! By default, momentum is advected vertically within the sponge, but !
! momentum is typically not damped within the sponge.                !

  filename = trim(inputdir)//trim(state_file)
  call log_param(param_file, mdl, "INPUTDIR/SPONGE_STATE_FILE", filename)
  if (.not.file_exists(filename, g%Domain)) &
    call mom_error(fatal, " initialize_sponges: Unable to open "//trim(filename))


!  The first call to set_up_sponge_field is for the interface height.!
  call read_data(filename, eta_var, eta(:,:,:), domain=g%Domain%mpp_domain)

  do j=js,je ; do i=is,ie
    eta(i,j,nz+1) = -g%bathyT(i,j)
  enddo ; enddo
  do k=nz,1,-1 ; do j=js,je ; do i=is,ie
    if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_z)) &
      eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_z
  enddo ; enddo ; enddo
! Set the inverse damping rates so that the model will know where to !
! apply the sponges, along with the interface heights.               !
  call initialize_sponge(idamp, eta, g, param_file, csp)

!   Now register all of the tracer fields which are damped in the    !
! sponge. By default, momentum is advected vertically within the     !
! sponge, but momentum is typically not damped within the sponge.    !

  if ( gv%nkml>0 ) then
!   This call to set_up_sponge_ML_density registers the target values of the
! mixed layer density, which is used in determining which layers can be
! inflated without causing static instabilities.
    do i=is-1,ie ; pres(i) = tv%P_Ref ; enddo

    call read_data(filename, potemp_var, tmp(:,:,:), domain=g%Domain%mpp_domain)
    call read_data(filename, salin_var, tmp2(:,:,:), domain=g%Domain%mpp_domain)

    do j=js,je
      call calculate_density(tmp(:,j,1), tmp2(:,j,1), pres, tmp_2d(:,j), &
                             is, ie-is+1, tv%eqn_of_state)
    enddo

    call set_up_sponge_ml_density(tmp_2d, g, csp)
  endif

!  The remaining calls to set_up_sponge_field can be in any order.   !
  if ( use_temperature ) then
    call read_data(filename, potemp_var, tmp(:,:,:), domain=g%Domain%mpp_domain)
    call set_up_sponge_field(tmp, tv%T, g, nz, csp)
    call read_data(filename, salin_var, tmp(:,:,:), domain=g%Domain%mpp_domain)
    call set_up_sponge_field(tmp, tv%S, g, nz, csp)
  endif

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(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(verticalgrid_type), intent(in)  GV, 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

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[out]	t	The potential temperature that is being
[out]	s	The salinity that is being initialized.
[in]	param_file	A structure to parse for run-time parameters.
	eqn_of_state	Integer that selects the equatio of state.
[in]	p_ref	The coordinate-density reference pressure in Pa.
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 1451 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

  type(ocean_grid_type),   intent(in)  :: g            !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)  :: gv           !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                           intent(out) :: t            !< The potential temperature that is being
                                                       !! initialized.
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                           intent(out) :: s            !< The salinity that is being initialized.
  type(param_file_type),   intent(in)  :: param_file   !< A structure to parse for run-time
                                                       !! parameters.
  type(eos_type),          pointer     :: eqn_of_state !< Integer that selects the equatio of state.
  real,                    intent(in)  :: p_ref        !< The coordinate-density reference pressure
                                                       !! in Pa.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                       !! only read parameters without changing h.
!  This function puts the initial layer temperatures and salinities  !
! into T(:,:,:) and S(:,:,:).                                        !

! Arguments: T - The potential temperature that is being initialized.
!  (out)     S - The salinity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      eqn_of_state - integer that selects the equatio of state
!  (in)      P_Ref - The coordinate-density reference pressure in Pa.
  real :: t0(szk_(g)), s0(szk_(g))
  real :: t_ref         ! Reference Temperature
  real :: s_ref         ! Reference Salinity
  real :: pres(szk_(g))      ! An array of the reference pressure in Pa.
  real :: drho_dt(szk_(g))   ! Derivative of density with temperature in kg m-3 K-1.                              !
  real :: drho_ds(szk_(g))   ! Derivative of density with salinity in kg m-3 PSU-1.                             !
  real :: rho_guess(szk_(g)) ! Potential density at T0 & S0 in kg m-3.
  logical :: fit_salin       ! If true, accept the prescribed temperature and fit the salinity.
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=40)  :: mdl = "initialize_temp_salt_fit" ! This subroutine's name.
  integer :: i, j, k, itt, nz
  nz = g%ke

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  call get_param(param_file, mdl, "T_REF", t_ref, &
                 "A reference temperature used in initialization.", &
                 units="degC", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "S_REF", s_ref, &
                 "A reference salinity used in initialization.", units="PSU", &
                 default=35.0, do_not_log=just_read)
  call get_param(param_file, mdl, "FIT_SALINITY", fit_salin, &
                 "If true, accept the prescribed temperature and fit the \n"//&
                 "salinity; otherwise take salinity and fit temperature.", &
                 default=.false., do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  do k=1,nz
    pres(k) = p_ref ; s0(k) = s_ref
    t0(k) = t_ref
  enddo

  call calculate_density(t0(1),s0(1),pres(1),rho_guess(1),eqn_of_state)
  call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,1,eqn_of_state)

  if (fit_salin) then
! A first guess of the layers' temperatures.
    do k=nz,1,-1
      s0(k) = max(0.0, s0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_ds(1))
    enddo
! Refine the guesses for each layer.
    do itt=1,6
      call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state)
      call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,eqn_of_state)
      do k=1,nz
        s0(k) = max(0.0, s0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_ds(k))
      enddo
    enddo
  else
! A first guess of the layers' temperatures.
    do k=nz,1,-1
      t0(k) = t0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_dt(1)
    enddo
    do itt=1,6
      call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state)
      call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,eqn_of_state)
      do k=1,nz
        t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
      enddo
    enddo
  endif

  do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
    t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
  enddo ; enddo ; enddo

  call calltree_leave(trim(mdl)//'()')

Here is the call graph for this function:

Here is the caller graph for this function:

initialize_temp_salt_from_file()

subroutine mom_state_initialization::initialize_temp_salt_from_file ( 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

Parameters

[in]	g	The ocean's grid structure
[out]	t	The potential temperature that is being initialized.
[out]	s	The salinity that is being initialized.
[in]	param_file	A structure to parse for run-time parameters
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 1332 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

  type(ocean_grid_type),                  intent(in)  :: g    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: t !< The potential temperature that is being initialized.
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: s !< The salinity that is being initialized.
  type(param_file_type),                  intent(in)  :: param_file !< A structure to parse for run-time parameters
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
!  This function puts the initial layer temperatures and salinities  !
! into T(:,:,:) and S(:,:,:).                                        !

! Arguments: T - The potential temperature that is being initialized.
!  (out)     S - The salinity that is being initialized.
!  (in)      from_file - .true. if the variables that are set here are to
!                        be read from a file; .false. to be set internally.
!  (in)      filename - The name of the file to read.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=200) :: filename, salt_filename ! Full paths to input files
  character(len=200) :: ts_file, salt_file, inputdir ! Strings for file/path
  character(len=40)  :: mdl = "initialize_temp_salt_from_file"
  character(len=64)  :: temp_var, salt_var ! Temperature and salinity names in files

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  call get_param(param_file, mdl, "TS_FILE", ts_file, &
                 "The initial condition file for temperature.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)

  filename = trim(inputdir)//trim(ts_file)
  if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/TS_FILE", filename)
  call get_param(param_file, mdl, "TEMP_IC_VAR", temp_var, &
                 "The initial condition variable for potential temperature.", &
                 default="PTEMP", do_not_log=just_read)
  call get_param(param_file, mdl, "SALT_IC_VAR", salt_var, &
                 "The initial condition variable for salinity.", &
                 default="SALT", do_not_log=just_read)
  call get_param(param_file, mdl, "SALT_FILE", salt_file, &
                 "The initial condition file for salinity.", &
                 default=trim(ts_file), do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
     " initialize_temp_salt_from_file: Unable to open "//trim(filename))

! Read the temperatures and salinities from netcdf files.           !
  call read_data(filename, temp_var, t(:,:,:), domain=g%Domain%mpp_domain)

  salt_filename = trim(inputdir)//trim(salt_file)
  if (.not.file_exists(salt_filename, g%Domain)) call mom_error(fatal, &
     " initialize_temp_salt_from_file: Unable to open "//trim(salt_filename))

  call read_data(salt_filename, salt_var, s(:,:,:), domain=g%Domain%mpp_domain)

  call calltree_leave(trim(mdl)//'()')

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

Parameters

[in]	g	The ocean's grid structure
[in]	param_file	A structure to parse for run-time parameters
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 1397 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

  type(ocean_grid_type),                  intent(in)  :: g    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: t, s
  type(param_file_type),                  intent(in)  :: param_file !< A structure to parse for run-time parameters
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
!  This function puts the initial layer temperatures and salinities  !
! into T(:,:,:) and S(:,:,:).                                        !

! Arguments: T - The potential temperature that is being initialized.
!  (out)     S - The salinity that is being initialized.
!  (in)      from_file - .true. if the variables that are set here are to
!                        be read from a file; .false. to be set internally.
!  (in)      filename - The name of the file to read.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
  real, dimension(SZK_(G)) :: t0, s0
  integer :: i, j, k
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=200) :: filename, ts_file, inputdir ! Strings for file/path
  character(len=40)  :: mdl = "initialize_temp_salt_from_profile"

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  call get_param(param_file, mdl, "TS_FILE", ts_file, &
                 "The file with the reference profiles for temperature \n"//&
                 "and salinity.", fail_if_missing=.not.just_read, do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
  filename = trim(inputdir)//trim(ts_file)
  call log_param(param_file, mdl, "INPUTDIR/TS_FILE", filename)
  if (.not.file_exists(filename)) call mom_error(fatal, &
     " initialize_temp_salt_from_profile: Unable to open "//trim(filename))

! Read the temperatures and salinities from a netcdf file.           !
  call read_data(filename,"PTEMP",t0(:),domain=g%Domain%mpp_domain)
  call read_data(filename,"SALT", s0(:),domain=g%Domain%mpp_domain)

  do k=1,g%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
    t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
  enddo ; enddo ; enddo

  call calltree_leave(trim(mdl)//'()')

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(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

Parameters

[in]	g	The ocean's grid structure
[in]	param_file	A structure to parse for run-time parameters
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 1551 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

  type(ocean_grid_type),                     intent(in)  :: g          !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: t, s
  type(param_file_type),                     intent(in)  :: param_file !< A structure to parse for
                                                                       !! run-time parameters
  logical,       optional,                   intent(in)  :: just_read_params !< If present and true,
                                                                       !! this call will only read
                                                                       !! parameters without
                                                                       !! changing h.

  ! This subroutine initializes linear profiles for T and S according to
  ! reference surface layer salinity and temperature and a specified range.
  ! Note that the linear distribution is set up with respect to the layer
  ! number, not the physical position).
  integer :: k;
  real  :: delta_s, delta_t
  real  :: s_top, t_top ! Reference salinity and temerature within surface layer
  real  :: s_range, t_range ! Range of salinities and temperatures over the vertical
  real  :: delta
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=40)  :: mdl = "initialize_temp_salt_linear" ! This subroutine's name.

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
  call get_param(param_file, mdl, "T_TOP", t_top, &
                 "Initial temperature of the top surface.", &
                 units="degC", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "T_RANGE", t_range, &
                 "Initial temperature difference (top-bottom).", &
                 units="degC", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "S_TOP", s_top, &
                 "Initial salinity of the top surface.", &
                 units="PSU", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "S_RANGE", s_range, &
                 "Initial salinity difference (top-bottom).", &
                 units="PSU", fail_if_missing=.not.just_read, do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

! ! Prescribe salinity
! delta_S = S_range / ( G%ke - 1.0 );
! S(:,:,1) = S_top;
! do k = 2,G%ke
!   S(:,:,k) = S(:,:,k-1) + delta_S;
! end do
  do k = 1,g%ke
    s(:,:,k) = s_top - s_range*((real(k)-0.5)/real(g%ke))
    t(:,:,k) = t_top - t_range*((real(k)-0.5)/real(g%ke))
  end do

! ! Prescribe temperature
! delta_T = T_range / ( G%ke - 1.0 );
! T(:,:,1) = T_top;
! do k = 2,G%ke
!   T(:,:,k) = T(:,:,k-1) + delta_T;
! end do
! delta = 1;
! T(:,:,G%ke/2 - (delta-1):G%ke/2 + delta) = 1.0;

  call calltree_leave(trim(mdl)//'()')

Here is the call graph for this function:

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(param_file_type), intent(in)  param_file, logical, intent(in)  file_has_thickness, logical, intent(in), optional  just_read_params  )

private

This subroutine reads the layer thicknesses or interface heights from a file.

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[out]	h	The thickness that is being initialized, in m.
[in]	param_file	A structure indicating the open file to parse for model parameter values.
[in]	file_has_thickness	If true, this file contains layer thicknesses; otherwise it contains interface heights.
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 580 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

  type(ocean_grid_type),   intent(in)  :: g    !< The ocean's grid structure
  type(verticalgrid_type), intent(in)  :: gv   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(out) :: h           !< The thickness that is being initialized, in m.
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                      !! to parse for model parameter values.
  logical,                 intent(in)  :: file_has_thickness !< If true, this file contains layer
                                                      !! thicknesses; otherwise it contains
                                                      !! interface heights.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      file_has_thickness - If true, this file contains thicknesses;
!                                 otherwise it contains interface heights.

!  This subroutine reads the layer thicknesses from file.
  real :: eta(szi_(g),szj_(g),szk_(g)+1)
  integer :: inconsistent = 0
  real :: dilate     ! The amount by which each layer is dilated to agree
                     ! with the bottom depth and free surface height, nondim.
  logical :: correct_thickness
  logical :: just_read    ! If true, just read parameters but set nothing.  character(len=20) :: verticalCoordinate
  character(len=40)  :: mdl = "initialize_thickness_from_file" ! This subroutine's name.
  character(len=200) :: filename, thickness_file, inputdir, mesg ! Strings for file/path
  integer :: i, j, k, is, ie, js, je, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) &
    call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".", do_not_log=just_read)
  inputdir = slasher(inputdir)
  call get_param(param_file, mdl, "THICKNESS_FILE", thickness_file, &
                 "The name of the thickness file.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)

  filename = trim(inputdir)//trim(thickness_file)
  if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/THICKNESS_FILE", filename)

  if ((.not.just_read) .and. (.not.file_exists(filename, g%Domain))) call mom_error(fatal, &
         " initialize_thickness_from_file: Unable to open "//trim(filename))

  if (file_has_thickness) then
    if (just_read) return ! All run-time parameters have been read, so return.
    call read_data(filename,"h",h(:,:,:),domain=g%Domain%mpp_domain)
  else
    call get_param(param_file, mdl, "ADJUST_THICKNESS", correct_thickness, &
                 "If true, all mass below the bottom removed if the \n"//&
                 "topography is shallower than the thickness input file \n"//&
                 "would indicate.", default=.false., do_not_log=just_read)
    if (just_read) return ! All run-time parameters have been read, so return.

    call read_data(filename,"eta",eta(:,:,:),domain=g%Domain%mpp_domain)

    if (correct_thickness) then
      call adjustetatofitbathymetry(g, gv, eta, h)
    else
      do k=nz,1,-1 ; do j=js,je ; do i=is,ie
        if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_z)) then
          eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_z
          h(i,j,k) = gv%Angstrom_z
        else
          h(i,j,k) = eta(i,j,k) - eta(i,j,k+1)
        endif
      enddo ; enddo ; enddo

      do j=js,je ; do i=is,ie
        if (abs(eta(i,j,nz+1) + g%bathyT(i,j)) > 1.0) &
          inconsistent = inconsistent + 1
      enddo ; enddo
      call sum_across_pes(inconsistent)

      if ((inconsistent > 0) .and. (is_root_pe())) then
        write(mesg,'("Thickness initial conditions are inconsistent ",'// &
                 '"with topography in ",I8," places.")') inconsistent
        call mom_error(warning, mesg)
      endif
    endif

  endif
  call calltree_leave(trim(mdl)//'()')

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

Definition at line 812 of file MOM_state_initialization.F90.

Referenced by mom_initialize_state().

! search density space for location of layers
  call mom_error(fatal,"  MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED")

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

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[out]	h	The thickness that is being initialized, in m.
[in]	param_file	A structure indicating the open file to parse for model parameter values.
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 749 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

  type(ocean_grid_type),   intent(in)  :: g           !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)  :: gv          !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(out) :: h           !< The thickness that is being initialized, in m.
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                      !! to parse for model parameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

! Arguments: h - The thickness that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.

!  This subroutine initializes the layer thicknesses to be uniform.
  character(len=40)  :: mdl = "initialize_thickness_uniform" ! This subroutine's name.
  real :: e0(szk_(g)+1)   ! The resting interface heights, in m, usually !
                          ! negative because it is positive upward.      !
  real :: eta1d(szk_(g)+1)! Interface height relative to the sea surface !
                          ! positive upward, in m.                       !
  logical :: just_read    ! If true, just read parameters but set nothing.  character(len=20) :: verticalCoordinate
  integer :: i, j, k, is, ie, js, je, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (just_read) return ! This subroutine has no run-time parameters.

  call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  if (g%max_depth<=0.) call mom_error(fatal,"initialize_thickness_uniform: "// &
      "MAXIMUM_DEPTH has a non-sensical value! Was it set?")

  do k=1,nz
    e0(k) = -g%max_depth * real(k-1) / real(nz)
  enddo

  do j=js,je ; do i=is,ie
!    This sets the initial thickness (in m) of the layers.  The      !
!  thicknesses are set to insure that: 1.  each layer is at least an !
!  Angstrom thick, and 2.  the interfaces are where they should be   !
!  based on the resting depths and interface height perturbations,   !
!  as long at this doesn't interfere with 1.                         !
    eta1d(nz+1) = -1.0*g%bathyT(i,j)
    do k=nz,1,-1
      eta1d(k) = e0(k)
      if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_z)) then
        eta1d(k) = eta1d(k+1) + gv%Angstrom_z
        h(i,j,k) = gv%Angstrom_z
      else
        h(i,j,k) = eta1d(k) - eta1d(k+1)
      endif
    enddo
  enddo ; enddo

  call calltree_leave(trim(mdl)//'()')

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(param_file_type), intent(in)  param_file, logical, intent(in), optional  just_read_params  )

private

Parameters

[in]	g	The ocean's grid structure
[out]	u	The zonal velocity that is being initialized, in m s-1
[out]	v	The meridional velocity that is being initialized, in m s-1
[in]	param_file	A structure indicating the open file to parse for modelparameter values.
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 1270 of file MOM_state_initialization.F90.

References my_psi().

Referenced by mom_initialize_state().

  type(ocean_grid_type),                   intent(in)  :: g    !< The ocean's grid structure
  real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u    !< The zonal velocity that is being initialized, in m s-1
  real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v    !< The meridional velocity that is being initialized, in m s-1
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file to
                                                      !! parse for modelparameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine sets the initial velocity components to be circular with
! no flow at edges of domain and center.
  character(len=200) :: mdl = "initialize_velocity_circular"
  real :: circular_max_u
  real :: dpi, psi1, psi2
  logical :: just_read    ! If true, just read parameters but set nothing.
  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  call get_param(param_file, mdl, "CIRCULAR_MAX_U", circular_max_u, &
                 "The amplitude of zonal flow from which to scale the\n"// &
                 "circular stream function (m/s).", &
                 units="m s-1", default=0., do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  dpi=acos(0.0)*2.0 ! pi

  do k=1,nz ; do j=js,je ; do i=isq,ieq
    psi1 = my_psi(i,j)
    psi2 = my_psi(i,j-1)
    u(i,j,k) = (psi1-psi2)/g%dy_Cu(i,j)! *(circular_max_u*G%len_lon/(2.0*dpi))
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    psi1 = my_psi(i,j)
    psi2 = my_psi(i-1,j)
    v(i,j,k) = (psi2-psi1)/g%dx_Cv(i,j)! *(circular_max_u*G%len_lon/(2.0*dpi))
  enddo ; enddo ; enddo

  contains

  real function my_psi(ig,jg) ! in-line function
    integer :: ig, jg
    real :: x, y, r
    x = 2.0*(g%geoLonBu(ig,jg)-g%west_lon)/g%len_lon-1.0  ! -1<x<1
    y = 2.0*(g%geoLatBu(ig,jg)-g%south_lat)/g%len_lat-1.0 ! -1<y<1
    r = sqrt( x**2 + y**2 ) ! Circulat stream fn nis fn of radius only
    r = min(1.0,r) ! Flatten stream function in corners of box
    my_psi = 0.5*(1.0 - cos(dpi*r))
    my_psi = my_psi * (circular_max_u*g%len_lon*1e3/dpi) ! len_lon is in km
  end function my_psi

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(szib_(g),szj_(g), szk_(g)), intent(out)  u, real, dimension(szi_(g),szjb_(g), szk_(g)), 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

Parameters

[in]	g	The ocean's grid structure
[out]	u	The zonal velocity that is being initialized, in m s-1
[out]	v	The meridional velocity that is being initialized, in m s-1
[in]	param_file	A structure indicating the open file to parse for modelparameter values.
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 1144 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

  type(ocean_grid_type),   intent(in)  :: g           !< The ocean's grid structure
  real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u    !< The zonal velocity that is being initialized, in m s-1
  real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v    !< The meridional velocity that is being initialized, in m s-1
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file to
                                                      !! parse for modelparameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine reads the initial velocity components from file
  character(len=40)  :: mdl = "initialize_velocity_from_file" ! This subroutine's name.
  character(len=200) :: filename,velocity_file,inputdir ! Strings for file/path
  logical :: just_read    ! If true, just read parameters but set nothing.

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  call get_param(param_file, mdl, "VELOCITY_FILE", velocity_file, &
                 "The name of the velocity initial condition file.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)

  if (just_read) return ! All run-time parameters have been read, so return.

  filename = trim(inputdir)//trim(velocity_file)
  call log_param(param_file, mdl, "INPUTDIR/VELOCITY_FILE", filename)

  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
         " initialize_velocity_from_file: Unable to open "//trim(filename))

  !  Read the velocities from a netcdf file.
  call read_data(filename,"u",u(:,:,:),domain=g%Domain%mpp_domain,position=east_face)
  call read_data(filename,"v",v(:,:,:),domain=g%Domain%mpp_domain,position=north_face)

  call calltree_leave(trim(mdl)//'()')

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(szib_(g),szj_(g), szk_(g)), intent(out)  u, real, dimension(szi_(g),szjb_(g), szk_(g)), 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

Parameters

[in]	g	The ocean's grid structure
[out]	u	The zonal velocity that is being initialized, in m s-1
[out]	v	The meridional velocity that is being initialized, in m s-1
[in]	param_file	A structure indicating the open file to parse for modelparameter values.
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 1227 of file MOM_state_initialization.F90.

Referenced by mom_initialize_state().

  type(ocean_grid_type),                   intent(in)  :: g    !< The ocean's grid structure
  real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u    !< The zonal velocity that is being initialized, in m s-1
  real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v    !< The meridional velocity that is being initialized, in m s-1
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file to
                                                      !! parse for modelparameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine sets the initial velocity components to uniform
  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
  real    :: initial_u_const, initial_v_const
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=200) :: mdl = "initialize_velocity_uniform" ! This subroutine's name.
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  call get_param(param_file, mdl, "INITIAL_U_CONST", initial_u_const, &
                 "A initial uniform value for the zonal flow.", &
                 units="m s-1", fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl, "INITIAL_V_CONST", initial_v_const, &
                 "A initial uniform value for the meridional flow.", &
                 units="m s-1", fail_if_missing=.not.just_read, do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  do k=1,nz ; do j=js,je ; do i=isq,ieq
    u(i,j,k) = initial_u_const
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    v(i,j,k) = initial_v_const
  enddo ; enddo ; enddo

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

Parameters

[in]	g	The ocean's grid structure
[out]	u	The zonal velocity that is being initialized, in m s-1
[out]	v	The meridional velocity that is being initialized, in m s-1
[in]	param_file	A structure indicating the open file to parse for modelparameter values.
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 1189 of file MOM_state_initialization.F90.

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

Referenced by mom_initialize_state().

  type(ocean_grid_type),                   intent(in)  :: g    !< The ocean's grid structure
  real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u    !< The zonal velocity that is being initialized, in m s-1
  real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v    !< The meridional velocity that is being initialized, in m s-1
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file to
                                                      !! parse for modelparameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
! Arguments: u - The zonal velocity that is being initialized.
!  (out)     v - The meridional velocity that is being initialized.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file -  parameter file type

!   This subroutine sets the initial velocity components to zero
  character(len=200) :: mdl = "initialize_velocity_zero" ! This subroutine's name.
  logical :: just_read    ! If true, just read parameters but set nothing.
  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")

  if (just_read) return ! All run-time parameters have been read, so return.

  do k=1,nz ; do j=js,je ; do i=isq,ieq
    u(i,j,k) = 0.0
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    v(i,j,k) = 0.0
  enddo ; enddo ; enddo

  call calltree_leave(trim(mdl)//'()')

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(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
	)

Parameters

[in,out]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[out]	u	The zonal velocity that is being initialized,
[out]	v	The meridional velocity that is being
[out]	h	Layer thicknesses, in H (usually m or
[in,out]	tv	A structure pointing to various thermodynamic variables
[in,out]	time	Time at the start of the run segment.
[in]	pf	A structure indicating the open file to parse for model parameter values.
[in]	dirs	A structure containing several relevant directory paths.
	restart_cs	A pointer to the restart control structure.
[in]	time_in	Time at the start of the run segment. Time_in overrides any value set for Time.

Definition at line 106 of file MOM_state_initialization.F90.

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(), dome2d_initialization::dome2d_initialize_sponges(), dome_initialization::dome_initialize_sponges(), 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_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(), mom_open_boundary::open_boundary_query(), phillips_initialization::phillips_initialize_sponges(), mom_restart::restore_state(), rossby_front_2d_initialization::rossby_front_initialize_velocity(), scm_cvmix_tests::scm_cvmix_tests_ts_init(), scm_idealized_hurricane::scm_idealized_hurricane_ts_init(), seamount_initialization::seamount_initialize_temperature_salinity(), mom_open_boundary::set_tracer_data(), sloshing_initialization::sloshing_initialize_temperature_salinity(), soliton_initialization::soliton_initialize_thickness(), supercritical_initialization::supercritical_set_obc_data(), trim_for_ice(), and user_initialization::user_initialize_sponges().

Referenced by mom::initialize_mom().

  type(ocean_grid_type),      intent(inout) :: g    !< The ocean's grid structure.
  type(verticalgrid_type),    intent(in)    :: gv   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                              intent(out)   :: u    !< The zonal velocity that is being initialized,
                                                    !! in m s-1
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                              intent(out)   :: v    !< The meridional velocity that is being
                                                    !! initialized, in m s-1
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                              intent(out)   :: h    !< Layer thicknesses, in H (usually m or
                                                    !! kg m-2)
  type(thermo_var_ptrs),      intent(inout) :: tv   !< A structure pointing to various thermodynamic
                                                    !! variables
  type(time_type),            intent(inout) :: time !< Time at the start of the run segment.
  type(param_file_type),      intent(in)    :: pf   !< A structure indicating the open file to parse
                                                    !! for model parameter values.
  type(directories),          intent(in)    :: dirs !< A structure containing several relevant
                                                    !! directory paths.
  type(mom_restart_cs),       pointer       :: restart_cs !< A pointer to the restart control
                                                    !! structure.
  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), optional,  intent(in)    :: time_in !< Time at the start of the run segment.
                                                     !! Time_in overrides any value set for Time.
! Arguments: u  - Zonal velocity, in m s-1.
!  (out)     v  - Meridional velocity, in m s-1.
!  (out)     h  - Layer thickness, in m.
!  (out)     tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (out)     Time    - Time at the start of the run segment.
!  (inout)   G       - The ocean's grid structure.
!  (in)      GV      - The ocean's vertical grid structure.
!  (in)      PF      - A structure indicating the open file to parse for
!                      model parameter values.
!  (in)      dirs    - A structure containing several relevant directory paths.
!  (inout)   restart_CS - A pointer to the restart control structure.
!  (inout)   CS      - A structure of pointers to be exchanged with MOM.F90.
!  (in)      Time_in - Time at the start of the run segment. Time_in overrides
!                      any value set for Time.

  character(len=200) :: filename   ! The name of an input file.
  character(len=200) :: filename2  ! The name of an input files.
  character(len=200) :: inputdir   ! The directory where NetCDF input files are.
  character(len=200) :: config
  logical :: from_z_file, useale
  logical :: new_sim
  integer :: write_geom
  logical :: use_temperature, use_sponge
  logical :: use_eos    ! If true, density is calculated from T & S using an
                        ! equation of state.
  logical :: depress_sfc ! If true, remove the mass that would be displaced
                         ! by a large surface pressure by squeezing the column.
  logical :: trim_ic_for_p_surf ! If true, remove the mass that would be displaced
                         ! by a large surface pressure, such as with an ice sheet.
  logical :: regrid_accelerate
  integer :: regrid_iterations
  logical :: analytic_fv_pgf, obsol_test
  logical :: convert
  logical :: just_read  ! If true, only read the parameters because this
                        ! is a run from a restart file; this option
                        ! allows the use of Fatal unused parameters.
  type(eos_type), pointer :: eos => null()
  logical :: debug    ! indicates whether to write debugging output
! This include declares and sets the variable "version".
#include "version_variable.h"
  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  call calltree_enter("MOM_initialize_state(), MOM_state_initialization.F90")
  call log_version(pf, mdl, version, "")
  call get_param(pf, mdl, "DEBUG", debug, default=.false.)

  new_sim = .false.
  if ((dirs%input_filename(1:1) == 'n') .and. &
      (len_trim(dirs%input_filename) == 1)) new_sim = .true.

  just_read = .not.new_sim

  call get_param(pf, mdl, "INPUTDIR", inputdir, &
         "The directory in which input files are found.", default=".")
  inputdir = slasher(inputdir)

  use_temperature = ASSOCIATED(tv%T)
  useale = associated(ale_csp)
  use_eos = associated(tv%eqn_of_state)
  if (use_eos) eos => tv%eqn_of_state

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

  if (new_sim) then
    call mom_mesg("Run initialized internally.", 3)

    if (present(time_in)) time = time_in
    ! Otherwise leave Time at its input value.

    ! This initialization should not be needed. Certainly restricting it
    ! to the computational domain helps detect possible uninitialized
    ! data in halos which should be covered by the pass_var(h) later.
    !do k = 1, nz; do j = js, je; do i = is, ie
    !  h(i,j,k) = 0.
    !enddo
  endif

  ! The remaining initialization calls are done, regardless of whether the
  ! fields are actually initialized here (if just_read=.false.) or whether it
  ! is just to make sure that all valid parameters are read to enable the
  ! detection of unused parameters.
  call get_param(pf, mdl, "INIT_LAYERS_FROM_Z_FILE", from_z_file, &
             "If true, intialize the layer thicknesses, temperatures, \n"//&
             "and salnities from a Z-space file on a latitude- \n"//&
             "longitude grid.", default=.false., do_not_log=just_read)

  if (from_z_file) then
!     Initialize thickness and T/S from z-coordinate data in a file.
    if (.NOT.use_temperature) call mom_error(fatal,"MOM_initialize_state : "//&
       "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true")

    call mom_temp_salt_initialize_from_z(h, tv, g, gv, pf, just_read_params=just_read)

  else
!     Initialize thickness, h.
    call get_param(pf, mdl, "THICKNESS_CONFIG", config, &
             "A string that determines how the initial layer \n"//&
             "thicknesses are specified for a new run: \n"//&
             " \t file - read interface heights from the file specified \n"//&
             " \t thickness_file - read thicknesses from the file specified \n"//&
             " \t\t by (THICKNESS_FILE).\n"//&
             " \t coord - determined by ALE coordinate.\n"//&
             " \t uniform - uniform thickness layers evenly distributed \n"//&
             " \t\t between the surface and MAXIMUM_DEPTH. \n"//&
             " \t DOME - use a slope and channel configuration for the \n"//&
             " \t\t DOME sill-overflow test case. \n"//&
             " \t ISOMIP - use a configuration for the \n"//&
             " \t\t ISOMIP test case. \n"//&
             " \t benchmark - use the benchmark test case thicknesses. \n"//&
             " \t search - search a density profile for the interface \n"//&
             " \t\t densities. This is not yet implemented. \n"//&
             " \t circle_obcs - the circle_obcs test case is used. \n"//&
             " \t DOME2D - 2D version of DOME initialization. \n"//&
             " \t adjustment2d - TBD AJA. \n"//&
             " \t sloshing - TBD AJA. \n"//&
             " \t seamount - TBD AJA. \n"//&
             " \t soliton - Equatorial Rossby soliton. \n"//&
             " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
             " \t USER - call a user modified routine.", &
             fail_if_missing=new_sim, do_not_log=just_read)
    select case (trim(config))
       case ("file"); call initialize_thickness_from_file(h, g, gv, pf, .false., just_read_params=just_read)
       case ("thickness_file"); call initialize_thickness_from_file(h, g, gv, pf, .true., just_read_params=just_read)
       case ("coord")
         if (new_sim .and. useale) then
           call ale_initthicknesstocoord( ale_csp, g, gv, h )
         elseif (new_sim) then
           call mom_error(fatal, "MOM_initialize_state: USE_REGRIDDING must be True "//&
                                 "for THICKNESS_CONFIG of 'coord'")
         endif
       case ("uniform"); call initialize_thickness_uniform(h, g, gv, pf, &
                                  just_read_params=just_read)
       case ("DOME"); call dome_initialize_thickness(h, g, gv, pf, &
                               just_read_params=just_read)
       case ("ISOMIP"); call isomip_initialize_thickness(h, g, gv, pf, tv, &
                                 just_read_params=just_read)
       case ("benchmark"); call benchmark_initialize_thickness(h, g, gv, pf, &
                                    tv%eqn_of_state, tv%P_Ref, just_read_params=just_read)
       case ("search"); call initialize_thickness_search
       case ("circle_obcs"); call circle_obcs_initialize_thickness(h, g, gv, pf, &
                                      just_read_params=just_read)
       case ("lock_exchange"); call lock_exchange_initialize_thickness(h, g, gv, &
                                        pf, just_read_params=just_read)
       case ("external_gwave"); call external_gwave_initialize_thickness(h, g, &
                                         pf, just_read_params=just_read)
       case ("DOME2D"); call dome2d_initialize_thickness(h, g, gv, pf, &
                                 just_read_params=just_read)
       case ("adjustment2d"); call adjustment_initialize_thickness(h, g, gv, &
                                       pf, just_read_params=just_read)
       case ("sloshing"); call sloshing_initialize_thickness(h, g, gv, pf, &
                                   just_read_params=just_read)
       case ("seamount"); call seamount_initialize_thickness(h, g, gv, pf, &
                                   just_read_params=just_read)
       case ("soliton"); call soliton_initialize_thickness(h, g)
       case ("phillips"); call phillips_initialize_thickness(h, g, gv, pf, &
                                   just_read_params=just_read)
       case ("rossby_front"); call rossby_front_initialize_thickness(h, g, gv, &
                                       pf, just_read_params=just_read)
       case ("USER"); call user_initialize_thickness(h, g, pf, tv%T, &
                               just_read_params=just_read)
       case default ; call mom_error(fatal,  "MOM_initialize_state: "//&
            "Unrecognized layer thickness configuration "//trim(config))
    end select

!     Initialize temperature and salinity (T and S).
    if ( use_temperature ) then
      call get_param(pf, mdl, "TS_CONFIG", config, &
             "A string that determines how the initial tempertures \n"//&
             "and salinities are specified for a new run: \n"//&
             " \t file - read velocities from the file specified \n"//&
             " \t\t by (TS_FILE). \n"//&
             " \t fit - find the temperatures that are consistent with \n"//&
             " \t\t the layer densities and salinity S_REF. \n"//&
             " \t TS_profile - use temperature and salinity profiles \n"//&
             " \t\t (read from TS_FILE) to set layer densities. \n"//&
             " \t benchmark - use the benchmark test case T & S. \n"//&
             " \t linear - linear in logical layer space. \n"//&
             " \t DOME2D - 2D DOME initialization. \n"//&
             " \t ISOMIP - ISOMIP initialization. \n"//&
             " \t adjustment2d - TBD AJA. \n"//&
             " \t sloshing - TBD AJA. \n"//&
             " \t seamount - TBD AJA. \n"//&
             " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
             " \t SCM_ideal_hurr - used in the SCM idealized hurricane test.\n"//&
             " \t SCM_CVmix_tests - used in the SCM CVmix tests.\n"//&
             " \t USER - call a user modified routine.", &
             fail_if_missing=new_sim, do_not_log=just_read)
!              " \t baroclinic_zone - an analytic baroclinic zone. \n"//&
      select case (trim(config))
        case ("fit"); call initialize_temp_salt_fit(tv%T, tv%S, g, gv, pf, &
                               eos, tv%P_Ref, just_read_params=just_read)
        case ("file"); call initialize_temp_salt_from_file(tv%T, tv%S, g, &
                                pf, just_read_params=just_read)
        case ("benchmark"); call benchmark_init_temperature_salinity(tv%T, tv%S, &
                                     g, gv, pf, eos, tv%P_Ref, just_read_params=just_read)
        case ("TS_profile") ; call initialize_temp_salt_from_profile(tv%T, tv%S, &
                                       g, pf, just_read_params=just_read)
        case ("linear"); call initialize_temp_salt_linear(tv%T, tv%S, g, pf, &
                                  just_read_params=just_read)
        case ("DOME2D"); call dome2d_initialize_temperature_salinity ( tv%T, &
                                  tv%S, h, g, pf, eos, just_read_params=just_read)
        case ("ISOMIP"); call isomip_initialize_temperature_salinity ( tv%T, &
                                  tv%S, h, g, gv, pf, eos, just_read_params=just_read)
        case ("adjustment2d"); call adjustment_initialize_temperature_salinity ( tv%T, &
                                        tv%S, h, g, pf, eos, just_read_params=just_read)
        case ("baroclinic_zone"); call baroclinic_zone_init_temperature_salinity( tv%T, &
                                           tv%S, h, g, pf, just_read_params=just_read)
        case ("sloshing"); call sloshing_initialize_temperature_salinity(tv%T, &
                                    tv%S, h, g, pf, eos, just_read_params=just_read)
        case ("seamount"); call seamount_initialize_temperature_salinity(tv%T, &
                                    tv%S, h, g, gv, pf, eos, just_read_params=just_read)
        case ("rossby_front"); call rossby_front_initialize_temperature_salinity ( tv%T, &
                                        tv%S, h, g, pf, eos, just_read_params=just_read)
        case ("SCM_ideal_hurr"); call scm_idealized_hurricane_ts_init ( tv%T, &
                                          tv%S, h, g, gv, pf, just_read_params=just_read)
        case ("SCM_CVmix_tests"); call scm_cvmix_tests_ts_init (tv%T, &
                                           tv%S, h, g, gv, pf, just_read_params=just_read)
        case ("dense"); call dense_water_initialize_ts(g, gv, pf, eos, tv%T, tv%S, &
                                 h, just_read_params=just_read)
        case ("USER"); call user_init_temperature_salinity(tv%T, tv%S, g, pf, eos, &
                                just_read_params=just_read)
        case default ; call mom_error(fatal,  "MOM_initialize_state: "//&
               "Unrecognized Temp & salt configuration "//trim(config))
      end select
    endif
  endif  ! not from_Z_file.

  ! The thicknesses in halo points might be needed to initialize the velocities.
  if (new_sim) call pass_var(h, g%Domain)

!   Initialize velocity components, u and v
  call get_param(pf, mdl, "VELOCITY_CONFIG", config, &
       "A string that determines how the initial velocities \n"//&
       "are specified for a new run: \n"//&
       " \t file - read velocities from the file specified \n"//&
       " \t\t by (VELOCITY_FILE). \n"//&
       " \t zero - the fluid is initially at rest. \n"//&
       " \t uniform - the flow is uniform (determined by\n"//&
       " \t\t parameters INITIAL_U_CONST and INITIAL_V_CONST).\n"//&
       " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
       " \t soliton - Equatorial Rossby soliton.\n"//&
       " \t USER - call a user modified routine.", default="zero", &
       do_not_log=just_read)
  select case (trim(config))
     case ("file"); call initialize_velocity_from_file(u, v, g, pf, &
                             just_read_params=just_read)
     case ("zero"); call initialize_velocity_zero(u, v, g, pf, &
                             just_read_params=just_read)
     case ("uniform"); call initialize_velocity_uniform(u, v, g, pf, &
                                just_read_params=just_read)
     case ("circular"); call initialize_velocity_circular(u, v, g, pf, &
                                 just_read_params=just_read)
     case ("phillips"); call phillips_initialize_velocity(u, v, g, gv, pf, &
                                 just_read_params=just_read)
     case ("rossby_front"); call rossby_front_initialize_velocity(u, v, h, &
                                     g, gv, pf, just_read_params=just_read)
     case ("soliton"); call soliton_initialize_velocity(u, v, h, g)
     case ("USER"); call user_initialize_velocity(u, v, g, pf, &
                             just_read_params=just_read)
     case default ; call mom_error(fatal,  "MOM_initialize_state: "//&
          "Unrecognized velocity configuration "//trim(config))
  end select

  if (new_sim) call pass_vector(u, v, g%Domain)
  if (debug .and. new_sim) then
    call uvchksum("MOM_initialize_state [uv]", u, v, g%HI, haloshift=1)
  endif

!   Optionally convert the thicknesses from m to kg m-2.  This is particularly
! useful in a non-Boussinesq model.
  call get_param(pf, mdl, "CONVERT_THICKNESS_UNITS", convert, &
               "If true,  convert the thickness initial conditions from \n"//&
               "units of m to kg m-2 or vice versa, depending on whether \n"//&
               "BOUSSINESQ is defined. This does not apply if a restart \n"//&
               "file is read.", default=.false., do_not_log=just_read)
  if (new_sim) then
    if (convert .and. .not. gv%Boussinesq) then
      ! Convert h from m to kg m-2 then to thickness units (H)
      call convert_thickness(h, g, gv, tv)
    elseif (gv%Boussinesq) then
      ! Convert h from m to thickness units (H)
      do k = 1, nz; do j = js, je; do i = is, ie
        h(i,j,k) = h(i,j,k)*gv%m_to_H
      enddo ; enddo ; enddo
    else
      do k = 1, nz; do j = js, je; do i = is, ie
        h(i,j,k) = h(i,j,k)*gv%kg_m2_to_H
      enddo ; enddo ; enddo
    endif
  endif

!  Remove the mass that would be displaced by an ice shelf or inverse barometer.
  call get_param(pf, mdl, "DEPRESS_INITIAL_SURFACE", depress_sfc, &
               "If true,  depress the initial surface to avoid huge \n"//&
               "tsunamis when a large surface pressure is applied.", &
               default=.false., do_not_log=just_read)
  call get_param(pf, mdl, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, &
               "If true, cuts way the top of the column for initial conditions\n"//&
               "at the depth where the hydrostatic presure matches the imposed\n"//&
               "surface pressure which is read from file.", default=.false., &
               do_not_log=just_read)
  if (depress_sfc .and. trim_ic_for_p_surf) call mom_error(fatal, "MOM_initialize_state: "//&
           "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True")
  if (depress_sfc) call depress_surface(h, g, gv, pf, tv, just_read_params=just_read)
  if (trim_ic_for_p_surf) call trim_for_ice(pf, g, gv, ale_csp, tv, h, just_read_params=just_read)

  ! Perhaps we want to run the regridding coordinate generator for multiple
  ! iterations here so the initial grid is consistent with the coordinate
  if (useale) then
    call get_param(pf, mdl, "REGRID_ACCELERATE_INIT", regrid_accelerate, &
         "If true, runs REGRID_ACCELERATE_ITERATIONS iterations of the regridding\n"//&
         "algorithm to push the initial grid to be consistent with the initial\n"//&
         "condition. Useful only for state-based and iterative coordinates.", &
         default=.false., do_not_log=just_read)
    if (regrid_accelerate) then
      call get_param(pf, mdl, "REGRID_ACCELERATE_ITERATIONS", regrid_iterations, &
           "The number of regridding iterations to perform to generate\n"//&
           "an initial grid that is consistent with the initial conditions.", &
           default=1, do_not_log=just_read)

      if (new_sim) &
        call ale_regrid_accelerated(ale_csp, g, gv, h, tv, regrid_iterations, h, u, v)
    endif
  endif
  ! This is the end of the block of code that might have initialized fields
  ! internally at the start of a new run.

  if (.not.new_sim) then ! This block restores the state from a restart file.
    !    This line calls a subroutine that reads the initial conditions  !
    !  from a previously generated file.                                 !
    call restore_state(dirs%input_filename, dirs%restart_input_dir, time, &
                       g, restart_cs)
    if (present(time_in)) time = time_in
  endif

  if ( use_temperature ) then
    call pass_var(tv%T, g%Domain, complete=.false.)
    call pass_var(tv%S, g%Domain, complete=.false.)
  endif
  call pass_var(h, g%Domain)

  if (debug) then
    call hchksum(h, "MOM_initialize_state: h ", g%HI, haloshift=1, scale=gv%H_to_m)
    if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", g%HI, haloshift=1)
    if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", g%HI, haloshift=1)
  endif

  call get_param(pf, mdl, "SPONGE", use_sponge, &
                 "If true, sponges may be applied anywhere in the domain. \n"//&
                 "The exact location and properties of those sponges are \n"//&
                 "specified via SPONGE_CONFIG.", default=.false.)
  if ( use_sponge ) then
    call get_param(pf, mdl, "SPONGE_CONFIG", config, &
                 "A string that sets how the sponges are configured: \n"//&
                 " \t file - read sponge properties from the file \n"//&
                 " \t\t specified by (SPONGE_FILE).\n"//&
                 " \t ISOMIP - apply ale sponge in the ISOMIP case \n"//&
                 " \t DOME - use a slope and channel configuration for the \n"//&
                 " \t\t DOME sill-overflow test case. \n"//&
                 " \t BFB - Sponge at the southern boundary of the domain\n"//&
                 " \t\t for buoyancy-forced basin case.\n"//&
                 " \t USER - call a user modified routine.", default="file")
    select case (trim(config))
      case ("DOME"); call dome_initialize_sponges(g, gv, tv, pf, sponge_csp)
      case ("DOME2D"); call dome2d_initialize_sponges(g, gv, tv, pf, useale, &
                                                      sponge_csp, ale_sponge_csp)
      case ("ISOMIP"); call isomip_initialize_sponges(g, gv, tv, pf, useale, &
                                                     sponge_csp, ale_sponge_csp)
      case ("USER"); call user_initialize_sponges(g, use_temperature, tv, &
                                               pf, sponge_csp, h)
      case ("BFB"); call bfb_initialize_sponges_southonly(g, use_temperature, tv, &
                                               pf, sponge_csp, h)
      case ("phillips"); call phillips_initialize_sponges(g, use_temperature, tv, &
                                               pf, sponge_csp, h)
      case ("dense"); call dense_water_initialize_sponges(g, gv, tv, pf, useale, &
           sponge_csp, ale_sponge_csp)
      case ("file"); call initialize_sponges_file(g, gv, use_temperature, tv, &
                                               pf, sponge_csp)
      case default ; call mom_error(fatal,  "MOM_initialize_state: "//&
             "Unrecognized sponge configuration "//trim(config))
    end select
  endif

  ! Reads OBC parameters not pertaining to the location of the boundaries
  call open_boundary_init(g, pf, obc)

  ! This controls user code for setting open boundary data
  if (associated(obc)) then
    call get_param(pf, mdl, "OBC_USER_CONFIG", config, &
                 "A string that sets how the user code is invoked to set open\n"//&
                 " boundary data: \n"//&
                 "   DOME - specified inflow on northern boundary\n"//&
                 "   tidal_bay - Flather with tidal forcing on eastern boundary\n"//&
                 "   supercritical - now only needed here for the allocations\n"//&
                 "   Kelvin - barotropic Kelvin wave forcing on the western boundary\n"//&
                 "   shelfwave - Flather with shelf wave forcing on western boundary\n"//&
                 "   USER - user specified", default="none")
    if (trim(config) == "DOME") then
      call dome_set_obc_data(obc, tv, g, gv, pf, tracer_reg)
    elseif (lowercase(trim(config)) == "supercritical") then
      call supercritical_set_obc_data(obc, g, pf)
    elseif (trim(config) == "tidal_bay") then
      obc%update_OBC = .true.
    elseif (trim(config) == "Kelvin") then
      obc%update_OBC = .true.
    elseif (trim(config) == "shelfwave") then
      obc%update_OBC = .true.
    elseif (trim(config) == "USER") then
      call user_set_obc_data(obc, tv, g, pf, tracer_reg)
    elseif (.not. trim(config) == "none") then
      call mom_error(fatal, "The open boundary conditions specified by "//&
              "OBC_USER_CONFIG = "//trim(config)//" have not been fully implemented.")
    endif
    if (open_boundary_query(obc, apply_open_obc=.true.)) then
      call set_tracer_data(obc, tv, h, g, pf, tracer_reg)
    endif
  endif
! if (open_boundary_query(OBC, apply_nudged_OBC=.true.)) then
!   call set_3D_OBC_data(OBC, tv, h, G, PF, tracer_Reg)
! endif
  ! Still need a way to specify the boundary values
  if (debug.and.associated(obc)) then
    call hchksum(g%mask2dT, 'MOM_initialize_state: mask2dT ', g%HI)
    call uvchksum('MOM_initialize_state: mask2dC[uv]', g%mask2dCu,  &
                  g%mask2dCv, g%HI)
    call qchksum(g%mask2dBu, 'MOM_initialize_state: mask2dBu ', g%HI)
  endif

  call calltree_leave('MOM_initialize_state()')

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(thermo_var_ptrs), intent(in)  tv  )

private

Run simple unit tests.

Parameters

[in,out]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	tv	Thermodynamics structure.

Definition at line 2252 of file MOM_state_initialization.F90.

References cut_off_column_top().

  type(ocean_grid_type),     intent(inout) :: g    !< The ocean's grid structure.
  type(verticalgrid_type),   intent(in)    :: gv   !< The ocean's vertical grid structure.
  type(thermo_var_ptrs),     intent(in)    :: tv   !< Thermodynamics structure.
  ! Local variables
  integer, parameter :: nk=5
  real, dimension(nk) :: t, t_t, t_b, s, s_t, s_b, rho, h, z
  real, dimension(nk+1) :: e
  integer :: k
  real :: p_tot, p_t, p_b, z_out
  type(remapping_cs), pointer :: remap_cs => null()

  do k = 1, nk
    h(k) = 100.
  enddo
  e(1) = 0.
  do k = 1, nk
    e(k+1) = e(k) - h(k)
  enddo
  p_tot = 0.
  do k = 1, nk
    z(k) = 0.5 * ( e(k) + e(k+1) )
    t_t(k) = 20.+(0./500.)*e(k)
    t(k)   = 20.+(0./500.)*z(k)
    t_b(k) = 20.+(0./500.)*e(k+1)
    s_t(k) = 35.-(0./500.)*e(k)
    s(k)   = 35.+(0./500.)*z(k)
    s_b(k) = 35.-(0./500.)*e(k+1)
    call calculate_density(0.5*(t_t(k)+t_b(k)), 0.5*(s_t(k)+s_b(k)), -gv%Rho0*gv%g_Earth*z(k), rho(k), tv%eqn_of_state)
    p_tot = p_tot + gv%g_Earth * rho(k) * h(k)
  enddo

  p_t = 0.
  do k = 1, nk
    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, gv%Rho0, gv%g_Earth, tv%eqn_of_state, p_b, z_out)
    write(0,*) k,p_t,p_b,0.5*p_tot,e(k),e(k+1),z_out
    p_t = p_b
  enddo
  write(0,*) p_b,p_tot

  write(0,*) ''
  write(0,*) ' ==================================================================== '
  write(0,*) ''
  write(0,*) h
  call cut_off_column_top(nk, tv, gv%Rho0, gv%g_Earth, -e(nk+1), gv%Angstrom, &
               t, t_t, t_b, s, s_t, s_b, 0.5*p_tot, h, remap_cs)
  write(0,*) h

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( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  h, type(thermo_var_ptrs), intent(inout)  tv, type(ocean_grid_type), intent(inout)  G, type(verticalgrid_type), intent(in)  GV, 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]	g	The ocean's grid structure
[out]	h	Layer thicknesses being initialized, in m
[in,out]	tv	A structure pointing to various thermodynamic variables including temperature and salinity
[in]	gv	The ocean's vertical grid structure
[in]	pf	A structure indicating the open file to parse for model parameter values.
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 1815 of file MOM_state_initialization.F90.

References adjustetatofitbathymetry(), mom_ale::ale_build_grid(), mom_ale::ale_remap_scalar(), mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_eos::convert_temp_salt_for_teos10(), midas_vertmap::determine_temperature(), mom_regridding::getcoordinateresolution(), and mom_regridding::set_regrid_params().

Referenced by mom_initialize_state().

! This subroutine was written by M. Harrison, with input from R. Hallberg & A. Adcroft.
!
! Arguments:
!  (out)     h  - Layer thickness, in m.
!  (out)     tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (inout)   G       - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      PF      - A structure indicating the open file to parse for
!                      model parameter values.

  type(ocean_grid_type),   intent(inout) :: g    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(out)   :: h    !< Layer thicknesses being initialized, in m
  type(thermo_var_ptrs),   intent(inout) :: tv   !< A structure pointing to various thermodynamic
                                                 !! variables including temperature and salinity
  type(verticalgrid_type), intent(in)    :: gv   !< The ocean's vertical grid structure
  type(param_file_type),   intent(in)    :: pf   !< A structure indicating the open file
                                                 !! to parse for model parameter values.
  logical,       optional, intent(in)    :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

  character(len=200) :: filename   ! The name of an input file containing temperature
                                   ! and salinity in z-space; also used for  ice shelf area.
  character(len=200) :: shelf_file ! The name of an input file used for  ice shelf area.
  character(len=200) :: inputdir ! The directory where NetCDF input filesare.
  character(len=200) :: mesg, area_varname, ice_shelf_file

  type(eos_type), pointer :: eos => null()

! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_initialize_layers_from_Z" ! This module's name.

  integer :: is, ie, js, je, nz ! compute domain indices
  integer :: isc,iec,jsc,jec    ! global compute domain indices
  integer :: isg, ieg, jsg, jeg ! global extent
  integer :: isd, ied, jsd, jed ! data domain indices

  integer :: i, j, k, ks, np, ni, nj
  integer :: idbg, jdbg
  integer :: nkml, nkbl         ! number of mixed and buffer layers

  integer :: kd, inconsistent
  real    :: pi_180             ! for conversion from degrees to radians

  real, dimension(:,:), pointer :: shelf_area
  real    :: min_depth
  real    :: dilate
  real    :: missing_value_temp, missing_value_salt
  logical :: correct_thickness
  character(len=40) :: potemp_var, salin_var
  character(len=8)  :: laynum

  integer, parameter :: niter=10   ! number of iterations for t/s adjustment to layer density
  logical :: just_read    ! If true, just read parameters but set nothing.
  logical            :: adjust_temperature = .true.  ! fit t/s to target densities
  real, parameter    :: missing_value = -1.e20
  real, parameter    :: temp_land_fill = 0.0, salt_land_fill = 35.0
  logical :: reentrant_x, tripolar_n,dbg
  logical :: debug = .false.  ! manually set this to true for verbose output

  !data arrays
  real, dimension(:), allocatable :: z_edges_in, z_in, rb
  real, dimension(:,:,:), allocatable, target :: temp_z, salt_z, mask_z
  real, dimension(:,:,:), allocatable :: rho_z
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: zi
  real, dimension(SZI_(G),SZJ_(G))  :: nlevs
  real, dimension(SZI_(G))   :: press


  ! Local variables for ALE remapping
  real, dimension(:), allocatable :: htarget
  real, dimension(:,:), allocatable :: area_shelf_h
  real, dimension(:,:), allocatable, target  :: frac_shelf_h
  real, dimension(:,:,:), allocatable :: tmpt1din, tmps1din, h1, tmp_mask_in
  real :: ztopofcell, zbottomofcell
  type(regridding_cs) :: regridcs ! Regridding parameters and work arrays
  type(remapping_cs) :: remapcs ! Remapping parameters and work arrays

  logical :: homogenize, usealeremapping, remap_full_column, remap_general, remap_old_alg
  logical :: use_ice_shelf
  character(len=10) :: remappingscheme
  real :: tempavg, saltavg
  integer :: npoints, ans
  integer :: id_clock_routine, id_clock_read, id_clock_interp, id_clock_fill, id_clock_ale

  id_clock_routine = cpu_clock_id('(Initialize from Z)', grain=clock_routine)
  id_clock_ale = cpu_clock_id('(Initialize from Z) ALE', grain=clock_loop)

  call cpu_clock_begin(id_clock_routine)

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg

  pi_180=atan(1.0)/45.

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
  if (.not.just_read) call log_version(pf, mdl, version, "")

  inputdir = "." ;  call get_param(pf, mdl, "INPUTDIR", inputdir)
  inputdir = slasher(inputdir)

  eos => tv%eqn_of_state

! call mpp_get_compute_domain(G%domain%mpp_domain,isc,iec,jsc,jec)

  reentrant_x = .false. ; call get_param(pf, mdl, "REENTRANT_X", reentrant_x,default=.true.)
  tripolar_n = .false. ;  call get_param(pf, mdl, "TRIPOLAR_N", tripolar_n, default=.false.)
  call get_param(pf, mdl, "MINIMUM_DEPTH", min_depth, default=0.0)

  call get_param(pf, mdl, "NKML",nkml,default=0)
  call get_param(pf, mdl, "NKBL",nkbl,default=0)

  call get_param(pf, mdl, "TEMP_SALT_Z_INIT_FILE",filename, &
                 "The name of the z-space input file used to initialize \n"//&
                 "the layer thicknesses, temperatures and salinities.", &
                 default="temp_salt_z.nc", do_not_log=just_read)
  filename = trim(inputdir)//trim(filename)
  call get_param(pf, mdl, "Z_INIT_FILE_PTEMP_VAR", potemp_var, &
                 "The name of the potential temperature variable in \n"//&
                 "TEMP_SALT_Z_INIT_FILE.", default="ptemp", do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_FILE_SALT_VAR", salin_var, &
                 "The name of the salinity variable in \n"//&
                 "TEMP_SALT_Z_INIT_FILE.", default="salt", do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_HOMOGENIZE", homogenize, &
                 "If True, then horizontally homogenize the interpolated \n"//&
                 "initial conditions.", default=.false., do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_ALE_REMAPPING", usealeremapping, &
                 "If True, then remap straight to model coordinate from file.",&
                 default=.false., do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_REMAPPING_SCHEME", remappingscheme, &
                 "The remapping scheme to use if using Z_INIT_ALE_REMAPPING\n"//&
                 "is True.", default="PPM_IH4", do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_REMAP_GENERAL", remap_general, &
                 "If false, only initializes to z* coordinates.\n"//&
                 "If true, allows initialization directly to general coordinates.",&
                 default=.false., do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_REMAP_FULL_COLUMN", remap_full_column, &
                 "If false, only reconstructs profiles for valid data points.\n"//&
                 "If true, inserts vanished layers below the valid data.",&
                 default=remap_general, do_not_log=just_read)
  call get_param(pf, mdl, "Z_INIT_REMAP_OLD_ALG", remap_old_alg, &
                 "If false, uses the preferred remapping algorithm for initialization.\n"//&
                 "If true, use an older, less robust algorithm for remapping.",&
                 default=.true., do_not_log=just_read)
  call get_param(pf, mdl, "ICE_SHELF", use_ice_shelf, default=.false.)
  if (use_ice_shelf) then
    call get_param(pf, mdl, "ICE_THICKNESS_FILE", ice_shelf_file, &
                 "The file from which the ice bathymetry and area are read.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
    shelf_file = trim(inputdir)//trim(ice_shelf_file)
    if (.not.just_read) call log_param(pf, mdl, "INPUTDIR/THICKNESS_FILE", shelf_file)
    call get_param(pf, mdl, "ICE_AREA_VARNAME", area_varname, &
                 "The name of the area variable in ICE_THICKNESS_FILE.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  endif
  if (.not.usealeremapping) then
    call get_param(pf, mdl, "ADJUST_THICKNESS", correct_thickness, &
                 "If true, all mass below the bottom removed if the \n"//&
                 "topography is shallower than the thickness input file \n"//&
                 "would indicate.", default=.false., do_not_log=just_read)

    call get_param(pf, mdl, "FIT_TO_TARGET_DENSITY_IC", adjust_temperature, &
                 "If true, all the interior layers are adjusted to \n"//&
                 "their target densities using mostly temperature \n"//&
                 "This approach can be problematic, particularly in the \n"//&
                 "high latitudes.", default=.true., do_not_log=just_read)
  endif
  if (just_read) then
    call cpu_clock_end(id_clock_routine)
    return ! All run-time parameters have been read, so return.
  endif

!   Read input grid coordinates for temperature and salinity field
!   in z-coordinate dataset. The file is REQUIRED to contain the
!   following:
!
!   dimension variables:
!            lon (degrees_E), lat (degrees_N), depth(meters)
!   variables:
!            ptemp(lon,lat,depth) : degC, potential temperature
!            salt (lon,lat,depth) : PSU, salinity
!
!   The first record will be read if there are multiple time levels.
!   The observation grid MUST tile the model grid. If the model grid extends
!   to the North/South Pole past the limits of the input data, they are extrapolated using the average
!   value at the northernmost/southernmost latitude.

  call horiz_interp_and_extrap_tracer(filename, potemp_var,1.0,1, &
       g, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, tripolar_n, homogenize)

  call horiz_interp_and_extrap_tracer(filename, salin_var,1.0,1, &
       g, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, tripolar_n, homogenize)

  kd = size(z_in,1)

  allocate(rho_z(isd:ied,jsd:jed,kd))
  allocate(area_shelf_h(isd:ied,jsd:jed))
  allocate(frac_shelf_h(isd:ied,jsd:jed))

  press(:)=tv%p_ref

  !Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 or NEMO
  call convert_temp_salt_for_teos10(temp_z,salt_z, press, g, kd, mask_z, eos)

  do k=1,kd
    do j=js,je
      call calculate_density(temp_z(:,j,k),salt_z(:,j,k), press, rho_z(:,j,k), is, ie, eos)
    enddo
  enddo ! kd

  call pass_var(temp_z,g%Domain)
  call pass_var(salt_z,g%Domain)
  call pass_var(mask_z,g%Domain)
  call pass_var(rho_z,g%Domain)

  ! This is needed for building an ALE grid under ice shelves
  if (use_ice_shelf) then
    if (.not.file_exists(shelf_file, g%Domain)) call mom_error(fatal, &
      "MOM_temp_salt_initialize_from_Z: Unable to open shelf file "//trim(shelf_file))

    call read_data(shelf_file,trim(area_varname),area_shelf_h,domain=g%Domain%mpp_domain)

    ! initialize frac_shelf_h with zeros (open water everywhere)
    frac_shelf_h(:,:) = 0.0
    ! compute fractional ice shelf coverage of h
    do j=jsd,jed ; do i=isd,ied
      if (g%areaT(i,j) > 0.0) &
        frac_shelf_h(i,j) = area_shelf_h(i,j) / g%areaT(i,j)
    enddo ; enddo
    ! pass to the pointer
    shelf_area => frac_shelf_h

  endif

! Done with horizontal interpolation.
! Now remap to model coordinates
  if (usealeremapping) then
    call cpu_clock_begin(id_clock_ale)
    ! The regridding tools (grid generation) are coded to work on model arrays of the same
    ! vertical shape. We need to re-write the regridding if the model has fewer layers
    ! than the data. -AJA
    if (kd>nz) call mom_error(fatal,"MOM_initialize_state, MOM_temp_salt_initialize_from_Z(): "//&
         "Data has more levels than the model - this has not been coded yet!")
    ! Build the source grid and copy data onto model-shaped arrays with vanished layers
    allocate( tmp_mask_in(isd:ied,jsd:jed,nz) ) ; tmp_mask_in(:,:,:) = 0.
    allocate( h1(isd:ied,jsd:jed,nz) ) ; h1(:,:,:) = 0.
    allocate( tmpt1din(isd:ied,jsd:jed,nz) ) ; tmpt1din(:,:,:) = 0.
    allocate( tmps1din(isd:ied,jsd:jed,nz) ) ; tmps1din(:,:,:) = 0.
    do j = js, je ; do i = is, ie
      if (g%mask2dT(i,j)>0.) then
        ztopofcell = 0. ; zbottomofcell = 0. ; npoints = 0
        tmp_mask_in(i,j,1:kd) = mask_z(i,j,:)
        do k = 1, nz
          if (tmp_mask_in(i,j,k)>0. .and. k<=kd) then
            zbottomofcell = -min( z_edges_in(k+1), g%bathyT(i,j) )
            tmpt1din(i,j,k) = temp_z(i,j,k)
            tmps1din(i,j,k) = salt_z(i,j,k)
          elseif (k>1) then
            zbottomofcell = -g%bathyT(i,j)
            tmpt1din(i,j,k) = tmpt1din(i,j,k-1)
            tmps1din(i,j,k) = tmps1din(i,j,k-1)
          else ! This next block should only ever be reached over land
            tmpt1din(i,j,k) = -99.9
            tmps1din(i,j,k) = -99.9
          endif
          h1(i,j,k) = ztopofcell - zbottomofcell
          if (h1(i,j,k)>0.) npoints = npoints + 1
          ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
        enddo
        h1(i,j,kd) = h1(i,j,kd) + ( ztopofcell + g%bathyT(i,j) ) ! In case data is deeper than model
      endif ! mask2dT
    enddo ; enddo
    deallocate( tmp_mask_in )
    call pass_var(h1, g%Domain)
    call pass_var(tmpt1din, g%Domain)
    call pass_var(tmps1din, g%Domain)

    ! Build the target grid (and set the model thickness to it)
    ! This call can be more general but is hard-coded for z* coordinates...  ????
    call ale_initregridding( gv, g%max_depth, pf, mdl, regridcs ) ! sets regridCS

    if (.not. remap_general) then
      ! This is the old way of initializing to z* coordinates only
      allocate( htarget(nz) )
      htarget = getcoordinateresolution( regridcs )
      do j = js, je ; do i = is, ie
        h(i,j,:) = 0.
        if (g%mask2dT(i,j)>0.) then
          ! Build the target grid combining hTarget and topography
          ztopofcell = 0. ; zbottomofcell = 0.
          do k = 1, nz
            zbottomofcell = max( ztopofcell - htarget(k), -g%bathyT(i,j) )
            h(i,j,k) = ztopofcell - zbottomofcell
            ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
          enddo
        else
          h(i,j,:) = 0.
        endif ! mask2dT
      enddo ; enddo
      call pass_var(h, g%Domain)
      deallocate( htarget )
    endif

    ! Now remap from source grid to target grid
    call initialize_remapping( remapcs, remappingscheme, boundary_extrapolation=.false. ) ! Reconstruction parameters
    if (remap_general) then
      call set_regrid_params( regridcs, min_thickness=0. )
      h(:,:,:) = h1(:,:,:) ; tv%T(:,:,:) = tmpt1din(:,:,:) ; tv%S(:,:,:) = tmps1din(:,:,:)
      do j = js, je ; do i = is, ie
        if (g%mask2dT(i,j)==0.) then ! Ensure there are no nonsense values on land
          h(i,j,:) = 0. ; tv%T(i,j,:) = 0. ; tv%S(i,j,:) = 0.
        endif
      enddo ; enddo
      call pass_var(h, g%Domain)    ! Regridding might eventually use spatial information and
      call pass_var(tv%T, g%Domain) ! thus needs to be up to date in the halo regions even though
      call pass_var(tv%S, g%Domain) ! ALE_build_grid() only updates h on the computational domain.

      if (use_ice_shelf) then
         call ale_build_grid( g, gv, regridcs, remapcs, h, tv, .true., shelf_area)
      else
         call ale_build_grid( g, gv, regridcs, remapcs, h, tv, .true. )
      endif
    endif
    call ale_remap_scalar( remapcs, g, gv, nz, h1, tmpt1din, h, tv%T, all_cells=remap_full_column, old_remap=remap_old_alg )
    call ale_remap_scalar( remapcs, g, gv, nz, h1, tmps1din, h, tv%S, all_cells=remap_full_column, old_remap=remap_old_alg )
    deallocate( h1 )
    deallocate( tmpt1din )
    deallocate( tmps1din )

    call cpu_clock_end(id_clock_ale)

  else ! remap to isopycnal layer space

! next find interface positions using local arrays
! nlevs contains the number of valid data points in each column
    nlevs = sum(mask_z,dim=3)

! Rb contains the layer interface densities
    allocate(rb(nz+1))
    do k=2,nz ; rb(k)=0.5*(gv%Rlay(k-1)+gv%Rlay(k)) ; enddo
    rb(1) = 0.0 ;  rb(nz+1) = 2.0*gv%Rlay(nz) - gv%Rlay(nz-1)

    zi(is:ie,js:je,:) = find_interfaces(rho_z(is:ie,js:je,:), z_in, rb, g%bathyT(is:ie,js:je), &
                         nlevs(is:ie,js:je), nkml, nkbl, min_depth)

    if (correct_thickness) then
      call adjustetatofitbathymetry(g, gv, zi, h)
    else
      do k=nz,1,-1 ; do j=js,je ; do i=is,ie
        if (zi(i,j,k) < (zi(i,j,k+1) + gv%Angstrom_z)) then
          zi(i,j,k) = zi(i,j,k+1) + gv%Angstrom_z
          h(i,j,k) = gv%Angstrom_z
        else
          h(i,j,k) = zi(i,j,k) - zi(i,j,k+1)
        endif
      enddo ; enddo ; enddo
      inconsistent=0
      do j=js,je ; do i=is,ie
        if (abs(zi(i,j,nz+1) + g%bathyT(i,j)) > 1.0) &
          inconsistent = inconsistent + 1
      enddo ; enddo
      call sum_across_pes(inconsistent)

      if ((inconsistent > 0) .and. (is_root_pe())) then
        write(mesg,'("Thickness initial conditions are inconsistent ",'// &
                 '"with topography in ",I5," places.")') inconsistent
        call mom_error(warning, mesg)
      endif
    endif

    tv%T(is:ie,js:je,:) = tracer_z_init(temp_z(is:ie,js:je,:),-1.0*z_edges_in,zi(is:ie,js:je,:), &
                                        nkml,nkbl,missing_value,g%mask2dT(is:ie,js:je),nz, &
                                        nlevs(is:ie,js:je),dbg,idbg,jdbg)
    tv%S(is:ie,js:je,:) = tracer_z_init(salt_z(is:ie,js:je,:),-1.0*z_edges_in,zi(is:ie,js:je,:), &
                                        nkml,nkbl,missing_value,g%mask2dT(is:ie,js:je),nz, &
                                        nlevs(is:ie,js:je))

    do k=1,nz
      npoints = 0 ; tempavg = 0. ; saltavg = 0.
      do j=js,je ; do i=is,ie ; if (g%mask2dT(i,j) >= 1.0) then
        npoints = npoints + 1
        tempavg = tempavg + tv%T(i,j,k)
        saltavg =saltavg + tv%S(i,j,k)
      endif ; enddo ; enddo

   ! Horizontally homogenize data to produce perfectly "flat" initial conditions
      if (homogenize) then
        call sum_across_pes(npoints)
        call sum_across_pes(tempavg)
        call sum_across_pes(saltavg)
        if (npoints>0) then
          tempavg = tempavg/real(npoints)
          saltavg = saltavg/real(npoints)
        endif
        tv%T(:,:,k) = tempavg
        tv%S(:,:,k) = saltavg
      endif
    enddo

  endif ! useALEremapping

! Fill land values
  do k=1,nz ; do j=js,je ; do i=is,ie
    if (tv%T(i,j,k) == missing_value) then
      tv%T(i,j,k)=temp_land_fill
      tv%S(i,j,k)=salt_land_fill
    endif
  enddo ; enddo ; enddo

! Finally adjust to target density
  ks=max(0,nkml)+max(0,nkbl)+1

  if (adjust_temperature .and. .not. usealeremapping) then
    call determine_temperature(tv%T(is:ie,js:je,:), tv%S(is:ie,js:je,:), &
            gv%Rlay(1:nz), tv%p_ref, niter, missing_value, h(is:ie,js:je,:), ks, eos)

  endif

  deallocate(z_in,z_edges_in,temp_z,salt_z,mask_z)

  call pass_var(h, g%Domain)
  call pass_var(tv%T, g%Domain)
  call pass_var(tv%S, g%Domain)

  call calltree_leave(trim(mdl)//'()')
  call cpu_clock_end(id_clock_routine)

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

Parameters

[in,out]

g

The ocean's grid structure

Definition at line 1759 of file MOM_state_initialization.F90.

  type(ocean_grid_type), intent(inout) :: g    !< The ocean's grid structure
  ! This subroutine sets the 4 bottom depths at velocity points to be the
  ! maximum of the adjacent depths.
  integer :: i, j

  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
    g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
    g%Dopen_u(i,j) = g%Dblock_u(i,j)
  enddo ; enddo
  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
    g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
    g%Dopen_v(i,j) = g%Dblock_v(i,j)
  enddo ; enddo

set_velocity_depth_min()

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

private

Parameters

[in,out]

g

The ocean's grid structure

Definition at line 1794 of file MOM_state_initialization.F90.

  type(ocean_grid_type), intent(inout) :: g    !< The ocean's grid structure
  ! This subroutine sets the 4 bottom depths at velocity points to be the
  ! minimum of the adjacent depths.
  integer :: i, j

  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
    g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
    g%Dopen_u(i,j) = g%Dblock_u(i,j)
  enddo ; enddo
  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
    g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
    g%Dopen_v(i,j) = g%Dblock_v(i,j)
  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(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]	pf	Parameter file structure
[in]	g	Ocean grid structure
[in]	gv	Vertical grid structure
	ale_csp	ALE control structure
[in,out]	tv	Thermodynamics structure
[in,out]	h	Layer thickness (H units, m or Pa)
[in]	just_read_params	If present and true, this call will only read parameters without changing h.

Definition at line 983 of file MOM_state_initialization.F90.

References cut_off_column_top().

Referenced by mom_initialize_state().

  type(param_file_type),   intent(in)    :: pf !< Parameter file structure
  type(ocean_grid_type),   intent(in)    :: g !< Ocean grid structure
  type(verticalgrid_type), intent(in)    :: gv !< Vertical grid structure
  type(ale_cs),            pointer       :: ale_csp !< ALE control structure
  type(thermo_var_ptrs),   intent(inout) :: tv !< Thermodynamics structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(inout) :: h !< Layer thickness (H units, m or Pa)
  logical,       optional, intent(in)    :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

  ! Local variables
  character(len=200) :: mdl = "trim_for_ice"
  real, dimension(SZI_(G),SZJ_(G)) :: p_surf ! Imposed pressure on ocean at surface (Pa)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: s_t, s_b, t_t, t_b ! Top and bottom edge values for reconstructions
                                                                 ! of salinity and temperature within each layer.
  character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path
  real :: scale_factor, min_thickness
  integer :: i, j, k
  logical :: just_read    ! If true, just read parameters but set nothing.
  logical :: use_remapping
  type(remapping_cs), pointer :: remap_cs => null()

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  call get_param(pf, mdl, "SURFACE_PRESSURE_FILE", p_surf_file, &
                 "The initial condition file for the surface height.", &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(pf, mdl, "SURFACE_PRESSURE_VAR", p_surf_var, &
                 "The initial condition variable for the surface height.", &
                 units="kg m-2", default="", do_not_log=just_read)
  call get_param(pf, mdl, "INPUTDIR", inputdir, default=".", do_not_log=.true.)
  filename = trim(slasher(inputdir))//trim(p_surf_file)
  if (.not.just_read) call log_param(pf,  mdl, "!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)

  call get_param(pf, mdl, "SURFACE_PRESSURE_SCALE", scale_factor, &
                 "A scaling factor to convert SURFACE_PRESSURE_VAR from\n"//&
                 "file SURFACE_PRESSURE_FILE into a surface pressure.", &
                 units="file dependent", default=1., do_not_log=just_read)
  call get_param(pf, mdl, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', &
                 units='m', default=1.e-3, do_not_log=just_read)
  call get_param(pf, mdl, "TRIMMING_USES_REMAPPING", use_remapping, &
                 'When trimming the column, also remap T and S.', &
                 default=.false., do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  call read_data(filename, p_surf_var, p_surf, domain=g%Domain%mpp_domain)
  if (scale_factor /= 1.) p_surf(:,:) = scale_factor * p_surf(:,:)

  if (use_remapping) then
    allocate(remap_cs)
    call initialize_remapping(remap_cs, 'PLM', boundary_extrapolation=.true.)
  endif

  ! Find edge values of T and S used in reconstructions
  if ( associated(ale_csp) ) then ! This should only be associated if we are in ALE mode
!   if ( PRScheme == PRESSURE_RECONSTRUCTION_PLM ) then
      call pressure_gradient_plm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h)
!   elseif ( PRScheme == PRESSURE_RECONSTRUCTION_PPM ) then
!     call pressure_gradient_ppm(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h)
!   endif
  else
!    call MOM_error(FATAL, "trim_for_ice: Does not work without ALE mode")
    do k=1,g%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
      t_t(i,j,k) = tv%T(i,j,k) ; t_b(i,j,k) = tv%T(i,j,k)
      s_t(i,j,k) = tv%S(i,j,k) ; s_b(i,j,k) = tv%S(i,j,k)
    enddo ; enddo ; enddo
  endif

  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    call cut_off_column_top(gv%ke, tv, gv%Rho0, gv%g_Earth, g%bathyT(i,j), min_thickness, &
               tv%T(i,j,:), t_t(i,j,:), t_b(i,j,:), tv%S(i,j,:), s_t(i,j,:), s_b(i,j,:), &
               p_surf(i,j), h(i,j,:), remap_cs)
  enddo ; enddo

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

Definition at line 98 of file MOM_state_initialization.F90.

Referenced by mom_initialize_state().

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