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

Last change on this file since 3115 was 3104, checked in by cetlod, 13 years ago

dev_LOCEAN_CMCC_INGV_MERCATOR_2011:add in changes dev_MERCATOR_INGV_2011_MERGE into the new branch
  • 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 eosbn2, ONLY : nn_eos
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   zdf_ric         ! called by step.F90
35   PUBLIC   zdf_ric_init    ! called by opa.F90
36
37   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .TRUE.   !: Richardson vertical mixing flag
38
39   !                                    !!* Namelist namzdf_ric : Richardson number dependent Kz *
40   INTEGER  ::   nn_ric   = 2            ! coefficient of the parameterization
41   REAL(wp) ::   rn_avmri = 100.e-4_wp   ! maximum value of the vertical eddy viscosity
42   REAL(wp) ::   rn_alp   =   5._wp      ! coefficient of the parameterization
43   REAL(wp) ::   rn_ekmfc =   0.7_wp     ! Ekman Factor Coeff
44   REAL(wp) ::   rn_mldmin=   1.0_wp     ! minimum mixed layer (ML) depth   
45   REAL(wp) ::   rn_mldmax=1000.0_wp     ! maximum mixed layer depth
46   REAL(wp) ::   rn_wtmix =  10.0_wp     ! Vertical eddy Diff. in the ML
47   REAL(wp) ::   rn_wvmix =  10.0_wp     ! Vertical eddy Visc. in the ML
48   LOGICAL  ::   ln_mldw  = .TRUE.       ! Use or not the MLD parameters
49
50   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tmric   !: coef. for the horizontal mean at t-point
51
52   !! * Substitutions
53#  include "domzgr_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
56   !! $Id$
57   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   INTEGER FUNCTION zdf_ric_alloc()
62      !!----------------------------------------------------------------------
63      !!                 ***  FUNCTION zdf_ric_alloc  ***
64      !!----------------------------------------------------------------------
65      ALLOCATE( tmric(jpi,jpj,jpk)   , STAT= zdf_ric_alloc )
66      !
67      IF( lk_mpp             )   CALL mpp_sum ( zdf_ric_alloc )
68      IF( zdf_ric_alloc /= 0 )   CALL ctl_warn('zdf_ric_alloc: failed to allocate arrays')
69   END FUNCTION zdf_ric_alloc
70
71
72   SUBROUTINE zdf_ric( kt )
73      !!----------------------------------------------------------------------
74      !!                 ***  ROUTINE zdfric  ***
75      !!                   
76      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
77      !!                a function of the local richardson number.
78      !!
79      !! ** Method  :   Local richardson number dependent formulation of the
80      !!                vertical eddy viscosity and diffusivity coefficients.
81      !!                The eddy coefficients are given by:
82      !!                    avm = avm0 + avmb
83      !!                    avt = avm0 / (1 + rn_alp*ri)
84      !!                with ri  = N^2 / dz(u)**2
85      !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]
86      !!                    avm0= rn_avmri / (1 + rn_alp*ri)**nn_ric
87      !!      Where ri is the before local Richardson number,
88      !!            rn_avmri is the maximum value reaches by avm and avt
89      !!            avmb and avtb are the background (or minimum) values
90      !!            and rn_alp, nn_ric are adjustable parameters.
91      !!      Typical values used are : avm0=1.e-2 m2/s, avmb=1.e-6 m2/s
92      !!      avtb=1.e-7 m2/s, rn_alp=5. and nn_ric=2.
93      !!      a numerical threshold is impose on the vertical shear (1.e-20)
94      !!      As second step compute Ekman depth from wind stress forcing
95      !!      and apply namelist provided vertical coeff within this depth.
96      !!      The Ekman depth is:
97      !!              Ustar = SQRT(Taum/rho0)
98      !!              ekd= rn_ekmfc * Ustar / f0
99      !!      Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation
100      !!      of the above equation indicates the value is somewhat arbitrary; therefore
101      !!      we allow the freedom to increase or decrease this value, if the
102      !!      Ekman depth estimate appears too shallow or too deep, respectively.
103      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the
104      !!      namelist
105      !!        N.B. the mask are required for implicit scheme, and surface
106      !!      and bottom value already set in zdfini.F90
107      !!
108      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
109      !!              PFJ Lermusiaux 2001.
110      !!----------------------------------------------------------------------
111      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
112      USE wrk_nemo, ONLY:   zwx => wrk_2d_1         ! 2D workspace
113      USE wrk_nemo, ONLY:   ekm_dep => wrk_2d_2     ! 2D workspace
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      !!----------------------------------------------------------------------
123
124      IF( wrk_in_use(2, 1,2) ) THEN
125         CALL ctl_stop('zdf_ric : requested workspace array unavailable')   ;   RETURN
126      ENDIF
127      !                                                ! ===============
128      DO jk = 2, jpkm1                                 ! Horizontal slab
129         !                                             ! ===============
130         ! Richardson number (put in zwx(ji,jj))
131         ! -----------------
132         DO jj = 2, jpjm1
133            DO ji = 2, jpim1
134               zcoef = 0.5 / fse3w(ji,jj,jk)
135               !                                            ! shear of horizontal velocity
136               zdku = zcoef * (  ub(ji-1,jj,jk-1) + ub(ji,jj,jk-1)   &
137                  &             -ub(ji-1,jj,jk  ) - ub(ji,jj,jk  )  )
138               zdkv = zcoef * (  vb(ji,jj-1,jk-1) + vb(ji,jj,jk-1)   &
139                  &             -vb(ji,jj-1,jk  ) - vb(ji,jj,jk  )  )
140               !                                            ! richardson number (minimum value set to zero)
141               zri = rn2(ji,jj,jk) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
142               zwx(ji,jj) = MAX( zri, 0.e0 )
143            END DO
144         END DO
145         CALL lbc_lnk( zwx, 'W', 1. )                       ! Boundary condition   (sign unchanged)
146
147         ! Vertical eddy viscosity and diffusivity coefficients
148         ! -------------------------------------------------------
149         z05alp = 0.5_wp * rn_alp
150         DO jj = 1, jpjm1                                   ! Eddy viscosity coefficients (avm)
151            DO ji = 1, jpim1
152               avmu(ji,jj,jk) = umask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji+1,jj)+zwx(ji,jj) ) )**nn_ric
153               avmv(ji,jj,jk) = vmask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji,jj+1)+zwx(ji,jj) ) )**nn_ric
154            END DO
155         END DO
156         DO jj = 2, jpjm1                                   ! Eddy diffusivity coefficients (avt)
157            DO ji = 2, jpim1
158               avt(ji,jj,jk) = tmric(ji,jj,jk) / ( 1._wp + rn_alp * zwx(ji,jj) )           &
159                  &                            * (  avmu(ji,jj,jk) + avmu(ji-1,jj,jk)      &
160                  &                               + avmv(ji,jj,jk) + avmv(ji,jj-1,jk)  )   &
161                  &          + avtb(jk) * tmask(ji,jj,jk)
162               !                                            ! Add the background coefficient on eddy viscosity
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, nn_eos - 1. ) ) * rau0
176      DO jj = 1, jpj
177         DO ji = 1, jpi
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(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 = 1, jpj
190         DO ji = 1, jpi
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 = 1, jpj
201            DO ji = 1, jpi
202               IF( fsdept(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 = 1, jpj
213            DO ji = 1, jpi
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      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('zdf_ric: failed to release workspace array')
227      !
228   END SUBROUTINE zdf_ric
229
230
231   SUBROUTINE zdf_ric_init
232      !!----------------------------------------------------------------------
233      !!                 ***  ROUTINE zdfbfr_init  ***
234      !!                   
235      !! ** Purpose :   Initialization of the vertical eddy diffusivity and
236      !!      viscosity coef. for the Richardson number dependent formulation.
237      !!
238      !! ** Method  :   Read the namzdf_ric namelist and check the parameter values
239      !!
240      !! ** input   :   Namelist namzdf_ric
241      !!
242      !! ** Action  :   increase by 1 the nstop flag is setting problem encounter
243      !!----------------------------------------------------------------------
244      INTEGER :: ji, jj, jk   ! dummy loop indices
245      !!
246      NAMELIST/namzdf_ric/ rn_avmri, rn_alp   , nn_ric  , rn_ekmfc,  &
247         &                rn_mldmin, rn_mldmax, rn_wtmix, rn_wvmix, ln_mldw
248      !!----------------------------------------------------------------------
249      !
250      REWIND( numnam )               ! Read Namelist namzdf_ric : richardson number dependent Kz
251      READ  ( numnam, namzdf_ric )
252      !
253      IF(lwp) THEN                   ! Control print
254         WRITE(numout,*)
255         WRITE(numout,*) 'zdf_ric : Ri depend vertical mixing scheme'
256         WRITE(numout,*) '~~~~~~~'
257         WRITE(numout,*) '   Namelist namzdf_ric : set Kz(Ri) parameters'
258         WRITE(numout,*) '      maximum vertical viscosity     rn_avmri  = ', rn_avmri
259         WRITE(numout,*) '      coefficient                    rn_alp    = ', rn_alp
260         WRITE(numout,*) '      coefficient                    nn_ric    = ', nn_ric
261         WRITE(numout,*) '      Ekman Factor Coeff             rn_ekmfc  = ', rn_ekmfc
262         WRITE(numout,*) '      minimum mixed layer depth      rn_mldmin = ', rn_mldmin
263         WRITE(numout,*) '      maximum mixed layer depth      rn_mldmax = ', rn_mldmax
264         WRITE(numout,*) '      Vertical eddy Diff. in the ML  rn_wtmix  = ', rn_wtmix
265         WRITE(numout,*) '      Vertical eddy Visc. in the ML  rn_wvmix  = ', rn_wvmix
266         WRITE(numout,*) '      Use the MLD parameterization   ln_mldw   = ', ln_mldw
267      ENDIF
268      !
269      !                              ! allocate zdfric arrays
270      IF( zdf_ric_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ric_init : unable to allocate arrays' )
271      !
272      DO jk = 1, jpk                 ! weighting mean array tmric for 4 T-points
273         DO jj = 2, jpj              ! which accounts for coastal boundary conditions           
274            DO ji = 2, jpi
275               tmric(ji,jj,jk) =  tmask(ji,jj,jk)                                  &
276                  &            / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
277                  &                      + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
278            END DO
279         END DO
280      END DO
281      tmric(:,1,:) = 0._wp
282      !
283      DO jk = 1, jpk                 ! Initialization of vertical eddy coef. to the background value
284         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
285         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
286         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
287      END DO
288      !
289   END SUBROUTINE zdf_ric_init
290
291#else
292   !!----------------------------------------------------------------------
293   !!   Dummy module :              NO Richardson dependent vertical mixing
294   !!----------------------------------------------------------------------
295   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .FALSE.   !: Richardson mixing flag
296CONTAINS
297   SUBROUTINE zdf_ric_init         ! Dummy routine
298   END SUBROUTINE zdf_ric_init
299   SUBROUTINE zdf_ric( kt )        ! Dummy routine
300      WRITE(*,*) 'zdf_ric: You should not have seen this print! error?', kt
301   END SUBROUTINE zdf_ric
302#endif
303
304   !!======================================================================
305END MODULE zdfric
Note: See TracBrowser for help on using the repository browser.