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 branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90 @ 8093

Last change on this file since 8093 was 8093, checked in by gm, 7 years ago

#1880 (HPC-09) - step-6: prepare some forthcoming evolutions (ZDF modules mainly)

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