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

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

#1880 (HPC-09): add zdfphy (the ZDF manager) + remove all key_...
  • Property svn:keywords set to Id
File size: 15.2 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   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!   zdf_ric       : update momentum and tracer Kz from the Richardson number
18   !!   zdf_ric_init  : initialization, namelist read, & parameters control
19   !!----------------------------------------------------------------------
20   USE oce            ! ocean dynamics and tracers variables
21   USE dom_oce        ! ocean space and time domain variables
22   USE zdf_oce        ! ocean vertical physics
23   USE in_out_manager ! I/O manager
24   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
25   USE lib_mpp        ! MPP library
26   USE wrk_nemo       ! work arrays
27   USE timing         ! Timing
28   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
29
30   USE eosbn2, ONLY : neos
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   zdf_ric         ! called by step.F90
36   PUBLIC   zdf_ric_init    ! called by opa.F90
37
38   !                           !!* Namelist namzdf_ric : Richardson number dependent Kz *
39   INTEGER  ::   nn_ric         ! coefficient of the parameterization
40   REAL(wp) ::   rn_avmri       ! maximum value of the vertical eddy viscosity
41   REAL(wp) ::   rn_alp         ! coefficient of the parameterization
42   REAL(wp) ::   rn_ekmfc       ! Ekman Factor Coeff
43   REAL(wp) ::   rn_mldmin      ! minimum mixed layer (ML) depth   
44   REAL(wp) ::   rn_mldmax      ! maximum mixed layer depth
45   REAL(wp) ::   rn_wtmix       ! Vertical eddy Diff. in the ML
46   REAL(wp) ::   rn_wvmix       ! Vertical eddy Visc. in the ML
47   LOGICAL  ::   ln_mldw        ! Use or not the MLD parameters
48
49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tmric   !: coef. for the horizontal mean at t-point
50
51   !! * Substitutions
52#  include "vectopt_loop_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
55   !! $Id$
56   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   INTEGER FUNCTION zdf_ric_alloc()
61      !!----------------------------------------------------------------------
62      !!                 ***  FUNCTION zdf_ric_alloc  ***
63      !!----------------------------------------------------------------------
64      ALLOCATE( tmric(jpi,jpj,jpk)   , STAT= zdf_ric_alloc )
65      !
66      IF( lk_mpp             )   CALL mpp_sum ( zdf_ric_alloc )
67      IF( zdf_ric_alloc /= 0 )   CALL ctl_warn('zdf_ric_alloc: failed to allocate arrays')
68   END FUNCTION zdf_ric_alloc
69
70
71   SUBROUTINE zdf_ric( kt )
72      !!----------------------------------------------------------------------
73      !!                 ***  ROUTINE zdfric  ***
74      !!                   
75      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
76      !!                a function of the local richardson number.
77      !!
78      !! ** Method  :   Local richardson number dependent formulation of the
79      !!                vertical eddy viscosity and diffusivity coefficients.
80      !!                The eddy coefficients are given by:
81      !!                    avm = avm0 + avmb
82      !!                    avt = avm0 / (1 + rn_alp*ri)
83      !!                with ri  = N^2 / dz(u)**2
84      !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]
85      !!                    avm0= rn_avmri / (1 + rn_alp*ri)**nn_ric
86      !!      Where ri is the before local Richardson number,
87      !!            rn_avmri is the maximum value reaches by avm and avt
88      !!            avmb and avtb are the background (or minimum) values
89      !!            and rn_alp, nn_ric are adjustable parameters.
90      !!      Typical values used are : avm0=1.e-2 m2/s, avmb=1.e-6 m2/s
91      !!      avtb=1.e-7 m2/s, rn_alp=5. and nn_ric=2.
92      !!      a numerical threshold is impose on the vertical shear (1.e-20)
93      !!      As second step compute Ekman depth from wind stress forcing
94      !!      and apply namelist provided vertical coeff within this depth.
95      !!      The Ekman depth is:
96      !!              Ustar = SQRT(Taum/rho0)
97      !!              ekd= rn_ekmfc * Ustar / f0
98      !!      Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation
99      !!      of the above equation indicates the value is somewhat arbitrary; therefore
100      !!      we allow the freedom to increase or decrease this value, if the
101      !!      Ekman depth estimate appears too shallow or too deep, respectively.
102      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the
103      !!      namelist
104      !!        N.B. the mask are required for implicit scheme, and surface
105      !!      and bottom value already set in zdfphy.F90
106      !!
107      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
108      !!              PFJ Lermusiaux 2001.
109      !!----------------------------------------------------------------------
110      USE phycst,   ONLY:   rsmall,rau0
111      USE sbc_oce,  ONLY:   taum
112      !!
113      INTEGER, INTENT( in ) ::   kt                           ! ocean time-step
114      !!
115      INTEGER  ::   ji, jj, jk                                ! dummy loop indices
116      REAL(wp) ::   zcoef, zdku, zdkv, zri, z05alp, zflageos  ! temporary scalars
117      REAL(wp) ::   zrhos, zustar
118      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, ekm_dep 
119      !!----------------------------------------------------------------------
120      !
121      IF( nn_timing == 1 )  CALL timing_start('zdf_ric')
122      !
123      CALL wrk_alloc( jpi,jpj, zwx, ekm_dep )
124      !                                                ! ===============
125      DO jk = 2, jpkm1                                 ! Horizontal slab
126         !                                             ! ===============
127         ! Richardson number (put in zwx(ji,jj))
128         ! -----------------
129         DO jj = 2, jpjm1
130            DO ji = fs_2, fs_jpim1
131               zcoef = 0.5 / e3w_n(ji,jj,jk)
132               !                                            ! shear of horizontal velocity
133               zdku = zcoef * (  ub(ji-1,jj,jk-1) + ub(ji,jj,jk-1)   &
134                  &             -ub(ji-1,jj,jk  ) - ub(ji,jj,jk  )  )
135               zdkv = zcoef * (  vb(ji,jj-1,jk-1) + vb(ji,jj,jk-1)   &
136                  &             -vb(ji,jj-1,jk  ) - vb(ji,jj,jk  )  )
137               !                                            ! richardson number (minimum value set to zero)
138               zri = rn2(ji,jj,jk) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
139               zwx(ji,jj) = MAX( zri, 0.e0 )
140            END DO
141         END DO
142         CALL lbc_lnk( zwx, 'W', 1. )                       ! Boundary condition   (sign unchanged)
143
144         ! Vertical eddy viscosity and diffusivity coefficients
145         ! -------------------------------------------------------
146         z05alp = 0.5_wp * rn_alp
147         DO jj = 1, jpjm1                                   ! Eddy viscosity coefficients (avm)
148            DO ji = 1, fs_jpim1
149               avmu(ji,jj,jk) = umask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji+1,jj)+zwx(ji,jj) ) )**nn_ric
150               avmv(ji,jj,jk) = vmask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji,jj+1)+zwx(ji,jj) ) )**nn_ric
151            END DO
152         END DO
153         DO jj = 2, jpjm1                                   ! Eddy diffusivity coefficients (avt)
154            DO ji = fs_2, fs_jpim1
155               avt(ji,jj,jk) = tmric(ji,jj,jk) / ( 1._wp + rn_alp * zwx(ji,jj) )           &
156                  &                            * (  avmu(ji,jj,jk) + avmu(ji-1,jj,jk)      &
157                  &                               + avmv(ji,jj,jk) + avmv(ji,jj-1,jk)  )   &
158                  &          + avtb(jk) * tmask(ji,jj,jk)
159            END DO
160         END DO
161         DO jj = 2, jpjm1                                   ! Add the background coefficient on eddy viscosity
162            DO ji = fs_2, fs_jpim1
163               avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk)
164               avmv(ji,jj,jk) = avmv(ji,jj,jk) + avmb(jk) * vmask(ji,jj,jk)
165            END DO
166         END DO
167         !                                             ! ===============
168      END DO                                           !   End of slab
169      !                                                ! ===============
170      !
171      IF( ln_mldw ) THEN
172
173      !  Compute Ekman depth from wind stress forcing.
174      ! -------------------------------------------------------
175      zflageos = ( 0.5 + SIGN( 0.5, neos - 1. ) ) * rau0
176      DO jj = 2, jpjm1
177            DO ji = fs_2, fs_jpim1
178            zrhos          = rhop(ji,jj,1) + zflageos * ( 1. - tmask(ji,jj,1) )
179            zustar         = SQRT( taum(ji,jj) / ( zrhos +  rsmall ) )
180            ekm_dep(ji,jj) = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )
181            ekm_dep(ji,jj) = MAX(ekm_dep(ji,jj),rn_mldmin) ! Minimun allowed
182            ekm_dep(ji,jj) = MIN(ekm_dep(ji,jj),rn_mldmax) ! Maximum allowed
183         END DO
184      END DO
185
186      ! In the first model level vertical diff/visc coeff.s
187      ! are always equal to the namelist values rn_wtmix/rn_wvmix
188      ! -------------------------------------------------------
189      DO jj = 2, jpjm1
190         DO ji = fs_2, fs_jpim1
191            avmv(ji,jj,1) = MAX( avmv(ji,jj,1), rn_wvmix )
192            avmu(ji,jj,1) = MAX( avmu(ji,jj,1), rn_wvmix )
193            avt( ji,jj,1) = MAX(  avt(ji,jj,1), rn_wtmix )
194         END DO
195      END DO
196
197      !  Force the vertical mixing coef within the Ekman depth
198      ! -------------------------------------------------------
199      DO jk = 2, jpkm1
200         DO jj = 2, jpjm1
201            DO ji = fs_2, fs_jpim1
202               IF( gdept_n(ji,jj,jk) < ekm_dep(ji,jj) ) THEN
203                  avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), rn_wvmix )
204                  avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), rn_wvmix )
205                  avt( ji,jj,jk) = MAX(  avt(ji,jj,jk), rn_wtmix )
206               ENDIF
207            END DO
208         END DO
209      END DO
210
211      DO jk = 1, jpkm1               
212         DO jj = 2, jpjm1
213            DO ji = fs_2, fs_jpim1
214               avmv(ji,jj,jk) = avmv(ji,jj,jk) * vmask(ji,jj,jk)
215               avmu(ji,jj,jk) = avmu(ji,jj,jk) * umask(ji,jj,jk)
216               avt( ji,jj,jk) = avt( ji,jj,jk) * tmask(ji,jj,jk)
217            END DO
218         END DO
219      END DO
220
221     ENDIF
222
223      CALL lbc_lnk( avt , 'W', 1. )                         ! Boundary conditions   (unchanged sign)
224      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
225      !
226      CALL wrk_dealloc( jpi,jpj, zwx, ekm_dep )
227      !
228      IF( nn_timing == 1 )  CALL timing_stop('zdf_ric')
229      !
230   END SUBROUTINE zdf_ric
231
232
233   SUBROUTINE zdf_ric_init
234      !!----------------------------------------------------------------------
235      !!                 ***  ROUTINE zdfbfr_init  ***
236      !!                   
237      !! ** Purpose :   Initialization of the vertical eddy diffusivity and
238      !!      viscosity coef. for the Richardson number dependent formulation.
239      !!
240      !! ** Method  :   Read the namzdf_ric namelist and check the parameter values
241      !!
242      !! ** input   :   Namelist namzdf_ric
243      !!
244      !! ** Action  :   increase by 1 the nstop flag is setting problem encounter
245      !!----------------------------------------------------------------------
246      INTEGER :: ji, jj, jk   ! dummy loop indices
247      INTEGER ::   ios        ! Local integer output status for namelist read
248      !!
249      NAMELIST/namzdf_ric/ rn_avmri, rn_alp   , nn_ric  , rn_ekmfc,  &
250         &                rn_mldmin, rn_mldmax, rn_wtmix, rn_wvmix, ln_mldw
251      !!----------------------------------------------------------------------
252      !
253      REWIND( numnam_ref )              ! Namelist namzdf_ric in reference namelist : Vertical diffusion Kz depends on Richardson number
254      READ  ( numnam_ref, namzdf_ric, IOSTAT = ios, ERR = 901)
255901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in reference namelist', lwp )
256
257      REWIND( numnam_cfg )              ! Namelist namzdf_ric in configuration namelist : Vertical diffusion Kz depends on Richardson number
258      READ  ( numnam_cfg, namzdf_ric, IOSTAT = ios, ERR = 902 )
259902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in configuration namelist', lwp )
260      IF(lwm) WRITE ( numond, namzdf_ric )
261      !
262      IF(lwp) THEN                   ! Control print
263         WRITE(numout,*)
264         WRITE(numout,*) 'zdf_ric : Ri depend vertical mixing scheme'
265         WRITE(numout,*) '~~~~~~~'
266         WRITE(numout,*) '   Namelist namzdf_ric : set Kz(Ri) parameters'
267         WRITE(numout,*) '      maximum vertical viscosity     rn_avmri  = ', rn_avmri
268         WRITE(numout,*) '      coefficient                    rn_alp    = ', rn_alp
269         WRITE(numout,*) '      coefficient                    nn_ric    = ', nn_ric
270         WRITE(numout,*) '      Ekman Factor Coeff             rn_ekmfc  = ', rn_ekmfc
271         WRITE(numout,*) '      minimum mixed layer depth      rn_mldmin = ', rn_mldmin
272         WRITE(numout,*) '      maximum mixed layer depth      rn_mldmax = ', rn_mldmax
273         WRITE(numout,*) '      Vertical eddy Diff. in the ML  rn_wtmix  = ', rn_wtmix
274         WRITE(numout,*) '      Vertical eddy Visc. in the ML  rn_wvmix  = ', rn_wvmix
275         WRITE(numout,*) '      Use the MLD parameterization   ln_mldw   = ', ln_mldw
276      ENDIF
277      !
278      !                              ! allocate zdfric arrays
279      IF( zdf_ric_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ric_init : unable to allocate arrays' )
280      !
281      DO jk = 1, jpk                 ! weighting mean array tmric for 4 T-points
282         DO jj = 2, jpj              ! which accounts for coastal boundary conditions           
283            DO ji = 2, jpi
284               tmric(ji,jj,jk) =  tmask(ji,jj,jk)                                  &
285                  &            / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
286                  &                      + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
287            END DO
288         END DO
289      END DO
290      tmric(:,1,:) = 0._wp
291      !
292      DO jk = 1, jpk                 ! Initialization of vertical eddy coef. to the background value
293         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
294         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
295         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
296      END DO
297      !
298   END SUBROUTINE zdf_ric_init
299
300   !!======================================================================
301END MODULE zdfric
Note: See TracBrowser for help on using the repository browser.