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.
zdfric.F90 in NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/ZDF – NEMO

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/ZDF/zdfric.F90 @ 11831

Last change on this file since 11831 was 11831, checked in by laurent, 4 years ago

Update the branch to r11830 of the trunk!

  • Property svn:keywords set to Id
File size: 12.8 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 "vectopt_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      REWIND( numnam_ref )              ! Namelist namzdf_ric in reference namelist : Vertical diffusion Kz depends on Richardson number
81      READ  ( numnam_ref, namzdf_ric, IOSTAT = ios, ERR = 901)
82901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in reference namelist' )
83
84      REWIND( numnam_cfg )              ! Namelist namzdf_ric in configuration namelist : Vertical diffusion Kz depends on Richardson number
85      READ  ( numnam_cfg, namzdf_ric, IOSTAT = ios, ERR = 902 )
86902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_ric in configuration namelist' )
87      IF(lwm) WRITE ( numond, namzdf_ric )
88      !
89      IF(lwp) THEN                   ! Control print
90         WRITE(numout,*)
91         WRITE(numout,*) 'zdf_ric_init : Ri depend vertical mixing scheme'
92         WRITE(numout,*) '~~~~~~~~~~~~'
93         WRITE(numout,*) '   Namelist namzdf_ric : set Kz=F(Ri) parameters'
94         WRITE(numout,*) '      maximum vertical viscosity        rn_avmri  = ', rn_avmri
95         WRITE(numout,*) '      coefficient                       rn_alp    = ', rn_alp
96         WRITE(numout,*) '      exponent                          nn_ric    = ', nn_ric
97         WRITE(numout,*) '      Ekman layer enhanced mixing       ln_mldw   = ', ln_mldw
98         WRITE(numout,*) '         Ekman Factor Coeff             rn_ekmfc  = ', rn_ekmfc
99         WRITE(numout,*) '         minimum mixed layer depth      rn_mldmin = ', rn_mldmin
100         WRITE(numout,*) '         maximum mixed layer depth      rn_mldmax = ', rn_mldmax
101         WRITE(numout,*) '         Vertical eddy Diff. in the ML  rn_wtmix  = ', rn_wtmix
102         WRITE(numout,*) '         Vertical eddy Visc. in the ML  rn_wvmix  = ', rn_wvmix
103      ENDIF
104      !
105      CALL ric_rst( nit000, 'READ' )  !* read or initialize all required files
106      !
107      IF( lwxios ) THEN
108         CALL iom_set_rstw_var_active('avt_k')
109         CALL iom_set_rstw_var_active('avm_k')
110      ENDIF
111   END SUBROUTINE zdf_ric_init
112
113
114   SUBROUTINE zdf_ric( kt, pdept, p_sh2, p_avm, p_avt )
115      !!----------------------------------------------------------------------
116      !!                 ***  ROUTINE zdfric  ***
117      !!                   
118      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
119      !!                a function of the local richardson number.
120      !!
121      !! ** Method  :   Local richardson number dependent formulation of the
122      !!                vertical eddy viscosity and diffusivity coefficients.
123      !!                The eddy coefficients are given by:
124      !!                    avm = avm0 + avmb
125      !!                    avt = avm0 / (1 + rn_alp*ri)
126      !!                with ri  = N^2 / dz(u)**2
127      !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]
128      !!                    avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric
129      !!                where ri is the before local Richardson number,
130      !!                rn_avmri is the maximum value reaches by avm and avt
131      !!                and rn_alp, nn_ric are adjustable parameters.
132      !!                Typical values : rn_alp=5. and nn_ric=2.
133      !!
134      !!      As second step compute Ekman depth from wind stress forcing
135      !!      and apply namelist provided vertical coeff within this depth.
136      !!      The Ekman depth is:
137      !!              Ustar = SQRT(Taum/rho0)
138      !!              ekd= rn_ekmfc * Ustar / f0
139      !!      Large et al. (1994, eq.24) suggest rn_ekmfc=0.7; however, the derivation
140      !!      of the above equation indicates the value is somewhat arbitrary; therefore
141      !!      we allow the freedom to increase or decrease this value, if the
142      !!      Ekman depth estimate appears too shallow or too deep, respectively.
143      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the
144      !!      namelist
145      !!        N.B. the mask are required for implicit scheme, and surface
146      !!      and bottom value already set in zdfphy.F90
147      !!
148      !! ** Action  :   avm, avt  mixing coeff (inner domain values only)
149      !!
150      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
151      !!              PFJ Lermusiaux 2001.
152      !!----------------------------------------------------------------------
153      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time-step
154      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdept          ! depth of t-point  [m]
155      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term
156      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   ! momentum and tracer Kz (w-points)
157      !!
158      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
159      REAL(wp) ::   zcfRi, zav, zustar, zhek    ! local scalars
160      REAL(wp), DIMENSION(jpi,jpj) ::   zh_ekm  ! 2D workspace
161      !!----------------------------------------------------------------------
162      !
163      !                       !==  avm and avt = F(Richardson number)  ==!
164      DO jk = 2, jpkm1
165         DO jj = 1, jpjm1
166            DO ji = 1, jpim1              ! coefficient = F(richardson number) (avm-weighted Ri)
167               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 ) )  )
168               zav   = rn_avmri * zcfRi**nn_ric
169               !                          ! avm and avt coefficients
170               p_avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk)
171               p_avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk)
172            END DO
173         END DO
174      END DO
175      !
176!!gm BUG <<<<====  This param can't work at low latitude
177!!gm               it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! )
178      !
179      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
180         !
181         DO jj = 2, jpjm1        !* Ekman depth
182            DO ji = 2, jpim1
183               zustar = SQRT( taum(ji,jj) * r1_rau0 )
184               zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth
185               zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed range
186            END DO
187         END DO
188         DO jk = 2, jpkm1        !* minimum mixing coeff. within the Ekman layer
189            DO jj = 2, jpjm1
190               DO ji = 2, jpim1
191                  IF( pdept(ji,jj,jk) < zh_ekm(ji,jj) ) THEN
192                     p_avm(ji,jj,jk) = MAX(  p_avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk)
193                     p_avt(ji,jj,jk) = MAX(  p_avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk)
194                  ENDIF
195               END DO
196            END DO
197         END DO
198      ENDIF
199      !
200   END SUBROUTINE zdf_ric
201
202
203   SUBROUTINE ric_rst( kt, cdrw )
204      !!---------------------------------------------------------------------
205      !!                   ***  ROUTINE ric_rst  ***
206      !!                     
207      !! ** Purpose :   Read or write TKE file (en) in restart file
208      !!
209      !! ** Method  :   use of IOM library
210      !!                if the restart does not contain TKE, en is either
211      !!                set to rn_emin or recomputed
212      !!----------------------------------------------------------------------
213      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
214      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
215      !
216      INTEGER ::   jit, jk    ! dummy loop indices
217      INTEGER ::   id1, id2   ! local integers
218      !!----------------------------------------------------------------------
219      !
220      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
221         !                                   ! ---------------
222         !           !* Read the restart file
223         IF( ln_rstart ) THEN
224            id1 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
225            id2 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
226            !
227            IF( MIN( id1, id2 ) > 0 ) THEN         ! restart exists => read it
228               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k, ldxios = lrxios )
229               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k, ldxios = lrxios )
230            ENDIF
231         ENDIF
232         !           !* otherwise Kz already set to the background value in zdf_phy_init
233         !
234      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
235         !                                   ! -------------------
236         IF(lwp) WRITE(numout,*) '---- ric-rst ----'
237         IF( lwxios ) CALL iom_swap(      cwxios_context          )
238         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k, ldxios = lwxios )
239         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k, ldxios = lwxios)
240         IF( lwxios ) CALL iom_swap(      cxios_context          )
241         !
242      ENDIF
243      !
244   END SUBROUTINE ric_rst
245
246   !!======================================================================
247END MODULE zdfric
Note: See TracBrowser for help on using the repository browser.