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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfphy.F90 @ 10954

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

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only

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