source: NEMO/trunk/src/OCE/ZDF/zdfric.F90 @ 12489

Last change on this file since 12489 was 12489, checked in by davestorkey, 8 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

  • Property svn:keywords set to Id
File size: 12.1 KB
Line 
1MODULE zdfric
2   !!======================================================================
3   !!                       ***  MODULE  zdfric  ***
4   !! Ocean physics:  vertical mixing coefficient compute from the local
5   !!                 Richardson number dependent formulation
6   !!======================================================================
7   !! History :  OPA  !  1987-09  (P. Andrich)  Original code
8   !!            4.0  !  1991-11  (G. Madec)
9   !!            7.0  !  1996-01  (G. Madec)  complete rewriting of multitasking suppression of common work arrays
10   !!            8.0  !  1997-06  (G. Madec)  complete rewriting of zdfmix
11   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
12   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
13   !!            3.3.1!  2011-09  (P. Oddo) Mixed layer depth parameterization
14   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   zdf_ric_init  : initialization, namelist read, & parameters control
19   !!   zdf_ric       : update momentum and tracer Kz from the Richardson number
20   !!   ric_rst       : read/write RIC restart in ocean restart file
21   !!----------------------------------------------------------------------
22   USE oce            ! ocean dynamics and tracers variables
23   USE dom_oce        ! ocean space and time domain variables
24   USE zdf_oce        ! vertical physics: variables
25   USE phycst         ! physical constants
26   USE sbc_oce,  ONLY :   taum
27   !
28   USE in_out_manager ! I/O manager
29   USE iom            ! I/O manager library
30   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
31
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   zdf_ric         ! called by zdfphy.F90
37   PUBLIC   ric_rst         ! called by zdfphy.F90
38   PUBLIC   zdf_ric_init    ! called by nemogcm.F90
39
40   !                        !!* Namelist namzdf_ric : Richardson number dependent Kz *
41   INTEGER  ::   nn_ric      ! coefficient of the parameterization
42   REAL(wp) ::   rn_avmri    ! maximum value of the vertical eddy viscosity
43   REAL(wp) ::   rn_alp      ! coefficient of the parameterization
44   REAL(wp) ::   rn_ekmfc    ! Ekman Factor Coeff
45   REAL(wp) ::   rn_mldmin   ! minimum mixed layer (ML) depth   
46   REAL(wp) ::   rn_mldmax   ! maximum mixed layer depth
47   REAL(wp) ::   rn_wtmix    ! Vertical eddy Diff. in the ML
48   REAL(wp) ::   rn_wvmix    ! Vertical eddy Visc. in the ML
49   LOGICAL  ::   ln_mldw     ! Use or not the MLD parameters
50
51   !! * Substitutions
52#  include "do_loop_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
55   !! $Id$
56   !! Software governed by the CeCILL license (see ./LICENSE)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   SUBROUTINE zdf_ric_init
61      !!----------------------------------------------------------------------
62      !!                 ***  ROUTINE zdf_ric_init  ***
63      !!                   
64      !! ** Purpose :   Initialization of the vertical eddy diffusivity and
65      !!      viscosity coef. for the Richardson number dependent formulation.
66      !!
67      !! ** Method  :   Read the namzdf_ric namelist and check the parameter values
68      !!
69      !! ** input   :   Namelist namzdf_ric
70      !!
71      !! ** Action  :   increase by 1 the nstop flag is setting problem encounter
72      !!----------------------------------------------------------------------
73      INTEGER ::   ji, jj, jk   ! dummy loop indices
74      INTEGER ::   ios          ! Local integer output status for namelist read
75      !!
76      NAMELIST/namzdf_ric/ rn_avmri, rn_alp   , nn_ric  , rn_ekmfc,  &
77         &                rn_mldmin, rn_mldmax, rn_wtmix, rn_wvmix, ln_mldw
78      !!----------------------------------------------------------------------
79      !
80      READ  ( numnam_ref, namzdf_ric, IOSTAT = ios, ERR = 901)
81901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in reference namelist' )
82
83      READ  ( numnam_cfg, namzdf_ric, IOSTAT = ios, ERR = 902 )
84902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_ric in configuration namelist' )
85      IF(lwm) WRITE ( numond, namzdf_ric )
86      !
87      IF(lwp) THEN                   ! Control print
88         WRITE(numout,*)
89         WRITE(numout,*) 'zdf_ric_init : Ri depend vertical mixing scheme'
90         WRITE(numout,*) '~~~~~~~~~~~~'
91         WRITE(numout,*) '   Namelist namzdf_ric : set Kz=F(Ri) parameters'
92         WRITE(numout,*) '      maximum vertical viscosity        rn_avmri  = ', rn_avmri
93         WRITE(numout,*) '      coefficient                       rn_alp    = ', rn_alp
94         WRITE(numout,*) '      exponent                          nn_ric    = ', nn_ric
95         WRITE(numout,*) '      Ekman layer enhanced mixing       ln_mldw   = ', ln_mldw
96         WRITE(numout,*) '         Ekman Factor Coeff             rn_ekmfc  = ', rn_ekmfc
97         WRITE(numout,*) '         minimum mixed layer depth      rn_mldmin = ', rn_mldmin
98         WRITE(numout,*) '         maximum mixed layer depth      rn_mldmax = ', rn_mldmax
99         WRITE(numout,*) '         Vertical eddy Diff. in the ML  rn_wtmix  = ', rn_wtmix
100         WRITE(numout,*) '         Vertical eddy Visc. in the ML  rn_wvmix  = ', rn_wvmix
101      ENDIF
102      !
103      CALL ric_rst( nit000, 'READ' )  !* read or initialize all required files
104      !
105      IF( lwxios ) THEN
106         CALL iom_set_rstw_var_active('avt_k')
107         CALL iom_set_rstw_var_active('avm_k')
108      ENDIF
109   END SUBROUTINE zdf_ric_init
110
111
112   SUBROUTINE zdf_ric( kt, Kmm, p_sh2, p_avm, p_avt )
113      !!----------------------------------------------------------------------
114      !!                 ***  ROUTINE zdfric  ***
115      !!                   
116      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
117      !!                a function of the local richardson number.
118      !!
119      !! ** Method  :   Local richardson number dependent formulation of the
120      !!                vertical eddy viscosity and diffusivity coefficients.
121      !!                The eddy coefficients are given by:
122      !!                    avm = avm0 + avmb
123      !!                    avt = avm0 / (1 + rn_alp*ri)
124      !!                with ri  = N^2 / dz(u)**2
125      !!                         = e3w**2 * rn2/[ mi( dk(uu(:,:,:,Kbb)) )+mj( dk(vv(:,:,:,Kbb)) ) ]
126      !!                    avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric
127      !!                where ri is the before local Richardson number,
128      !!                rn_avmri is the maximum value reaches by avm and avt
129      !!                and rn_alp, nn_ric are adjustable parameters.
130      !!                Typical values : rn_alp=5. and nn_ric=2.
131      !!
132      !!      As second step compute Ekman depth from wind stress forcing
133      !!      and apply namelist provided vertical coeff within this depth.
134      !!      The Ekman depth is:
135      !!              Ustar = SQRT(Taum/rho0)
136      !!              ekd= rn_ekmfc * Ustar / f0
137      !!      Large et al. (1994, eq.24) suggest rn_ekmfc=0.7; however, the derivation
138      !!      of the above equation indicates the value is somewhat arbitrary; therefore
139      !!      we allow the freedom to increase or decrease this value, if the
140      !!      Ekman depth estimate appears too shallow or too deep, respectively.
141      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the
142      !!      namelist
143      !!        N.B. the mask are required for implicit scheme, and surface
144      !!      and bottom value already set in zdfphy.F90
145      !!
146      !! ** Action  :   avm, avt  mixing coeff (inner domain values only)
147      !!
148      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
149      !!              PFJ Lermusiaux 2001.
150      !!----------------------------------------------------------------------
151      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time-step
152      INTEGER                   , INTENT(in   ) ::   Kmm            ! ocean time level index
153      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term
154      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   ! momentum and tracer Kz (w-points)
155      !!
156      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
157      REAL(wp) ::   zcfRi, zav, zustar, zhek    ! local scalars
158      REAL(wp), DIMENSION(jpi,jpj) ::   zh_ekm  ! 2D workspace
159      !!----------------------------------------------------------------------
160      !
161      !                       !==  avm and avt = F(Richardson number)  ==!
162      DO_3D_10_10( 2, jpkm1 )
163         zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) )  )
164         zav   = rn_avmri * zcfRi**nn_ric
165         !                          ! avm and avt coefficients
166         p_avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk)
167         p_avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk)
168      END_3D
169      !
170!!gm BUG <<<<====  This param can't work at low latitude
171!!gm               it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! )
172      !
173      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
174         !
175         DO_2D_00_00
176            zustar = SQRT( taum(ji,jj) * r1_rho0 )
177            zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth
178            zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed range
179         END_2D
180         DO_3D_00_00( 2, jpkm1 )
181            IF( gdept(ji,jj,jk,Kmm) < zh_ekm(ji,jj) ) THEN
182               p_avm(ji,jj,jk) = MAX(  p_avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk)
183               p_avt(ji,jj,jk) = MAX(  p_avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk)
184            ENDIF
185         END_3D
186      ENDIF
187      !
188   END SUBROUTINE zdf_ric
189
190
191   SUBROUTINE ric_rst( kt, cdrw )
192      !!---------------------------------------------------------------------
193      !!                   ***  ROUTINE ric_rst  ***
194      !!                     
195      !! ** Purpose :   Read or write TKE file (en) in restart file
196      !!
197      !! ** Method  :   use of IOM library
198      !!                if the restart does not contain TKE, en is either
199      !!                set to rn_emin or recomputed
200      !!----------------------------------------------------------------------
201      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
202      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
203      !
204      INTEGER ::   jit, jk    ! dummy loop indices
205      INTEGER ::   id1, id2   ! local integers
206      !!----------------------------------------------------------------------
207      !
208      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
209         !                                   ! ---------------
210         !           !* Read the restart file
211         IF( ln_rstart ) THEN
212            id1 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
213            id2 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
214            !
215            IF( MIN( id1, id2 ) > 0 ) THEN         ! restart exists => read it
216               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k, ldxios = lrxios )
217               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k, ldxios = lrxios )
218            ENDIF
219         ENDIF
220         !           !* otherwise Kz already set to the background value in zdf_phy_init
221         !
222      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
223         !                                   ! -------------------
224         IF(lwp) WRITE(numout,*) '---- ric-rst ----'
225         IF( lwxios ) CALL iom_swap(      cwxios_context          )
226         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k, ldxios = lwxios )
227         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k, ldxios = lwxios)
228         IF( lwxios ) CALL iom_swap(      cxios_context          )
229         !
230      ENDIF
231      !
232   END SUBROUTINE ric_rst
233
234   !!======================================================================
235END MODULE zdfric
Note: See TracBrowser for help on using the repository browser.