source: NEMO/trunk/src/OCE/ZDF/zdfiwm.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 10 months 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.

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