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/trunk/src/OCE/ZDF – NEMO

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

Last change on this file since 13970 was 13970, checked in by andmirek, 3 years ago

Ticket #2462 into the trunk

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