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/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago

The Dr Hook changes from my perl code.

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