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 @ 1537

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

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