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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90 @ 4460

Last change on this file since 4460 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 12.1 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)  complet 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   !!----------------------------------------------------------------------
14#if defined key_zdfric   ||   defined key_esopa
15   !!----------------------------------------------------------------------
16   !!   'key_zdfric'                                             Kz = f(Ri)
17   !!----------------------------------------------------------------------
18   !!   zdf_ric      : update momentum and tracer Kz from the Richardson
19   !!                  number computation
20   !!   zdf_ric_init : initialization, namelist read, & parameters control
21   !!----------------------------------------------------------------------
22   USE oce             ! ocean dynamics and tracers variables
23   USE dom_oce         ! ocean space and time domain variables
24   USE zdf_oce         ! ocean vertical physics
25   USE in_out_manager  ! I/O manager
26   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
27   USE lib_mpp         ! MPP library
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   zdf_ric         ! called by step.F90
33   PUBLIC   zdf_ric_init    ! called by opa.F90
34
35   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .TRUE.   !: Richardson vertical mixing flag
36
37   !                                    !!* Namelist namzdf_ric : Richardson number dependent Kz *
38   INTEGER  ::   nn_ric   = 2            ! coefficient of the parameterization
39   REAL(wp) ::   rn_avmri = 100.e-4_wp   ! maximum value of the vertical eddy viscosity
40   REAL(wp) ::   rn_alp   =   5._wp      ! coefficient of the parameterization
41
42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tmric   !: coef. for the horizontal mean at t-point
43
44   !! * Control permutation of array indices
45#  include "oce_ftrans.h90"
46#  include "dom_oce_ftrans.h90"
47#  include "zdf_oce_ftrans.h90"
48!FTRANS tmric :I :I :z
49
50   !! * Substitutions
51#  include "domzgr_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
54   !! $Id$
55   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   INTEGER FUNCTION zdf_ric_alloc()
60      !!----------------------------------------------------------------------
61      !!                 ***  FUNCTION zdf_ric_alloc  ***
62      !!----------------------------------------------------------------------
63      ALLOCATE( tmric(jpi,jpj,jpk)   , STAT= zdf_ric_alloc )
64      !
65      IF( lk_mpp             )   CALL mpp_sum ( zdf_ric_alloc )
66      IF( zdf_ric_alloc /= 0 )   CALL ctl_warn('zdf_ric_alloc: failed to allocate arrays')
67   END FUNCTION zdf_ric_alloc
68
69
70   SUBROUTINE zdf_ric( kt )
71      !!----------------------------------------------------------------------
72      !!                 ***  ROUTINE zdfric  ***
73      !!                   
74      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
75      !!              a function of the local richardson number.
76      !!
77      !! ** Method  :   Local richardson number dependent formulation of the
78      !!              vertical eddy viscosity and diffusivity coefficients.
79      !!                The eddy coefficients are given by:
80      !!                    avm = avm0 + avmb
81      !!                    avt = avm0 / (1 + rn_alp*ri)
82      !!              with ri  = N^2 / dz(u)**2
83      !!                       = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]
84      !!                   avm0= rn_avmri / (1 + rn_alp*ri)**nn_ric
85      !!      Where ri is the before local Richardson number,
86      !!            rn_avmri is the maximum value reaches by avm and avt
87      !!            avmb and avtb are the background (or minimum) values
88      !!            and rn_alp, nn_ric are adjustable parameters.
89      !!      Typical values used are : avm0=1.e-2 m2/s, avmb=1.e-6 m2/s
90      !!      avtb=1.e-7 m2/s, rn_alp=5. and nn_ric=2.
91      !!      a numerical threshold is impose on the vertical shear (1.e-20)
92      !!        N.B. the mask are required for implicit scheme, and surface
93      !!      and bottom value already set in zdfini.F90
94      !!
95      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
96      !!----------------------------------------------------------------------
97      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
98      USE wrk_nemo, ONLY:   zwx => wrk_2d_1     ! 2D workspace
99      !!
100      INTEGER, INTENT( in ) ::   kt         ! ocean time-step indexocean time step
101      !!
102      INTEGER  ::   ji, jj, jk               ! dummy loop indices
103      REAL(wp) ::   zcoef, zdku, zdkv, zri, z05alp     ! temporary scalars
104      !!----------------------------------------------------------------------
105
106      IF( wrk_in_use(2, 1) ) THEN
107         CALL ctl_stop('zdf_ric : requested workspace array unavailable')   ;   RETURN
108      ENDIF
109
110      !! DCSE_NEMO: To optimise this loop for z_first indexing, make zwx 3-dimensional
111
112      !                                                ! ===============
113      DO jk = 2, jpkm1                                 ! Horizontal slab
114         !                                             ! ===============
115         ! Richardson number (put in zwx(ji,jj))
116         ! -----------------
117         DO jj = 2, jpjm1
118            DO ji = 2, jpim1
119               zcoef = 0.5 / fse3w(ji,jj,jk)
120               !                                            ! shear of horizontal velocity
121               zdku = zcoef * (  ub(ji-1,jj,jk-1) + ub(ji,jj,jk-1)   &
122                  &             -ub(ji-1,jj,jk  ) - ub(ji,jj,jk  )  )
123               zdkv = zcoef * (  vb(ji,jj-1,jk-1) + vb(ji,jj,jk-1)   &
124                  &             -vb(ji,jj-1,jk  ) - vb(ji,jj,jk  )  )
125               !                                            ! richardson number (minimum value set to zero)
126               zri = rn2(ji,jj,jk) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
127               zwx(ji,jj) = MAX( zri, 0.e0 )
128            END DO
129         END DO
130         CALL lbc_lnk( zwx, 'W', 1. )                       ! Boundary condition   (sign unchanged)
131
132         ! Vertical eddy viscosity and diffusivity coefficients
133         ! -------------------------------------------------------
134         z05alp = 0.5_wp * rn_alp
135         DO jj = 1, jpjm1                                   ! Eddy viscosity coefficients (avm)
136            DO ji = 1, jpim1
137               avmu(ji,jj,jk) = umask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji+1,jj)+zwx(ji,jj) ) )**nn_ric
138               avmv(ji,jj,jk) = vmask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji,jj+1)+zwx(ji,jj) ) )**nn_ric
139            END DO
140         END DO
141         DO jj = 2, jpjm1                                   ! Eddy diffusivity coefficients (avt)
142            DO ji = 2, jpim1
143               avt(ji,jj,jk) = tmric(ji,jj,jk) / ( 1._wp + rn_alp * zwx(ji,jj) )           &
144                  &                            * (  avmu(ji,jj,jk) + avmu(ji-1,jj,jk)      &
145                  &                               + avmv(ji,jj,jk) + avmv(ji,jj-1,jk)  )   &
146                  &          + avtb(jk) * tmask(ji,jj,jk)
147               !                                            ! Add the background coefficient on eddy viscosity
148               avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk)
149               avmv(ji,jj,jk) = avmv(ji,jj,jk) + avmb(jk) * vmask(ji,jj,jk)
150            END DO
151         END DO
152         !                                             ! ===============
153      END DO                                           !   End of slab
154      !                                                ! ===============
155      !
156      CALL lbc_lnk( avt , 'W', 1. )                         ! Boundary conditions   (unchanged sign)
157      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
158      !
159      IF( wrk_not_released(2, 1) )   CALL ctl_stop('zdf_ric: failed to release workspace array')
160      !
161   END SUBROUTINE zdf_ric
162
163
164   SUBROUTINE zdf_ric_init
165      !!----------------------------------------------------------------------
166      !!                 ***  ROUTINE zdfbfr_init  ***
167      !!                   
168      !! ** Purpose :   Initialization of the vertical eddy diffusivity and
169      !!      viscosity coef. for the Richardson number dependent formulation.
170      !!
171      !! ** Method  :   Read the namzdf_ric namelist and check the parameter values
172      !!
173      !! ** input   :   Namelist namzdf_ric
174      !!
175      !! ** Action  :   increase by 1 the nstop flag is setting problem encounter
176      !!----------------------------------------------------------------------
177      INTEGER :: ji, jj, jk   ! dummy loop indices
178      !!
179      NAMELIST/namzdf_ric/ rn_avmri, rn_alp, nn_ric
180      !!----------------------------------------------------------------------
181      !
182      REWIND( numnam )               ! Read Namelist namzdf_ric : richardson number dependent Kz
183      READ  ( numnam, namzdf_ric )
184      !
185      IF(lwp) THEN                   ! Control print
186         WRITE(numout,*)
187         WRITE(numout,*) 'zdf_ric : Ri depend vertical mixing scheme'
188         WRITE(numout,*) '~~~~~~~'
189         WRITE(numout,*) '   Namelist namzdf_ric : set Kz(Ri) parameters'
190         WRITE(numout,*) '      maximum vertical viscosity     rn_avmri = ', rn_avmri
191         WRITE(numout,*) '      coefficient                    rn_alp   = ', rn_alp
192         WRITE(numout,*) '      coefficient                    nn_ric   = ', nn_ric
193      ENDIF
194      !
195      !                              ! allocate zdfric arrays
196      IF( zdf_ric_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ric_init : unable to allocate arrays' )
197      !
198#if defined key_z_first
199      DO jj = 2, jpj 
200         DO ji = 2, jpi
201            DO jk = 1, jpk
202#else
203      DO jk = 1, jpk                 ! weighting mean array tmric for 4 T-points
204         DO jj = 2, jpj              ! which accounts for coastal boundary conditions           
205            DO ji = 2, jpi
206#endif
207               tmric(ji,jj,jk) =  tmask(ji,jj,jk)                                  &
208                  &            / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
209                  &                      + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
210            END DO
211         END DO
212      END DO
213      tmric(:,1,:) = 0._wp
214      !
215#if defined key_z_first
216      DO jj = 1, jpj
217         DO ji = 1, jpi
218            DO jk = 1, jpk           ! Initialization of vertical eddy coef. to the background value
219               avt (ji,jj,jk) = avtb(jk) * tmask(ji,jj,jk)
220               avmu(ji,jj,jk) = avmb(jk) * umask(ji,jj,jk)
221               avmv(ji,jj,jk) = avmb(jk) * vmask(ji,jj,jk)
222            END DO
223         END DO
224      END DO
225#else
226      DO jk = 1, jpk                 ! Initialization of vertical eddy coef. to the background value
227         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
228         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
229         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
230      END DO
231#endif
232      !
233   END SUBROUTINE zdf_ric_init
234
235#else
236   !!----------------------------------------------------------------------
237   !!   Dummy module :              NO Richardson dependent vertical mixing
238   !!----------------------------------------------------------------------
239   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .FALSE.   !: Richardson mixing flag
240CONTAINS
241   SUBROUTINE zdf_ric_init         ! Dummy routine
242   END SUBROUTINE zdf_ric_init
243   SUBROUTINE zdf_ric( kt )        ! Dummy routine
244      WRITE(*,*) 'zdf_ric: You should not have seen this print! error?', kt
245   END SUBROUTINE zdf_ric
246#endif
247
248   !!======================================================================
249END MODULE zdfric
Note: See TracBrowser for help on using the repository browser.