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/trunk/src/OCE/ZDF – NEMO

source: NEMO/trunk/src/OCE/ZDF/zdfddm.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 8.6 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   !!----------------------------------------------------------------------
34   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
35   !! $Id$
36   !! Software governed by the CeCILL license (see ./LICENSE)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE zdf_ddm( kt, Kmm, p_avm, p_avt, p_avs )
41      !!----------------------------------------------------------------------
42      !!                  ***  ROUTINE zdf_ddm  ***
43      !!                   
44      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
45      !!              effect of salt fingering and diffusive convection.
46      !!
47      !! ** Method  :   Diapycnal mixing is increased in case of double
48      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
49      !!      following Merryfield et al. (1999). The rate of double diffusive
50      !!      mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]):
51      !!         * salt fingering (Schmitt 1981):
52      !!      for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 )
53      !!      for R > 1 and rn2 > 0 : zavfs = O
54      !!      otherwise                : zavft = 0.7 zavs / R
55      !!         * diffusive layering (Federov 1988):
56      !!      for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) )
57      !!      otherwise                   : zavdt = 0
58      !!      for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85)
59      !!      for  0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R     
60      !!      otherwise                     : zavds = 0
61      !!         * update the eddy diffusivity:
62      !!      avt = avt + zavft + zavdt
63      !!      avs = avs + zavfs + zavds
64      !!      avm is required to remain at least above avt and avs.
65      !!     
66      !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
67      !!
68      !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
69      !!----------------------------------------------------------------------
70      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
71      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
72      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points)
73      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points)
74      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avs   !  Kz on salinity    (w-points)
75      !
76      INTEGER  ::   ji, jj , jk     ! dummy loop indices
77      REAL(wp) ::   zaw, zbw, zrw   ! local scalars
78      REAL(wp) ::   zdt, zds
79      REAL(wp) ::   zinr, zrr       !   -      -
80      REAL(wp) ::   zavft, zavfs    !   -      -
81      REAL(wp) ::   zavdt, zavds    !   -      -
82      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
83      !!----------------------------------------------------------------------
84      !
85      !                                                ! ===============
86      DO jk = 2, jpkm1                                 ! Horizontal slab
87         !                                             ! ===============
88         ! Define the mask
89         ! ---------------
90!!gm  WORK to be done:   change the code from vector optimisation to scalar one.
91!!gm                     ==>>>  test in the loop instead of use of mask arrays
92!!gm                            and many acces in memory
93         
94         DO_2D_11_11
95            zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
96!!gm please, use e3w(:,:,:,Kmm) below
97               &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
98            !
99            zaw = (  rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw  )  &
100                &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
101            zbw = (  rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw  )  &
102                &    * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
103            !
104            zdt = zaw * ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) )
105            zds = zbw * ( ts(ji,jj,jk-1,jp_sal,Kmm) - ts(ji,jj,jk,jp_sal,Kmm) ) 
106            IF( ABS( zds) <= 1.e-20_wp )   zds = 1.e-20_wp
107            zrau(ji,jj) = MAX(  1.e-20, zdt / zds  )    ! only retains positive value of zrau
108         END_2D
109
110         DO_2D_11_11
111            ! stability indicator: msks=1 if rn2>0; 0 elsewhere
112            IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
113            ELSE                                       ;   zmsks(ji,jj) = 1._wp
114            ENDIF
115            ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere           
116            IF( zrau(ji,jj) <= 1.             ) THEN   ;   zmskf(ji,jj) = 0._wp
117            ELSE                                       ;   zmskf(ji,jj) = 1._wp
118            ENDIF
119            ! diffusive layering indicators:
120            !     ! mskdl1=1 if 0< R <1; 0 elsewhere
121            IF( zrau(ji,jj) >= 1.             ) THEN   ;   zmskd1(ji,jj) = 0._wp
122            ELSE                                       ;   zmskd1(ji,jj) = 1._wp
123            ENDIF
124            !     ! mskdl2=1 if 0< R <0.5; 0 elsewhere
125            IF( zrau(ji,jj) >= 0.5            ) THEN   ;   zmskd2(ji,jj) = 0._wp
126            ELSE                                       ;   zmskd2(ji,jj) = 1._wp
127            ENDIF
128            !   mskdl3=1 if 0.5< R <1; 0 elsewhere
129            IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN   ;   zmskd3(ji,jj) = 0._wp
130            ELSE                                                   ;   zmskd3(ji,jj) = 1._wp
131            ENDIF
132         END_2D
133         ! mask zmsk in order to have avt and avs masked
134         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
135
136
137         ! Update avt and avs
138         ! ------------------
139         ! Constant eddy coefficient: reset to the background value
140         DO_2D_11_11
141            zinr = 1._wp / zrau(ji,jj)
142            ! salt fingering
143            zrr = zrau(ji,jj) / rn_hsbfr
144            zrr = zrr * zrr
145            zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj)
146            zavft = 0.7 * zavfs * zinr
147            ! diffusive layering
148            zavdt = 1.3635e-6 * EXP(  4.6 * EXP( -0.54*(zinr-1.) )  ) * zmsks(ji,jj) * zmskd1(ji,jj)
149            zavds = zavdt * zmsks(ji,jj) * (  ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj)   &
150               &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
151            ! add to the eddy viscosity coef. previously computed
152            p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds
153            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt
154            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
155         END_2D
156         !                                                ! ===============
157      END DO                                              !   End of slab
158      !                                                   ! ===============
159      !
160      IF(sn_cfctl%l_prtctl) THEN
161         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', kdim=jpk)
162      ENDIF
163      !
164   END SUBROUTINE zdf_ddm
165   
166   !!======================================================================
167END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.