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.
zdftmx.F90 in trunk/NEMO/OPA_SRC/ZDF – NEMO

source: trunk/NEMO/OPA_SRC/ZDF/zdftmx.F90 @ 1437

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

add zdftmx.F90 module linked to the calculation of vertical viscosity/diffusion coefficients due to internal tidal mixing, see ticket: #425
File size: 22.5 KB
Line 
1MODULE zdftmx
2   !!========================================================================
3   !!                       ***  MODULE  zdftmx  ***
4   !! Ocean physics: vertical tidal mixing coefficient
5   !!========================================================================
6   !! History :  1.0  !  2004-04  (L. Bessieres, G. Madec)  Original code
7   !!             -   !  2006-08  (A. Koch-Larrouy) Indonesian strait
8   !!----------------------------------------------------------------------
9#if defined key_zdftmx
10   !!----------------------------------------------------------------------
11   !!   'key_zdftmx'                                  Tidal vertical mixing
12   !!----------------------------------------------------------------------
13   !!   zdf_tmx      : global     momentum & tracer Kz with tidal induced Kz
14   !!   tmx_itf      : Indonesian momentum & tracer Kz with tidal induced Kz
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers variables
17   USE dom_oce         ! ocean space and time domain variables
18   USE zdf_oce         ! ocean vertical physics variables
19   USE in_out_manager  ! I/O manager
20   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
21   USE eosbn2          ! ocean equation of state
22   USE phycst          ! physical constants
23   USE prtctl          ! Print control
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC zdf_tmx          ! called in step module
29
30   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: tidal mixing flag
31
32   !                                                    !!* Namelist  namtmx *
33   REAL(wp) ::  rn_htmx  = 500.                          ! vertical decay scale for turbulence (meters)
34   REAL(wp) ::  rn_n2min = 1.e-8                         ! threshold of the Brunt-Vaisala frequency (s-1)
35   REAL(wp) ::  rn_tfe   = 1./3.                         ! tidal dissipation efficiency (St Laurent et al. 2002)
36   REAL(wp) ::  rn_tfe_itf  = 1.                         ! ITF tidal dissipation efficiency (St Laurent et al. 2002)
37   REAL(wp) ::  rn_me    = 0.2                           ! mixing efficiency (Osborn 1980)
38
39   REAL(wp), DIMENSION(jpi,jpj) ::   en_tmx              ! energy available for tidal mixing (W/m2)
40   REAL(wp), DIMENSION(jpi,jpj) ::   mask_itf            ! mask to use over Indonesian area
41
42   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   az_tmx          ! coefficient used to evaluate the tidal induced Kz
43   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   az_tmx_itf      ! coefficient used to evaluate the tidal induced Kz
44   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   az_tsm          ! avt(tidal) = avt0 + az_tmx / n2
45   !                                                     ! vertical structure function taking into account
46   !                                                     ! of the subscale topographic effect
47   REAL(wp), DIMENSION(jpk) ::   zempba                  ! zempba profil vertical d'energie
48   REAL(wp), DIMENSION(jpi,jpj) ::  zhdep                ! Ocean depth
49   REAL(wp), DIMENSION(jpi,jpj) ::  zsum1, zsum2, zsum   
50   REAL(wp), DIMENSION(jpi,jpj) ::  z1sum1, z1sum2, z1sum 
51
52   REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zempba_3d_1, zempba_3d_2
53   REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zempba_3d, drn2dz                              !
54
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58   !!----------------------------------------------------------------------
59   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
60   !! $Id:$
61   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
62   !!----------------------------------------------------------------------
63
64CONTAINS
65
66   SUBROUTINE zdf_tmx( kt )
67      !!----------------------------------------------------------------------
68      !!                  ***  ROUTINE zdf_tmx  ***
69      !!                   
70      !! ** Purpose :   add to the vertical mixing coefficients the effect of
71      !!      tidal mixing (Simmons et al 2004).
72      !!
73      !! ** Method  : - tidal-induced vertical mixing is given by:
74      !!      Kz_tides = az_tmx / rn2, where az_tmx is a 2D coefficient set in
75      !!      zdf_tmx_init. This mixing is added to avt, avmu and avmv as
76      !!      follows:  avt  = avt  +    az_tmx
77      !!                avmu = avmu + mi(az_tmx)
78      !!                avmv = avmv + mj(az_tmx)
79      !!
80      !! ** Action  :   Update avt, avmu and avmv
81      !!
82      !! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
83      !!              St Laurent et al. 2002, GRL, 29, 2106, 10.1029/2002GL015 633.
84      !!              Osborn 1982, JPO, 10, 83-89.
85      !!----------------------------------------------------------------------
86      INTEGER, INTENT(in) ::   kt   ! ocean time-step
87      !!
88      INTEGER  ::   ji, jj, jk   ! dummy loop indices
89      REAL(wp) ::   ztpc         ! scalar workspace
90      REAL(wp), DIMENSION(jpi,jpj) ::   zkz   ! temporary 2D workspace
91      !!----------------------------------------------------------------------
92
93      !                             ! Initialization (first time-step only)
94      IF( kt == nit000  )   CALL zdf_tmx_init
95
96      !                             ! First estimation (with n2 bound by rn_n2min) bounded by 60 cm2/s
97      av_tide(:,:,:) = MIN(  60.e-4, az_tmx(:,:,:) / MAX( rn_n2min, rn2(:,:,:) )  )
98
99      zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column
100      DO jk = 2, jpkm1
101         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * av_tide(:,:,jk)* tmask(:,:,jk)
102      END DO
103
104      DO jj = 1, jpj                ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
105         DO ji = 1, jpi
106            IF( zkz(ji,jj) /= 0.e0 )   zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
107         END DO
108      END DO
109
110      DO jk = 2, jpkm1              ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> av_tide bound by 300 cm2/s
111         av_tide(:,:,jk) = av_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s
112      END DO
113
114      IF( kt == nit000 ) THEN       ! diagnose the nergy consumed by av_tide
115         ztpc = 0.e0
116         DO jk= 1, jpk
117            DO jj= 1, jpj
118               DO ji= 1, jpi
119                  ztpc = ztpc + fse3w(ji,jj,jk)*e1t(ji,jj)*e2t(ji,jj)   &
120                     &         *MAX(0.e0,rn2(ji,jj,jk))*av_tide(ji,jj,jk)*tmask(ji,jj,jk)*tmask_i(ji,jj)
121               END DO
122            END DO
123         END DO
124         ztpc= rau0 * 1/(rn_tfe * rn_me) * ztpc
125         IF(lwp) WRITE(numout,*) 
126         IF(lwp) WRITE(numout,*) '          N Total power consumption by av_tide    : ztpc = ', ztpc * 1.e-12 ,'TW'
127      ENDIF
128
129      DO jk = 2, jpkm1              ! update momentum & tracer diffusivity with tidal mixing
130         av_tide(:,:,jk) = av_tide(:,:,jk) * ( 1. - mask_itf(:,:) )
131         avt    (:,:,jk) = avt    (:,:,jk) + av_tide(:,:,jk)
132         DO jj = 2, jpjm1
133            DO ji = fs_2, fs_jpim1  ! vector opt.
134               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( av_tide(ji,jj,jk) + av_tide(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
135               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( av_tide(ji,jj,jk) + av_tide(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
136            END DO
137         END DO
138      END DO
139      CALL lbc_lnk( avmu, 'U', 1. )   
140      CALL lbc_lnk( avmv, 'V', 1. )
141       
142      CALL tmx_itf( kt )            ! Indonesian throughflow tidal vertical mixing
143
144      IF(ln_ctl) THEN
145         CALL prt_ctl_info( clinfo1='tmx  - t: ', ivar1=INT(SUM( avt (:,:,:) ) ) )
146         CALL prt_ctl_info( clinfo1='tmx  - u: ', ivar1=INT(SUM( avmu(:,:,:) ) ),   &
147            &               clinfo2='     - v: ', ivar2=INT(SUM( avmv(:,:,:) ) ) )
148      ENDIF
149      !
150   END SUBROUTINE zdf_tmx
151
152
153   SUBROUTINE tmx_itf( kt )
154      !!----------------------------------------------------------------------
155      !!                  ***  ROUTINE tmx_itf  ***
156      !!                   
157      !! ** Purpose :   add to the vertical mixing coefficients the effect of
158      !!      tidal mixing (st laurent et al. 2004) + specificities for the ITF
159      !!      region (vertical profil of energy (P. Bouruet Aubertot + Theo Gerkema)
160      !!      , q=1 all the energy remains confined)
161      !!
162      !! ** Method  : - tidal-induced vertical mixing is given by:
163      !!      Kz_tides = az_tmx / rn2, where az_tmx is a 2D coefficient set in
164      !!      zdf_tmx_init. This mixing is added to avt, avmu and avmv as
165      !!      follows:  avt  = avt  +    az_tmx
166      !!                avmu = avmu + mi(az_tmx)
167      !!                avmv = avmv + mj(az_tmx)
168      !!
169      !! ** Action  :   Update avt, avmu and avmv
170      !!
171      !! References :  Koch-Larrouy et al. 2007, GRL
172      !!----------------------------------------------------------------------
173      INTEGER, INTENT(in) ::   kt   ! ocean time-step
174      !!
175      INTEGER  ::   ji, jj, jk   ! dummy loop indices
176      REAL(wp) ::   ztpc         ! total power consumption
177      REAL(wp), DIMENSION(jpi,jpj) ::   zkz   ! temporary 2D workspace
178      !!----------------------------------------------------------------------
179
180      !                             ! compute the form function using N2 at each time step
181      zempba_3d_1(:,:,:) = 0.e0
182      zempba_3d_2(:,:,:) = 0.e0
183   
184      DO jk = 1, jpk-1              ! Vertical profile dN2/dz
185         DO jj = 1, jpj               
186            DO ji = 1, jpi
187               drn2dz(ji,jj,jk) =(rn2(ji,jj,jk) - rn2(ji,jj,jk+1) )
188               IF( rn2(ji,jj,jk) > 0 ) THEN
189                  zempba_3d_1(ji,jj,jk)=SQRT(rn2(ji,jj,jk))*tmask(ji,jj,jk)
190                  zempba_3d_2(ji,jj,jk)=rn2(ji,jj,jk)*tmask(ji,jj,jk)
191               ENDIF
192            END DO
193         END DO
194      END DO
195
196      zsum1(:,:)=0.e0
197      zsum2(:,:)=0.e0
198      DO jk= 2, jpk
199         DO jj = 1, jpj
200            DO ji = 1, jpi
201               zsum1(ji,jj)=zsum1(ji,jj)+zempba_3d_1(ji,jj,jk)*fse3w(ji,jj,jk)
202               zsum2(ji,jj)=zsum2(ji,jj)+zempba_3d_2(ji,jj,jk)*fse3w(ji,jj,jk)               
203            END DO
204         END DO
205      END DO
206
207      z1sum1(:,:)=0.e0
208      z1sum2(:,:)=0.e0
209      DO jj = 1, jpj
210        DO ji = 1, jpi
211         IF (zsum1(ji,jj) /= 0) z1sum1(ji,jj)=1.0/zsum1(ji,jj)
212         IF (zsum2(ji,jj) /= 0) z1sum2(ji,jj)=1.0/zsum2(ji,jj)               
213        END DO
214      END DO
215
216      zsum (:,:)=0.e0
217      z1sum(:,:)=0.e0
218      DO jk= 1, jpk
219         DO jj = 1, jpj
220            DO ji = 1, jpi
221               zempba_3d(ji,jj,jk)=zempba_3d_1(ji,jj,jk)*z1sum1(ji,jj)*0.5*(1-SIGN(1.0,drn2dz(ji,jj,jk)))       &
222                  &               +zempba_3d_2(ji,jj,jk)*z1sum2(ji,jj)*0.5*(1+SIGN(1.0,drn2dz(ji,jj,jk)))
223               zsum(ji,jj)=zsum(ji,jj)+zempba_3d(ji,jj,jk)*fse3w(ji,jj,jk)
224            END DO
225         END DO
226       END DO
227
228       DO jj = 1, jpj
229          DO ji = 1, jpi
230             IF( zsum(ji,jj) > 0 ) z1sum(ji,jj)=1.0/zsum(ji,jj)               
231          END DO
232       END DO
233
234       DO jk= 1, jpk
235          DO jj = 1, jpj
236             DO ji = 1, jpi
237                az_tmx_itf(ji,jj,jk) = en_tmx(ji,jj) * rn_tfe_itf / rn_tfe / rau0    &
238                    &                  * zempba_3d(ji,jj,jk) * z1sum(ji,jj) * tmask(ji,jj,jk) 
239             END DO
240          END DO
241       END DO           
242
243      !                             ! first estimation (with n2 bound by rn_n2min) bounded by 10 cm2/s
244      av_tide_itf(:,:,:) = MIN(  10.e-4, az_tmx_itf(:,:,:) / MAX( rn_n2min, rn2(:,:,:) )  )
245
246      zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column
247      DO jk = 2, jpkm1
248         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * av_tide_itf(:,:,jk)* tmask(:,:,jk)
249      END DO
250
251      DO jj = 1, jpj                ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
252         DO ji = 1, jpi
253            IF( zkz(ji,jj) /= 0.e0 )   zkz(ji,jj) = en_tmx(ji,jj) * rn_tfe_itf / rn_tfe / zkz(ji,jj)
254         END DO
255      END DO
256
257      DO jk = 2, jpkm1              ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> av_tide_itf bound by 300 cm2/s
258         av_tide_itf(:,:,jk) = av_tide_itf(:,:,jk) * MIN( zkz(:,:), 120./10. )   ! kz max = 120 cm2/s
259      END DO
260
261      IF( kt == nit000 ) THEN       ! diagnose the nergy consumed by av_tide_itf
262         ztpc = 0.e0
263         DO jk= 1, jpk
264            DO jj= 1, jpj
265               DO ji= 1, jpi
266                  ztpc = ztpc + fse3w(ji,jj,jk)*e1t(ji,jj)*e2t(ji,jj) * 1/(rn_tfe_itf)  &
267                     &         *MAX(0.e0,rn2(ji,jj,jk))*av_tide_itf(ji,jj,jk)*tmask(ji,jj,jk)*tmask_i(ji,jj)
268               END DO
269            END DO
270         END DO
271         ztpc= rau0 * 1/( rn_me) * ztpc
272         IF(lwp) WRITE(numout,*) '          N Total power consumption by av_tide_itf: ztpc = ', ztpc * 1.e-12 ,'TW'
273      ENDIF
274
275      DO jk = 2, jpkm1
276         av_tide_itf(:,:,jk) = av_tide_itf(:,:,jk) * mask_itf(:,:) 
277         avt(:,:,jk) = avt(:,:,jk) + av_tide_itf(:,:,jk)
278         !
279         DO jj = 2, jpjm1
280            DO ji = fs_2, fs_jpim1  ! vector opt.
281               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( av_tide_itf(ji,jj,jk) + av_tide_itf(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
282               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( av_tide_itf(ji,jj,jk) + av_tide_itf(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
283            END DO
284         END DO
285      END DO
286      CALL lbc_lnk( avmu, 'U', 1. )   
287      CALL lbc_lnk( avmv, 'V', 1. )
288
289      !
290   END SUBROUTINE tmx_itf
291
292
293   SUBROUTINE zdf_tmx_init
294      !!----------------------------------------------------------------------
295      !!                  ***  ROUTINE zdf_tmx_init  ***
296      !!                     
297      !! ** Purpose :   Initialization of the vertical tidal mixing, Reading
298      !!      of M2 and K1 tidal energy in nc files
299      !!
300      !! ** Method  :   Read the namtmx namelist and check the parameters
301      !!      called at the first timestep (nit000),
302      !!              - Read the M2 and K1 tidal energy in NetCDF files.
303      !!              - Compute az_tmx, a 3D coefficient that allows to compute
304      !!      the tidal-induced vertical mixing as follows:
305      !!         Kz_tides = avt0 + az_tmx / max( rn_n2min, N^2 )
306      !!      where az_tmx is given by.....
307      !!
308      !! ** input   :   - Namlist namtmx
309      !!                - NetCDF file  : M2_ORCA2.nc and K1_ORCA2.nc
310      !!
311      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter
312      !!              - defined az_tmx used to compute tidal-induced mixing
313      !!----------------------------------------------------------------------
314      USE iom                                           ! to read input files
315
316      INTEGER ::   ji, jj, jk                           ! dummy loop indices
317      INTEGER ::   inum                                 ! temporary logical unit
318      REAL(wp) ::   ztpc, ze_z                          ! total power consumption
319      REAL(wp), DIMENSION(jpi,jpj) ::  zem2, zek1       ! read M2 and K1 tidal energy
320      REAL(wp), DIMENSION(jpi,jpj) ::  zkz              ! total M2, K1 and S2 tidal energy
321      REAL(wp), DIMENSION(jpi,jpj) ::  zfact            ! used for vertical structure function
322      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zpc          ! power consumption
323
324      !                                   ! * tidal mixing namelist * (namtmx)
325      NAMELIST/namtmx/ rn_htmx, rn_n2min, rn_tfe, rn_tfe_itf, rn_me
326      !!----------------------------------------------------------------------
327
328      IF(lwp) THEN
329         WRITE(numout,*)
330         WRITE(numout,*) 'zdf_tmx_init : tidal mixing'
331         WRITE(numout,*) '~~~~~~~~~~~~'
332      ENDIF
333
334      ! Read Namelist namtmx : Tidal Mixing
335      ! --------------------
336      REWIND ( numnam )
337      READ   ( numnam, namtmx )
338
339      ! Control print
340      IF(lwp) THEN
341         WRITE(numout,*) '          Namelist namtmx : set tidal mixing parameters'
342         WRITE(numout,*) '             Vertical decay scale for turbulence   = ', rn_htmx 
343         WRITE(numout,*) '             Brunt-Vaisala frequency threshold     = ', rn_n2min
344         WRITE(numout,*) '             Tidal dissipation efficiency          = ', rn_tfe
345         WRITE(numout,*) '             ITF tidal dissipation efficiency      = ', rn_tfe_itf
346         WRITE(numout,*) '             Mixing efficiency                     = ', rn_me
347         WRITE(numout,*)
348      ENDIF
349
350      ! read mask kz banda
351      CALL iom_open('mask_itf',inum)
352      CALL iom_get (inum, jpdom_data, 'tmaskitf',mask_itf,1) !
353      CALL iom_close(inum)
354
355      ! read M2 tidal energy flux : W/m2  ( zem2 < 0 )
356      CALL iom_open('M2rowdrg',inum)
357      CALL iom_get (inum, jpdom_data, 'field',zem2,1) !
358      CALL iom_close(inum)
359
360      ! read K1 tidal energy flux : W/m2  ( zek1 < 0 )
361      CALL iom_open('K1rowdrg',inum)
362      CALL iom_get (inum, jpdom_data, 'field',zek1,1) !
363      CALL iom_close(inum)
364 
365      ! Total tidal energy ( M2, S2 and K1  with S2=(1/2)^2 * M2 )
366      ! only the energy available for mixing is taken into account,
367      ! (mixing efficiency tidal dissipation efficiency)
368      en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * tmask(:,:,1)
369
370      ! Bound total energy required to be <= 1 W/m2
371      ! Vertical structure
372      ! part independent of the level
373      DO jj = 1, jpj
374         DO ji = 1, jpi
375            zhdep(ji,jj) = fsdepw(ji,jj,mbathy(ji,jj))       ! depth of the ocean
376            zfact(ji,jj) = rau0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) )
377            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj)
378         END DO
379      END DO
380      DO jk= 1, jpk                 ! complete with the level-dependent part
381         DO jj = 1, jpj
382            DO ji = 1, jpi
383               az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-fsdepw(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk)
384            END DO
385         END DO
386      END DO
387
388      IF( nprint == 1 .AND. lwp ) THEN
389         ! Control print
390         ! Total power consumption due to vertical mixing
391         ! zpc = rau0 * 1/rn_me * rn2 * av_tide
392         av_tide(:,:,:) = 0.e0
393         DO jk = 2, jpkm1
394            av_tide(:,:,jk) = az_tmx(:,:,jk) / MAX( rn_n2min, rn2(:,:,jk) )
395         END DO
396
397         ztpc = 0.e0
398         zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * av_tide(:,:,:)
399         DO jk= 2, jpkm1
400            DO jj = 1, jpj
401               DO ji = 1, jpi
402                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
403               END DO
404            END DO
405         END DO
406         ztpc= rau0 * 1/(rn_tfe * rn_me) * ztpc
407
408         WRITE(numout,*) 
409         WRITE(numout,*) '          Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
410
411
412         ! control print 2
413         av_tide(:,:,:) = MIN( av_tide(:,:,:), 60.e-4 )   
414         zkz(:,:) = 0.e0
415         DO jk = 2, jpkm1
416         DO jj = 1, jpj
417            DO ji = 1, jpi
418               zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * av_tide(ji,jj,jk)* tmask(ji,jj,jk)
419            END DO
420         END DO
421         END DO
422         ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz
423         DO jj = 1, jpj
424            DO ji = 1, jpi
425               IF( zkz(ji,jj) /= 0.e0 )   THEN
426                   zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
427               ENDIF
428            END DO
429         END DO
430         ztpc = 1.e50
431         DO jj = 1, jpj
432            DO ji = 1, jpi
433               IF( zkz(ji,jj) /= 0.e0 )   THEN
434                   ztpc = Min( zkz(ji,jj), ztpc)
435               ENDIF
436            END DO
437         END DO
438         WRITE(numout,*) '          Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) )
439
440         DO jk = 2, jpkm1
441            av_tide(:,:,jk) = av_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s
442         END DO
443         ztpc = 0.e0
444         zpc(:,:,:) = Max(0.e0,rn2(:,:,:)) * av_tide(:,:,:)
445         DO jk= 1, jpk
446            DO jj = 1, jpj
447               DO ji = 1, jpi
448                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
449               END DO
450            END DO
451         END DO
452         ztpc= rau0 * 1/(rn_tfe * rn_me) * ztpc
453         WRITE(numout,*) '          2 Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
454
455         DO jk = 1, jpk
456            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * av_tide(:,:,jk)     * tmask_i(:,:) )   &
457               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
458            ztpc = 1.E50
459            DO jj = 1, jpj
460               DO ji = 1, jpi
461                  IF( av_tide(ji,jj,jk) /= 0.e0 )   ztpc =Min( ztpc, av_tide(ji,jj,jk) )
462               END DO
463            END DO
464            WRITE(numout,*) '            N2 min - jk= ', jk,'   ', ze_z * 1.e4,' cm2/s min= ',ztpc*1.e4,   &
465               &       'max= ', MAXVAL(av_tide(:,:,jk) )*1.e4, ' cm2/s'
466         END DO
467
468         WRITE(numout,*) '          e_tide : ', SUM( e1t*e2t*en_tmx ) / ( rn_tfe * rn_me ) * 1.e-12, 'TW'
469         WRITE(numout,*) 
470         WRITE(numout,*) '          Initial profile of tidal vertical mixing'
471         DO jk = 1, jpk
472            DO jj = 1,jpj
473               DO ji = 1,jpi
474                  zkz(ji,jj) = az_tmx(ji,jj,jk) /MAX( rn_n2min, rn2(ji,jj,jk) )
475               END DO
476            END DO
477            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   &
478               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
479            WRITE(numout,*) '                jk= ', jk,'   ', ze_z * 1.e4,' cm2/s'
480         END DO
481         DO jk = 1, jpk
482            zkz(:,:) = az_tmx(:,:,jk) /rn_n2min
483            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   &
484               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
485            WRITE(numout,*) 
486            WRITE(numout,*) '          N2 min - jk= ', jk,'   ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4,   &
487               &       'max= ', MAXVAL(zkz)*1.e4, ' cm2/s'
488         END DO
489         !
490      ENDIF
491
492   END SUBROUTINE zdf_tmx_init
493
494#else
495   !!----------------------------------------------------------------------
496   !!   Default option          Dummy module                NO Tidal MiXing
497   !!----------------------------------------------------------------------
498   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .FALSE.   !: tidal mixing flag
499CONTAINS
500   SUBROUTINE zdf_tmx( kt )          ! Empty routine
501      WRITE(*,*) 'zdf_tmx: You should not have seen this print! error?', kt
502   END SUBROUTINE zdf_tmx
503#endif
504
505   !!======================================================================
506END MODULE zdftmx
Note: See TracBrowser for help on using the repository browser.