New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdfddm.F90 in NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF – NEMO

source: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/ZDF/zdfddm.F90 @ 14601

Last change on this file since 14601 was 14601, checked in by francesca, 3 years ago

[comm_cleanup: ZDF files] - ticket #2607
  • Property svn:keywords set to Id
File size: 9.1 KB
RevLine 
[3]1MODULE zdfddm
2   !!======================================================================
3   !!                       ***  MODULE  zdfddm  ***
4   !! Ocean physics : double diffusion mixing parameterization
5   !!======================================================================
[1601]6   !! History :  OPA  ! 2000-08  (G. Madec)  double diffusive mixing
7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
[4990]8   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
9   !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta
[9019]10   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
[1601]11   !!----------------------------------------------------------------------
[9019]12
[3]13   !!----------------------------------------------------------------------
[9019]14   !!   zdf_ddm       : compute the Kz for salinity
[3]15   !!----------------------------------------------------------------------
[9019]16   USE oce            ! ocean dynamics and tracers variables
17   USE dom_oce        ! ocean space and time domain variables
18   USE zdf_oce        ! ocean vertical physics variables
[4990]19   USE eosbn2         ! equation of state
20   !
[9019]21   USE in_out_manager ! I/O manager
22   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
23   USE prtctl         ! Print control
24   USE lib_mpp        ! MPP library
[3]25
26   IMPLICIT NONE
27   PRIVATE
28
[2528]29   PUBLIC   zdf_ddm       ! called by step.F90
[3]30
31   !! * Substitutions
[12377]32#  include "do_loop_substitute.h90"
[14053]33#  include "domzgr_substitute.h90"
[3]34   !!----------------------------------------------------------------------
[9598]35   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2528]36   !! $Id$
[10068]37   !! Software governed by the CeCILL license (see ./LICENSE)
[3]38   !!----------------------------------------------------------------------
39CONTAINS
40
[12377]41   SUBROUTINE zdf_ddm( kt, Kmm, p_avm, p_avt, p_avs )
[2715]42      !!----------------------------------------------------------------------
[3]43      !!                  ***  ROUTINE zdf_ddm  ***
44      !!                   
45      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
[1601]46      !!              effect of salt fingering and diffusive convection.
[3]47      !!
48      !! ** Method  :   Diapycnal mixing is increased in case of double
49      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
50      !!      following Merryfield et al. (1999). The rate of double diffusive
[4990]51      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
[3]52      !!         * salt fingering (Schmitt 1981):
[4990]53      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 )
54      !!      for R > 1 and rn2 > 0 : zavfs = O
55      !!      otherwise                : zavft = 0.7 zavs / R
[3]56      !!         * diffusive layering (Federov 1988):
[4990]57      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
[3]58      !!      otherwise                   : zavdt = 0
[4990]59      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85)
60      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R     
[3]61      !!      otherwise                     : zavds = 0
62      !!         * update the eddy diffusivity:
63      !!      avt = avt + zavft + zavdt
64      !!      avs = avs + zavfs + zavds
[9019]65      !!      avm is required to remain at least above avt and avs.
[3]66      !!     
[1601]67      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
[3]68      !!
[1601]69      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
[3]70      !!----------------------------------------------------------------------
[12377]71      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
72      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
[9019]73      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points)
74      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points)
75      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avs   !  Kz on salinity    (w-points)
[2715]76      !
[1601]77      INTEGER  ::   ji, jj , jk     ! dummy loop indices
[4990]78      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
79      REAL(wp) ::   zdt, zds
[13226]80      REAL(wp) ::   zinr            !   -      -
81      REAL(dp) ::         zrr       !   -      -
82      REAL(wp) ::   zavft           !   -      -
83      REAL(dp) ::          zavfs    !   -      -
[4990]84      REAL(wp) ::   zavdt, zavds    !   -      -
[9019]85      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
[3]86      !!----------------------------------------------------------------------
[3294]87      !
[3]88      !                                                ! ===============
89      DO jk = 2, jpkm1                                 ! Horizontal slab
90         !                                             ! ===============
91         ! Define the mask
92         ! ---------------
[9019]93!!gm  WORK to be done:   change the code from vector optimisation to scalar one.
94!!gm                     ==>>>  test in the loop instead of use of mask arrays
95!!gm                            and many acces in memory
96         
[14601]97         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )           !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
98         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )           !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
[12377]99            zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
[13237]100!!gm please, use e3w at Kmm below
[12377]101               &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
102            !
103            zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
104                &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
105            zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
106                &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
107            !
108            zdt = zaw * ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) )
109            zds = zbw * ( ts(ji,jj,jk-1,jp_sal,Kmm) - ts(ji,jj,jk,jp_sal,Kmm) ) 
110            IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
111            zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
112         END_2D
[3]113
[14601]114         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )           !==  indicators  ==!
115         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )           !==  indicators  ==!
[12377]116            ! stability indicator: msks=1 if rn2>0; 0 elsewhere
117            IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
118            ELSE                                       ;   zmsks(ji,jj) = 1._wp
119            ENDIF
120            ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
121            IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
122            ELSE                                       ;   zmskf(ji,jj) = 1._wp
123            ENDIF
124            ! diffusive layering indicators:
125            !     ! mskdl1=1 if 0< R <1; 0 elsewhere
126            IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
127            ELSE                                       ;   zmskd1(ji,jj) = 1._wp
128            ENDIF
129            !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
130            IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
131            ELSE                                       ;   zmskd2(ji,jj) = 1._wp
132            ENDIF
133            !   mskdl3=1 if 0.5< R <1; 0 elsewhere
134            IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
135            ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
136            ENDIF
137         END_2D
[3]138         ! mask zmsk in order to have avt and avs masked
[7753]139         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
[3]140
141
142         ! Update avt and avs
143         ! ------------------
144         ! Constant eddy coefficient: reset to the background value
[14601]145         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
146         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[12377]147            zinr = 1._wp / zrau(ji,jj)
148            ! salt fingering
149            zrr = zrau(ji,jj) / rn_hsbfr
150            zrr = zrr * zrr
151            zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
152            zavft = 0.7 * zavfs * zinr
153            ! diffusive layering
154            zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
155            zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
156               &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
157            ! add to the eddy viscosity coef. previously computed
158            p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds
159            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt
160            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
161         END_2D
[3]162         !                                                ! ===============
163      END DO                                              !   End of slab
164      !                                                   ! ===============
[1601]165      !
[12377]166      IF(sn_cfctl%l_prtctl) THEN
[9440]167         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', kdim=jpk)
[49]168      ENDIF
[1601]169      !
[3]170   END SUBROUTINE zdf_ddm
171   
172   !!======================================================================
173END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.