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

source: trunk/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 1527

Last change on this file since 1527 was 1527, checked in by ctlod, 15 years ago

tke2: minor bug correction on the TKE vertical diffusion term, see ticket: #484
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 11.3 KB
Line 
1MODULE zdfddm
2   !!======================================================================
3   !!                       ***  MODULE  zdfddm  ***
4   !! Ocean physics : double diffusion mixing parameterization
5   !!======================================================================
6#if defined key_zdfddm   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_zdfddm' :                                     double diffusion
9   !!----------------------------------------------------------------------
10   !!   zdf_ddm       : compute the Ks for salinity
11   !!   zdf_ddm_init  : read namelist and control the parameters
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce             ! ocean dynamics and tracers variables
15   USE dom_oce         ! ocean space and time domain variables
16   USE zdf_oce         ! ocean vertical physics variables
17   USE in_out_manager  ! I/O manager
18   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
19   USE prtctl          ! Print control
20
21   IMPLICIT NONE
22   PRIVATE
23
24   !! * Routine accessibility
25   PUBLIC zdf_ddm     ! called by step.F90
26
27   !! * Shared module variables
28   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.    !: double diffusive mixing flag
29   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
30      avs ,               &  !: salinity vertical diffusivity coeff. at w-point
31      rrau                   !: heat/salt buoyancy flux ratio
32
33   !! * Module variables
34   REAL(wp) ::            & !!! * double diffusive mixing namelist *
35      avts  = 1.e-4_wp ,  &  ! maximum value of avs for salt fingering
36      hsbfr = 1.6_wp         ! heat/salt buoyancy flux ratio
37
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !!   OPA 9.0 , LOCEAN-IPSL (2005)
42   !! $Id$
43   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
44   !!----------------------------------------------------------------------
45
46CONTAINS
47
48   SUBROUTINE zdf_ddm( kt )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE zdf_ddm  ***
51      !!                   
52      !! ** Purpose :   Add to the vertical eddy diffusivity coefficient the
53      !!      effect of salt fingering and diffusive convection.
54      !!
55      !! ** Method  :   Diapycnal mixing is increased in case of double
56      !!      diffusive mixing (i.e. salt fingering and diffusive layering)
57      !!      following Merryfield et al. (1999). The rate of double diffusive
58      !!      mixing depend on the buoyancy ratio: Rrau=alpha/beta dk[T]/dk[S]
59      !!      which is computed in rn2.F
60      !!         * salt fingering (Schmitt 1981):
61      !!      for Rrau > 1 and rn2 > 0 : zavfs = avts / ( 1 + (Rrau/hsbfr)^6 )
62      !!      for Rrau > 1 and rn2 > 0 : zavfs = O
63      !!      otherwise                : zavft = 0.7 zavs / Rrau
64      !!         * diffusive layering (Federov 1988):
65      !!      for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6 
66      !!                                 * exp( 4.6 exp(-0.54 (1/Rrau-1) ) )
67      !!      otherwise                   : zavdt = 0
68      !!      for .5 < Rrau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau -0.85)
69      !!      for  0 < Rrau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau     
70      !!      otherwise                     : zavds = 0
71      !!         * update the eddy diffusivity:
72      !!      avt = avt + zavft + zavdt
73      !!      avs = avs + zavfs + zavds
74      !!      avmu, avmv are required to remain at least above avt and avs.
75      !!     
76      !! ** Action  :   avt, avs : update vertical eddy diffusivity coef.
77      !!                           for temperature and salinity
78      !!
79      !! References :
80      !!      Merryfield et al., JPO, 29, 1124-1142, 1999.
81      !! History :
82      !!        !  00-08  (G. Madec)  double diffusive mixing
83      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
84      !!----------------------------------------------------------------------
85      !! * Arguments
86      INTEGER, INTENT( in ) ::   kt         ! ocean time-step indexocean time step
87
88      !! * Local declarations
89      INTEGER ::   ji, jj , jk              ! dummy loop indices
90      REAL(wp), DIMENSION(jpi,jpj) ::   &
91         zmsks, zmskf,                    & ! temporary workspace
92         zmskd1, zmskd2, zmskd3             !    "           "
93      REAL(wp) ::   &
94         zinr, zrr,                       & ! temporary scalars
95         zavft, zavfs,                    & !    "         "
96         zavdt, zavds                       !    "         "
97      !!----------------------------------------------------------------------
98
99
100      IF ( kt == nit000 )   CALL zdf_ddm_init          ! Initialization (first time-step only)
101
102
103      ! Compute avs
104      ! -----------
105      !                                                ! ===============
106      DO jk = 2, jpkm1                                 ! Horizontal slab
107         !                                             ! ===============
108         ! Define the mask
109         ! ---------------
110         ! only retains positive value of rrau
111         rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )
112
113         ! indicators:
114         DO jj = 1, jpj
115            DO ji = 1, jpi
116               ! stability indicator: msks=1 if rn2>0; 0 elsewhere
117               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN
118                  zmsks(ji,jj) = 0.e0
119               ELSE
120                  zmsks(ji,jj) = 1.e0
121               ENDIF
122               ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere           
123               IF( rrau(ji,jj,jk) <= 1. ) THEN
124                  zmskf(ji,jj) = 0.e0
125               ELSE
126                  zmskf(ji,jj) = 1.e0
127               ENDIF
128               ! diffusive layering indicators:
129               !   mskdl1=1 if 0<rrau<1; 0 elsewhere
130               IF( rrau(ji,jj,jk) >= 1. ) THEN
131                  zmskd1(ji,jj) = 0.e0
132               ELSE
133                  zmskd1(ji,jj) = 1.e0
134               ENDIF
135               !   mskdl2=1 if 0<rrau<0.5; 0 elsewhere
136               IF( rrau(ji,jj,jk) >= 0.5 ) THEN
137                  zmskd2(ji,jj) = 0.e0
138               ELSE
139                  zmskd2(ji,jj) = 1.e0
140               ENDIF
141               !   mskdl3=1 if 0.5<rrau<1; 0 elsewhere
142               IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN
143                  zmskd3(ji,jj) = 0.e0
144               ELSE
145                  zmskd3(ji,jj) = 1.e0
146               ENDIF
147            END DO
148         END DO
149         ! mask zmsk in order to have avt and avs masked
150         zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
151
152
153         ! Update avt and avs
154         ! ------------------
155         ! Constant eddy coefficient: reset to the background value
156!CDIR NOVERRCHK
157         DO jj = 1, jpj
158!CDIR NOVERRCHK
159            DO ji = 1, jpi
160               zinr = 1./rrau(ji,jj,jk)
161               ! salt fingering
162               zrr = rrau(ji,jj,jk)/hsbfr
163               zrr = zrr * zrr
164               zavfs = avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) *zmskf(ji,jj)
165               zavft = 0.7 * zavfs * zinr
166               ! diffusive layering
167               zavdt = 1.3635e-6 * EXP(4.6*EXP(-0.54*(zinr-1.) ) )   &
168                                 * zmsks(ji,jj) * zmskd1(ji,jj)
169               zavds = zavdt * zmsks(ji,jj)   &
170                     * ( (1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj)   &
171                       +  0.15 * rrau(ji,jj,jk)          * zmskd2(ji,jj)  )
172               ! add to the eddy viscosity coef. previously computed
173               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
174               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
175               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
176            END DO
177         END DO
178
179
180         ! Increase avmu, avmv if necessary
181         ! --------------------------------
182!!gm to be changed following the definition of avm.
183         DO jj = 1, jpjm1
184            DO ji = 1, fs_jpim1   ! vector opt.
185               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
186                                     avt(ji,jj,jk), avt(ji+1,jj,jk),   &
187                                     avs(ji,jj,jk), avs(ji+1,jj,jk) )   &
188                              * umask(ji,jj,jk)
189               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
190                                     avt(ji,jj,jk), avt(ji,jj+1,jk),   &
191                                     avs(ji,jj,jk), avs(ji,jj+1,jk) )   &
192                              * vmask(ji,jj,jk)
193            END DO
194         END DO
195         !                                                ! ===============
196      END DO                                              !   End of slab
197      !                                                   ! ===============
198     
199      ! Lateral boundary conditions on ( avt, avs, avmu, avmv )   (unchanged sign)
200      ! -------------------------------========================
201      CALL lbc_lnk( avt , 'W', 1. )
202      CALL lbc_lnk( avs , 'W', 1. )
203      CALL lbc_lnk( avm , 'W', 1. )
204      CALL lbc_lnk( avmu, 'U', 1. ) 
205      CALL lbc_lnk( avmv, 'V', 1. )
206
207      IF(ln_ctl) THEN
208         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
209         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
210            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
211      ENDIF
212     
213   END SUBROUTINE zdf_ddm
214   
215   
216   SUBROUTINE zdf_ddm_init
217      !!----------------------------------------------------------------------
218      !!                  ***  ROUTINE zdf_ddm_init  ***
219      !!
220      !! ** Purpose :   Initialization of double diffusion mixing scheme
221      !!
222      !! ** Method  :   Read the nammbf namelist and check the parameter values
223      !!      called by zdf_ddm at the first timestep (nit000)
224      !!
225      !! History :
226      !!   8.5  !  02-08  (G. Madec)  Original code
227      !!----------------------------------------------------------------------
228      NAMELIST/namddm/ avts, hsbfr
229      !!----------------------------------------------------------------------
230
231      ! Read Namelist namddm : double diffusion mixing scheme
232      ! --------------------
233      REWIND ( numnam )
234      READ   ( numnam, namddm )
235
236
237      ! Parameter control and print
238      ! ---------------------------
239      IF(lwp) THEN
240         WRITE(numout,*)
241         WRITE(numout,*) 'zdf_ddm : double diffusive mixing'
242         WRITE(numout,*) '~~~~~~~'
243         WRITE(numout,*) '          Namelist namddm : set dd mixing parameter'
244         WRITE(numout,*) '             maximum avs for dd mixing      avts   = ', avts
245         WRITE(numout,*) '             heat/salt buoyancy flux ratio  hsbfr  = ', hsbfr
246         WRITE(numout,*)
247      ENDIF
248
249   END SUBROUTINE zdf_ddm_init
250
251#else
252   !!----------------------------------------------------------------------
253   !!   Default option :          Dummy module          No double diffusion
254   !!----------------------------------------------------------------------
255   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
256CONTAINS
257   SUBROUTINE zdf_ddm( kt )           ! Dummy routine
258      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
259   END SUBROUTINE zdf_ddm
260#endif
261
262   !!======================================================================
263END MODULE zdfddm
Note: See TracBrowser for help on using the repository browser.