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

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90 @ 4839

Last change on this file since 4839 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.

  • Property svn:keywords set to Id
File size: 26.6 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   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
9   !!----------------------------------------------------------------------
10#if defined key_zdftmx   ||   defined key_esopa
11   !!----------------------------------------------------------------------
12   !!   'key_zdftmx'                                  Tidal vertical mixing
13   !!----------------------------------------------------------------------
14   !!   zdf_tmx       : global     momentum & tracer Kz with tidal induced Kz
15   !!   tmx_itf       : Indonesian momentum & tracer Kz with tidal induced Kz
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 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   USE in_out_manager ! I/O manager
25   USE iom            ! I/O Manager
26   USE lib_mpp        ! MPP library
27   USE wrk_nemo       ! work arrays
28   USE timing         ! Timing
29   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   zdf_tmx         ! called in step module
35   PUBLIC   zdf_tmx_init    ! called in opa module
36   PUBLIC   zdf_tmx_alloc   ! called in nemogcm module
37
38   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: tidal mixing flag
39
40   !                       !!* Namelist  namzdf_tmx : tidal mixing *
41   REAL(wp) ::  rn_htmx     ! vertical decay scale for turbulence (meters)
42   REAL(wp) ::  rn_n2min    ! threshold of the Brunt-Vaisala frequency (s-1)
43   REAL(wp) ::  rn_tfe      ! tidal dissipation efficiency (St Laurent et al. 2002)
44   REAL(wp) ::  rn_me       ! mixing efficiency (Osborn 1980)
45   LOGICAL  ::  ln_tmx_itf  ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization
46   REAL(wp) ::  rn_tfe_itf  ! ITF tidal dissipation efficiency (St Laurent et al. 2002)
47
48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   en_tmx     ! energy available for tidal mixing (W/m2)
49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mask_itf   ! mask to use over Indonesian area
50   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   az_tmx     ! coefficient used to evaluate the tidal induced Kz
51
52   !! * Substitutions
53#  include "domzgr_substitute.h90"
54#  include "vectopt_loop_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
57   !! $Id$
58   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   INTEGER FUNCTION zdf_tmx_alloc()
63      !!----------------------------------------------------------------------
64      !!                ***  FUNCTION zdf_tmx_alloc  ***
65      !!----------------------------------------------------------------------
66      ALLOCATE(en_tmx(jpi,jpj), mask_itf(jpi,jpj), az_tmx(jpi,jpj,jpk), STAT=zdf_tmx_alloc )
67      !
68      IF( lk_mpp             )   CALL mpp_sum ( zdf_tmx_alloc )
69      IF( zdf_tmx_alloc /= 0 )   CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays')
70   END FUNCTION zdf_tmx_alloc
71
72
73   SUBROUTINE zdf_tmx( kt )
74      !!----------------------------------------------------------------------
75      !!                  ***  ROUTINE zdf_tmx  ***
76      !!                   
77      !! ** Purpose :   add to the vertical mixing coefficients the effect of
78      !!              tidal mixing (Simmons et al 2004).
79      !!
80      !! ** Method  : - tidal-induced vertical mixing is given by:
81      !!                  Kz_tides = az_tmx / max( rn_n2min, N^2 )
82      !!              where az_tmx is a coefficient that specified the 3D space
83      !!              distribution of the faction of tidal energy taht is used
84      !!              for mixing. Its expression is set in zdf_tmx_init routine,
85      !!              following Simmons et al. 2004.
86      !!                NB: a specific bounding procedure is performed on av_tide
87      !!              so that the input tidal energy is actually almost used. The
88      !!              basic maximum value is 60 cm2/s, but values of 300 cm2/s
89      !!              can be reached in area where bottom stratification is too
90      !!              weak.
91      !!
92      !!              - update av_tide in the Indonesian Through Flow area
93      !!              following Koch-Larrouy et al. (2007) parameterisation
94      !!              (see tmx_itf routine).
95      !!
96      !!              - update the model vertical eddy viscosity and diffusivity:
97      !!                     avt  = avt  +    av_tides
98      !!                     avm  = avm  +    av_tides
99      !!                     avmu = avmu + mi(av_tides)
100      !!                     avmv = avmv + mj(av_tides)
101      !!
102      !! ** Action  :   avt, avm, avmu, avmv   increased by tidal mixing
103      !!
104      !! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
105      !!              Koch-Larrouy et al. 2007, GRL.
106      !!----------------------------------------------------------------------
107      USE oce, zav_tide  =>   ua    ! use ua as workspace
108      !!
109      INTEGER, INTENT(in) ::   kt   ! ocean time-step
110      !!
111      INTEGER  ::   ji, jj, jk   ! dummy loop indices
112      REAL(wp) ::   ztpc         ! scalar workspace
113      REAL(wp), POINTER, DIMENSION(:,:) ::   zkz
114      !!----------------------------------------------------------------------
115      !
116      IF( nn_timing == 1 )  CALL timing_start('zdf_tmx')
117      !
118      CALL wrk_alloc( jpi,jpj, zkz )
119
120      !                          ! ----------------------- !
121      !                          !  Standard tidal mixing  !  (compute zav_tide)
122      !                          ! ----------------------- !
123      !                             !* First estimation (with n2 bound by rn_n2min) bounded by 60 cm2/s
124      zav_tide(:,:,:) = MIN(  60.e-4, az_tmx(:,:,:) / MAX( rn_n2min, rn2(:,:,:) )  )
125
126      zkz(:,:) = 0.e0               !* Associated potential energy consummed over the whole water column
127      DO jk = 2, jpkm1
128         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk)
129      END DO
130
131      DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
132         DO ji = 1, jpi
133            IF( zkz(ji,jj) /= 0.e0 )   zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
134         END DO
135      END DO
136
137      DO jk = 2, jpkm1              !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
138         zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s
139      END DO
140
141      IF( kt == nit000 ) THEN       !* check at first time-step: diagnose the energy consumed by zav_tide
142         ztpc = 0.e0
143         DO jk= 1, jpk
144            DO jj= 1, jpj
145               DO ji= 1, jpi
146                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj)   &
147                     &         * MAX( 0.e0, rn2(ji,jj,jk) ) * zav_tide(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
148               END DO
149            END DO
150         END DO
151         ztpc= rau0 / ( rn_tfe * rn_me ) * ztpc
152         IF(lwp) WRITE(numout,*) 
153         IF(lwp) WRITE(numout,*) '          N Total power consumption by av_tide    : ztpc = ', ztpc * 1.e-12 ,'TW'
154      ENDIF
155       
156      !                          ! ----------------------- !
157      !                          !    ITF  tidal mixing    !  (update zav_tide)
158      !                          ! ----------------------- !
159      IF( ln_tmx_itf )   CALL tmx_itf( kt, zav_tide )
160
161      !                          ! ----------------------- !
162      !                          !   Update  mixing coefs  !                         
163      !                          ! ----------------------- !
164      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
165         avt(:,:,jk) = avt(:,:,jk) + zav_tide(:,:,jk)
166         avm(:,:,jk) = avm(:,:,jk) + zav_tide(:,:,jk)
167         DO jj = 2, jpjm1
168            DO ji = fs_2, fs_jpim1  ! vector opt.
169               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
170               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
171            END DO
172         END DO
173      END DO
174      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! lateral boundary condition
175
176      !                             !* output tidal mixing coefficient
177      CALL iom_put( "av_tide", zav_tide )
178
179      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_tide , clinfo1=' tmx - av_tide: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
180      !
181      CALL wrk_dealloc( jpi,jpj, zkz )
182      !
183      IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx')
184      !
185   END SUBROUTINE zdf_tmx
186
187
188   SUBROUTINE tmx_itf( kt, pav )
189      !!----------------------------------------------------------------------
190      !!                  ***  ROUTINE tmx_itf  ***
191      !!                   
192      !! ** Purpose :   modify the vertical eddy diffusivity coefficients
193      !!              (pav) in the Indonesian Through Flow area (ITF).
194      !!
195      !! ** Method  : - Following Koch-Larrouy et al. (2007), in the ITF defined
196      !!                by msk_itf (read in a file, see tmx_init), the tidal
197      !!                mixing coefficient is computed with :
198      !!                  * q=1 (i.e. all the tidal energy remains trapped in
199      !!                         the area and thus is used for mixing)
200      !!                  * the vertical distribution of the tifal energy is a
201      !!                    proportional to N above the thermocline (d(N^2)/dz > 0)
202      !!                    and to N^2 below the thermocline (d(N^2)/dz < 0)
203      !!
204      !! ** Action  :   av_tide   updated in the ITF area (msk_itf)
205      !!
206      !! References :  Koch-Larrouy et al. 2007, GRL
207      !!----------------------------------------------------------------------
208      INTEGER , INTENT(in   )                         ::   kt   ! ocean time-step
209      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pav  ! Tidal mixing coef.
210      !!
211      INTEGER  ::   ji, jj, jk    ! dummy loop indices
212      REAL(wp) ::   zcoef, ztpc   ! temporary scalar
213      REAL(wp), DIMENSION(:,:)  , POINTER ::   zkz                        ! 2D workspace
214      REAL(wp), DIMENSION(:,:)  , POINTER ::   zsum1 , zsum2 , zsum       !  -      -
215      REAL(wp), DIMENSION(:,:,:), POINTER ::   zempba_3d_1, zempba_3d_2   ! 3D workspace
216      REAL(wp), DIMENSION(:,:,:), POINTER ::   zempba_3d  , zdn2dz        !  -      -
217      REAL(wp), DIMENSION(:,:,:), POINTER ::   zavt_itf                   !  -      -
218      !!----------------------------------------------------------------------
219      !
220      IF( nn_timing == 1 )  CALL timing_start('tmx_itf')
221      !
222      CALL wrk_alloc( jpi,jpj, zkz, zsum1 , zsum2 , zsum )
223      CALL wrk_alloc( jpi,jpj,jpk, zempba_3d_1, zempba_3d_2, zempba_3d, zdn2dz, zavt_itf )
224
225      !                             ! compute the form function using N2 at each time step
226      zempba_3d_1(:,:,jpk) = 0.e0
227      zempba_3d_2(:,:,jpk) = 0.e0
228      DO jk = 1, jpkm1             
229         zdn2dz     (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1)           ! Vertical profile of dN2/dz
230!CDIR NOVERRCHK
231         zempba_3d_1(:,:,jk) = SQRT(  MAX( 0.e0, rn2(:,:,jk) )  )    !    -        -    of N
232         zempba_3d_2(:,:,jk) =        MAX( 0.e0, rn2(:,:,jk) )       !    -        -    of N^2
233      END DO
234      !
235      zsum (:,:) = 0.e0
236      zsum1(:,:) = 0.e0
237      zsum2(:,:) = 0.e0
238      DO jk= 2, jpk
239         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk)
240         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk)               
241      END DO
242      DO jj = 1, jpj
243         DO ji = 1, jpi
244            IF( zsum1(ji,jj) /= 0.e0 )   zsum1(ji,jj) = 1.e0 / zsum1(ji,jj)
245            IF( zsum2(ji,jj) /= 0.e0 )   zsum2(ji,jj) = 1.e0 / zsum2(ji,jj)               
246         END DO
247      END DO
248
249      DO jk= 1, jpk
250         DO jj = 1, jpj
251            DO ji = 1, jpi
252               zcoef = 0.5 - SIGN( 0.5, zdn2dz(ji,jj,jk) )       ! =0 if dN2/dz > 0, =1 otherwise
253               ztpc  = zempba_3d_1(ji,jj,jk) * zsum1(ji,jj) *        zcoef     &
254                  &  + zempba_3d_2(ji,jj,jk) * zsum2(ji,jj) * ( 1. - zcoef )
255               !
256               zempba_3d(ji,jj,jk) =               ztpc 
257               zsum     (ji,jj)    = zsum(ji,jj) + ztpc * fse3w(ji,jj,jk)
258            END DO
259         END DO
260       END DO
261       DO jj = 1, jpj
262          DO ji = 1, jpi
263             IF( zsum(ji,jj) > 0.e0 )   zsum(ji,jj) = 1.e0 / zsum(ji,jj)               
264          END DO
265       END DO
266
267      !                             ! first estimation bounded by 10 cm2/s (with n2 bounded by rn_n2min)
268      zcoef = rn_tfe_itf / ( rn_tfe * rau0 )
269      DO jk = 1, jpk
270         zavt_itf(:,:,jk) = MIN(  10.e-4, zcoef * en_tmx(:,:) * zsum(:,:) * zempba_3d(:,:,jk)   &
271            &                                      / MAX( rn_n2min, rn2(:,:,jk) ) * tmask(:,:,jk)  )
272      END DO           
273
274      zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column
275      DO jk = 2, jpkm1
276         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk)
277      END DO
278
279      DO jj = 1, jpj                ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
280         DO ji = 1, jpi
281            IF( zkz(ji,jj) /= 0.e0 )   zkz(ji,jj) = en_tmx(ji,jj) * rn_tfe_itf / rn_tfe / zkz(ji,jj)
282         END DO
283      END DO
284
285      DO jk = 2, jpkm1              ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s
286         zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. )   ! kz max = 120 cm2/s
287      END DO
288
289      IF( kt == nit000 ) THEN       ! diagnose the nergy consumed by zavt_itf
290         ztpc = 0.e0
291         DO jk= 1, jpk
292            DO jj= 1, jpj
293               DO ji= 1, jpi
294                  ztpc = ztpc + e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) )   &
295                     &                     * zavt_itf(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
296               END DO
297            END DO
298         END DO
299         ztpc= rau0 * ztpc / ( rn_me * rn_tfe_itf )
300         IF(lwp) WRITE(numout,*) '          N Total power consumption by zavt_itf: ztpc = ', ztpc * 1.e-12 ,'TW'
301      ENDIF
302
303      !                             ! Update pav with the ITF mixing coefficient
304      DO jk = 2, jpkm1
305         pav(:,:,jk) = pav     (:,:,jk) * ( 1.e0 - mask_itf(:,:) )   &
306            &        + zavt_itf(:,:,jk) *          mask_itf(:,:) 
307      END DO
308      !
309      CALL wrk_dealloc( jpi,jpj, zkz, zsum1 , zsum2 , zsum )
310      CALL wrk_dealloc( jpi,jpj,jpk, zempba_3d_1, zempba_3d_2, zempba_3d, zdn2dz, zavt_itf )
311      !
312      IF( nn_timing == 1 )  CALL timing_stop('tmx_itf')
313      !
314   END SUBROUTINE tmx_itf
315
316
317   SUBROUTINE zdf_tmx_init
318      !!----------------------------------------------------------------------
319      !!                  ***  ROUTINE zdf_tmx_init  ***
320      !!                     
321      !! ** Purpose :   Initialization of the vertical tidal mixing, Reading
322      !!              of M2 and K1 tidal energy in nc files
323      !!
324      !! ** Method  : - Read the namtmx namelist and check the parameters
325      !!
326      !!              - Read the input data in NetCDF files :
327      !!              M2 and K1 tidal energy. The total tidal energy, en_tmx,
328      !!              is the sum of M2, K1 and S2 energy where S2 is assumed
329      !!              to be: S2=(1/2)^2 * M2
330      !!              mask_itf, a mask array that determine where substituing
331      !!              the standard Simmons et al. (2005) formulation with the
332      !!              one of Koch_Larrouy et al. (2007).
333      !!
334      !!              - Compute az_tmx, a 3D coefficient that allows to compute
335      !!             the standard tidal-induced vertical mixing as follows:
336      !!                  Kz_tides = az_tmx / max( rn_n2min, N^2 )
337      !!             with az_tmx a bottom intensified coefficient is given by:
338      !!                 az_tmx(z) = en_tmx / ( rau0 * rn_htmx ) * EXP( -(H-z)/rn_htmx )
339      !!                                                  / ( 1. - EXP( - H   /rn_htmx ) )
340      !!             where rn_htmx the characteristic length scale of the bottom
341      !!             intensification, en_tmx the tidal energy, and H the ocean depth
342      !!
343      !! ** input   :   - Namlist namtmx
344      !!                - NetCDF file : M2_ORCA2.nc, K1_ORCA2.nc, and mask_itf.nc
345      !!
346      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter
347      !!              - defined az_tmx used to compute tidal-induced mixing
348      !!
349      !! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
350      !!              Koch-Larrouy et al. 2007, GRL.
351      !!----------------------------------------------------------------------
352      USE oce     ,         zav_tide =>  ua         ! ua used as workspace
353      !!
354      INTEGER  ::   ji, jj, jk   ! dummy loop indices
355      INTEGER  ::   inum         ! local integer
356      INTEGER  ::   ios
357      REAL(wp) ::   ztpc, ze_z   ! local scalars
358      REAL(wp), DIMENSION(:,:)  , POINTER ::  zem2, zek1   ! read M2 and K1 tidal energy
359      REAL(wp), DIMENSION(:,:)  , POINTER ::  zkz          ! total M2, K1 and S2 tidal energy
360      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact        ! used for vertical structure function
361      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep        ! Ocean depth
362      REAL(wp), DIMENSION(:,:,:), POINTER ::  zpc      ! power consumption
363      !!
364      NAMELIST/namzdf_tmx/ rn_htmx, rn_n2min, rn_tfe, rn_me, ln_tmx_itf, rn_tfe_itf
365      !!----------------------------------------------------------------------
366      !
367      IF( nn_timing == 1 )  CALL timing_start('zdf_tmx_init')
368      !
369      CALL wrk_alloc( jpi,jpj, zem2, zek1, zkz, zfact, zhdep )
370      CALL wrk_alloc( jpi,jpj,jpk, zpc )
371     
372      REWIND( numnam_ref )              ! Namelist namzdf_tmx in reference namelist : Tidal Mixing
373      READ  ( numnam_ref, namzdf_tmx, IOSTAT = ios, ERR = 901)
374901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp )
375
376      REWIND( numnam_cfg )              ! Namelist namzdf_tmx in configuration namelist : Tidal Mixing
377      READ  ( numnam_cfg, namzdf_tmx, IOSTAT = ios, ERR = 902 )
378902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp )
379      IF(lwm) WRITE ( numond, namzdf_tmx )
380
381      IF(lwp) THEN                   ! Control print
382         WRITE(numout,*)
383         WRITE(numout,*) 'zdf_tmx_init : tidal mixing'
384         WRITE(numout,*) '~~~~~~~~~~~~'
385         WRITE(numout,*) '   Namelist namzdf_tmx : set tidal mixing parameters'
386         WRITE(numout,*) '      Vertical decay scale for turbulence   = ', rn_htmx 
387         WRITE(numout,*) '      Brunt-Vaisala frequency threshold     = ', rn_n2min
388         WRITE(numout,*) '      Tidal dissipation efficiency          = ', rn_tfe
389         WRITE(numout,*) '      Mixing efficiency                     = ', rn_me
390         WRITE(numout,*) '      ITF specific parameterisation         = ', ln_tmx_itf
391         WRITE(numout,*) '      ITF tidal dissipation efficiency      = ', rn_tfe_itf
392      ENDIF
393
394      !                              ! allocate tmx arrays
395      IF( zdf_tmx_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' )
396
397      IF( ln_tmx_itf ) THEN          ! read the Indonesian Through Flow mask
398         CALL iom_open('mask_itf',inum)
399         CALL iom_get (inum, jpdom_data, 'tmaskitf',mask_itf,1) !
400         CALL iom_close(inum)
401      ENDIF
402
403      ! read M2 tidal energy flux : W/m2  ( zem2 < 0 )
404      CALL iom_open('M2rowdrg',inum)
405      CALL iom_get (inum, jpdom_data, 'field',zem2,1) !
406      CALL iom_close(inum)
407
408      ! read K1 tidal energy flux : W/m2  ( zek1 < 0 )
409      CALL iom_open('K1rowdrg',inum)
410      CALL iom_get (inum, jpdom_data, 'field',zek1,1) !
411      CALL iom_close(inum)
412 
413      ! Total tidal energy ( M2, S2 and K1  with S2=(1/2)^2 * M2 )
414      ! only the energy available for mixing is taken into account,
415      ! (mixing efficiency tidal dissipation efficiency)
416      en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * tmask(:,:,1)
417
418      ! Vertical structure (az_tmx)
419      DO jj = 1, jpj                ! part independent of the level
420         DO ji = 1, jpi
421            zhdep(ji,jj) = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
422            zfact(ji,jj) = rau0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) )
423            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj)
424         END DO
425      END DO
426      DO jk= 1, jpk                 ! complete with the level-dependent part
427         DO jj = 1, jpj
428            DO ji = 1, jpi
429               az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-fsdepw(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk)
430            END DO
431         END DO
432      END DO
433
434      IF( nprint == 1 .AND. lwp ) THEN
435         ! Control print
436         ! Total power consumption due to vertical mixing
437         ! zpc = rau0 * 1/rn_me * rn2 * zav_tide
438         zav_tide(:,:,:) = 0.e0
439         DO jk = 2, jpkm1
440            zav_tide(:,:,jk) = az_tmx(:,:,jk) / MAX( rn_n2min, rn2(:,:,jk) )
441         END DO
442
443         ztpc = 0.e0
444         zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:)
445         DO jk= 2, jpkm1
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
454         WRITE(numout,*) 
455         WRITE(numout,*) '          Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
456
457
458         ! control print 2
459         zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )   
460         zkz(:,:) = 0.e0
461         DO jk = 2, jpkm1
462         DO jj = 1, jpj
463            DO ji = 1, jpi
464               zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* tmask(ji,jj,jk)
465            END DO
466         END DO
467         END DO
468         ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz
469         DO jj = 1, jpj
470            DO ji = 1, jpi
471               IF( zkz(ji,jj) /= 0.e0 )   THEN
472                   zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
473               ENDIF
474            END DO
475         END DO
476         ztpc = 1.e50
477         DO jj = 1, jpj
478            DO ji = 1, jpi
479               IF( zkz(ji,jj) /= 0.e0 )   THEN
480                   ztpc = Min( zkz(ji,jj), ztpc)
481               ENDIF
482            END DO
483         END DO
484         WRITE(numout,*) '          Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) )
485
486         DO jk = 2, jpkm1
487            zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s
488         END DO
489         ztpc = 0.e0
490         zpc(:,:,:) = Max(0.e0,rn2(:,:,:)) * zav_tide(:,:,:)
491         DO jk= 1, jpk
492            DO jj = 1, jpj
493               DO ji = 1, jpi
494                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
495               END DO
496            END DO
497         END DO
498         ztpc= rau0 * 1/(rn_tfe * rn_me) * ztpc
499         WRITE(numout,*) '          2 Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
500
501         DO jk = 1, jpk
502            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zav_tide(:,:,jk)     * tmask_i(:,:) )   &
503               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
504            ztpc = 1.E50
505            DO jj = 1, jpj
506               DO ji = 1, jpi
507                  IF( zav_tide(ji,jj,jk) /= 0.e0 )   ztpc =Min( ztpc, zav_tide(ji,jj,jk) )
508               END DO
509            END DO
510            WRITE(numout,*) '            N2 min - jk= ', jk,'   ', ze_z * 1.e4,' cm2/s min= ',ztpc*1.e4,   &
511               &       'max= ', MAXVAL(zav_tide(:,:,jk) )*1.e4, ' cm2/s'
512         END DO
513
514         WRITE(numout,*) '          e_tide : ', SUM( e1t*e2t*en_tmx ) / ( rn_tfe * rn_me ) * 1.e-12, 'TW'
515         WRITE(numout,*) 
516         WRITE(numout,*) '          Initial profile of tidal vertical mixing'
517         DO jk = 1, jpk
518            DO jj = 1,jpj
519               DO ji = 1,jpi
520                  zkz(ji,jj) = az_tmx(ji,jj,jk) /MAX( rn_n2min, rn2(ji,jj,jk) )
521               END DO
522            END DO
523            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   &
524               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
525            WRITE(numout,*) '                jk= ', jk,'   ', ze_z * 1.e4,' cm2/s'
526         END DO
527         DO jk = 1, jpk
528            zkz(:,:) = az_tmx(:,:,jk) /rn_n2min
529            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   &
530               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
531            WRITE(numout,*) 
532            WRITE(numout,*) '          N2 min - jk= ', jk,'   ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4,   &
533               &       'max= ', MAXVAL(zkz)*1.e4, ' cm2/s'
534         END DO
535         !
536      ENDIF
537      !
538      CALL wrk_dealloc( jpi,jpj, zem2, zek1, zkz, zfact, zhdep )
539      CALL wrk_dealloc( jpi,jpj,jpk, zpc )
540      !
541      IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx_init')
542      !
543   END SUBROUTINE zdf_tmx_init
544
545#else
546   !!----------------------------------------------------------------------
547   !!   Default option          Dummy module                NO Tidal MiXing
548   !!----------------------------------------------------------------------
549   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .FALSE.   !: tidal mixing flag
550CONTAINS
551   SUBROUTINE zdf_tmx_init           ! Dummy routine
552      WRITE(*,*) 'zdf_tmx: You should not have seen this print! error?'
553   END SUBROUTINE zdf_tmx_init
554   SUBROUTINE zdf_tmx( kt )          ! Dummy routine
555      WRITE(*,*) 'zdf_tmx: You should not have seen this print! error?', kt
556   END SUBROUTINE zdf_tmx
557#endif
558
559   !!======================================================================
560END MODULE zdftmx
Note: See TracBrowser for help on using the repository browser.