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
Line 
1MODULE zdfddm
2   !!======================================================================
3   !!                       ***  MODULE  zdfddm  ***
4   !! Ocean physics : double diffusion mixing parameterization
5   !!======================================================================
6   !! History :  OPA  ! 2000-08  (G. Madec)  double diffusive mixing
7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
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
10   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   zdf_ddm       : compute the Kz for salinity
15   !!----------------------------------------------------------------------
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
19   USE eosbn2         ! equation of state
20   !
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
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   zdf_ddm       ! called by step.F90
30
31   !! * Substitutions
32#  include "do_loop_substitute.h90"
33#  include "domzgr_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
36   !! $Id$
37   !! Software governed by the CeCILL license (see ./LICENSE)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE zdf_ddm( kt, Kmm, p_avm, p_avt, p_avs )
42      !!----------------------------------------------------------------------
43      !!                  ***  ROUTINE zdf_ddm  ***
44      !!                   
45      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
46      !!              effect of salt fingering and diffusive convection.
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
51      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
52      !!         * salt fingering (Schmitt 1981):
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
56      !!         * diffusive layering (Federov 1988):
57      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
58      !!      otherwise                   : zavdt = 0
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     
61      !!      otherwise                     : zavds = 0
62      !!         * update the eddy diffusivity:
63      !!      avt = avt + zavft + zavdt
64      !!      avs = avs + zavfs + zavds
65      !!      avm is required to remain at least above avt and avs.
66      !!     
67      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
68      !!
69      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
70      !!----------------------------------------------------------------------
71      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
72      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
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)
76      !
77      INTEGER  ::   ji, jj , jk     ! dummy loop indices
78      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
79      REAL(wp) ::   zdt, zds
80      REAL(wp) ::   zinr            !   -      -
81      REAL(dp) ::         zrr       !   -      -
82      REAL(wp) ::   zavft           !   -      -
83      REAL(dp) ::          zavfs    !   -      -
84      REAL(wp) ::   zavdt, zavds    !   -      -
85      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
86      !!----------------------------------------------------------------------
87      !
88      !                                                ! ===============
89      DO jk = 2, jpkm1                                 ! Horizontal slab
90         !                                             ! ===============
91         ! Define the mask
92         ! ---------------
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         
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])  ==!
99            zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
100!!gm please, use e3w at Kmm below
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
113
114         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )           !==  indicators  ==!
115         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )           !==  indicators  ==!
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
138         ! mask zmsk in order to have avt and avs masked
139         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
140
141
142         ! Update avt and avs
143         ! ------------------
144         ! Constant eddy coefficient: reset to the background value
145         ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 )
146         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
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
162         !                                                ! ===============
163      END DO                                              !   End of slab
164      !                                                   ! ===============
165      !
166      IF(sn_cfctl%l_prtctl) THEN
167         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', kdim=jpk)
168      ENDIF
169      !
170   END SUBROUTINE zdf_ddm
171   
172   !!======================================================================
173END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.