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 branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 7037

Last change on this file since 7037 was 7037, checked in by mocavero, 8 years ago

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