Skip to content
Snippets Groups Projects
zdfgls.F90 66.3 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
MODULE zdfgls
   !!======================================================================
   !!                       ***  MODULE  zdfgls  ***
   !! Ocean physics:  vertical mixing coefficient computed from the gls
   !!                 turbulent closure parameterization
   !!======================================================================
   !! History :  3.0  !  2009-09  (G. Reffray)  Original code
   !!            3.3  !  2010-10  (C. Bricaud)  Add in the reference
   !!            4.0  !  2017-04  (G. Madec)  remove CPP keys & avm at t-point only
   !!             -   !  2017-05  (G. Madec)  add top friction as boundary condition
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   zdf_gls       : update momentum and tracer Kz from a gls scheme
   !!   zdf_gls_init  : initialization, namelist read, and parameters control
   !!   gls_rst       : read/write gls restart in ocean restart file
   !!----------------------------------------------------------------------
   USE oce            ! ocean dynamics and active tracers
   USE dom_oce        ! ocean space and time domain
   USE domvvl         ! ocean space and time domain : variable volume layer
   USE zdfdrg  , ONLY : ln_drg_OFF            ! top/bottom free-slip flag
   USE zdfdrg  , ONLY : r_z0_top , r_z0_bot   ! top/bottom roughness
   USE zdfdrg  , ONLY : rCdU_top , rCdU_bot   ! top/bottom friction
   USE sbc_oce        ! surface boundary condition: ocean
   USE phycst         ! physical constants
   USE zdfmxl         ! mixed layer
   USE sbcwave , ONLY : hsw   ! significant wave height
#if defined key_si3
   USE ice, ONLY: hm_i, h_i
#endif
#if defined key_cice
   USE sbc_ice, ONLY: h_i
#endif
   !
   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
   USE lib_mpp        ! MPP manager
   USE prtctl         ! Print control
   USE in_out_manager ! I/O manager
   USE iom            ! I/O manager library
   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)

   IMPLICIT NONE
   PRIVATE

   PUBLIC   zdf_gls        ! called in zdfphy
   PUBLIC   zdf_gls_init   ! called in zdfphy
   PUBLIC   gls_rst        ! called in zdfphy

   !
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hmxl_n    !: now mixing length
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_surf !: Squared surface velocity scale at T-points
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_top  !: Squared top     velocity scale at T-points
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_bot  !: Squared bottom  velocity scale at T-points

   !                              !! ** Namelist  namzdf_gls  **
   LOGICAL  ::   ln_length_lim     ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988)
   LOGICAL  ::   ln_sigpsi         ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing
   INTEGER  ::   nn_mxlice         ! type of scaling under sea-ice (=0/1/2/3)
   INTEGER  ::   nn_bc_surf        ! surface boundary condition (=0/1)
   INTEGER  ::   nn_bc_bot         ! bottom boundary condition (=0/1)
   INTEGER  ::   nn_z0_met         ! Method for surface roughness computation
   INTEGER  ::   nn_z0_ice         ! Roughness accounting for sea ice
   INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2)
   INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen
   REAL(wp) ::   rn_clim_galp      ! Holt 2008 value for k-eps: 0.267
   REAL(wp) ::   rn_epsmin         ! minimum value of dissipation (m2/s3)
   REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2)
   REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
   REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing
   REAL(wp) ::   rn_hsro           ! Minimum surface roughness
   REAL(wp) ::   rn_hsri           ! Ice ocean roughness
   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)

   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1
   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn
   REAL(wp) ::   rghmin        = -0.28_wp
   REAL(wp) ::   rgh0          =  0.0329_wp
   REAL(wp) ::   rghcri        =  0.03_wp
   REAL(wp) ::   ra1           =  0.92_wp
   REAL(wp) ::   ra2           =  0.74_wp
   REAL(wp) ::   rb1           = 16.60_wp
   REAL(wp) ::   rb2           = 10.10_wp
   REAL(wp) ::   re2           =  1.33_wp
   REAL(wp) ::   rl1           =  0.107_wp
   REAL(wp) ::   rl2           =  0.0032_wp
   REAL(wp) ::   rl3           =  0.0864_wp
   REAL(wp) ::   rl4           =  0.12_wp
   REAL(wp) ::   rl5           = 11.9_wp
   REAL(wp) ::   rl6           =  0.4_wp
   REAL(wp) ::   rl7           =  0.0_wp
   REAL(wp) ::   rl8           =  0.48_wp
   REAL(wp) ::   rm1           =  0.127_wp
   REAL(wp) ::   rm2           =  0.00336_wp
   REAL(wp) ::   rm3           =  0.0906_wp
   REAL(wp) ::   rm4           =  0.101_wp
   REAL(wp) ::   rm5           = 11.2_wp
   REAL(wp) ::   rm6           =  0.4_wp
   REAL(wp) ::   rm7           =  0.0_wp
   REAL(wp) ::   rm8           =  0.318_wp
   REAL(wp) ::   rtrans        =  0.1_wp
   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters
   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        -
   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rfact_psi                  !     -           -           -        -
   REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        -
   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        -
   REAL(wp) ::   rs0, rs1, rs2, rs4, rs5, rs6                     !     -           -           -        -
   REAL(wp) ::   rd0, rd1, rd2, rd3, rd4, rd5                     !     -           -           -        -
   REAL(wp) ::   rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0  !     -           -           -        -
   REAL(wp) ::   rpsi3m, rpsi3p, rpp, rmm, rnn                    !     -           -           -        -
   !
   REAL(wp) ::   r2_3 = 2._wp/3._wp   ! constant=2/3

   !! * Substitutions
#  include "do_loop_substitute.h90"
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: zdfgls.F90 15145 2021-07-26 16:16:45Z smasson $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

   INTEGER FUNCTION zdf_gls_alloc()
      !!----------------------------------------------------------------------
      !!                ***  FUNCTION zdf_gls_alloc  ***
      !!----------------------------------------------------------------------
      ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustar2_surf(jpi,jpj) ,                     &
         &      zwall (jpi,jpj,jpk) , ustar2_top (jpi,jpj) , ustar2_bot(jpi,jpj) , STAT= zdf_gls_alloc )
         !
      CALL mpp_sum ( 'zdfgls', zdf_gls_alloc )
      IF( zdf_gls_alloc /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_alloc: failed to allocate arrays' )
   END FUNCTION zdf_gls_alloc


   SUBROUTINE zdf_gls( kt, Kbb, Kmm, p_sh2, p_avm, p_avt )
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE zdf_gls  ***
      !!
      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
      !!              coefficients using the GLS turbulent closure scheme.
      !!----------------------------------------------------------------------
      USE zdf_oce , ONLY : en, avtb, avmb   ! ocean vertical physics
      !!
      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time step
      INTEGER                             , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   p_sh2          ! shear production term
      REAL(wp), DIMENSION(:,:,:)          , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
      !
      INTEGER  ::   ji, jj, jk    ! dummy loop arguments
      INTEGER  ::   ibot, ibotm1  ! local integers
      INTEGER  ::   itop, itopp1  !   -       -
      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1 , zex2  ! local scalars
      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof, zdir  !   -      -
      REAL(wp) ::   zratio, zrn2, zflxb, sh     , z_en  !   -      -
      REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      -
      REAL(wp) ::   gh, gm, shr, dif, zsqen, zavt, zavm !   -      -
      REAL(wp) ::   zmsku, zmskv                        !   -      -
      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zdep
      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zkar
      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zflxs                 ! Turbulence fluxed induced by internal waves
      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zhsro                 ! Surface roughness (surface waves)
      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zice_fra              ! Tapering of wave breaking under sea ice
      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   eb                    ! tke at time before
      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   hmxl_b                ! mixing length at time before
      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   eps                   ! dissipation rate
      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwall_psi             ! Wall function use in the wb case (ln_sigpsi)
      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   psi                   ! psi at time now
      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zd_lw, zd_up, zdiag   ! lower, upper  and diagonal of the matrix
      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zstt, zstm            ! stability function on tracer and momentum
      !!--------------------------------------------------------------------
      !
      ! Preliminary computing
      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
         ustar2_surf(ji,jj) = 0._wp   ;   ustar2_top(ji,jj) = 0._wp   ;   ustar2_bot(ji,jj) = 0._wp
      END_2D
      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
         psi(ji,jj,jk) = 0._wp   ;   zwall_psi(ji,jj,jk) = 0._wp
      END_3D

      SELECT CASE ( nn_z0_ice )
      CASE( 0 )   ;   zice_fra(:,:) = 0._wp
      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(A2D(nn_hls)) * 10._wp )
      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(A2D(nn_hls))
      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp )
      END SELECT

      ! Compute surface, top and bottom friction at T-points
      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          !==  surface ocean friction
         ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1)   ! surface friction
      END_2D
      !
      !!gm Rq we may add here r_ke0(_top/_bot) ?  ==>> think about that...
      !
      IF( .NOT.ln_drg_OFF ) THEN     !== top/bottom friction   (explicit before friction)
         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          ! bottom friction (explicit before friction)
            zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
            zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
            ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  &
Loading
Loading full blame...