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

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

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