DOME_initialization.F90

Go to the documentation of this file.

module dome_initialization
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

use mom_sponge, only : sponge_cs, set_up_sponge_field, initialize_sponge
use mom_dyn_horgrid, only : dyn_horgrid_type
use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_get_input, only : directories
use mom_grid, only : ocean_grid_type
use mom_open_boundary, only : ocean_obc_type, obc_none, obc_simple
use mom_open_boundary, only : obc_segment_type
use mom_tracer_registry, only : tracer_registry_type, add_tracer_obc_values
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density, calculate_density_derivs, eos_type

implicit none ; private

#include <MOM_memory.h>

public dome_initialize_topography
public dome_initialize_thickness
public dome_initialize_sponges
public dome_set_obc_data

contains

! -----------------------------------------------------------------------------
!> This subroutine sets up the DOME topography
subroutine dome_initialize_topography(D, G, param_file, max_depth)
  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

end subroutine dome_initialize_topography
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> This subroutine initializes layer thicknesses for the DOME experiment
subroutine dome_initialize_thickness(h, G, GV, param_file, just_read_params)
  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

end subroutine dome_initialize_thickness
! -----------------------------------------------------------------------------

! -----------------------------------------------------------------------------
!> 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.                                        !
subroutine dome_initialize_sponges(G, GV, tv, PF, CSp)
  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

end subroutine dome_initialize_sponges

!> 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.
subroutine dome_set_obc_data(OBC, tv, G, GV, param_file, tr_Reg)
  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

end subroutine dome_set_obc_data

!> \namespace dome_initialization
!!
!! The module configures the model for the "DOME" experiment.
!! DOME = Dynamics of Overflows and Mixing Experiment
end module dome_initialization