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/UKMO/dev_r10171_test_crs_AMM7/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/dev_r10171_test_crs_AMM7/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90 @ 10207

Last change on this file since 10207 was 10207, checked in by cmao, 6 years ago

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