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.
zdfiwm.F90 in NEMO/trunk/tests/BENCH/MY_SRC – NEMO

source: NEMO/trunk/tests/BENCH/MY_SRC/zdfiwm.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
File size: 25.2 KB
Line 
1MODULE zdfiwm
2   !!========================================================================
3   !!                       ***  MODULE  zdfiwm  ***
4   !! Ocean physics: Internal gravity wave-driven vertical mixing
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   !!            3.6  !  2016-03  (C. de Lavergne)  New param: internal wave-driven mixing
10   !!            4.0  !  2017-04  (G. Madec)  renamed module, remove the old param. and the CPP keys
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   zdf_iwm       : global     momentum & tracer Kz with wave induced Kz
15   !!   zdf_iwm_init  : global     momentum & tracer Kz with wave 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 zdfddm         ! ocean vertical physics: double diffusive mixing
21   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
22   USE eosbn2         ! ocean equation of state
23   USE phycst         ! physical constants
24   !
25   USE prtctl         ! Print control
26   USE in_out_manager ! I/O manager
27   USE iom            ! I/O Manager
28   USE lib_mpp        ! MPP library
29   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   zdf_iwm        ! called in step module
35   PUBLIC   zdf_iwm_init   ! called in nemogcm module
36
37   !                      !!* Namelist  namzdf_iwm : internal wave-driven mixing *
38   INTEGER ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2)
39   LOGICAL ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency
40   LOGICAL ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F)
41
42   REAL(wp)::  r1_6 = 1._wp / 6._wp
43
44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ebot_iwm   ! power available from high-mode wave breaking (W/m2)
45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   epyc_iwm   ! power available from low-mode, pycnocline-intensified wave breaking (W/m2)
46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ecri_iwm   ! power available from low-mode, critical slope wave breaking (W/m2)
47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbot_iwm   ! WKB decay scale for high-mode energy dissipation (m)
48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hcri_iwm   ! decay scale for low-mode critical slope dissipation (m)
49
50   !!----------------------------------------------------------------------
51   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
52   !! $Id: zdfiwm.F90 8093 2017-05-30 08:13:14Z gm $
53   !! Software governed by the CeCILL licence     (./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   INTEGER FUNCTION zdf_iwm_alloc()
58      !!----------------------------------------------------------------------
59      !!                ***  FUNCTION zdf_iwm_alloc  ***
60      !!----------------------------------------------------------------------
61      ALLOCATE( ebot_iwm(jpi,jpj),  epyc_iwm(jpi,jpj),  ecri_iwm(jpi,jpj) ,     &
62      &         hbot_iwm(jpi,jpj),  hcri_iwm(jpi,jpj)                     , STAT=zdf_iwm_alloc )
63      !
64      CALL mpp_sum ( 'zdfiwm', zdf_iwm_alloc )
65      IF( zdf_iwm_alloc /= 0 )   CALL ctl_stop( 'STOP', 'zdf_iwm_alloc: failed to allocate arrays' )
66   END FUNCTION zdf_iwm_alloc
67
68
69   SUBROUTINE zdf_iwm( kt, Kmm, p_avm, p_avt, p_avs )
70      !!----------------------------------------------------------------------
71      !!                  ***  ROUTINE zdf_iwm  ***
72      !!                   
73      !! ** Purpose :   add to the vertical mixing coefficients the effect of
74      !!              breaking internal waves.
75      !!
76      !! ** Method  : - internal wave-driven vertical mixing is given by:
77      !!                  Kz_wave = min(  100 cm2/s, f(  Reb = zemx_iwm /( Nu * N^2 )  )
78      !!              where zemx_iwm is the 3D space distribution of the wave-breaking
79      !!              energy and Nu the molecular kinematic viscosity.
80      !!              The function f(Reb) is linear (constant mixing efficiency)
81      !!              if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T.
82      !!
83      !!              - Compute zemx_iwm, the 3D power density that allows to compute
84      !!              Reb and therefrom the wave-induced vertical diffusivity.
85      !!              This is divided into three components:
86      !!                 1. Bottom-intensified low-mode dissipation at critical slopes
87      !!                     zemx_iwm(z) = ( ecri_iwm / rau0 ) * EXP( -(H-z)/hcri_iwm )
88      !!                                   / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm
89      !!              where hcri_iwm is the characteristic length scale of the bottom
90      !!              intensification, ecri_iwm a map of available power, and H the ocean depth.
91      !!                 2. Pycnocline-intensified low-mode dissipation
92      !!                     zemx_iwm(z) = ( epyc_iwm / rau0 ) * ( sqrt(rn2(z))^nn_zpyc )
93      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) )
94      !!              where epyc_iwm is a map of available power, and nn_zpyc
95      !!              is the chosen stratification-dependence of the internal wave
96      !!              energy dissipation.
97      !!                 3. WKB-height dependent high mode dissipation
98      !!                     zemx_iwm(z) = ( ebot_iwm / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm)
99      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) )
100      !!              where hbot_iwm is the characteristic length scale of the WKB bottom
101      !!              intensification, ebot_iwm is a map of available power, and z_wkb is the
102      !!              WKB-stretched height above bottom defined as
103      !!                    z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) )
104      !!                                 / SUM( sqrt(rn2(z'))    * e3w(z')    )
105      !!
106      !!              - update the model vertical eddy viscosity and diffusivity:
107      !!                     avt  = avt  +    av_wave
108      !!                     avm  = avm  +    av_wave
109      !!
110      !!              - if namelist parameter ln_tsdiff = T, account for differential mixing:
111      !!                     avs  = avt  +    av_wave * diffusivity_ratio(Reb)
112      !!
113      !! ** Action  : - avt, avs, avm, increased by tide internal wave-driven mixing   
114      !!
115      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep.
116      !!----------------------------------------------------------------------
117      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step
118      INTEGER                    , INTENT(in   ) ::   Kmm            ! time level index
119      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm          ! momentum Kz (w-points)
120      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avt, p_avs   ! tracer   Kz (w-points)
121      !
122      INTEGER  ::   ji, jj, jk   ! dummy loop indices
123      REAL(wp) ::   zztmp        ! scalar workspace
124      REAL(wp), DIMENSION(jpi,jpj)     ::   zfact       ! Used for vertical structure
125      REAL(wp), DIMENSION(jpi,jpj)     ::   zhdep       ! Ocean depth
126      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwkb        ! WKB-stretched height above bottom
127      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zweight     ! Weight for high mode vertical distribution
128      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znu_t       ! Molecular kinematic viscosity (T grid)
129      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znu_w       ! Molecular kinematic viscosity (W grid)
130      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zReb        ! Turbulence intensity parameter
131      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zemx_iwm    ! local energy density available for mixing (W/kg)
132      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zav_ratio   ! S/T diffusivity ratio (only for ln_tsdiff=T)
133      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zav_wave    ! Internal wave-induced diffusivity
134      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z3d  ! 3D workspace used for iom_put
135      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   z2d  ! 2D     -      -    -     -
136      !!----------------------------------------------------------------------
137      !
138      !                       !* Set to zero the 1st and last vertical levels of appropriate variables
139      zemx_iwm (:,:,1) = 0._wp   ;   zemx_iwm (:,:,jpk) = 0._wp
140      zav_ratio(:,:,1) = 0._wp   ;   zav_ratio(:,:,jpk) = 0._wp
141      zav_wave (:,:,1) = 0._wp   ;   zav_wave (:,:,jpk) = 0._wp
142      !
143      !                       ! ----------------------------- !
144      !                       !  Internal wave-driven mixing  !  (compute zav_wave)
145      !                       ! ----------------------------- !
146      !                             
147      !                       !* Critical slope mixing: distribute energy over the time-varying ocean depth,
148      !                                                 using an exponential decay from the seafloor.
149      DO jj = 1, jpj                ! part independent of the level
150         DO ji = 1, jpi
151            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
152            zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  )
153            IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj)
154         END DO
155      END DO
156!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(Kmm) - ssh(Kmm)
157      DO jk = 2, jpkm1              ! complete with the level-dependent part
158         zemx_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      &
159            &                             - EXP( ( gde3w(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) )  ) * wmask(:,:,jk)   &
160            &                          / ( gde3w(:,:,jk) - gde3w(:,:,jk-1) )
161
162!!gm delta(gde3w) = e3t(Kmm)  !!  Please verify the grid-point position w versus t-point
163!!gm it seems to me that only 1/hcri_iwm  is used ==>  compute it one for all
164
165      END DO
166
167      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying
168      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc
169      !                                          ! (NB: N2 is masked, so no use of wmask here)
170      SELECT CASE ( nn_zpyc )
171      !
172      CASE ( 1 )               ! Dissipation scales as N (recommended)
173         !
174         zfact(:,:) = 0._wp
175         DO jk = 2, jpkm1              ! part independent of the level
176            zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
177         END DO
178         !
179         DO jj = 1, jpj
180            DO ji = 1, jpi
181               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
182            END DO
183         END DO
184         !
185         DO jk = 2, jpkm1              ! complete with the level-dependent part
186            zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
187         END DO
188         !
189      CASE ( 2 )               ! Dissipation scales as N^2
190         !
191         zfact(:,:) = 0._wp
192         DO jk = 2, jpkm1              ! part independent of the level
193            zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
194         END DO
195         !
196         DO jj= 1, jpj
197            DO ji = 1, jpi
198               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
199            END DO
200         END DO
201         !
202         DO jk = 2, jpkm1              ! complete with the level-dependent part
203            zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
204         END DO
205         !
206      END SELECT
207
208      !                        !* WKB-height dependent mixing: distribute energy over the time-varying
209      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot)
210      !
211      zwkb (:,:,:) = 0._wp
212      zfact(:,:)   = 0._wp
213      DO jk = 2, jpkm1
214         zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
215         zwkb(:,:,jk) = zfact(:,:)
216      END DO
217!!gm even better:
218!      DO jk = 2, jpkm1
219!         zwkb(:,:) = zwkb(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  )
220!      END DO
221!      zfact(:,:) = zwkb(:,:,jpkm1)
222!!gm or just use zwkb(k=jpk-1) instead of zfact...
223!!gm
224      !
225      DO jk = 2, jpkm1
226         DO jj = 1, jpj
227            DO ji = 1, jpi
228               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
229                  &                                     * wmask(ji,jj,jk) / zfact(ji,jj)
230            END DO
231         END DO
232      END DO
233      zwkb(:,:,1) = zhdep(:,:) * wmask(:,:,1)
234      !
235      zweight(:,:,:) = 0._wp
236      DO jk = 2, jpkm1
237         zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_iwm(:,:) * wmask(:,:,jk)                    &
238            &   * (  EXP( -zwkb(:,:,jk) / hbot_iwm(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_iwm(:,:) )  )
239      END DO
240      !
241      zfact(:,:) = 0._wp
242      DO jk = 2, jpkm1              ! part independent of the level
243         zfact(:,:) = zfact(:,:) + zweight(:,:,jk)
244      END DO
245      !
246      DO jj = 1, jpj
247         DO ji = 1, jpi
248            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
249         END DO
250      END DO
251      !
252      DO jk = 2, jpkm1              ! complete with the level-dependent part
253         zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   &
254            &                                / ( gde3w(:,:,jk) - gde3w(:,:,jk-1) )
255!!gm  use of e3t(:,:,:,Kmm) just above?
256      END DO
257      !
258!!gm  this is to be replaced by just a constant value znu=1.e-6 m2/s
259      ! Calculate molecular kinematic viscosity
260      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * ts(:,:,:,jp_tem,Kmm) + 0.00694_wp * ts(:,:,:,jp_tem,Kmm) * ts(:,:,:,jp_tem,Kmm)  &
261         &                                  + 0.02305_wp * ts(:,:,:,jp_sal,Kmm)  ) * tmask(:,:,:) * r1_rau0
262      DO jk = 2, jpkm1
263         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk)
264      END DO
265!!gm end
266      !
267      ! Calculate turbulence intensity parameter Reb
268      DO jk = 2, jpkm1
269         zReb(:,:,jk) = zemx_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) )
270      END DO
271      !
272      ! Define internal wave-induced diffusivity
273      DO jk = 2, jpkm1
274         zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6
275      END DO
276      !
277      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the
278         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
279            DO jj = 1, jpj
280               DO ji = 1, jpi
281                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
282                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
283                  ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
284                     zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
285                  ENDIF
286               END DO
287            END DO
288         END DO
289      ENDIF
290      !
291      DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s
292         zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk)
293      END DO
294      !
295      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave
296         zztmp = 0._wp
297!!gm used of glosum 3D....
298         DO jk = 2, jpkm1
299            DO jj = 1, jpj
300               DO ji = 1, jpi
301                  zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj)   &
302                     &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
303               END DO
304            END DO
305         END DO
306         CALL mpp_sum( 'zdfiwm', zztmp )
307         zztmp = rau0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing
308         !
309         IF(lwp) THEN
310            WRITE(numout,*)
311            WRITE(numout,*) 'zdf_iwm : Internal wave-driven mixing (iwm)'
312            WRITE(numout,*) '~~~~~~~ '
313            WRITE(numout,*)
314            WRITE(numout,*) '      Total power consumption by av_wave =  ', zztmp * 1.e-12_wp, 'TW'
315         ENDIF
316      ENDIF
317
318      !                          ! ----------------------- !
319      !                          !   Update  mixing coefs  !                         
320      !                          ! ----------------------- !
321      !     
322      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature
323         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb
324            DO jj = 1, jpj
325               DO ji = 1, jpi
326                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  &
327                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   &
328                      &                 ) * wmask(ji,jj,jk)
329               END DO
330            END DO
331         END DO
332         CALL iom_put( "av_ratio", zav_ratio )
333         DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing
334            p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk)
335            p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk)
336            p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk)
337         END DO
338         !
339      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing
340         DO jk = 2, jpkm1
341            p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk)
342            p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk)
343            p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk)
344         END DO
345      ENDIF
346
347      !                             !* output internal wave-driven mixing coefficient
348      CALL iom_put( "av_wave", zav_wave )
349                                    !* output useful diagnostics: Kz*N^2 ,
350!!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5)
351                                    !  vertical integral of rau0 * Kz * N^2 , energy density (zemx_iwm)
352      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
353         ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) )
354         z3d(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:)
355         z2d(:,:) = 0._wp
356         DO jk = 2, jpkm1
357            z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)
358         END DO
359         z2d(:,:) = rau0 * z2d(:,:)
360         CALL iom_put( "bflx_iwm", z3d )
361         CALL iom_put( "pcmap_iwm", z2d )
362         DEALLOCATE( z2d , z3d )
363      ENDIF
364      CALL iom_put( "emix_iwm", zemx_iwm )
365     
366      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk)
367      !
368   END SUBROUTINE zdf_iwm
369
370
371   SUBROUTINE zdf_iwm_init
372      !!----------------------------------------------------------------------
373      !!                  ***  ROUTINE zdf_iwm_init  ***
374      !!                     
375      !! ** Purpose :   Initialization of the wave-driven vertical mixing, reading
376      !!              of input power maps and decay length scales in netcdf files.
377      !!
378      !! ** Method  : - Read the namzdf_iwm namelist and check the parameters
379      !!
380      !!              - Read the input data in NetCDF files :
381      !!              power available from high-mode wave breaking (mixing_power_bot.nc)
382      !!              power available from pycnocline-intensified wave-breaking (mixing_power_pyc.nc)
383      !!              power available from critical slope wave-breaking (mixing_power_cri.nc)
384      !!              WKB decay scale for high-mode wave-breaking (decay_scale_bot.nc)
385      !!              decay scale for critical slope wave-breaking (decay_scale_cri.nc)
386      !!
387      !! ** input   : - Namlist namzdf_iwm
388      !!              - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc,
389      !!              decay_scale_bot.nc decay_scale_cri.nc
390      !!
391      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter
392      !!              - Define ebot_iwm, epyc_iwm, ecri_iwm, hbot_iwm, hcri_iwm
393      !!
394      !! References : de Lavergne et al. JPO, 2015 ; de Lavergne PhD 2016
395      !!              de Lavergne et al. in prep., 2017
396      !!----------------------------------------------------------------------
397      INTEGER  ::   ji, jj, jk   ! dummy loop indices
398      INTEGER  ::   inum         ! local integer
399      INTEGER  ::   ios
400      REAL(wp) ::   zbot, zpyc, zcri   ! local scalars
401      !!
402      NAMELIST/namzdf_iwm/ nn_zpyc, ln_mevar, ln_tsdiff
403      !!----------------------------------------------------------------------
404      !
405      READ  ( numnam_ref, namzdf_iwm, IOSTAT = ios, ERR = 901)
406901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf_iwm in reference namelist' )
407      !
408      READ  ( numnam_cfg, namzdf_iwm, IOSTAT = ios, ERR = 902 )
409902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf_iwm in configuration namelist' )
410      IF(lwm) WRITE ( numond, namzdf_iwm )
411      !
412      IF(lwp) THEN                  ! Control print
413         WRITE(numout,*)
414         WRITE(numout,*) 'zdf_iwm_init : internal wave-driven mixing'
415         WRITE(numout,*) '~~~~~~~~~~~~'
416         WRITE(numout,*) '   Namelist namzdf_iwm : set wave-driven mixing parameters'
417         WRITE(numout,*) '      Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc
418         WRITE(numout,*) '      Variable (T) or constant (F) mixing efficiency            = ', ln_mevar
419         WRITE(numout,*) '      Differential internal wave-driven mixing (T) or not (F)   = ', ln_tsdiff
420      ENDIF
421     
422      ! The new wave-driven mixing parameterization elevates avt and avm in the interior, and
423      ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should
424      ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6).
425      avmb(:) = 1.4e-6_wp        ! viscous molecular value
426      avtb(:) = 1.e-10_wp        ! very small diffusive minimum (background avt is specified in zdf_iwm)   
427      avtb_2d(:,:) = 1.e0_wp     ! uniform
428      IF(lwp) THEN                  ! Control print
429         WRITE(numout,*)
430         WRITE(numout,*) '   Force the background value applied to avm & avt in TKE to be everywhere ',   &
431            &               'the viscous molecular value & a very small diffusive value, resp.'
432      ENDIF
433           
434      !                             ! allocate iwm arrays
435      IF( zdf_iwm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_iwm_init : unable to allocate iwm arrays' )
436      !
437      !                             ! read necessary fields
438!!$      CALL iom_open('mixing_power_bot',inum)       ! energy flux for high-mode wave breaking [W/m2]
439!!$      CALL iom_get  (inum, jpdom_data, 'field', ebot_iwm, 1 )
440!!$      CALL iom_close(inum)
441      ebot_iwm(:,:) = 1.e-6
442      !
443!!$      CALL iom_open('mixing_power_pyc',inum)       ! energy flux for pynocline-intensified wave breaking [W/m2]
444!!$      CALL iom_get  (inum, jpdom_data, 'field', epyc_iwm, 1 )
445!!$      CALL iom_close(inum)
446      epyc_iwm(:,:) = 1.e-6
447      !
448!!$      CALL iom_open('mixing_power_cri',inum)       ! energy flux for critical slope wave breaking [W/m2]
449!!$      CALL iom_get  (inum, jpdom_data, 'field', ecri_iwm, 1 )
450!!$      CALL iom_close(inum)
451      ecri_iwm(:,:) = 1.e-10
452      !
453!!$      CALL iom_open('decay_scale_bot',inum)        ! spatially variable decay scale for high-mode wave breaking [m]
454!!$      CALL iom_get  (inum, jpdom_data, 'field', hbot_iwm, 1 )
455!!$      CALL iom_close(inum)
456      hbot_iwm(:,:) = 100.
457      !
458!!$      CALL iom_open('decay_scale_cri',inum)        ! spatially variable decay scale for critical slope wave breaking [m]
459!!$      CALL iom_get  (inum, jpdom_data, 'field', hcri_iwm, 1 )
460!!$      CALL iom_close(inum)
461      hcri_iwm(:,:) = 100.
462
463      ebot_iwm(:,:) = ebot_iwm(:,:) * ssmask(:,:)
464      epyc_iwm(:,:) = epyc_iwm(:,:) * ssmask(:,:)
465      ecri_iwm(:,:) = ecri_iwm(:,:) * ssmask(:,:)
466
467      zbot = glob_sum( 'zdfiwm', e1e2t(:,:) * ebot_iwm(:,:) )
468      zpyc = glob_sum( 'zdfiwm', e1e2t(:,:) * epyc_iwm(:,:) )
469      zcri = glob_sum( 'zdfiwm', e1e2t(:,:) * ecri_iwm(:,:) )
470      IF(lwp) THEN
471         WRITE(numout,*) '      High-mode wave-breaking energy:             ', zbot * 1.e-12_wp, 'TW'
472         WRITE(numout,*) '      Pycnocline-intensifed wave-breaking energy: ', zpyc * 1.e-12_wp, 'TW'
473         WRITE(numout,*) '      Critical slope wave-breaking energy:        ', zcri * 1.e-12_wp, 'TW'
474      ENDIF
475      !
476   END SUBROUTINE zdf_iwm_init
477
478   !!======================================================================
479END MODULE zdfiwm
Note: See TracBrowser for help on using the repository browser.