benchmark_initialization Module Reference

Functions/Subroutines
subroutine, public	benchmark_initialize_topography (D, G, param_file, max_depth)
	This subroutine sets up the benchmark test case topography. More...
subroutine, public	benchmark_initialize_thickness (h, G, GV, param_file, eqn_of_state, P_ref, just_read_params)
	This subroutine initializes layer thicknesses for the benchmark test case, by finding the depths of interfaces in a specified latitude-dependent temperature profile with an exponentially decaying thermocline on top of a linear stratification. More...
subroutine, public	benchmark_init_temperature_salinity (T, S, G, GV, param_file, eqn_of_state, P_Ref, just_read_params)
	This function puts the initial layer temperatures and salinities into T(:,:,:) and S(:,:,:). More...

Function/Subroutine Documentation

benchmark_init_temperature_salinity()

subroutine, public benchmark_initialization::benchmark_init_temperature_salinity	(	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
	)

This function puts the initial layer temperatures and salinities into T(:,:,:) and S(:,:,:).

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[out]	t	The potential temperature that is being initialized.
[out]	s	The salinity that is being initialized.
[in]	param_file	A structure indicating the open file to parse for model parameter values.
	eqn_of_state	integer that selects the equation 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 220 of file benchmark_initialization.F90.

References mom_eos::calculate_density_derivs().

Referenced by mom_state_initialization::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 indicating the
                                                                   !! open file to parse for
                                                                   !! model parameter values.
  type(eos_type),                      pointer     :: eqn_of_state !< integer that selects the
                                                                   !! equation 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.

  real :: t0(szk_(g)), s0(szk_(g))
  real :: pres(szk_(g))      ! Reference pressure in kg m-3.             !
  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.   !
  real :: pi        ! 3.1415926... calculated as 4*atan(1)
  real :: sst       !  The initial sea surface temperature, in deg C.
  real :: lat
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=40)  :: mdl = "benchmark_init_temperature_salinity" ! This subroutine's name.
  integer :: i, j, k, k1, is, ie, js, je, nz, itt

  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 ! All run-time parameters have been read, so return.

  k1 = gv%nk_rho_varies + 1

  do k=1,nz
    pres(k) = p_ref ; s0(k) = 35.0
  enddo

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

! A first guess of the layers' temperatures.                         !
  do k=1,nz
    t0(k) = t0(k1) + (gv%Rlay(k) - rho_guess(k1)) / drho_dt(k1)
  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
      t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
    enddo
  enddo

  do k=1,nz ; do i=is,ie ; do j=js,je
    t(i,j,k) = t0(k)
    s(i,j,k) = s0(k)
  enddo ; enddo ; enddo
  pi = 4.0*atan(1.0)
  do i=is,ie ; do j=js,je
    sst = 0.5*(t0(k1)+t0(nz)) - 0.9*0.5*(t0(k1)-t0(nz)) * &
                               cos(pi*(g%geoLatT(i,j)-g%south_lat)/(g%len_lat))
    do k=1,k1-1
      t(i,j,k) = sst
    enddo
  enddo ; enddo

Here is the call graph for this function:

Here is the caller graph for this function:

benchmark_initialize_thickness()

subroutine, public benchmark_initialization::benchmark_initialize_thickness	(	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,
		type(eos_type), pointer	eqn_of_state,
		real, intent(in)	P_ref,
		logical, intent(in), optional	just_read_params
	)

This subroutine initializes layer thicknesses for the benchmark test case, by finding the depths of interfaces in a specified latitude-dependent temperature profile with an exponentially decaying thermocline on top of a linear stratification.

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.
	eqn_of_state	integer that selects the equation 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 96 of file benchmark_initialization.F90.

References mom_eos::calculate_density_derivs(), and mom_error_handler::mom_mesg().

  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.
  type(eos_type),          pointer     :: eqn_of_state !< integer that selects the
                                                      !! equation 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.

  real :: e0(szk_(gv)+1)     ! The resting interface heights, in m, usually !
                             ! negative because it is positive upward.      !
  real :: e_pert(szk_(gv)+1) ! Interface height perturbations, positive     !
                             ! upward, in m.                                !
  real :: eta1d(szk_(gv)+1)  ! Interface height relative to the sea surface !
                             ! positive upward, in m.                       !
  real :: sst       !  The initial sea surface temperature, in deg C.
  real :: t_int     !  The initial temperature of an interface, in deg C.
  real :: ml_depth  !  The specified initial mixed layer depth, in m.
  real :: thermocline_scale ! The e-folding scale of the thermocline, in m.
  real, dimension(SZK_(GV)) :: t0, pres, s0, rho_guess, drho, drho_dt, drho_ds
  real :: a_exp      ! The fraction of the overall stratification that is exponential.
  real :: i_ts, i_md ! Inverse lengthscales in m-1.
  real :: t_frac     ! A ratio of the interface temperature to the range
                     ! between SST and the bottom temperature.
  real :: err, derr_dz  ! The error between the profile's temperature and the
                     ! interface temperature for a given z and its derivative.
  real :: pi, z
  logical :: just_read
  character(len=40)  :: mdl = "benchmark_initialize_thickness" ! This subroutine's name.
  integer :: i, j, k, k1, is, ie, js, je, nz, itt

  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 mom_mesg("  benchmark_initialization.F90, benchmark_initialize_thickness: setting thickness", 5)

  k1 = gv%nk_rho_varies + 1

  ml_depth = 50.0
  thermocline_scale = 500.0
  a_exp = 0.9

! This block calculates T0(k) for the purpose of diagnosing where the
! interfaces will be found.
  do k=1,nz
    pres(k) = p_ref ; s0(k) = 35.0
  enddo
  t0(k1) = 29.0
  call calculate_density(t0(k1),s0(k1),pres(k1),rho_guess(k1),eqn_of_state)
  call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,k1,1,eqn_of_state)

! A first guess of the layers' temperatures.
  do k=1,nz
    t0(k) = t0(k1) + (gv%Rlay(k) - rho_guess(k1)) / drho_dt(k1)
  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
      t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
    enddo
  enddo

  pi = 4.0*atan(1.0)
  i_ts = 1.0 / thermocline_scale
  i_md = 1.0 / g%max_depth
  do j=js,je ; do i=is,ie
    sst = 0.5*(t0(k1)+t0(nz)) - 0.9*0.5*(t0(k1)-t0(nz)) * &
                               cos(pi*(g%geoLatT(i,j)-g%south_lat)/(g%len_lat))

    do k=1,nz ; e_pert(k) = 0.0 ; enddo

!  The remainder of this subroutine should not be changed.           !

!    This sets the initial thickness (in m) of the layers.  The      !
!  thicknesses are set to insure that: 1.  each layer is at least    !
!  Gv%Angstrom_z 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,2,-1
      t_int = 0.5*(t0(k) + t0(k-1))
      t_frac = (t_int - t0(nz)) / (sst - t0(nz))
      ! Find the z such that T_frac = a exp(z/thermocline_scale) + (1-a) (z+D)/D
      z = 0.0
      do itt=1,6
        err = a_exp * exp(z*i_ts) + (1.0 - a_exp) * (z*i_md + 1.0) - t_frac
        derr_dz = a_exp * i_ts * exp(z*i_ts) + (1.0 - a_exp) * i_md
        z = z - err / derr_dz
      enddo
      e0(k) = z
!       e0(K) = -ML_depth + thermocline_scale * log((T_int - T0(nz)) / (SST - T0(nz)))

      eta1d(k) = e0(k) + e_pert(k)

      if (eta1d(k) > -ml_depth) eta1d(k) = -ml_depth

      if (eta1d(k) < eta1d(k+1) + gv%Angstrom_z) &
        eta1d(k) = eta1d(k+1) + gv%Angstrom_z

      h(i,j,k) = max(eta1d(k) - eta1d(k+1), gv%Angstrom_z)
    enddo
    h(i,j,1) = max(0.0 - eta1d(2), gv%Angstrom_z)

  enddo ; enddo

Here is the call graph for this function:

benchmark_initialize_topography()

subroutine, public benchmark_initialization::benchmark_initialize_topography	(	real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(out)	D,
		type(dyn_horgrid_type), intent(in)	G,
		type(param_file_type), intent(in)	param_file,
		real, intent(in)	max_depth
	)

This subroutine sets up the benchmark test case topography.

Parameters

[in]	g	The dynamic horizontal grid type
[out]	d	Ocean bottom depth in m
[in]	param_file	Parameter file structure
[in]	max_depth	Maximum depth of model in m

Definition at line 46 of file benchmark_initialization.F90.

References mom_error_handler::mom_mesg().

Referenced by mom_fixed_initialization::mom_initialize_topography().

  type(dyn_horgrid_type),             intent(in)  :: g !< The dynamic horizontal grid type
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                                      intent(out) :: d !< Ocean bottom depth in m
  type(param_file_type),              intent(in)  :: param_file !< Parameter file structure
  real,                               intent(in)  :: max_depth  !< Maximum depth of model in m

! This subroutine sets up the benchmark test case topography
  real :: min_depth            ! The minimum and maximum depths in m.
  real :: pi                   ! 3.1415926... calculated as 4*atan(1)
  real :: d0                   ! A constant to make the maximum     !
                               ! basin depth MAXIMUM_DEPTH.         !
  real :: x, y
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "benchmark_initialize_topography" ! This subroutine's name.
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  call mom_mesg("  benchmark_initialization.F90, benchmark_initialize_topography: setting topography", 5)

  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
                 "The minimum depth of the ocean.", units="m", default=0.0)

  pi = 4.0*atan(1.0)
  d0 = max_depth / 0.5;

!  Calculate the depth of the bottom.
  do i=is,ie ; do j=js,je
    x=(g%geoLonT(i,j)-g%west_lon)/g%len_lon
    y=(g%geoLatT(i,j)-g%south_lat)/g%len_lat
!  This sets topography that has a reentrant channel to the south.
    d(i,j) = -d0 * ( y*(1.0 + 0.6*cos(4.0*pi*x)) &
                   + 0.75*exp(-6.0*y) &
                   + 0.05*cos(10.0*pi*x) - 0.7 )
    if (d(i,j) > max_depth) d(i,j) = max_depth
    if (d(i,j) < min_depth) d(i,j) = 0.
  enddo ; enddo

Here is the call graph for this function:

Here is the caller graph for this function: