source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 4147

Last change on this file since 4147 was 4147, checked in by cetlod, 7 years ago

merge in dev_LOCEAN_2013, the 1st development branch dev_r3853_CNRS9_Confsetting, from its starting point ( r3853 ) on the trunk: see ticket #1169

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