RGC_initialization.F90

module rgc_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                      *
!* By Elizabeth Yankovsky, May 2018                                    *
!***********************************************************************

use mom_ale_sponge, only : ale_sponge_cs, set_up_ale_sponge_field, initialize_ale_sponge
use mom_ale_sponge, only : set_up_ale_sponge_vel_field
use mom_domains, only : pass_var
use mom_dyn_horgrid, only : dyn_horgrid_type
use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe, warning
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_io, only : file_exists, read_data
use mom_io, only : slasher
use mom_sponge, only : sponge_cs, set_up_sponge_field, initialize_sponge
use mom_sponge, only : set_up_sponge_ml_density
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density, calculate_density_derivs, eos_type, eos_domain
implicit none ; private

#include <MOM_memory.h>

character(len=40) :: mod = "RGC_initialization" ! This module's name.
public rgc_initialize_sponges

contains

!> Sets up the the inverse restoration time, and the values towards which the interface heights,
!! velocities and tracers should be restored within the sponges for the RGC test case.
subroutine rgc_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp)
  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, potential temperature and
                                            !! salinity or mixed layer density.
                                            !! Absent fields have NULL ptrs.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 target, intent(in) :: u    !< Array with the u velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                 target, intent(in) :: v    !< Array with the v velocity [L T-1 ~> m s-1]
  type(param_file_type), intent(in) :: PF   !< A structure indicating the
                                            !! open file to parse for model
                                            !! parameter values.
  logical, intent(in) :: use_ALE            !< If true, indicates model is in ALE mode
  type(sponge_CS),   pointer    :: CSp      !< Layer-mode sponge structure
  type(ALE_sponge_CS),   pointer    :: ACSp !< ALE-mode sponge structure

! Local variables
  real :: T(SZI_(G),SZJ_(G),SZK_(G))  ! A temporary array for temp
  real :: S(SZI_(G),SZJ_(G),SZK_(G))  ! A temporary array for salt
  real :: U1(SZIB_(G),SZJ_(G),SZK_(G))  ! A temporary array for u [L T-1 ~> m s-1]
  real :: V1(SZI_(G),SZJB_(G),SZK_(G))  ! A temporary array for v [L T-1 ~> m s-1]
  real :: RHO(SZI_(G),SZJ_(G),SZK_(G))  ! A temporary array for RHO
  real :: tmp(SZI_(G),SZJ_(G))        ! A temporary array for tracers.
  real :: h(SZI_(G),SZJ_(G),SZK_(G))  ! A temporary array for thickness at h points
  real :: Idamp(SZI_(G),SZJ_(G))    ! The inverse damping rate at h points [T-1 ~> s-1].
  real :: TNUDG                     ! Nudging time scale [T ~> s]
  real :: pres(SZI_(G))             ! An array of the reference pressure [R L2 T-2 ~> Pa]
  real :: e0(SZK_(G)+1)               ! The resting interface heights, in m, usually !
                                    ! negative because it is positive upward.      !
  real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1) ! A temporary array for eta.
                                    ! positive upward, in m.
  logical :: sponge_uv              ! Nudge velocities (u and v) towards zero
  real :: min_depth, dummy1, z, delta_h
  real :: rho_dummy, min_thickness, rho_tmp, xi0
  real :: lenlat, lenlon, lensponge
  character(len=40) :: filename, state_file
  character(len=40) :: temp_var, salt_var, eta_var, inputdir, h_var

  character(len=40)  :: mod = "RGC_initialize_sponges" ! This subroutine's name.
  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, iscB, iecB, jscB, jecB

  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
  iscb = g%iscB ; iecb = g%iecB; jscb = g%jscB ; jecb = g%jecB

  call get_param(pf,mod,"MIN_THICKNESS",min_thickness,'Minimum layer thickness',units='m',default=1.e-3)

  call get_param(pf, mod, "RGC_TNUDG", tnudg, 'Nudging time scale for sponge layers (days)', &
                 default=0.0, scale=86400.0*us%s_to_T)

  call get_param(pf, mod, "LENLAT", lenlat, &
                  "The latitudinal or y-direction length of the domain", &
                 fail_if_missing=.true., do_not_log=.true.)

  call get_param(pf, mod, "LENLON", lenlon, &
                  "The longitudinal or x-direction length of the domain", &
                 fail_if_missing=.true., do_not_log=.true.)

  call get_param(pf, mod, "LENSPONGE", lensponge, &
                 "The length of the sponge layer (km).", &
                 default=10.0)

  call get_param(pf, mod, "SPONGE_UV", sponge_uv, &
                 "Nudge velocities (u and v) towards zero in the sponge layer.", &
                 default=.false., do_not_log=.true.)

  t(:,:,:) = 0.0 ; s(:,:,:) = 0.0 ; idamp(:,:) = 0.0; rho(:,:,:) = 0.0

  call get_param(pf, mod, "MINIMUM_DEPTH", min_depth, &
                 "The minimum depth of the ocean.", units="m", default=0.0)

  if (associated(csp)) call mom_error(fatal, &
          "RGC_initialize_sponges called with an associated control structure.")
  if (associated(acsp)) call mom_error(fatal, &
          "RGC_initialize_sponges called with an associated ALE-sponge control structure.")

  !  Here the inverse damping time [T-1 ~> 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.

  do i=is,ie ; do j=js,je
    if ((g%bathyT(i,j) <= min_depth) .or. (g%geoLonT(i,j) <= lensponge)) then
      idamp(i,j) = 0.0
    elseif (g%geoLonT(i,j) >= (lenlon - lensponge) .AND. g%geoLonT(i,j) <= lenlon) then
      dummy1 = (g%geoLonT(i,j)-(lenlon - lensponge))/(lensponge)
      idamp(i,j) = (1.0/tnudg) * max(0.0,dummy1)
    else
      idamp(i,j) = 0.0
    endif
  enddo ; enddo


  ! 1) Read eta, salt and temp from IC file
  call get_param(pf, mod, "INPUTDIR", inputdir, default=".")
  inputdir = slasher(inputdir)
   ! GM: get two different files, one with temp and one with salt values
   ! this is work around to avoid having wrong values near the surface
   ! because of the FIT_SALINITY option. To get salt values right in the
   ! sponge, FIT_SALINITY=False. The oposite is true for temp. One can
   ! combined the *correct* temp and salt values in one file instead.
  call get_param(pf, mod, "RGC_SPONGE_FILE", state_file, &
              "The name of the file with temps., salts. and interfaces to \n"// &
              " damp toward.", fail_if_missing=.true.)
  call get_param(pf, mod, "SPONGE_PTEMP_VAR", temp_var, &
              "The name of the potential temperature variable in \n"//&
              "SPONGE_STATE_FILE.", default="Temp")
  call get_param(pf, mod, "SPONGE_SALT_VAR", salt_var, &
              "The name of the salinity variable in \n"//&
              "SPONGE_STATE_FILE.", default="Salt")
  call get_param(pf, mod, "SPONGE_ETA_VAR", eta_var, &
              "The name of the interface height variable in \n"//&
              "SPONGE_STATE_FILE.", default="eta")
  call get_param(pf, mod, "SPONGE_H_VAR", h_var, &
              "The name of the layer thickness variable in \n"//&
              "SPONGE_STATE_FILE.", default="h")

  !read temp and eta
  filename = trim(inputdir)//trim(state_file)
  if (.not.file_exists(filename, g%Domain)) &
      call mom_error(fatal, " RGC_initialize_sponges: Unable to open "//trim(filename))
  call read_data(filename,temp_var,t(:,:,:), domain=g%Domain%mpp_domain)
  call read_data(filename,salt_var,s(:,:,:), domain=g%Domain%mpp_domain)

  if (use_ale) then

    call read_data(filename,h_var,h(:,:,:), domain=g%Domain%mpp_domain)
    call pass_var(h, g%domain)

    !call initialize_ALE_sponge(Idamp, h, nz, G, PF, ACSp)
    call initialize_ale_sponge(idamp, g, pf, acsp, h, nz)

    !  The remaining calls to set_up_sponge_field can be in any order. !
    if ( associated(tv%T) ) then
      call set_up_ale_sponge_field(t,g,tv%T,acsp)
    endif
    if ( associated(tv%S) ) then
      call set_up_ale_sponge_field(s,g,tv%S,acsp)
    endif

    if (sponge_uv) then
      u1(:,:,:) = 0.0; v1(:,:,:) = 0.0
      call set_up_ale_sponge_vel_field(u1,v1,g,u,v,acsp)
    endif


  else ! layer mode

    !read eta
    call read_data(filename,eta_var,eta(:,:,:), domain=g%Domain%mpp_domain)

    ! 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, pf, csp, gv)

    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
      eosdom(:) = eos_domain(g%HI)
      do j=js,je
        call calculate_density(t(:,j,1), s(:,j,1), pres, tmp(:,j), tv%eqn_of_state, eosdom)
      enddo

      call set_up_sponge_ml_density(tmp, g, csp)
    endif

    ! Apply sponge in tracer fields
    call set_up_sponge_field(t, tv%T, g, nz, csp)
    call set_up_sponge_field(s, tv%S, g, nz, csp)

  endif

end subroutine rgc_initialize_sponges

end module rgc_initialization