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

source: trunk/NEMO/OPA_SRC/ZDF/zdfric.F90 @ 1537

Last change on this file since 1537 was 1537, checked in by ctlod, 15 years ago

ensure the restartability of the 2nd order advection scheme,see ticket: 489

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 10.3 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#if defined key_zdfric   ||   defined key_esopa
8   !!----------------------------------------------------------------------
9   !!   'key_zdfric'                                             Kz = f(Ri)
10   !!----------------------------------------------------------------------
11   !!   zdf_ric      : update momentum and tracer Kz from the Richardson
12   !!                  number computation
13   !!   zdf_ric_init : initialization, namelist read, & parameters control
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers variables
16   USE dom_oce         ! ocean space and time domain variables
17   USE zdf_oce         ! ocean vertical physics
18   USE in_out_manager  ! I/O manager
19   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
20
21   IMPLICIT NONE
22   PRIVATE
23
24   PUBLIC   zdf_ric    ! called by step.F90
25
26   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .TRUE.   !: Richardson vertical mixing flag
27
28   !                                    !!* Namelist nam_ric : Richardson number dependent Kz *
29   INTEGER  ::   nn_ric  = 2             ! coefficient of the parameterization
30   REAL(wp) ::   rn_avmri = 100.e-4_wp   ! maximum value of the vertical eddy viscosity
31   REAL(wp) ::   rn_alp   =   5._wp      ! coefficient of the parameterization
32
33   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   tmric                    ! coef. for the horizontal mean at t-point
34
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
39   !! $Id$
40   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE zdf_ric( kt )
45      !!----------------------------------------------------------------------
46      !!                 ***  ROUTINE zdfric  ***
47      !!                   
48      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
49      !!      a function of the local richardson number.
50      !!
51      !! ** Method  :   Local richardson number dependent formulation of the
52      !!      vertical eddy viscosity and diffusivity coefficients. the eddy
53      !!      coefficients are given by:
54      !!              avm = avm0 + avmb
55      !!              avt = avm0 / (1 + rn_alp*ri)
56      !!      with    ri  = N^2 / dz(u)**2
57      !!                  = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]
58      !!              avm0= rn_avmri / (1 + rn_alp*ri)**nn_ric
59      !!      Where ri is the before local Richardson number, rn_avmri the maximum
60      !!      value reaches by the vertical eddy coefficients, avmb and avtb
61      !!      the background (or minimum) values of these coefficients for
62      !!      momemtum and tracers, and rn_alp, nn_ric are adjustable parameters.
63      !!      typical values used are : avm0=1.e-2 m2/s, avmb=1.e-6 m2/s
64      !!      avtb=1.e-7 m2/s, rn_alp=5. and nn_ric=2.
65      !!      this formulation needs ri>=0 : ri is set to zero if dz(rau)<0.
66      !!      a numerical threshold is impose on the vertical shear (1.e-20)
67      !!        N.B. the mask are required for implicit scheme, and surface
68      !!      and bottom value already set in inimix.F
69      !!
70      !! References :
71      !!      pacanowski & philander 1981, j. phys. oceanogr., 1441-1451.
72      !! History :
73      !!        !  87-09  (P. Andrich)  Original code
74      !!        !  91-11  (G. Madec)
75      !!        !  93-03  (M. Guyon)  symetrical conditions
76      !!        !  96-01  (G. Madec)  complet rewriting of multitasking
77      !!                                  suppression of common work arrays
78      !!        !  97-06 (G. Madec)  complete rewriting of zdfmix
79      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
80      !!----------------------------------------------------------------------
81      INTEGER, INTENT( in ) ::   kt         ! ocean time-step indexocean time step
82      !!
83      INTEGER  ::   ji, jj, jk               ! dummy loop indices
84      REAL(wp) ::   zcoef, zdku, zdkv, zri, z05alp     ! temporary scalars
85      REAL(wp), DIMENSION(jpi,jpj) ::   zwx ! temporary workspace
86      !!----------------------------------------------------------------------
87
88      IF( kt == nit000  ) CALL zdf_ric_init            ! Initialization (first time-step only)
89
90      !                                                ! ===============
91      DO jk = 2, jpkm1                                 ! Horizontal slab
92         !                                             ! ===============
93         ! Richardson number (put in zwx(ji,jj))
94         ! -----------------
95         ! minimum value set to zero
96         DO jj = 2, jpjm1
97            DO ji = 2, jpim1
98               zcoef = 0.5 / fse3w(ji,jj,jk)
99               ! shear of horizontal velocity
100               zdku = zcoef * (  ub(ji-1,jj,jk-1) + ub(ji,jj,jk-1)   &
101                  &             -ub(ji-1,jj,jk  ) - ub(ji,jj,jk  )  )
102               zdkv = zcoef * (  vb(ji,jj-1,jk-1) + vb(ji,jj,jk-1)   &
103                  &             -vb(ji,jj-1,jk  ) - vb(ji,jj,jk  )  )
104               ! richardson number (minimum value set to zero)
105               zri = rn2(ji,jj,jk) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
106               zwx(ji,jj) = MAX( zri, 0.e0 )
107            END DO
108         END DO
109
110         ! Boundary condition on zwx   (sign unchanged)
111         CALL lbc_lnk( zwx, 'W', 1. )
112
113
114         ! Vertical eddy viscosity and diffusivity coefficients
115         ! -------------------------------------------------------
116         ! Eddy viscosity coefficients
117         z05alp = 0.5 * rn_alp
118         DO jj = 1, jpjm1
119            DO ji = 1, jpim1
120               avmu(ji,jj,jk) = umask(ji,jj,jk)   &
121                  &           * rn_avmri / ( 1. + z05alp*( zwx(ji+1,jj)+zwx(ji,jj) ) )**nn_ric
122               avmv(ji,jj,jk) = vmask(ji,jj,jk)   &
123                  &           * rn_avmri / ( 1. + z05alp*( zwx(ji,jj+1)+zwx(ji,jj) ) )**nn_ric
124            END DO
125         END DO
126
127         ! Eddy diffusivity coefficients
128         DO jj = 2, jpjm1
129            DO ji = 2, jpim1
130               avt(ji,jj,jk) = tmric(ji,jj,jk) / ( 1. + rn_alp * zwx(ji,jj) )   &
131                  &          * (  avmu(ji,jj,jk) + avmu(ji-1, jj ,jk)        &
132                  &             + avmv(ji,jj,jk) + avmv( ji ,jj-1,jk)  )     &
133                  &          + avtb(jk) * tmask(ji,jj,jk)
134            END DO
135         END DO
136
137         ! Add the background coefficient on eddy viscosity
138         DO jj = 2, jpjm1
139            DO ji = 2, jpim1
140               avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk)
141               avmv(ji,jj,jk) = avmv(ji,jj,jk) + avmb(jk) * vmask(ji,jj,jk)
142            END DO
143         END DO
144         !                                             ! ===============
145      END DO                                           !   End of slab
146      !                                                ! ===============
147
148      ! Boundary conditions on (avt,avmu,avmv)   (unchanged sign)
149      ! -----------------------===============
150      CALL lbc_lnk( avt , 'W', 1. )
151      CALL lbc_lnk( avmu, 'U', 1. )
152      CALL lbc_lnk( avmv, 'V', 1. )
153
154   END SUBROUTINE zdf_ric
155
156
157   SUBROUTINE zdf_ric_init
158      !!----------------------------------------------------------------------
159      !!                 ***  ROUTINE zdfbfr_init  ***
160      !!                   
161      !! ** Purpose :   Initialization of the vertical eddy diffusivity and
162      !!      viscosity coef. for the Richardson number dependent formulation.
163      !!
164      !! ** Method  :   Read the namric namelist and check the parameter values
165      !!
166      !! ** input   :   Namelist namric
167      !!
168      !! ** Action  :   increase by 1 the nstop flag is setting problem encounter
169      !!
170      !! history :
171      !!  8.5  !  02-06  (G. Madec)  original code
172      !!----------------------------------------------------------------------
173      INTEGER :: ji, jj, jk        ! dummy loop indices
174      !!
175      NAMELIST/nam_ric/ rn_avmri, rn_alp, nn_ric
176      !!----------------------------------------------------------------------
177
178      REWIND ( numnam )               ! Read Namelist nam_ric : richardson number dependent Kz
179      READ   ( numnam, nam_ric )
180
181      IF(lwp) THEN                    ! Control print
182         WRITE(numout,*)
183         WRITE(numout,*) 'zdf_ric : Ri depend vertical mixing scheme'
184         WRITE(numout,*) '~~~~~~~'
185         WRITE(numout,*) '   Namelist namric : set Kz(Ri) parameters'
186         WRITE(numout,*) '      maximum vertical viscosity     rn_avmri = ', rn_avmri
187         WRITE(numout,*) '      coefficient                    rn_alp   = ', rn_alp
188         WRITE(numout,*) '      coefficient                    nn_ric   = ', nn_ric
189      ENDIF
190
191      ! weighting mean array tmric for 4 T-points which accounts for coastal boundary conditions.
192      DO jk = 1, jpk
193         DO jj = 2, jpj
194            DO ji = 2, jpi
195               tmric(ji,jj,jk) =  tmask(ji,jj,jk)                                  &
196                               / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
197                                         + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
198            END DO
199         END DO
200      END DO
201      tmric(:,1,:) = 0.e0
202
203      ! Initialization of vertical eddy coef. to the background value
204      DO jk = 1, jpk
205         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
206         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
207         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
208      END DO
209
210   END SUBROUTINE zdf_ric_init
211
212#else
213   !!----------------------------------------------------------------------
214   !!   Dummy module :              NO Richardson dependent vertical mixing
215   !!----------------------------------------------------------------------
216   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .FALSE.   !: Richardson mixing flag
217CONTAINS
218   SUBROUTINE zdf_ric( kt )        ! Dummy routine
219      WRITE(*,*) 'zdf_ric: You should not have seen this print! error?', kt
220   END SUBROUTINE zdf_ric
221#endif
222
223   !!======================================================================
224END MODULE zdfric
Note: See TracBrowser for help on using the repository browser.