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.
zdfphy.F90 in NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/ZDF – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/ZDF/zdfphy.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.

  • Property svn:keywords set to Id
File size: 19.0 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
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 ::   jk            ! dummy loop indices
73      INTEGER ::   ioptio, ios   ! local integers
74      !!
75      NAMELIST/namzdf/ ln_zdfcst, ln_zdfric, ln_zdftke, ln_zdfgls,   &     ! type of closure scheme
76         &             ln_zdfosm,                                    &     ! type of closure scheme
77         &             ln_zdfevd, nn_evdm, rn_evd ,                  &     ! convection : evd
78         &             ln_zdfnpc, nn_npc , nn_npcp,                  &     ! convection : npc
79         &             ln_zdfddm, rn_avts, rn_hsbfr,                 &     ! double diffusion
80         &             ln_zdfswm,                                    &     ! surface  wave-induced mixing
81         &             ln_zdfiwm,                                    &     ! internal  -      -      -
82         &             ln_zad_Aimp,                                  &     ! apdative-implicit vertical advection
83         &             rn_avm0, rn_avt0, nn_avb, nn_havtb                  ! coefficients
84      !!----------------------------------------------------------------------
85      !
86      IF(lwp) THEN
87         WRITE(numout,*)
88         WRITE(numout,*) 'zdf_phy_init: ocean vertical physics'
89         WRITE(numout,*) '~~~~~~~~~~~~'
90      ENDIF
91      !
92      !                           !==  Namelist  ==!
93      READ  ( numnam_ref, namzdf, IOSTAT = ios, ERR = 901)
94901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf in reference namelist' )
95      !
96      READ  ( numnam_cfg, namzdf, IOSTAT = ios, ERR = 902 )
97902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf in configuration namelist' )
98      IF(lwm)   WRITE ( numond, namzdf )
99      !
100      IF(lwp) THEN                      ! Parameter print
101         WRITE(numout,*) '   Namelist namzdf : set vertical mixing mixing parameters'
102         WRITE(numout,*) '      adaptive-implicit vertical advection'
103         WRITE(numout,*) '         Courant number targeted application   ln_zad_Aimp = ', ln_zad_Aimp
104         WRITE(numout,*) '      vertical closure scheme'
105         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfcst = ', ln_zdfcst
106         WRITE(numout,*) '         Richardson number dependent closure     ln_zdfric = ', ln_zdfric
107         WRITE(numout,*) '         Turbulent Kinetic Energy closure (TKE)  ln_zdftke = ', ln_zdftke
108         WRITE(numout,*) '         Generic Length Scale closure (GLS)      ln_zdfgls = ', ln_zdfgls
109         WRITE(numout,*) '         OSMOSIS-OBL closure (OSM)               ln_zdfosm = ', ln_zdfosm
110         WRITE(numout,*) '      convection: '
111         WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd
112         WRITE(numout,*) '            applied on momentum (=1/0)             nn_evdm = ', nn_evdm
113         WRITE(numout,*) '            vertical coefficient for evd           rn_evd  = ', rn_evd
114         WRITE(numout,*) '         non-penetrative convection (npc)        ln_zdfnpc = ', ln_zdfnpc
115         WRITE(numout,*) '            npc call  frequency                    nn_npc  = ', nn_npc
116         WRITE(numout,*) '            npc print frequency                    nn_npcp = ', nn_npcp
117         WRITE(numout,*) '      double diffusive mixing                    ln_zdfddm = ', ln_zdfddm
118         WRITE(numout,*) '         maximum avs for dd mixing                 rn_avts = ', rn_avts
119         WRITE(numout,*) '         heat/salt buoyancy flux ratio             rn_hsbfr= ', rn_hsbfr
120         WRITE(numout,*) '      gravity wave-induced mixing'
121         WRITE(numout,*) '         surface  wave (Qiao et al 2010)         ln_zdfswm = ', ln_zdfswm                                          ! surface wave induced mixing
122         WRITE(numout,*) '         internal wave (de Lavergne et al 2017)  ln_zdfiwm = ', ln_zdfiwm
123         WRITE(numout,*) '      coefficients : '
124         WRITE(numout,*) '         vertical eddy viscosity                 rn_avm0   = ', rn_avm0
125         WRITE(numout,*) '         vertical eddy diffusivity               rn_avt0   = ', rn_avt0
126         WRITE(numout,*) '         constant background or profile          nn_avb    = ', nn_avb
127         WRITE(numout,*) '         horizontal variation for avtb           nn_havtb  = ', nn_havtb
128      ENDIF
129
130      IF( ln_zad_Aimp ) THEN
131         IF( zdf_phy_alloc() /= 0 )   &
132            &       CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' )
133         Cu_adv(:,:,:) = 0._wp
134         wi    (:,:,:) = 0._wp
135      ENDIF
136      !                          !==  Background eddy viscosity and diffusivity  ==!
137      IF( nn_avb == 0 ) THEN             ! Define avmb, avtb from namelist parameter
138         avmb(:) = rn_avm0
139         avtb(:) = rn_avt0                     
140      ELSE                               ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
141         avmb(:) = rn_avm0
142         avtb(:) = rn_avt0 + ( 3.e-4_wp - 2._wp * rn_avt0 ) * 1.e-4_wp * gdepw_1d(:)   ! m2/s
143         IF(ln_sco .AND. lwp)   CALL ctl_warn( 'avtb profile not valid in sco' )
144      ENDIF
145      !                                  ! 2D shape of the avtb
146      avtb_2d(:,:) = 1._wp                   ! uniform
147      !
148      IF( nn_havtb == 1 ) THEN               ! decrease avtb by a factor of ten in the equatorial band
149           !                                 !   -15S -5S : linear decrease from avt0 to avt0/10.
150           !                                 !   -5S  +5N : cst value avt0/10.
151           !                                 !    5N  15N : linear increase from avt0/10, to avt0
152           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
153           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
154           WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
155      ENDIF
156      !
157      DO jk = 1, jpk                      ! set turbulent closure Kz to the background value (avt_k, avm_k)
158         avt_k(:,:,jk) = avtb_2d(:,:) * avtb(jk) * wmask (:,:,jk)
159         avm_k(:,:,jk) =                avmb(jk) * wmask (:,:,jk)
160      END DO
161!!gm  to be tested only the 1st & last levels
162!      avt  (:,:, 1 ) = 0._wp   ;   avs(:,:, 1 ) = 0._wp   ;   avm  (:,:, 1 ) = 0._wp
163!      avt  (:,:,jpk) = 0._wp   ;   avs(:,:,jpk) = 0._wp   ;   avm  (:,:,jpk) = 0._wp
164!!gm
165      avt  (:,:,:) = 0._wp   ;   avs(:,:,:) = 0._wp   ;   avm  (:,:,:) = 0._wp
166
167      !                          !==  Convection  ==!
168      !
169      IF( ln_zdfnpc .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfnpc and ln_zdfevd' )
170      IF( ln_zdfosm .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfosm and ln_zdfevd' )
171      IF( lk_top    .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' )
172      IF( lk_top    .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: osmosis scheme is not working with key_top' )
173      IF(lwp) THEN
174         WRITE(numout,*)
175         IF    ( ln_zdfnpc ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use non penetrative convective scheme'
176         ELSEIF( ln_zdfevd ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use enhanced vertical diffusion scheme'
177         ELSE                      ;   WRITE(numout,*) '   ==>>>   convection: no specific scheme used'
178         ENDIF
179      ENDIF
180
181      IF(lwp) THEN               !==  Double Diffusion Mixing parameterization  ==!   (ddm)
182         WRITE(numout,*)
183         IF( ln_zdfddm ) THEN   ;   WRITE(numout,*) '   ==>>>   use double diffusive mixing: avs /= avt'
184         ELSE                   ;   WRITE(numout,*) '   ==>>>   No  double diffusive mixing: avs = avt'
185         ENDIF
186      ENDIF
187
188      !                          !==  type of vertical turbulent closure  ==!    (set nzdf_phy)
189      ioptio = 0 
190      IF( ln_zdfcst ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_CST   ;   ENDIF
191      IF( ln_zdfric ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_RIC   ;   CALL zdf_ric_init   ;   ENDIF
192      IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init   ;   ENDIF
193      IF( ln_zdfgls ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_GLS   ;   CALL zdf_gls_init   ;   ENDIF
194      IF( ln_zdfosm ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_OSM   ;   CALL zdf_osm_init   ;   ENDIF
195      !
196      IF( ioptio /= 1 )    CALL ctl_stop( 'zdf_phy_init: one and only one vertical diffusion option has to be defined ' )
197      IF( ln_isfcav ) THEN
198      IF( ln_zdfric .OR. ln_zdfgls )    CALL ctl_stop( 'zdf_phy_init: zdfric and zdfgls never tested with ice shelves cavities ' )
199      ENDIF
200      !                                ! shear production term flag
201      IF( ln_zdfcst ) THEN   ;   l_zdfsh2 = .FALSE.
202      ELSE                   ;   l_zdfsh2 = .TRUE.
203      ENDIF
204
205      !                          !== gravity wave-driven mixing  ==!
206      IF( ln_zdfiwm )   CALL zdf_iwm_init       ! internal wave-driven mixing
207      IF( ln_zdfswm )   CALL zdf_swm_init       ! surface  wave-driven mixing
208
209      !                          !== top/bottom friction  ==!
210      CALL zdf_drg_init
211      !
212      !                          !== time-stepping  ==!
213      ! Check/update of time stepping done in dynzdf_init/trazdf_init
214      !!gm move it here ?
215      !
216   END SUBROUTINE zdf_phy_init
217
218
219   SUBROUTINE zdf_phy( kt )
220      !!----------------------------------------------------------------------
221      !!                     ***  ROUTINE zdf_phy  ***
222      !!
223      !! ** Purpose :  Update ocean physics at each time-step
224      !!
225      !! ** Method  :
226      !!
227      !! ** Action  :   avm, avt vertical eddy viscosity and diffusivity at w-points
228      !!                nmld ??? mixed layer depth in level and meters   <<<<====verifier !
229      !!                bottom stress.....                               <<<<====verifier !
230      !!----------------------------------------------------------------------
231      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
232      !
233      INTEGER ::   ji, jj, jk   ! dummy loop indice
234      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2   ! shear production
235      !! ---------------------------------------------------------------------
236      !
237      IF( ln_timing )   CALL timing_start('zdf_phy')
238      !
239      IF( l_zdfdrg ) THEN     !==  update top/bottom drag  ==!   (non-linear cases)
240         !
241         !                       !* bottom drag
242         CALL zdf_drg( kt, mbkt    , r_Cdmin_bot, r_Cdmax_bot,   &   ! <<== in
243            &              r_z0_bot,   r_ke0_bot,    rCd0_bot,   &
244            &                                        rCdU_bot  )     ! ==>> out : bottom drag [m/s]
245         IF( ln_isfcav ) THEN    !* top drag   (ocean cavities)
246            CALL zdf_drg( kt, mikt    , r_Cdmin_top, r_Cdmax_top,   &   ! <<== in
247               &              r_z0_top,   r_ke0_top,    rCd0_top,   &
248               &                                        rCdU_top  )     ! ==>> out : bottom drag [m/s]
249         ENDIF
250      ENDIF
251      !
252      !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k)
253      !
254      IF( l_zdfsh2 )   &         !* shear production at w-points (energy conserving form)
255         CALL zdf_sh2( ub, vb, un, vn, avm_k,   &     ! <<== in
256            &                           zsh2    )     ! ==>> out : shear production
257      !
258      SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points
259      CASE( np_RIC )   ;   CALL zdf_ric( kt, gdept_n, zsh2, avm_k, avt_k )    ! Richardson number dependent Kz
260      CASE( np_TKE )   ;   CALL zdf_tke( kt         , zsh2, avm_k, avt_k )    ! TKE closure scheme for Kz
261      CASE( np_GLS )   ;   CALL zdf_gls( kt         , zsh2, avm_k, avt_k )    ! GLS closure scheme for Kz
262      CASE( np_OSM )   ;   CALL zdf_osm( kt               , avm_k, avt_k )    ! OSMOSIS closure scheme for Kz
263!     CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value)
264!         ! avt_k and avm_k set one for all at initialisation phase
265!!gm         avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
266!!gm         avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
267      END SELECT
268     
269      !                          !==  ocean Kz  ==!   (avt, avs, avm)
270      !
271      !                                         !* start from turbulent closure values
272      avt(:,:,2:jpkm1) = avt_k(:,:,2:jpkm1)
273      avm(:,:,2:jpkm1) = avm_k(:,:,2:jpkm1)
274      !
275      IF( ln_rnf_mouth ) THEN                   !* increase diffusivity at rivers mouths
276         DO jk = 2, nkrnf
277            avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * wmask(:,:,jk)
278         END DO
279      ENDIF
280      !
281      IF( ln_zdfevd )   CALL zdf_evd( kt, avm, avt )  !* convection: enhanced vertical eddy diffusivity
282      !
283      !                                         !* double diffusive mixing
284      IF( ln_zdfddm ) THEN                            ! update avt and compute avs
285                        CALL zdf_ddm( kt, avm, avt, avs )
286      ELSE                                            ! same mixing on all tracers
287         avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1)
288      ENDIF
289      !
290      !                                         !* wave-induced mixing
291      IF( ln_zdfswm )   CALL zdf_swm( kt, avm, avt, avs )   ! surface  wave (Qiao et al. 2004)
292      IF( ln_zdfiwm )   CALL zdf_iwm( kt, avm, avt, avs )   ! internal wave (de Lavergne et al 2017)
293
294#if defined key_agrif 
295      ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2)
296      IF( l_zdfsh2 )   CALL Agrif_avm
297#endif
298
299      !                                         !* Lateral boundary conditions (sign unchanged)
300      IF( l_zdfsh2 ) THEN
301         CALL lbc_lnk_multi( 'zdfphy', avm_k, 'W', 1. , avt_k, 'W', 1.,   &
302            &                avm  , 'W', 1. , avt  , 'W', 1. , avs , 'W', 1. )
303      ELSE
304         CALL lbc_lnk_multi( 'zdfphy', avm  , 'W', 1. , avt  , 'W', 1. , avs , 'W', 1. )
305      ENDIF
306      !
307      IF( l_zdfdrg ) THEN     ! drag  have been updated (non-linear cases)
308         IF( ln_isfcav ) THEN   ;  CALL lbc_lnk_multi( 'zdfphy', rCdU_top, 'T', 1. , rCdU_bot, 'T', 1. )   ! top & bot drag
309         ELSE                   ;  CALL lbc_lnk      ( 'zdfphy', rCdU_bot, 'T', 1. )                       ! bottom drag only
310         ENDIF
311      ENDIF
312      !
313      CALL zdf_mxl( kt )                        !* mixed layer depth, and level
314      !
315      IF( lrst_oce ) THEN                       !* write TKE, GLS or RIC fields in the restart file
316         IF( ln_zdftke )   CALL tke_rst( kt, 'WRITE' )
317         IF( ln_zdfgls )   CALL gls_rst( kt, 'WRITE' )
318         IF( ln_zdfric )   CALL ric_rst( kt, 'WRITE' ) 
319         ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after wn has been updated
320      ENDIF
321      !
322      IF( ln_timing )   CALL timing_stop('zdf_phy')
323      !
324   END SUBROUTINE zdf_phy
325   INTEGER FUNCTION zdf_phy_alloc()
326      !!----------------------------------------------------------------------
327      !!                 ***  FUNCTION zdf_phy_alloc  ***
328      !!----------------------------------------------------------------------
329     ! Allocate wi array (declared in oce.F90) for use with the adaptive-implicit vertical velocity option
330     ALLOCATE(     wi(jpi,jpj,jpk), Cu_adv(jpi,jpj,jpk),  STAT= zdf_phy_alloc )
331     IF( zdf_phy_alloc /= 0 )   CALL ctl_warn('zdf_phy_alloc: failed to allocate ln_zad_Aimp=T required arrays')
332     CALL mpp_sum ( 'zdfphy', zdf_phy_alloc )
333   END FUNCTION zdf_phy_alloc
334
335   !!======================================================================
336END MODULE zdfphy
Note: See TracBrowser for help on using the repository browser.