dome_initialization Module Reference

Detailed Description

The module configures the model for the "DOME" experiment. DOME = Dynamics of Overflows and Mixing Experiment.

Functions/Subroutines
subroutine, public	dome_initialize_topography (D, G, param_file, max_depth)
	This subroutine sets up the DOME topography. More...
subroutine, public	dome_initialize_thickness (h, G, GV, param_file, just_read_params)
	This subroutine initializes layer thicknesses for the DOME experiment. More...
subroutine, public	dome_initialize_sponges (G, GV, tv, PF, 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, public	dome_set_obc_data (OBC, tv, G, GV, param_file, tr_Reg)
	This subroutine sets the properties of flow at open boundary conditions. This particular example is for the DOME inflow describe in Legg et al. 2006. More...

Function/Subroutine Documentation

dome_initialize_sponges()

subroutine, public dome_initialization::dome_initialize_sponges	(	type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(thermo_var_ptrs), intent(in)	tv,
		type(param_file_type), intent(in)	PF,
		type(sponge_cs), pointer	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. !

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[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]	pf	A structure indicating the open file to parse for model parameter values.
	csp	A pointer that is set to point to the control structure for this module.

Definition at line 152 of file DOME_initialization.F90.

References mom_sponge::initialize_sponge(), mom_error_handler::mom_error(), and mom_sponge::set_up_sponge_field().

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.
  type(thermo_var_ptrs), intent(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.
  type(param_file_type), intent(in) :: pf   !< A structure indicating the open file to
                                            !! parse for model parameter values.
  type(sponge_cs),       pointer    :: csp  !< A pointer that is set to point to the control
                                            !! structure for this module.

  real :: eta(szi_(g),szj_(g),szk_(g)+1) ! A temporary array for eta.
  real :: temp(szi_(g),szj_(g),szk_(g))  ! A temporary array for other variables. !
  real :: idamp(szi_(g),szj_(g))    ! The inverse damping rate, in s-1.

  real :: h0(szk_(g))
  real :: min_depth
  real :: damp, e_dense, damp_new
  character(len=40)  :: mdl = "DOME_initialize_sponges" ! This subroutine's name.
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz

  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

  eta(:,:,:) = 0.0 ; temp(:,:,:) = 0.0 ; idamp(:,:) = 0.0

!  Here the inverse damping time, in s-1, is set. Set Idamp to 0     !
!  wherever there is no sponge, and the subroutines that are called  !
!  will automatically set up the sponges only where Idamp is positive!
!  and mask2dT is 1.                                                   !

!   Set up sponges for DOME configuration
  call get_param(pf, mdl, "MINIMUM_DEPTH", min_depth, &
                 "The minimum depth of the ocean.", units="m", default=0.0)

  h0(1) = 0.0
  do k=2,nz ; h0(k) = -(real(k-1)-0.5)*g%max_depth/real(nz-1) ; enddo
  do i=is,ie; do j=js,je
    if (g%geoLonT(i,j) < 100.0) then ; damp = 10.0
    elseif (g%geoLonT(i,j) < 200.0) then
      damp = 10.0*(200.0-g%geoLonT(i,j))/100.0
    else ; damp=0.0
    endif

    if (g%geoLonT(i,j) > 1400.0) then ; damp_new = 10.0
    elseif (g%geoLonT(i,j) > 1300.0) then
       damp_new = 10.0*(g%geoLonT(i,j)-1300.0)/100.0
    else ; damp_new = 0.0
    endif

    if (damp <= damp_new) damp=damp_new

    ! These will be stretched inside of apply_sponge, so they can be in
    ! depth space for Boussinesq or non-Boussinesq models.
    eta(i,j,1) = 0.0
    do k=2,nz
!     eta(i,j,K)=max(H0(k), -G%bathyT(i,j), GV%Angstrom_z*(nz-k+1)-G%bathyT(i,j))
      e_dense = -g%bathyT(i,j)
      if (e_dense >= h0(k)) then ; eta(i,j,k) = e_dense
      else ; eta(i,j,k) = h0(k) ; endif
      if (eta(i,j,k) < gv%Angstrom_z*(nz-k+1)-g%bathyT(i,j)) &
          eta(i,j,k) = gv%Angstrom_z*(nz-k+1)-g%bathyT(i,j)
    enddo
    eta(i,j,nz+1) = -g%bathyT(i,j)

    if (g%bathyT(i,j) > min_depth) then
      idamp(i,j) = damp/86400.0
    else ; idamp(i,j) = 0.0 ; endif
  enddo ; enddo

!  This call sets up the damping rates and interface heights.
!  This sets the inverse damping timescale fields in the sponges.    !
  call initialize_sponge(idamp, eta, g, pf, csp)

!   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.                !

! At this point, the DOME configuration is done. The following are here as a
! template for other configurations.

!  The remaining calls to set_up_sponge_field can be in any order. !
  if ( associated(tv%T) ) then
    call mom_error(fatal,"DOME_initialize_sponges is not set up for use with"//&
                         " a temperatures defined.")
    ! This should use the target values of T in temp.
    call set_up_sponge_field(temp, tv%T, g, nz, csp)
    ! This should use the target values of S in temp.
    call set_up_sponge_field(temp, tv%S, g, nz, csp)
  endif

Here is the call graph for this function:

Here is the caller graph for this function:

dome_initialize_thickness()

subroutine, public dome_initialization::dome_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,
		logical, intent(in), optional	just_read_params
	)

This subroutine initializes layer thicknesses for the DOME experiment.

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 94 of file DOME_initialization.F90.

References 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.
  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 :: eta1d(szk_(gv)+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=40)  :: mdl = "DOME_initialize_thickness" ! This subroutine's name.
  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 mom_mesg("  DOME_initialization.F90, DOME_initialize_thickness: setting thickness", 5)

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

  do j=g%jsc,g%jec ; do i=g%isc,g%iec
!    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

Here is the call graph for this function:

dome_initialize_topography()

subroutine, public dome_initialization::dome_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 DOME 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 49 of file DOME_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

  real :: min_depth ! The minimum and maximum depths in m.
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "DOME_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("  DOME_initialization.F90, DOME_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)

  do j=js,je ; do i=is,ie
    if (g%geoLatT(i,j) < 600.0) then
      if (g%geoLatT(i,j) < 300.0) then
        d(i,j)=max_depth
      else
        d(i,j)=max_depth-10.0*(g%geoLatT(i,j)-300.0)
      endif
    else
      if ((g%geoLonT(i,j) > 1000.0).AND.(g%geoLonT(i,j) < 1100.0)) then
        d(i,j)=600.0
      else
        d(i,j)=0.5*min_depth
      endif
    endif

    if (d(i,j) > max_depth) d(i,j) = max_depth
    if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
  enddo ; enddo

Here is the call graph for this function:

Here is the caller graph for this function:

dome_set_obc_data()

subroutine, public dome_initialization::dome_set_obc_data	(	type(ocean_obc_type), pointer	OBC,
		type(thermo_var_ptrs), intent(in)	tv,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(param_file_type), intent(in)	param_file,
		type(tracer_registry_type), pointer	tr_Reg
	)

This subroutine sets the properties of flow at open boundary conditions. This particular example is for the DOME inflow describe in Legg et al. 2006.

Parameters

	obc	This open boundary condition type specifies whether, where, and what open boundary conditions are used.
[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]	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.
	tr_reg	Tracer registry.

Definition at line 247 of file DOME_initialization.F90.

References mom_tracer_registry::add_tracer_obc_values(), and mom_eos::calculate_density_derivs().

  type(ocean_obc_type),       pointer    :: obc !< This open boundary condition type specifies
                                                !! whether, where, and what open boundary
                                                !! conditions are used.
  type(thermo_var_ptrs),      intent(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.
  type(ocean_grid_type),      intent(in) :: g   !< The ocean's grid structure.
  type(verticalgrid_type),    intent(in) :: gv  !< The ocean's vertical grid structure.
  type(param_file_type),      intent(in) :: param_file !< A structure indicating the open file
                              !! to parse for model parameter values.
  type(tracer_registry_type), pointer    :: tr_reg !< Tracer registry.

  real, pointer, dimension(:,:,:) :: &
    obc_t_v => null(), &    ! specify the values of T and S that should come
    obc_s_v => null()       ! boundary conditions, in C and psu.
  ! The following variables are used to set the target temperature and salinity.
  real :: t0(szk_(g)), s0(szk_(g))
  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.
  ! The following variables are used to set up the transport in the DOME example.
  real :: tr_0, y1, y2, tr_k, rst, rsb, rc, v_k, lon_im1
  real :: d_edge            ! The thickness in m of the dense fluid at the
                            ! inner edge of the inflow.
  real :: g_prime_tot       ! The reduced gravity across all layers, m s-2.
  real :: def_rad           ! The deformation radius, based on fluid of
                            ! thickness D_edge, in the same units as lat.
  real :: ri_trans          ! The shear Richardson number in the transition
                            ! region of the specified shear profile.
  character(len=40)  :: mdl = "DOME_set_OBC_data" ! This subroutine's name.
  integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz
  integer :: isdb, iedb, jsdb, jedb
  type(obc_segment_type), pointer :: segment

  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
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  ! The following variables should be transformed into runtime parameters.
  d_edge = 300.0  ! The thickness of dense fluid in the inflow.
  ri_trans = 1.0/3.0 ! The shear Richardson number in the transition region
                     ! region of the specified shear profile.

  if (.not.associated(obc)) return

  g_prime_tot = (gv%g_Earth/gv%Rho0)*2.0
  def_rad = sqrt(d_edge*g_prime_tot) / (1.0e-4*1000.0)
  tr_0 = (-d_edge*sqrt(d_edge*g_prime_tot)*0.5e3*def_rad) * gv%m_to_H

  if (obc%number_of_segments .ne. 1) then
    print *, 'Error in DOME OBC segment setup'
    return   !!! Need a better error message here
  endif
  segment => obc%segment(1)
  if (.not. segment%on_pe) return

  do k=1,nz
    rst = -1.0
    if (k>1) rst = -1.0 + (real(k-1)-0.5)/real(nz-1)

    rsb = 0.0
    if (k<nz) rsb = -1.0 + (real(k-1)+0.5)/real(nz-1)
    rc = -1.0 + real(k-1)/real(nz-1)

    ! These come from assuming geostrophy and a constant Ri profile.
    y1 = (2.0*ri_trans*rst + ri_trans + 2.0)/(2.0 - ri_trans)
    y2 = (2.0*ri_trans*rsb + ri_trans + 2.0)/(2.0 - ri_trans)
    tr_k = tr_0 * (2.0/(ri_trans*(2.0-ri_trans))) * &
           ((log(y1)+1.0)/y1 - (log(y2)+1.0)/y2)
    v_k = -sqrt(d_edge*g_prime_tot)*log((2.0 + ri_trans*(1.0 + 2.0*rc)) / &
                                        (2.0 - ri_trans))
    if (k == nz)  tr_k = tr_k + tr_0 * (2.0/(ri_trans*(2.0+ri_trans))) * &
                                       log((2.0+ri_trans)/(2.0-ri_trans))
    ! New way
    isd = segment%HI%isd ; ied = segment%HI%ied
    jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
    do j=jsdb,jedb ; do i=isd,ied
      lon_im1 = 2.0*g%geoLonCv(i,j) - g%geoLonBu(i,j)
      segment%normal_trans(i,j,k) = tr_k * (exp(-2.0*(lon_im1 - 1000.0)/def_rad) -&
                                  exp(-2.0*(g%geoLonBu(i,j) - 1000.0)/def_rad))
      segment%normal_vel(i,j,k) = v_k * exp(-2.0*(g%geoLonCv(i,j) - 1000.0)/def_rad)
    enddo ; enddo
  enddo

  !   The inflow values of temperature and salinity also need to be set here if
  ! these variables are used.  The following code is just a naive example.
  if (associated(tv%S)) then
    ! In this example, all S inflows have values of 35 psu.
    call add_tracer_obc_values("S", tr_reg, obc_inflow=35.0)
  endif
  if (associated(tv%T)) then
    ! In this example, the T values are set to be consistent with the layer
    ! target density and a salinity of 35 psu.  This code is taken from
    ! USER_initialize_temp_sal.
    pres(:) = tv%P_Ref ; s0(:) = 35.0 ; t0(1) = 25.0
    call calculate_density(t0(1),s0(1),pres(1),rho_guess(1),tv%eqn_of_state)
    call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,1,tv%eqn_of_state)

    do k=1,nz ; 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,tv%eqn_of_state)
      call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,tv%eqn_of_state)
      do k=1,nz ; t0(k) = t0(k) + (gv%Rlay(k)-rho_guess(k)) / drho_dt(k) ; enddo
    enddo

    allocate(obc_t_v(isd:ied,jsdb:jedb,nz))
    do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
      obc_t_v(i,j,k) = t0(k)
    enddo ; enddo ; enddo
    call add_tracer_obc_values("T", tr_reg, obc_in_v=obc_t_v)
  endif

Here is the call graph for this function: