source: NEMO/trunk/src/OCE/ZDF/zdfphy.F90 @ 12702

Last change on this file since 12702 was 12377, checked in by acc, 13 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: 19.3 KB
Line 
1MODULE zdfphy
2   !!======================================================================
3   !!                      ***  MODULE  zdfphy  ***
4   !! Vertical ocean physics :   manager of all vertical physics packages
5   !!======================================================================
6   !! History :  4.0  !  2017-04  (G. Madec)  original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   zdf_phy_init  : initialization of all vertical physics packages
11   !!   zdf_phy       : upadate at each time-step the vertical mixing coeff.
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and tracers variables
14   USE zdf_oce        ! vertical physics: shared variables         
15   USE zdfdrg         ! vertical physics: top/bottom drag coef.
16   USE zdfsh2         ! vertical physics: shear production term of TKE
17   USE zdfric         ! vertical physics: RIChardson dependent vertical mixing   
18   USE zdftke         ! vertical physics: TKE vertical mixing
19   USE zdfgls         ! vertical physics: GLS vertical mixing
20   USE zdfosm         ! vertical physics: OSMOSIS vertical mixing
21   USE zdfddm         ! vertical physics: double diffusion mixing     
22   USE zdfevd         ! vertical physics: convection via enhanced vertical diffusion 
23   USE zdfiwm         ! vertical physics: internal wave-induced mixing 
24   USE zdfswm         ! vertical physics: surface  wave-induced mixing
25   USE zdfmxl         ! vertical physics: mixed layer
26   USE tranpc         ! convection: non penetrative adjustment
27   USE trc_oce        ! variables shared between passive tracer & ocean           
28   USE sbc_oce        ! surface module (only for nn_isf in the option compatibility test)
29   USE sbcrnf         ! surface boundary condition: runoff variables
30#if defined key_agrif
31   USE agrif_oce_interp   ! interpavm
32#endif
33   !
34   USE in_out_manager ! I/O manager
35   USE iom            ! IOM library
36   USE lbclnk         ! lateral boundary conditions
37   USE lib_mpp        ! distribued memory computing
38   USE timing         ! Timing
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   zdf_phy_init  ! called by nemogcm.F90
44   PUBLIC   zdf_phy       ! called by step.F90
45
46   INTEGER ::   nzdf_phy   ! type of vertical closure used
47   !                       ! associated indicators
48   INTEGER, PARAMETER ::   np_CST = 1   ! Constant Kz
49   INTEGER, PARAMETER ::   np_RIC = 2   ! Richardson number dependent Kz
50   INTEGER, PARAMETER ::   np_TKE = 3   ! Turbulente Kinetic Eenergy closure scheme for Kz
51   INTEGER, PARAMETER ::   np_GLS = 4   ! Generic Length Scale closure scheme for Kz
52   INTEGER, PARAMETER ::   np_OSM = 5   ! OSMOSIS-OBL closure scheme for Kz
53
54   LOGICAL ::   l_zdfsh2   ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC))
55
56   !!----------------------------------------------------------------------
57   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
58   !! $Id$
59   !! Software governed by the CeCILL license (see ./LICENSE)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE zdf_phy_init( Kmm )
64      !!----------------------------------------------------------------------
65      !!                  ***  ROUTINE zdf_phy_init  ***
66      !!
67      !! ** Purpose :   initializations of the vertical ocean physics
68      !!
69      !! ** Method  :   Read namelist namzdf, control logicals
70      !!                set horizontal shape and vertical profile of background mixing coef.
71      !!----------------------------------------------------------------------
72      INTEGER, INTENT(in)    :: Kmm ! time level index (middle)
73      !
74      INTEGER ::   jk            ! dummy loop indices
75      INTEGER ::   ioptio, ios   ! local integers
76      !!
77      NAMELIST/namzdf/ ln_zdfcst, ln_zdfric, ln_zdftke, ln_zdfgls,   &     ! type of closure scheme
78         &             ln_zdfosm,                                    &     ! type of closure scheme
79         &             ln_zdfevd, nn_evdm, rn_evd ,                  &     ! convection : evd
80         &             ln_zdfnpc, nn_npc , nn_npcp,                  &     ! convection : npc
81         &             ln_zdfddm, rn_avts, rn_hsbfr,                 &     ! double diffusion
82         &             ln_zdfswm,                                    &     ! surface  wave-induced mixing
83         &             ln_zdfiwm,                                    &     ! internal  -      -      -
84         &             ln_zad_Aimp,                                  &     ! apdative-implicit vertical advection
85         &             rn_avm0, rn_avt0, nn_avb, nn_havtb                  ! coefficients
86      !!----------------------------------------------------------------------
87      !
88      IF(lwp) THEN
89         WRITE(numout,*)
90         WRITE(numout,*) 'zdf_phy_init: ocean vertical physics'
91         WRITE(numout,*) '~~~~~~~~~~~~'
92      ENDIF
93      !
94      !                           !==  Namelist  ==!
95      READ  ( numnam_ref, namzdf, IOSTAT = ios, ERR = 901)
96901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf in reference namelist' )
97      !
98      READ  ( numnam_cfg, namzdf, IOSTAT = ios, ERR = 902 )
99902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf in configuration namelist' )
100      IF(lwm)   WRITE ( numond, namzdf )
101      !
102      IF(lwp) THEN                      ! Parameter print
103         WRITE(numout,*) '   Namelist namzdf : set vertical mixing mixing parameters'
104         WRITE(numout,*) '      adaptive-implicit vertical advection'
105         WRITE(numout,*) '         Courant number targeted application   ln_zad_Aimp = ', ln_zad_Aimp
106         WRITE(numout,*) '      vertical closure scheme'
107         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfcst = ', ln_zdfcst
108         WRITE(numout,*) '         Richardson number dependent closure     ln_zdfric = ', ln_zdfric
109         WRITE(numout,*) '         Turbulent Kinetic Energy closure (TKE)  ln_zdftke = ', ln_zdftke
110         WRITE(numout,*) '         Generic Length Scale closure (GLS)      ln_zdfgls = ', ln_zdfgls
111         WRITE(numout,*) '         OSMOSIS-OBL closure (OSM)               ln_zdfosm = ', ln_zdfosm
112         WRITE(numout,*) '      convection: '
113         WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd
114         WRITE(numout,*) '            applied on momentum (=1/0)             nn_evdm = ', nn_evdm
115         WRITE(numout,*) '            vertical coefficient for evd           rn_evd  = ', rn_evd
116         WRITE(numout,*) '         non-penetrative convection (npc)        ln_zdfnpc = ', ln_zdfnpc
117         WRITE(numout,*) '            npc call  frequency                    nn_npc  = ', nn_npc
118         WRITE(numout,*) '            npc print frequency                    nn_npcp = ', nn_npcp
119         WRITE(numout,*) '      double diffusive mixing                    ln_zdfddm = ', ln_zdfddm
120         WRITE(numout,*) '         maximum avs for dd mixing                 rn_avts = ', rn_avts
121         WRITE(numout,*) '         heat/salt buoyancy flux ratio             rn_hsbfr= ', rn_hsbfr
122         WRITE(numout,*) '      gravity wave-induced mixing'
123         WRITE(numout,*) '         surface  wave (Qiao et al 2010)         ln_zdfswm = ', ln_zdfswm                                          ! surface wave induced mixing
124         WRITE(numout,*) '         internal wave (de Lavergne et al 2017)  ln_zdfiwm = ', ln_zdfiwm
125         WRITE(numout,*) '      coefficients : '
126         WRITE(numout,*) '         vertical eddy viscosity                 rn_avm0   = ', rn_avm0
127         WRITE(numout,*) '         vertical eddy diffusivity               rn_avt0   = ', rn_avt0
128         WRITE(numout,*) '         constant background or profile          nn_avb    = ', nn_avb
129         WRITE(numout,*) '         horizontal variation for avtb           nn_havtb  = ', nn_havtb
130      ENDIF
131
132      IF( ln_zad_Aimp ) THEN
133         IF( zdf_phy_alloc() /= 0 )   &
134            &       CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' )
135         Cu_adv(:,:,:) = 0._wp
136         wi    (:,:,:) = 0._wp
137      ENDIF
138      !                          !==  Background eddy viscosity and diffusivity  ==!
139      IF( nn_avb == 0 ) THEN             ! Define avmb, avtb from namelist parameter
140         avmb(:) = rn_avm0
141         avtb(:) = rn_avt0                     
142      ELSE                               ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
143         avmb(:) = rn_avm0
144         avtb(:) = rn_avt0 + ( 3.e-4_wp - 2._wp * rn_avt0 ) * 1.e-4_wp * gdepw_1d(:)   ! m2/s
145         IF(ln_sco .AND. lwp)   CALL ctl_warn( 'avtb profile not valid in sco' )
146      ENDIF
147      !                                  ! 2D shape of the avtb
148      avtb_2d(:,:) = 1._wp                   ! uniform
149      !
150      IF( nn_havtb == 1 ) THEN               ! decrease avtb by a factor of ten in the equatorial band
151           !                                 !   -15S -5S : linear decrease from avt0 to avt0/10.
152           !                                 !   -5S  +5N : cst value avt0/10.
153           !                                 !    5N  15N : linear increase from avt0/10, to avt0
154           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
155           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
156           WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
157      ENDIF
158      !
159      DO jk = 1, jpk                      ! set turbulent closure Kz to the background value (avt_k, avm_k)
160         avt_k(:,:,jk) = avtb_2d(:,:) * avtb(jk) * wmask (:,:,jk)
161         avm_k(:,:,jk) =                avmb(jk) * wmask (:,:,jk)
162      END DO
163!!gm  to be tested only the 1st & last levels
164!      avt  (:,:, 1 ) = 0._wp   ;   avs(:,:, 1 ) = 0._wp   ;   avm  (:,:, 1 ) = 0._wp
165!      avt  (:,:,jpk) = 0._wp   ;   avs(:,:,jpk) = 0._wp   ;   avm  (:,:,jpk) = 0._wp
166!!gm
167      avt  (:,:,:) = 0._wp   ;   avs(:,:,:) = 0._wp   ;   avm  (:,:,:) = 0._wp
168
169      !                          !==  Convection  ==!
170      !
171      IF( ln_zdfnpc .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfnpc and ln_zdfevd' )
172      IF( ln_zdfosm .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfosm and ln_zdfevd' )
173      IF( lk_top    .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' )
174      IF( lk_top    .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: osmosis scheme is not working with key_top' )
175      IF(lwp) THEN
176         WRITE(numout,*)
177         IF    ( ln_zdfnpc ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use non penetrative convective scheme'
178         ELSEIF( ln_zdfevd ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use enhanced vertical diffusion scheme'
179         ELSE                      ;   WRITE(numout,*) '   ==>>>   convection: no specific scheme used'
180         ENDIF
181      ENDIF
182
183      IF(lwp) THEN               !==  Double Diffusion Mixing parameterization  ==!   (ddm)
184         WRITE(numout,*)
185         IF( ln_zdfddm ) THEN   ;   WRITE(numout,*) '   ==>>>   use double diffusive mixing: avs /= avt'
186         ELSE                   ;   WRITE(numout,*) '   ==>>>   No  double diffusive mixing: avs = avt'
187         ENDIF
188      ENDIF
189
190      !                          !==  type of vertical turbulent closure  ==!    (set nzdf_phy)
191      ioptio = 0 
192      IF( ln_zdfcst ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_CST   ;   ENDIF
193      IF( ln_zdfric ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_RIC   ;   CALL zdf_ric_init          ;   ENDIF
194      IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init( Kmm )   ;   ENDIF
195      IF( ln_zdfgls ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_GLS   ;   CALL zdf_gls_init          ;   ENDIF
196      IF( ln_zdfosm ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_OSM   ;   CALL zdf_osm_init( Kmm )   ;   ENDIF
197      !
198      IF( ioptio /= 1 )    CALL ctl_stop( 'zdf_phy_init: one and only one vertical diffusion option has to be defined ' )
199      IF( ln_isfcav ) THEN
200      IF( ln_zdfric .OR. ln_zdfgls )    CALL ctl_stop( 'zdf_phy_init: zdfric and zdfgls never tested with ice shelves cavities ' )
201      ENDIF
202      !                                ! shear production term flag
203      IF( ln_zdfcst ) THEN   ;   l_zdfsh2 = .FALSE.
204      ELSE                   ;   l_zdfsh2 = .TRUE.
205      ENDIF
206
207      !                          !== gravity wave-driven mixing  ==!
208      IF( ln_zdfiwm )   CALL zdf_iwm_init       ! internal wave-driven mixing
209      IF( ln_zdfswm )   CALL zdf_swm_init       ! surface  wave-driven mixing
210
211      !                          !== top/bottom friction  ==!
212      CALL zdf_drg_init
213      !
214      !                          !== time-stepping  ==!
215      ! Check/update of time stepping done in dynzdf_init/trazdf_init
216      !!gm move it here ?
217      !
218   END SUBROUTINE zdf_phy_init
219
220
221   SUBROUTINE zdf_phy( kt, Kbb, Kmm, Krhs )
222      !!----------------------------------------------------------------------
223      !!                     ***  ROUTINE zdf_phy  ***
224      !!
225      !! ** Purpose :  Update ocean physics at each time-step
226      !!
227      !! ** Method  :
228      !!
229      !! ** Action  :   avm, avt vertical eddy viscosity and diffusivity at w-points
230      !!                nmld ??? mixed layer depth in level and meters   <<<<====verifier !
231      !!                bottom stress.....                               <<<<====verifier !
232      !!----------------------------------------------------------------------
233      INTEGER, INTENT(in) ::   kt         ! ocean time-step index
234      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs   ! ocean time level indices
235      !
236      INTEGER ::   ji, jj, jk   ! dummy loop indice
237      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2   ! shear production
238      !! ---------------------------------------------------------------------
239      !
240      IF( ln_timing )   CALL timing_start('zdf_phy')
241      !
242      IF( l_zdfdrg ) THEN     !==  update top/bottom drag  ==!   (non-linear cases)
243         !
244         !                       !* bottom drag
245         CALL zdf_drg( kt, Kmm, mbkt , r_Cdmin_bot, r_Cdmax_bot,   &   ! <<== in
246            &              r_z0_bot,   r_ke0_bot,    rCd0_bot,   &
247            &                                        rCdU_bot  )     ! ==>> out : bottom drag [m/s]
248         IF( ln_isfcav ) THEN    !* top drag   (ocean cavities)
249            CALL zdf_drg( kt, Kmm, mikt , r_Cdmin_top, r_Cdmax_top,   &   ! <<== in
250               &              r_z0_top,   r_ke0_top,    rCd0_top,   &
251               &                                        rCdU_top  )     ! ==>> out : bottom drag [m/s]
252         ENDIF
253      ENDIF
254      !
255      !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k)
256      !
257      IF( l_zdfsh2 )   &         !* shear production at w-points (energy conserving form)
258         CALL zdf_sh2( Kbb, Kmm, avm_k,   &     ! <<== in
259            &                     zsh2    )     ! ==>> out : shear production
260      !
261      SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points
262      CASE( np_RIC )   ;   CALL zdf_ric( kt,      Kmm, zsh2, avm_k, avt_k )    ! Richardson number dependent Kz
263      CASE( np_TKE )   ;   CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k )    ! TKE closure scheme for Kz
264      CASE( np_GLS )   ;   CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k )    ! GLS closure scheme for Kz
265      CASE( np_OSM )   ;   CALL zdf_osm( kt, Kbb, Kmm, Krhs, avm_k, avt_k )    ! OSMOSIS closure scheme for Kz
266!     CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value)
267!         ! avt_k and avm_k set one for all at initialisation phase
268!!gm         avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
269!!gm         avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
270      END SELECT
271     
272      !                          !==  ocean Kz  ==!   (avt, avs, avm)
273      !
274      !                                         !* start from turbulent closure values
275      avt(:,:,2:jpkm1) = avt_k(:,:,2:jpkm1)
276      avm(:,:,2:jpkm1) = avm_k(:,:,2:jpkm1)
277      !
278      IF( ln_rnf_mouth ) THEN                   !* increase diffusivity at rivers mouths
279         DO jk = 2, nkrnf
280            avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * wmask(:,:,jk)
281         END DO
282      ENDIF
283      !
284      IF( ln_zdfevd )   CALL zdf_evd( kt, Kmm, Krhs, avm, avt )  !* convection: enhanced vertical eddy diffusivity
285      !
286      !                                         !* double diffusive mixing
287      IF( ln_zdfddm ) THEN                            ! update avt and compute avs
288                        CALL zdf_ddm( kt, Kmm,  avm, avt, avs )
289      ELSE                                            ! same mixing on all tracers
290         avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1)
291      ENDIF
292      !
293      !                                         !* wave-induced mixing
294      IF( ln_zdfswm )   CALL zdf_swm( kt, Kmm, avm, avt, avs )   ! surface  wave (Qiao et al. 2004)
295      IF( ln_zdfiwm )   CALL zdf_iwm( kt, Kmm, avm, avt, avs )   ! internal wave (de Lavergne et al 2017)
296
297#if defined key_agrif 
298      ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2)
299      IF( l_zdfsh2 )   CALL Agrif_avm
300#endif
301
302      !                                         !* Lateral boundary conditions (sign unchanged)
303      IF( l_zdfsh2 ) THEN
304         CALL lbc_lnk_multi( 'zdfphy', avm_k, 'W', 1. , avt_k, 'W', 1.,   &
305            &                avm  , 'W', 1. , avt  , 'W', 1. , avs , 'W', 1. )
306      ELSE
307         CALL lbc_lnk_multi( 'zdfphy', avm  , 'W', 1. , avt  , 'W', 1. , avs , 'W', 1. )
308      ENDIF
309      !
310      IF( l_zdfdrg ) THEN     ! drag  have been updated (non-linear cases)
311         IF( ln_isfcav ) THEN   ;  CALL lbc_lnk_multi( 'zdfphy', rCdU_top, 'T', 1. , rCdU_bot, 'T', 1. )   ! top & bot drag
312         ELSE                   ;  CALL lbc_lnk      ( 'zdfphy', rCdU_bot, 'T', 1. )                       ! bottom drag only
313         ENDIF
314      ENDIF
315      !
316      CALL zdf_mxl( kt, Kmm )                        !* mixed layer depth, and level
317      !
318      IF( lrst_oce ) THEN                       !* write TKE, GLS or RIC fields in the restart file
319         IF( ln_zdftke )   CALL tke_rst( kt, 'WRITE' )
320         IF( ln_zdfgls )   CALL gls_rst( kt, 'WRITE' )
321         IF( ln_zdfric )   CALL ric_rst( kt, 'WRITE' ) 
322         ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after ww has been updated
323      ENDIF
324      !
325      IF( ln_timing )   CALL timing_stop('zdf_phy')
326      !
327   END SUBROUTINE zdf_phy
328   INTEGER FUNCTION zdf_phy_alloc()
329      !!----------------------------------------------------------------------
330      !!                 ***  FUNCTION zdf_phy_alloc  ***
331      !!----------------------------------------------------------------------
332     ! Allocate wi array (declared in oce.F90) for use with the adaptive-implicit vertical velocity option
333     ALLOCATE(     wi(jpi,jpj,jpk), Cu_adv(jpi,jpj,jpk),  STAT= zdf_phy_alloc )
334     IF( zdf_phy_alloc /= 0 )   CALL ctl_warn('zdf_phy_alloc: failed to allocate ln_zad_Aimp=T required arrays')
335     CALL mpp_sum ( 'zdfphy', zdf_phy_alloc )
336   END FUNCTION zdf_phy_alloc
337
338   !!======================================================================
339END MODULE zdfphy
Note: See TracBrowser for help on using the repository browser.