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

Last change on this file since 106 was 106, checked in by opalod, 20 years ago

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