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/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF/zdfddm.F90 @ 13998

Last change on this file since 13998 was 13998, checked in by techene, 3 years ago

branch updated with trunk 13787
  • Property svn:keywords set to Id
File size: 8.8 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         DO_2D( 1, 1, 1, 1 )           !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
98            zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
99!!gm please, use e3w at Kmm below
100               &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
101            !
102            zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
103                &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
104            zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
105                &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
106            !
107            zdt = zaw * ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) )
108            zds = zbw * ( ts(ji,jj,jk-1,jp_sal,Kmm) - ts(ji,jj,jk,jp_sal,Kmm) ) 
109            IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
110            zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
111         END_2D
112
113         DO_2D( 1, 1, 1, 1 )           !==  indicators  ==!
114            ! stability indicator: msks=1 if rn2>0; 0 elsewhere
115            IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
116            ELSE                                       ;   zmsks(ji,jj) = 1._wp
117            ENDIF
118            ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
119            IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
120            ELSE                                       ;   zmskf(ji,jj) = 1._wp
121            ENDIF
122            ! diffusive layering indicators:
123            !     ! mskdl1=1 if 0< R <1; 0 elsewhere
124            IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
125            ELSE                                       ;   zmskd1(ji,jj) = 1._wp
126            ENDIF
127            !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
128            IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
129            ELSE                                       ;   zmskd2(ji,jj) = 1._wp
130            ENDIF
131            !   mskdl3=1 if 0.5< R <1; 0 elsewhere
132            IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
133            ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
134            ENDIF
135         END_2D
136         ! mask zmsk in order to have avt and avs masked
137         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
138
139
140         ! Update avt and avs
141         ! ------------------
142         ! Constant eddy coefficient: reset to the background value
143         DO_2D( 1, 1, 1, 1 )
144            zinr = 1._wp / zrau(ji,jj)
145            ! salt fingering
146            zrr = zrau(ji,jj) / rn_hsbfr
147            zrr = zrr * zrr
148            zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
149            zavft = 0.7 * zavfs * zinr
150            ! diffusive layering
151            zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
152            zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
153               &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
154            ! add to the eddy viscosity coef. previously computed
155            p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds
156            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt
157            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
158         END_2D
159         !                                                ! ===============
160      END DO                                              !   End of slab
161      !                                                   ! ===============
162      !
163      IF(sn_cfctl%l_prtctl) THEN
164         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', kdim=jpk)
165      ENDIF
166      !
167   END SUBROUTINE zdf_ddm
168   
169   !!======================================================================
170END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.