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

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago

#2050 first version. Compiled OK in moci test suite

File size: 16.7 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)
[3294]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
[1601]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
[3294]13   !!            3.3.1! 2011-09  (P. Oddo) Mixed layer depth parameterization
[1601]14   !!----------------------------------------------------------------------
[3]15#if defined key_zdfric   ||   defined key_esopa
16   !!----------------------------------------------------------------------
17   !!   'key_zdfric'                                             Kz = f(Ri)
18   !!----------------------------------------------------------------------
[3625]19   !!   zdf_ric       : update momentum and tracer Kz from the Richardson
[3]20   !!                  number computation
[3625]21   !!   zdf_ric_init  : initialization, namelist read, & parameters control
[3]22   !!----------------------------------------------------------------------
[3625]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 timing         ! Timing
31   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[3]32
[3294]33   USE eosbn2, ONLY : nn_eos
34
[3]35   IMPLICIT NONE
36   PRIVATE
37
[2528]38   PUBLIC   zdf_ric         ! called by step.F90
39   PUBLIC   zdf_ric_init    ! called by opa.F90
[9366]40   PRIVATE  ric_namelist
[3]41
[1537]42   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .TRUE.   !: Richardson vertical mixing flag
[3]43
[4147]44   !                           !!* Namelist namzdf_ric : Richardson number dependent Kz *
45   INTEGER  ::   nn_ric         ! coefficient of the parameterization
46   REAL(wp) ::   rn_avmri       ! maximum value of the vertical eddy viscosity
47   REAL(wp) ::   rn_alp         ! coefficient of the parameterization
48   REAL(wp) ::   rn_ekmfc       ! Ekman Factor Coeff
49   REAL(wp) ::   rn_mldmin      ! minimum mixed layer (ML) depth   
50   REAL(wp) ::   rn_mldmax      ! maximum mixed layer depth
51   REAL(wp) ::   rn_wtmix       ! Vertical eddy Diff. in the ML
52   REAL(wp) ::   rn_wvmix       ! Vertical eddy Visc. in the ML
53   LOGICAL  ::   ln_mldw        ! Use or not the MLD parameters
[3]54
[2715]55   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tmric   !: coef. for the horizontal mean at t-point
[1537]56
[3]57   !! * Substitutions
58#  include "domzgr_substitute.h90"
59   !!----------------------------------------------------------------------
[2715]60   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
[2528]61   !! $Id$
[2715]62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[247]63   !!----------------------------------------------------------------------
[3]64CONTAINS
65
[2715]66   INTEGER FUNCTION zdf_ric_alloc()
67      !!----------------------------------------------------------------------
68      !!                 ***  FUNCTION zdf_ric_alloc  ***
69      !!----------------------------------------------------------------------
70      ALLOCATE( tmric(jpi,jpj,jpk)   , STAT= zdf_ric_alloc )
71      !
72      IF( lk_mpp             )   CALL mpp_sum ( zdf_ric_alloc )
73      IF( zdf_ric_alloc /= 0 )   CALL ctl_warn('zdf_ric_alloc: failed to allocate arrays')
74   END FUNCTION zdf_ric_alloc
75
76
[3]77   SUBROUTINE zdf_ric( kt )
78      !!----------------------------------------------------------------------
79      !!                 ***  ROUTINE zdfric  ***
80      !!                   
81      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
[3294]82      !!                a function of the local richardson number.
[3]83      !!
84      !! ** Method  :   Local richardson number dependent formulation of the
[3294]85      !!                vertical eddy viscosity and diffusivity coefficients.
[1601]86      !!                The eddy coefficients are given by:
87      !!                    avm = avm0 + avmb
88      !!                    avt = avm0 / (1 + rn_alp*ri)
[3294]89      !!                with ri  = N^2 / dz(u)**2
90      !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]
91      !!                    avm0= rn_avmri / (1 + rn_alp*ri)**nn_ric
[1601]92      !!      Where ri is the before local Richardson number,
93      !!            rn_avmri is the maximum value reaches by avm and avt
94      !!            avmb and avtb are the background (or minimum) values
95      !!            and rn_alp, nn_ric are adjustable parameters.
96      !!      Typical values used are : avm0=1.e-2 m2/s, avmb=1.e-6 m2/s
[1537]97      !!      avtb=1.e-7 m2/s, rn_alp=5. and nn_ric=2.
[3]98      !!      a numerical threshold is impose on the vertical shear (1.e-20)
[3294]99      !!      As second step compute Ekman depth from wind stress forcing
100      !!      and apply namelist provided vertical coeff within this depth.
101      !!      The Ekman depth is:
102      !!              Ustar = SQRT(Taum/rho0)
103      !!              ekd= rn_ekmfc * Ustar / f0
104      !!      Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation
105      !!      of the above equation indicates the value is somewhat arbitrary; therefore
106      !!      we allow the freedom to increase or decrease this value, if the
107      !!      Ekman depth estimate appears too shallow or too deep, respectively.
108      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the
109      !!      namelist
[3]110      !!        N.B. the mask are required for implicit scheme, and surface
[1601]111      !!      and bottom value already set in zdfini.F90
[3]112      !!
[1601]113      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
[3294]114      !!              PFJ Lermusiaux 2001.
[3]115      !!----------------------------------------------------------------------
[3294]116      USE phycst,   ONLY:   rsmall,rau0
117      USE sbc_oce,  ONLY:   taum
[2715]118      !!
[3294]119      INTEGER, INTENT( in ) ::   kt                           ! ocean time-step
[1537]120      !!
[3294]121      INTEGER  ::   ji, jj, jk                                ! dummy loop indices
122      REAL(wp) ::   zcoef, zdku, zdkv, zri, z05alp, zflageos  ! temporary scalars
123      REAL(wp) ::   zrhos, zustar
124      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, ekm_dep 
[1537]125      !!----------------------------------------------------------------------
[3294]126      !
127      IF( nn_timing == 1 )  CALL timing_start('zdf_ric')
128      !
129      CALL wrk_alloc( jpi,jpj, zwx, ekm_dep )
[3]130      !                                                ! ===============
131      DO jk = 2, jpkm1                                 ! Horizontal slab
132         !                                             ! ===============
133         ! Richardson number (put in zwx(ji,jj))
134         ! -----------------
135         DO jj = 2, jpjm1
136            DO ji = 2, jpim1
137               zcoef = 0.5 / fse3w(ji,jj,jk)
[1601]138               !                                            ! shear of horizontal velocity
[3]139               zdku = zcoef * (  ub(ji-1,jj,jk-1) + ub(ji,jj,jk-1)   &
[1537]140                  &             -ub(ji-1,jj,jk  ) - ub(ji,jj,jk  )  )
[3]141               zdkv = zcoef * (  vb(ji,jj-1,jk-1) + vb(ji,jj,jk-1)   &
[1537]142                  &             -vb(ji,jj-1,jk  ) - vb(ji,jj,jk  )  )
[1601]143               !                                            ! richardson number (minimum value set to zero)
[3]144               zri = rn2(ji,jj,jk) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
145               zwx(ji,jj) = MAX( zri, 0.e0 )
146            END DO
147         END DO
[1601]148         CALL lbc_lnk( zwx, 'W', 1. )                       ! Boundary condition   (sign unchanged)
[3]149
150         ! Vertical eddy viscosity and diffusivity coefficients
151         ! -------------------------------------------------------
[2715]152         z05alp = 0.5_wp * rn_alp
[1601]153         DO jj = 1, jpjm1                                   ! Eddy viscosity coefficients (avm)
[3]154            DO ji = 1, jpim1
[2715]155               avmu(ji,jj,jk) = umask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji+1,jj)+zwx(ji,jj) ) )**nn_ric
156               avmv(ji,jj,jk) = vmask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji,jj+1)+zwx(ji,jj) ) )**nn_ric
[3]157            END DO
158         END DO
[1601]159         DO jj = 2, jpjm1                                   ! Eddy diffusivity coefficients (avt)
[3]160            DO ji = 2, jpim1
[2715]161               avt(ji,jj,jk) = tmric(ji,jj,jk) / ( 1._wp + rn_alp * zwx(ji,jj) )           &
162                  &                            * (  avmu(ji,jj,jk) + avmu(ji-1,jj,jk)      &
163                  &                               + avmv(ji,jj,jk) + avmv(ji,jj-1,jk)  )   &
[1537]164                  &          + avtb(jk) * tmask(ji,jj,jk)
[1601]165               !                                            ! Add the background coefficient on eddy viscosity
[3]166               avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk)
167               avmv(ji,jj,jk) = avmv(ji,jj,jk) + avmb(jk) * vmask(ji,jj,jk)
168            END DO
169         END DO
170         !                                             ! ===============
171      END DO                                           !   End of slab
172      !                                                ! ===============
[1601]173      !
[3294]174      IF( ln_mldw ) THEN
175
176      !  Compute Ekman depth from wind stress forcing.
177      ! -------------------------------------------------------
178      zflageos = ( 0.5 + SIGN( 0.5, nn_eos - 1. ) ) * rau0
179      DO jj = 1, jpj
180         DO ji = 1, jpi
181            zrhos          = rhop(ji,jj,1) + zflageos * ( 1. - tmask(ji,jj,1) )
182            zustar         = SQRT( taum(ji,jj) / ( zrhos +  rsmall ) )
183            ekm_dep(ji,jj) = rn_ekmfc * zustar / ( ABS( ff(ji,jj) ) + rsmall )
184            ekm_dep(ji,jj) = MAX(ekm_dep(ji,jj),rn_mldmin) ! Minimun allowed
185            ekm_dep(ji,jj) = MIN(ekm_dep(ji,jj),rn_mldmax) ! Maximum allowed
186         END DO
187      END DO
188
189      ! In the first model level vertical diff/visc coeff.s
190      ! are always equal to the namelist values rn_wtmix/rn_wvmix
191      ! -------------------------------------------------------
192      DO jj = 1, jpj
193         DO ji = 1, jpi
194            avmv(ji,jj,1) = MAX( avmv(ji,jj,1), rn_wvmix )
195            avmu(ji,jj,1) = MAX( avmu(ji,jj,1), rn_wvmix )
196            avt( ji,jj,1) = MAX(  avt(ji,jj,1), rn_wtmix )
197         END DO
198      END DO
199
200      !  Force the vertical mixing coef within the Ekman depth
201      ! -------------------------------------------------------
202      DO jk = 2, jpkm1
203         DO jj = 1, jpj
204            DO ji = 1, jpi
205               IF( fsdept(ji,jj,jk) < ekm_dep(ji,jj) ) THEN
206                  avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), rn_wvmix )
207                  avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), rn_wvmix )
208                  avt( ji,jj,jk) = MAX(  avt(ji,jj,jk), rn_wtmix )
209               ENDIF
210            END DO
211         END DO
212      END DO
213
214      DO jk = 1, jpkm1               
215         DO jj = 1, jpj
216            DO ji = 1, jpi
217               avmv(ji,jj,jk) = avmv(ji,jj,jk) * vmask(ji,jj,jk)
218               avmu(ji,jj,jk) = avmu(ji,jj,jk) * umask(ji,jj,jk)
219               avt( ji,jj,jk) = avt( ji,jj,jk) * tmask(ji,jj,jk)
220            END DO
221         END DO
222      END DO
223
224     ENDIF
225
[1601]226      CALL lbc_lnk( avt , 'W', 1. )                         ! Boundary conditions   (unchanged sign)
227      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
228      !
[3294]229      CALL wrk_dealloc( jpi,jpj, zwx, ekm_dep )
[2715]230      !
[3294]231      IF( nn_timing == 1 )  CALL timing_stop('zdf_ric')
232      !
[3]233   END SUBROUTINE zdf_ric
234
235
236   SUBROUTINE zdf_ric_init
237      !!----------------------------------------------------------------------
238      !!                 ***  ROUTINE zdfbfr_init  ***
239      !!                   
240      !! ** Purpose :   Initialization of the vertical eddy diffusivity and
241      !!      viscosity coef. for the Richardson number dependent formulation.
242      !!
[1601]243      !! ** Method  :   Read the namzdf_ric namelist and check the parameter values
[3]244      !!
[1601]245      !! ** input   :   Namelist namzdf_ric
[3]246      !!
247      !! ** Action  :   increase by 1 the nstop flag is setting problem encounter
248      !!----------------------------------------------------------------------
[2715]249      INTEGER :: ji, jj, jk   ! dummy loop indices
[4147]250      INTEGER ::   ios        ! Local integer output status for namelist read
[1537]251      !!
[3294]252      NAMELIST/namzdf_ric/ rn_avmri, rn_alp   , nn_ric  , rn_ekmfc,  &
253         &                rn_mldmin, rn_mldmax, rn_wtmix, rn_wvmix, ln_mldw
[3]254      !!----------------------------------------------------------------------
[1601]255      !
[9366]256      IF(lwm) THEN
257         REWIND( numnam_ref )              ! Namelist namzdf_ric in reference namelist : Vertical diffusion Kz depends on Richardson number
258         READ  ( numnam_ref, namzdf_ric, IOSTAT = ios, ERR = 901)
259901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in reference namelist', lwm )
260         REWIND( numnam_cfg )              ! Namelist namzdf_ric in configuration namelist : Vertical diffusion Kz depends on Richardson number
261         READ  ( numnam_cfg, namzdf_ric, IOSTAT = ios, ERR = 902 )
262902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in configuration namelist', lwm )
263      ENDIF
[4147]264
[4624]265      IF(lwm) WRITE ( numond, namzdf_ric )
[1601]266      !
[9366]267      CALL ric_namelist()
268
[1601]269      IF(lwp) THEN                   ! Control print
[3]270         WRITE(numout,*)
[1537]271         WRITE(numout,*) 'zdf_ric : Ri depend vertical mixing scheme'
272         WRITE(numout,*) '~~~~~~~'
[1601]273         WRITE(numout,*) '   Namelist namzdf_ric : set Kz(Ri) parameters'
[3294]274         WRITE(numout,*) '      maximum vertical viscosity     rn_avmri  = ', rn_avmri
275         WRITE(numout,*) '      coefficient                    rn_alp    = ', rn_alp
276         WRITE(numout,*) '      coefficient                    nn_ric    = ', nn_ric
277         WRITE(numout,*) '      Ekman Factor Coeff             rn_ekmfc  = ', rn_ekmfc
278         WRITE(numout,*) '      minimum mixed layer depth      rn_mldmin = ', rn_mldmin
279         WRITE(numout,*) '      maximum mixed layer depth      rn_mldmax = ', rn_mldmax
280         WRITE(numout,*) '      Vertical eddy Diff. in the ML  rn_wtmix  = ', rn_wtmix
281         WRITE(numout,*) '      Vertical eddy Visc. in the ML  rn_wvmix  = ', rn_wvmix
282         WRITE(numout,*) '      Use the MLD parameterization   ln_mldw   = ', ln_mldw
[3]283      ENDIF
[1601]284      !
[2715]285      !                              ! allocate zdfric arrays
286      IF( zdf_ric_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ric_init : unable to allocate arrays' )
287      !
288      DO jk = 1, jpk                 ! weighting mean array tmric for 4 T-points
289         DO jj = 2, jpj              ! which accounts for coastal boundary conditions           
[3]290            DO ji = 2, jpi
291               tmric(ji,jj,jk) =  tmask(ji,jj,jk)                                  &
[1601]292                  &            / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
293                  &                      + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
[3]294            END DO
295         END DO
296      END DO
[2715]297      tmric(:,1,:) = 0._wp
[1601]298      !
299      DO jk = 1, jpk                 ! Initialization of vertical eddy coef. to the background value
[3]300         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
301         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
302         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
303      END DO
[1601]304      !
[3]305   END SUBROUTINE zdf_ric_init
306
[9366]307   SUBROUTINE ric_namelist()
308     !!---------------------------------------------------------------------
309     !!                   ***  ROUTINE ric_namelist  ***
310     !!                     
311     !! ** Purpose :   Broadcast namelist variables read by procesor lwm
312     !!
313     !! ** Method  :   use lib_mpp
314     !!----------------------------------------------------------------------
315#if defined key_mpp_mpi
316     CALL mpp_bcast(rn_avmri)
317     CALL mpp_bcast(rn_alp)
318     CALL mpp_bcast(nn_ric)
319     CALL mpp_bcast(rn_ekmfc)
320     CALL mpp_bcast(rn_mldmin)
321     CALL mpp_bcast(rn_mldmax)
322     CALL mpp_bcast(rn_wtmix)
323     CALL mpp_bcast(rn_wvmix)
324     CALL mpp_bcast(ln_mldw)
325#endif
326    END SUBROUTINE ric_namelist
[3]327#else
328   !!----------------------------------------------------------------------
329   !!   Dummy module :              NO Richardson dependent vertical mixing
330   !!----------------------------------------------------------------------
331   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .FALSE.   !: Richardson mixing flag
332CONTAINS
[2528]333   SUBROUTINE zdf_ric_init         ! Dummy routine
334   END SUBROUTINE zdf_ric_init
[3]335   SUBROUTINE zdf_ric( kt )        ! Dummy routine
[16]336      WRITE(*,*) 'zdf_ric: You should not have seen this print! error?', kt
[3]337   END SUBROUTINE zdf_ric
338#endif
339
340   !!======================================================================
341END MODULE zdfric
Note: See TracBrowser for help on using the repository browser.