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

Last change on this file since 3042 was 2715, checked in by rblod, 13 years ago

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