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

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90 @ 3186

Last change on this file since 3186 was 3186, checked in by smasson, 12 years ago

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