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

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90 @ 8143

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

#1880 (HPC-09) - step-7: top/bottom drag computed at T-points, zdfbfr.F90 replaced by zdfdrg.F90 + changes in namelist

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