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_sponge, only : sponge_cs, set_up_sponge_field, initialize_sponge
use mom_sponge, only : set_up_sponge_ml_density
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_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density, calculate_density_derivs, eos_type
use mom_domains, only: pass_var
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 (Idamp), and
! the values towards which the interface heights and an arbitrary
! number of tracers should be restored within each sponge.
subroutine rgc_initialize_sponges(G, GV, 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(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)), intent(in)   :: u !< u velocity.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)   :: v !< v velocity.
  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
  real :: V1(SZI_(G), SZJB_(G), SZK_(G))  ! A temporary array for v
  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, in s-1.
  real :: TNUDG                     ! Nudging time scale, days
  real :: pres(SZI_(G))             ! An array of the reference pressure, in 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 :: damp, 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 :: 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)
 
  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, 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.
 
   do i=is,ie; do j=js,je
      if (g%geoLonT(i,j) <= lensponge) then
        dummy1 = -(g%geoLonT(i,j))/lensponge + 1.0
        !damp = 1.0/TNUDG * max(0.0,dummy1)
        damp = 0.0
        !write(*,*)'1st, G%geoLonT(i,j), damp',G%geoLonT(i,j), damp
 
      elseif (g%geoLonT(i,j) >= (lenlon - lensponge) .AND. g%geoLonT(i,j) <= lenlon) then
 
  ! 1 / day
        dummy1=(g%geoLonT(i,j)-(lenlon - lensponge))/(lensponge)
        damp = (1.0/tnudg) * max(0.0,dummy1)
 
      else ; damp=0.0
      endif
 
  ! convert to 1 / seconds
      if (g%bathyT(i,j) > min_depth) then
          idamp(i,j) = damp/86400.0
      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
 
          do j=js,je
            call calculate_density(t(:,j,1), s(:,j,1), pres, tmp(:,j), &
                             is, ie-is+1, tv%eqn_of_state)
          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