dome_initialization Module Reference

Detailed Description

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, US)
	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, US, 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, US, 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(unit_scale_type), intent(in)	US,
		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]	us	A dimensional unit scaling type
[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 148 of file DOME_initialization.F90.

  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(unit_scale_type),   intent(in) :: US !< A dimensional unit scaling type
  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 [Z ~> m].
  real :: temp(SZI_(G),SZJ_(G),SZK_(G))  ! A temporary array for other variables. !
  real :: Idamp(SZI_(G),SZJ_(G))    ! The inverse damping rate [T-1 ~> s-1].
 
  real :: H0(SZK_(G))  ! Interface heights [Z ~> m].
  real :: min_depth    ! The minimum depth at which to apply damping [Z ~> m]
  real :: damp, damp_new ! Damping rates in the sponge [days]
  real :: e_dense      ! The depth of the densest interfaces [Z ~> m]
  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 [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, scale=us%m_to_Z)
 
  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
    damp = us%T_to_s*damp
 
    ! 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, gv)
 
!   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
 

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 [H ~> m or kg m-2].
[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 90 of file DOME_initialization.F90.

  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 [H ~> m or kg m-2].
  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 [Z ~> m], usually
                            ! negative because it is positive upward [Z ~> m].
  real :: eta1D(SZK_(GV)+1) ! Interface height relative to the sea surface
                            ! positive upward [Z ~> 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) = -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_H
      else
        h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
      endif
    enddo
  enddo ; enddo
 

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,
		type(unit_scale_type), intent(in), optional	US
	)

This subroutine sets up the DOME topography.

Parameters

[in]	g	The dynamic horizontal grid type
[out]	d	Ocean bottom depth in m or Z if US is present
[in]	param_file	Parameter file structure
[in]	max_depth	Maximum model depth in the units of D
[in]	us	A dimensional unit scaling type

Definition at line 40 of file DOME_initialization.F90.

  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 or Z if US is present
  type(param_file_type),           intent(in)  :: param_file !< Parameter file structure
  real,                            intent(in)  :: max_depth !< Maximum model depth in the units of D
  type(unit_scale_type), optional, intent(in)  :: US !< A dimensional unit scaling type
 
  ! Local variables
  real :: m_to_Z  ! A dimensional rescaling factor.
  real :: min_depth ! The minimum and maximum depths [Z ~> 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)
 
  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z
 
  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, scale=m_to_z)
 
  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*m_to_z * (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*m_to_z
      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
 

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(unit_scale_type), intent(in)	US,
		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]	us	A dimensional unit scaling type
[in]	param_file	A structure indicating the open file to parse for model parameter values.
	tr_reg	Tracer registry.

Definition at line 246 of file DOME_initialization.F90.

  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(unit_scale_type),      intent(in) :: US  !< A dimensional unit scaling type
  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.
 
! Local variables
  ! 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 [R L2 T-2 ~> Pa].
  real :: drho_dT(SZK_(G))   ! Derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
  real :: drho_dS(SZK_(G))   ! Derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
  real :: rho_guess(SZK_(G)) ! Potential density at T0 & S0 [R ~> 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 [Z ~> m], of the dense fluid at the
                            ! inner edge of the inflow.
  real :: g_prime_tot       ! The reduced gravity across all layers [L2 Z-1 T-2 ~> m s-2].
  real :: Def_Rad           ! The deformation radius, based on fluid of
                            ! thickness D_edge, in the same units as lat [m].
  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.
  character(len=32)  :: name
  integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, m, nz, NTR
  integer :: IsdB, IedB, JsdB, JedB
  type(OBC_segment_type), pointer :: segment => null()
  type(tracer_type), pointer      :: tr_ptr => null()
 
  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*us%m_to_Z  ! 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*us%kg_m3_to_R
  def_rad = us%L_to_m*sqrt(d_edge*g_prime_tot) / (1.0e-4*us%T_to_s * 1000.0)
  tr_0 = (-d_edge*sqrt(d_edge*g_prime_tot)*0.5e3*us%m_to_L*def_rad) * gv%Z_to_H
 
  if (obc%number_of_segments /= 1) then
    call mom_error(warning, 'Error in DOME OBC segment setup', .true.)
    return   !!! Need a better error message here
  endif
 
  ntr = tr_reg%NTR
 
  ! Stash this information away for the messy tracer restarts.
  obc%ntr = ntr
  if (.not. associated(obc%tracer_x_reservoirs_used)) then
    allocate(obc%tracer_x_reservoirs_used(ntr))
    allocate(obc%tracer_y_reservoirs_used(ntr))
    obc%tracer_x_reservoirs_used(:) = .false.
    obc%tracer_y_reservoirs_used(:) = .false.
    obc%tracer_y_reservoirs_used(1) = .true.
  endif
 
  segment => obc%segment(1)
  if (.not. segment%on_pe) return
 
  allocate(segment%field(ntr))
 
  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.
    name = 'salt'
    call tracer_name_lookup(tr_reg, tr_ptr, name)
    call register_segment_tracer(tr_ptr, param_file, gv, segment, obc_scalar=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, tv%eqn_of_state, (/1,1/) )
 
    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, tv%eqn_of_state)
      call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, tv%eqn_of_state)
      do k=1,nz ; t0(k) = t0(k) + (gv%Rlay(k)-rho_guess(k)) / drho_dt(k) ; enddo
    enddo
 
    ! Temperature on tracer 1???
    allocate(segment%field(1)%buffer_src(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB,nz))
    do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
      segment%field(1)%buffer_src(i,j,k) = t0(k)
    enddo ; enddo ; enddo
    name = 'temp'
    call tracer_name_lookup(tr_reg, tr_ptr, name)
    call register_segment_tracer(tr_ptr, param_file, gv, segment, obc_array=.true.)
  endif
 
  ! Dye tracers - fight with T,S???
  ! First dye - only one with OBC values
  ! This field(1) requires tr_D1 to be the first tracer.
  allocate(segment%field(1)%buffer_src(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB,nz))
  do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed ; do i=segment%HI%isd,segment%HI%ied
    if (k < nz/2) then ; segment%field(1)%buffer_src(i,j,k) = 0.0
    else ; segment%field(1)%buffer_src(i,j,k) = 1.0 ; endif
  enddo ; enddo ; enddo
  name = 'tr_D1'
  call tracer_name_lookup(tr_reg, tr_ptr, name)
  call register_segment_tracer(tr_ptr, param_file, gv, &
                               obc%segment(1), obc_array=.true.)
 
  ! All tracers but the first have 0 concentration in their inflows. As this
  ! is the default value, the following calls are unnecessary.
  do m=2,ntr
    if (m < 10) then ; write(name,'("tr_D",I1.1)') m
    else ; write(name,'("tr_D",I2.2)') m ; endif
    call tracer_name_lookup(tr_reg, tr_ptr, name)
    call register_segment_tracer(tr_ptr, param_file, gv, &
                                 obc%segment(1), obc_scalar=0.0)
  enddo