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 @ 14982

Last change on this file since 14982 was 14982, checked in by hadcv, 20 months ago

#2665: Various fixes for code enabled with key_qco/key_linssh

  • Property svn:keywords set to Id
File size: 12.0 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#  include "domzgr_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
56   !! $Id$
57   !! Software governed by the CeCILL license (see ./LICENSE)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE zdf_ric_init
62      !!----------------------------------------------------------------------
63      !!                 ***  ROUTINE zdf_ric_init  ***
64      !!
65      !! ** Purpose :   Initialization of the vertical eddy diffusivity and
66      !!      viscosity coef. for the Richardson number dependent formulation.
67      !!
68      !! ** Method  :   Read the namzdf_ric namelist and check the parameter values
69      !!
70      !! ** input   :   Namelist namzdf_ric
71      !!
72      !! ** Action  :   increase by 1 the nstop flag is setting problem encounter
73      !!----------------------------------------------------------------------
74      INTEGER ::   ji, jj, jk   ! dummy loop indices
75      INTEGER ::   ios          ! Local integer output status for namelist read
76      !!
77      NAMELIST/namzdf_ric/ rn_avmri, rn_alp   , nn_ric  , rn_ekmfc,  &
78         &                rn_mldmin, rn_mldmax, rn_wtmix, rn_wvmix, ln_mldw
79      !!----------------------------------------------------------------------
80      !
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      READ  ( numnam_cfg, namzdf_ric, IOSTAT = ios, ERR = 902 )
85902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_ric in configuration namelist' )
86      IF(lwm) WRITE ( numond, namzdf_ric )
87      !
88      IF(lwp) THEN                   ! Control print
89         WRITE(numout,*)
90         WRITE(numout,*) 'zdf_ric_init : Ri depend vertical mixing scheme'
91         WRITE(numout,*) '~~~~~~~~~~~~'
92         WRITE(numout,*) '   Namelist namzdf_ric : set Kz=F(Ri) parameters'
93         WRITE(numout,*) '      maximum vertical viscosity        rn_avmri  = ', rn_avmri
94         WRITE(numout,*) '      coefficient                       rn_alp    = ', rn_alp
95         WRITE(numout,*) '      exponent                          nn_ric    = ', nn_ric
96         WRITE(numout,*) '      Ekman layer enhanced mixing       ln_mldw   = ', ln_mldw
97         WRITE(numout,*) '         Ekman Factor Coeff             rn_ekmfc  = ', rn_ekmfc
98         WRITE(numout,*) '         minimum mixed layer depth      rn_mldmin = ', rn_mldmin
99         WRITE(numout,*) '         maximum mixed layer depth      rn_mldmax = ', rn_mldmax
100         WRITE(numout,*) '         Vertical eddy Diff. in the ML  rn_wtmix  = ', rn_wtmix
101         WRITE(numout,*) '         Vertical eddy Visc. in the ML  rn_wvmix  = ', rn_wvmix
102      ENDIF
103      !
104      CALL ric_rst( nit000, 'READ' )  !* read or initialize all required files
105      !
106   END SUBROUTINE zdf_ric_init
107
108
109   SUBROUTINE zdf_ric( kt, Kmm, p_sh2, p_avm, p_avt )
110      !!----------------------------------------------------------------------
111      !!                 ***  ROUTINE zdfric  ***
112      !!
113      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
114      !!                a function of the local richardson number.
115      !!
116      !! ** Method  :   Local richardson number dependent formulation of the
117      !!                vertical eddy viscosity and diffusivity coefficients.
118      !!                The eddy coefficients are given by:
119      !!                    avm = avm0 + avmb
120      !!                    avt = avm0 / (1 + rn_alp*ri)
121      !!                with ri  = N^2 / dz(u)**2
122      !!                         = e3w**2 * rn2/[ mi( dk(uu(:,:,:,Kbb)) )+mj( dk(vv(:,:,:,Kbb)) ) ]
123      !!                    avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric
124      !!                where ri is the before local Richardson number,
125      !!                rn_avmri is the maximum value reaches by avm and avt
126      !!                and rn_alp, nn_ric are adjustable parameters.
127      !!                Typical values : rn_alp=5. and nn_ric=2.
128      !!
129      !!      As second step compute Ekman depth from wind stress forcing
130      !!      and apply namelist provided vertical coeff within this depth.
131      !!      The Ekman depth is:
132      !!              Ustar = SQRT(Taum/rho0)
133      !!              ekd= rn_ekmfc * Ustar / f0
134      !!      Large et al. (1994, eq.24) suggest rn_ekmfc=0.7; however, the derivation
135      !!      of the above equation indicates the value is somewhat arbitrary; therefore
136      !!      we allow the freedom to increase or decrease this value, if the
137      !!      Ekman depth estimate appears too shallow or too deep, respectively.
138      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the
139      !!      namelist
140      !!        N.B. the mask are required for implicit scheme, and surface
141      !!      and bottom value already set in zdfphy.F90
142      !!
143      !! ** Action  :   avm, avt  mixing coeff (inner domain values only)
144      !!
145      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
146      !!              PFJ Lermusiaux 2001.
147      !!----------------------------------------------------------------------
148      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time-step
149      INTEGER                             , INTENT(in   ) ::   Kmm            ! ocean time level index
150      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   p_sh2          ! shear production term
151      REAL(wp), DIMENSION(:,:,:)          , INTENT(inout) ::   p_avm, p_avt   ! momentum and tracer Kz (w-points)
152      !!
153      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
154      REAL(wp) ::   zcfRi, zav, zustar, zhek    ! local scalars
155      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zh_ekm  ! 2D workspace
156      !!----------------------------------------------------------------------
157      !
158      !                       !==  avm and avt = F(Richardson number)  ==!
159      DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, jpkm1 )       ! coefficient = F(richardson number) (avm-weighted Ri)
160         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 ) )  )
161         zav   = rn_avmri * zcfRi**nn_ric
162         !                          ! avm and avt coefficients
163         p_avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk)
164         p_avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk)
165      END_3D
166      !
167!!gm BUG <<<<====  This param can't work at low latitude
168!!gm               it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! )
169      !
170      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
171         !
172         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
173            zustar = SQRT( taum(ji,jj) * r1_rho0 )
174            zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth
175            zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed range
176         END_2D
177         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )   !* minimum mixing coeff. within the Ekman layer
178            IF( gdept(ji,jj,jk,Kmm) < zh_ekm(ji,jj) ) THEN
179               p_avm(ji,jj,jk) = MAX(  p_avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk)
180               p_avt(ji,jj,jk) = MAX(  p_avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk)
181            ENDIF
182         END_3D
183      ENDIF
184      !
185   END SUBROUTINE zdf_ric
186
187
188   SUBROUTINE ric_rst( kt, cdrw )
189      !!---------------------------------------------------------------------
190      !!                   ***  ROUTINE ric_rst  ***
191      !!
192      !! ** Purpose :   Read or write TKE file (en) in restart file
193      !!
194      !! ** Method  :   use of IOM library
195      !!                if the restart does not contain TKE, en is either
196      !!                set to rn_emin or recomputed
197      !!----------------------------------------------------------------------
198      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
199      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
200      !
201      INTEGER ::   jit, jk    ! dummy loop indices
202      INTEGER ::   id1, id2   ! local integers
203      !!----------------------------------------------------------------------
204      !
205      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
206         !                                   ! ---------------
207         !           !* Read the restart file
208         IF( ln_rstart ) THEN
209            id1 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
210            id2 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
211            !
212            IF( MIN( id1, id2 ) > 0 ) THEN         ! restart exists => read it
213               CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k )
214               CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k )
215            ENDIF
216         ENDIF
217         !           !* otherwise Kz already set to the background value in zdf_phy_init
218         !
219      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
220         !                                   ! -------------------
221         IF(lwp) WRITE(numout,*) '---- ric-rst ----'
222         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
223         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k)
224         !
225      ENDIF
226      !
227   END SUBROUTINE ric_rst
228
229   !!======================================================================
230END MODULE zdfric
Note: See TracBrowser for help on using the repository browser.