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.
ldfdyn.F90 in NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90 @ 12688

Last change on this file since 12688 was 12688, checked in by davestorkey, 4 years ago

UKMO/NEMO_4.0.1_biharmonic_GM : Switch to bhm coefficient at W points.

File size: 41.7 KB
RevLine 
[3]1MODULE ldfdyn
2   !!======================================================================
3   !!                       ***  MODULE  ldfdyn  ***
4   !! Ocean physics:  lateral viscosity coefficient
5   !!=====================================================================
[1601]6   !! History :  OPA  ! 1997-07  (G. Madec)  multi dimensional coefficients
7   !!   NEMO     1.0  ! 2002-09  (G. Madec)  F90: Free form and module
[5836]8   !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
9   !!                 !                                  add velocity dependent coefficient and optional read in file
[12535]10   !!            4.0  ! 2020-02  (C. Wilson, ...) add bhm coefficient for bi-Laplacian GM implementation via momentum     
[1601]11   !!----------------------------------------------------------------------
[3]12
13   !!----------------------------------------------------------------------
[5836]14   !!   ldf_dyn_init  : initialization, namelist read, and parameters control
15   !!   ldf_dyn       : update lateral eddy viscosity coefficients at each time step
[3]16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers   
[12535]18   USE dom_oce         ! ocean space and time domain
[3]19   USE phycst          ! physical constants
[9490]20   USE ldfslp          ! lateral diffusion: slopes of mixing orientation
[5836]21   USE ldfc1d_c2d      ! lateral diffusion: 1D and 2D cases
22   !
[3]23   USE in_out_manager  ! I/O manager
[5836]24   USE iom             ! I/O module for ehanced bottom friction file
25   USE timing          ! Timing
[3]26   USE lib_mpp         ! distribued memory computing library
27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
28
29   IMPLICIT NONE
30   PRIVATE
31
[5836]32   PUBLIC   ldf_dyn_init   ! called by nemogcm.F90
33   PUBLIC   ldf_dyn        ! called by step.F90
[3]34
[9490]35   !                                    !!* Namelist namdyn_ldf : lateral mixing on momentum *
[9526]36   LOGICAL , PUBLIC ::   ln_dynldf_OFF   !: No operator (i.e. no explicit diffusion)
[5836]37   LOGICAL , PUBLIC ::   ln_dynldf_lap   !: laplacian operator
[12535]38   LOGICAL , PUBLIC ::   ln_dynldf_bgm   !: bi-Laplacian GM implementation via momentum
[5836]39   LOGICAL , PUBLIC ::   ln_dynldf_blp   !: bilaplacian operator
40   LOGICAL , PUBLIC ::   ln_dynldf_lev   !: iso-level direction
41   LOGICAL , PUBLIC ::   ln_dynldf_hor   !: horizontal (geopotential) direction
[9490]42!  LOGICAL , PUBLIC ::   ln_dynldf_iso   !: iso-neutral direction                        (see ldfslp)
[5836]43   INTEGER , PUBLIC ::   nn_ahm_ijk_t    !: choice of time & space variations of the lateral eddy viscosity coef.
[9490]44   !                                        !  time invariant coefficients:  aht = 1/2  Ud*Ld   (lap case)
45   !                                           !                             bht = 1/12 Ud*Ld^3 (blp case)
[12535]46   INTEGER , PUBLIC ::   nn_bhm_ijk_t    !: choice of time & space variations of the lateral bi-Laplacian GM coef.
[9490]47   REAL(wp), PUBLIC ::   rn_Uv                 !: lateral viscous velocity  [m/s]
48   REAL(wp), PUBLIC ::   rn_Lv                 !: lateral viscous length    [m]
49   !                                        ! Smagorinsky viscosity  (nn_ahm_ijk_t = 32)
50   REAL(wp), PUBLIC ::   rn_csmc               !: Smagorinsky constant of proportionality
51   REAL(wp), PUBLIC ::   rn_minfac             !: Multiplicative factor of theorectical minimum Smagorinsky viscosity
52   REAL(wp), PUBLIC ::   rn_maxfac             !: Multiplicative factor of theorectical maximum Smagorinsky viscosity
53   !                                        ! iso-neutral laplacian (ln_dynldf_lap=ln_dynldf_iso=T)
54   REAL(wp), PUBLIC ::   rn_ahm_b              !: lateral laplacian background eddy viscosity  [m2/s]
[12535]55   REAL(wp), PUBLIC ::   rn_bhm_b              !: lateral bi-Laplacian GM background eddy viscosity  [m4/s]
56   REAL(wp), PUBLIC ::   rn_bgmzcrit           !: critical depth (>=0) for linear tapering of bi-Lap. GM coef. to zero at surface [m]
[9490]57   !                                    !!* Parameter to control the type of lateral viscous operator
58   INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10                       !: error in setting the operator
59   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00                       !: without operator (i.e. no lateral viscous trend)
60   !                          !!      laplacian     !    bilaplacian    !
61   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  !: iso-level operator
62   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       !: iso-neutral or geopotential operator
[12535]63   INTEGER, PARAMETER, PUBLIC ::   np_bgm    = 30                       !: iso-level operator, bi-Laplacian GM
[9490]64   INTEGER           , PUBLIC ::   nldf_dyn         !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals)
65   LOGICAL           , PUBLIC ::   l_ldfdyn_time    !: flag for time variation of the lateral eddy viscosity coef.
[5836]66
[9490]67   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy viscosity coef. at T- and F-points [m2/s or m4/s]
[12688]68   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bhm          !: eddy viscosity coef. for bi-Laplacian GM at W-points (cf. avm) [m4/s]
[10784]69   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dtensq       !: horizontal tension squared         (Smagorinsky only)
70   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dshesq       !: horizontal shearing strain squared (Smagorinsky only)
[7646]71   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   esqt, esqf   !: Square of the local gridscale (e1e2/(e1+e2))**2           
[5836]72
[9490]73   REAL(wp) ::   r1_2    = 0.5_wp            ! =1/2
[5836]74   REAL(wp) ::   r1_4    = 0.25_wp           ! =1/4
[7646]75   REAL(wp) ::   r1_8    = 0.125_wp          ! =1/8
[9490]76   REAL(wp) ::   r1_12   = 1._wp / 12._wp    ! =1/12
[5836]77   REAL(wp) ::   r1_288  = 1._wp / 288._wp   ! =1/( 12^2 * 2 )
78
[3]79   !! * Substitutions
[5836]80#  include "vectopt_loop_substitute.h90"
[3]81   !!----------------------------------------------------------------------
[9598]82   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[12535]83   !! $Id: ldfdyn.F90 11653 2019-10-04 12:34:02Z davestorkey $
[10068]84   !! Software governed by the CeCILL license (see ./LICENSE)
[3]85   !!----------------------------------------------------------------------
86CONTAINS
87
88   SUBROUTINE ldf_dyn_init
89      !!----------------------------------------------------------------------
90      !!                  ***  ROUTINE ldf_dyn_init  ***
91      !!                   
92      !! ** Purpose :   set the horizontal ocean dynamics physics
93      !!
[5836]94      !! ** Method  :   the eddy viscosity coef. specification depends on:
95      !!              - the operator:
96      !!             ln_dynldf_lap = T     laplacian operator
97      !!             ln_dynldf_blp = T   bilaplacian operator
[12535]98      !!             ln_dynldf_bgm = T   bi-Laplacian GM operator
[5836]99      !!              - the parameter nn_ahm_ijk_t:
100      !!    nn_ahm_ijk_t  =  0 => = constant
101      !!                  = 10 => = F(z) :     = constant with a reduction of 1/4 with depth
102      !!                  =-20 => = F(i,j)     = shape read in 'eddy_viscosity.nc' file
103      !!                  = 20    = F(i,j)     = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
104      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_viscosity.nc'  file
105      !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10)
106      !!                  = 31    = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator
107      !!                                                           or |u|e^3/12 bilaplacian operator )
[7646]108      !!                  = 32    = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 
109      !!                                                           (   L^2|D|      laplacian operator
110      !!                                                           or  L^4|D|/8  bilaplacian operator )
[12535]111      !!                - the parameter nn_bhm_ijk_t:
112      !!    nn_bhm_ijk_t  = 11 => = F(z) :     = constant, with BCs set to zero at top and bottom
113      !!    nn_bhm_ijk_t  = 12 => = F(z) :     = constant in interior, linear taper to zero at surface,
114      !!                                         for depths < rn_bgmzcrit, with BCs set to zero at top and bottom
115      !!    nn_bhm_ijk_t  = 13 => = F(z) :     = bhm0 X dx^2 X dz^2, with BCs set to zero at top and bottom
[3]116      !!----------------------------------------------------------------------
[9490]117      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
118      INTEGER  ::   ioptio, ierr, inum, ios, inn   ! local integer
[12688]119      INTEGER  ::   ikw                            ! local integer
[9490]120      REAL(wp) ::   zah0, zah_max, zUfac           ! local scalar
[12535]121      REAL(wp) ::   bhm0, lhscrit, lhscritmax      ! local scalar
[9490]122      CHARACTER(len=5) ::   cl_Units               ! units (m2/s or m4/s)
123      !!
[9526]124      NAMELIST/namdyn_ldf/ ln_dynldf_OFF, ln_dynldf_lap, ln_dynldf_blp,   &   ! type of operator
[12535]125         &                 ln_dynldf_bgm,                                 &   ! type of operator
[9526]126         &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso,   &   ! acting direction of the operator
127         &                 nn_ahm_ijk_t , rn_Uv    , rn_Lv,   rn_ahm_b,   &   ! lateral eddy coefficient
[12535]128         &                 nn_bhm_ijk_t, rn_bhm_b, rn_bgmzcrit,           &   ! lateral bi-Lap. GM coefficient
[9526]129         &                 rn_csmc      , rn_minfac    , rn_maxfac            ! Smagorinsky settings
[5836]130      !!----------------------------------------------------------------------
131      !
[4147]132      REWIND( numnam_ref )              ! Namelist namdyn_ldf in reference namelist : Lateral physics
133      READ  ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901)
[11536]134901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist' )
[3]135
[4147]136      REWIND( numnam_cfg )              ! Namelist namdyn_ldf in configuration namelist : Lateral physics
137      READ  ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 )
[11536]138902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist' )
[4624]139      IF(lwm) WRITE ( numond, namdyn_ldf )
[4147]140
[1601]141      IF(lwp) THEN                      ! Parameter print
[3]142         WRITE(numout,*)
143         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
144         WRITE(numout,*) '~~~~~~~'
[4147]145         WRITE(numout,*) '   Namelist namdyn_ldf : set lateral mixing parameters'
[5836]146         !
147         WRITE(numout,*) '      type :'
[9526]148         WRITE(numout,*) '         no explicit diffusion                ln_dynldf_OFF = ', ln_dynldf_OFF
[5836]149         WRITE(numout,*) '         laplacian operator                   ln_dynldf_lap = ', ln_dynldf_lap
150         WRITE(numout,*) '         bilaplacian operator                 ln_dynldf_blp = ', ln_dynldf_blp
[12535]151         WRITE(numout,*) '         bi-Laplacian GM operator             ln_dynldf_bgm = ', ln_dynldf_bgm
152         WRITE(numout,*) '         critical depth for taper (if used)   rn_bgmzcrit   = ', rn_bgmzcrit
[5836]153         !
154         WRITE(numout,*) '      direction of action :'
155         WRITE(numout,*) '         iso-level                            ln_dynldf_lev = ', ln_dynldf_lev
156         WRITE(numout,*) '         horizontal (geopotential)            ln_dynldf_hor = ', ln_dynldf_hor
157         WRITE(numout,*) '         iso-neutral                          ln_dynldf_iso = ', ln_dynldf_iso
158         !
159         WRITE(numout,*) '      coefficients :'
160         WRITE(numout,*) '         type of time-space variation         nn_ahm_ijk_t  = ', nn_ahm_ijk_t
[9490]161         WRITE(numout,*) '         lateral viscous velocity  (if cst)      rn_Uv      = ', rn_Uv, ' m/s'
162         WRITE(numout,*) '         lateral viscous length    (if cst)      rn_Lv      = ', rn_Lv, ' m'
163         WRITE(numout,*) '         background viscosity (iso-lap case)     rn_ahm_b   = ', rn_ahm_b, ' m2/s'
[12535]164         WRITE(numout,*) '         bi-Laplacian GM background viscosity    rn_bhm_b   = ', rn_bhm_b, ' m4/s'
[9019]165         WRITE(numout,*) '      Smagorinsky settings (nn_ahm_ijk_t  = 32) :'
[7646]166         WRITE(numout,*) '         Smagorinsky coefficient              rn_csmc       = ', rn_csmc
[9490]167         WRITE(numout,*) '         factor multiplier for eddy visc.'
168         WRITE(numout,*) '            lower limit (default 1.0)         rn_minfac    = ', rn_minfac
169         WRITE(numout,*) '            upper limit (default 1.0)         rn_maxfac    = ', rn_maxfac
[3]170      ENDIF
171
[9490]172      !
173      !           !==  type of lateral operator used  ==!   (set nldf_dyn)
174      !           !=====================================!
175      !
176      nldf_dyn = np_ERROR
177      ioptio = 0
[12535]178      !CW exclude combination of lap+blp (as in existing code), but allow other combinations
[9526]179      IF( ln_dynldf_OFF ) THEN   ;   nldf_dyn = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF
[12535]180      IF( ln_dynldf_lap ) THEN   ;                              ioptio = ioptio + 5   ;   ENDIF
181      IF( ln_dynldf_blp ) THEN   ;                              ioptio = ioptio + 5   ;   ENDIF
182      IF( ln_dynldf_bgm ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF
[9490]183      !
[12535]184      IF( (ioptio < 1).OR.(ioptio > 6) )   CALL ctl_stop( 'dyn_ldf_init: use at least ONE of the 4 operator options (NONE/lap/blp/bgm) but not (lap+blp)' )
185      !
186
[9526]187      IF(.NOT.ln_dynldf_OFF ) THEN     !==  direction ==>> type of operator  ==!
[9490]188         ioptio = 0
189         IF( ln_dynldf_lev )   ioptio = ioptio + 1
190         IF( ln_dynldf_hor )   ioptio = ioptio + 1
191         IF( ln_dynldf_iso )   ioptio = ioptio + 1
192         IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' )
193         !
194         !                             ! Set nldf_dyn, the type of lateral diffusion, from ln_dynldf_... logicals
195         ierr = 0
196         IF( ln_dynldf_lap ) THEN         ! laplacian operator
197            IF( ln_zco ) THEN                   ! z-coordinate
198               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
199               IF ( ln_dynldf_hor )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
200               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
201            ENDIF
202            IF( ln_zps ) THEN                   ! z-coordinate with partial step
203               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level              (no rotation)
204               IF ( ln_dynldf_hor )   nldf_dyn = np_lap     ! iso-level              (no rotation)
205               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
206            ENDIF
207            IF( ln_sco ) THEN                   ! s-coordinate
208               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
209               IF ( ln_dynldf_hor )   nldf_dyn = np_lap_i   ! horizontal             (   rotation)
210               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
211            ENDIF
212         ENDIF
213         !
214         IF( ln_dynldf_blp ) THEN         ! bilaplacian operator
215            IF( ln_zco ) THEN                   ! z-coordinate
216               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level = horizontal (no rotation)
217               IF( ln_dynldf_hor )   nldf_dyn = np_blp   ! iso-level = horizontal (no rotation)
218               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
219            ENDIF
220            IF( ln_zps ) THEN                   ! z-coordinate with partial step
221               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level              (no rotation)
222               IF( ln_dynldf_hor )   nldf_dyn = np_blp   ! iso-level              (no rotation)
223               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
224            ENDIF
225            IF( ln_sco ) THEN                   ! s-coordinate
226               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level              (no rotation)
227               IF( ln_dynldf_hor )   ierr = 2            ! horizontal             (   rotation)
228               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
229            ENDIF
230         ENDIF
231         !
[12535]232         IF( ln_dynldf_bgm ) THEN         ! bi-Laplacian GM operator
233            IF( ln_zco ) THEN                   ! z-coordinate
234               IF( ln_dynldf_lev )   nldf_dyn = np_bgm   ! iso-level = horizontal (no rotation)
235               IF( ln_dynldf_hor )   nldf_dyn = np_bgm   ! iso-level = horizontal (no rotation)
236               IF( ln_dynldf_iso )   ierr = 3            ! iso-neutral            (   rotation)
237            ENDIF
238            IF( ln_zps ) THEN                   ! z-coordinate with partial step
239               IF( ln_dynldf_lev )   nldf_dyn = np_bgm   ! iso-level              (no rotation)
240               IF( ln_dynldf_hor )   nldf_dyn = np_bgm   ! iso-level              (no rotation)
241               IF( ln_dynldf_iso )   ierr = 3            ! iso-neutral            (   rotation)
242            ENDIF
243            IF( ln_sco ) THEN                   ! s-coordinate
244               IF( ln_dynldf_lev )   nldf_dyn = np_bgm   ! iso-level              (no rotation)
245               IF( ln_dynldf_hor )   ierr = 3            ! horizontal             (   rotation)
246               IF( ln_dynldf_iso )   ierr = 3            ! iso-neutral            (   rotation)
247            ENDIF
248         ENDIF
249         !         
250
[9490]251         IF( ierr == 2 )   CALL ctl_stop( 'rotated bi-laplacian operator does not exist' )
252         !
[12535]253         IF( ierr == 3 )   CALL ctl_stop( 'rotated bi-Laplacian GM operator does not exist' )
254         !
255
[9490]256         IF( nldf_dyn == np_lap_i )   l_ldfslp = .TRUE.  ! rotation require the computation of the slopes
257         !
[3]258      ENDIF
[5836]259      !
[9490]260      IF(lwp) THEN
261         WRITE(numout,*)
262         SELECT CASE( nldf_dyn )
263         CASE( np_no_ldf )   ;   WRITE(numout,*) '   ==>>>   NO lateral viscosity'
264         CASE( np_lap    )   ;   WRITE(numout,*) '   ==>>>   iso-level laplacian operator'
265         CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background'
266         CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator'
[12535]267         CASE( np_bgm    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-Laplacian GM operator'
[9490]268         END SELECT
269         WRITE(numout,*)
[3]270      ENDIF
[9490]271     
[1601]272      !
[9490]273      !           !==  Space/time variation of eddy coefficients  ==!
274      !           !=================================================!
[5836]275      !
[9490]276      l_ldfdyn_time = .FALSE.                ! no time variation except in case defined below
[1601]277      !
[9526]278      IF( ln_dynldf_OFF ) THEN
[12535]279         !CW
280         !IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt and ahmf are not allocated'
[12688]281         IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt, ahmf, bhm are not allocated'
[12535]282         !
[9490]283         RETURN
[5836]284         !
[9490]285      ELSE                                   !==  a lateral diffusion operator is used  ==!
[5836]286         !
[9490]287         !                                         ! allocate the ahm arrays
[12688]288         ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , bhm(jpi,jpj,jpk) , STAT=ierr )
[12535]289         !
[9490]290         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays')
291         !
[10784]292         ahmt(:,:,:) = 0._wp                       ! init to 0 needed
293         ahmf(:,:,:) = 0._wp
[12688]294         bhm(:,:,:)  = 0._wp                       ! init to 0 needed
[12535]295         !                                        ! value of lap/blp eddy mixing coef.
[9490]296         IF(     ln_dynldf_lap ) THEN   ;   zUfac = r1_2 *rn_Uv   ;   inn = 1   ;   cl_Units = ' m2/s'   !   laplacian
297         ELSEIF( ln_dynldf_blp ) THEN   ;   zUfac = r1_12*rn_Uv   ;   inn = 3   ;   cl_Units = ' m4/s'   ! bilaplacian
[12535]298         ELSEIF( ln_dynldf_bgm ) THEN   ;   bhm0 = rn_bhm_b       ;             ;   cl_Units = ' m4/s'   ! bi-Laplacian GM
[9490]299         ENDIF
300         zah0    = zUfac *    rn_Lv**inn              ! mixing coefficient
301         zah_max = zUfac * (ra*rad)**inn              ! maximum reachable coefficient (value at the Equator)
302         !
303         SELECT CASE(  nn_ahm_ijk_t  )             !* Specification of space-time variations of ahmt, ahmf
304         !
[5836]305         CASE(   0  )      !==  constant  ==!
[9490]306            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity. = constant = ', zah0, cl_Units
307            ahmt(:,:,1:jpkm1) = zah0
308            ahmf(:,:,1:jpkm1) = zah0
[5836]309            !
310         CASE(  10  )      !==  fixed profile  ==!
[9490]311            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( depth )'
312            IF(lwp) WRITE(numout,*) '           surface viscous coef. = constant = ', zah0, cl_Units
313            ahmt(:,:,1) = zah0                        ! constant surface value
314            ahmf(:,:,1) = zah0
315            CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
[12535]316            !       
[5836]317         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
[9490]318            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j) read in eddy_viscosity.nc file'
[5836]319            CALL iom_open( 'eddy_viscosity_2D.nc', inum )
320            CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) )
321            CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) )
322            CALL iom_close( inum )
323            DO jk = 2, jpkm1
[9490]324               ahmt(:,:,jk) = ahmt(:,:,1)
325               ahmf(:,:,jk) = ahmf(:,:,1)
[5836]326            END DO
327            !
328         CASE(  20  )      !== fixed horizontal shape  ==!
[9490]329            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)'
330            IF(lwp) WRITE(numout,*) '           using a fixed viscous velocity = ', rn_Uv  ,' m/s   and   Lv = Max(e1,e2)'
331            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)'
332            CALL ldf_c2d( 'DYN', zUfac      , inn        , ahmt, ahmf )         ! surface value proportional to scale factor^inn
[5836]333            !
334         CASE( -30  )      !== fixed 3D shape read in file  ==!
[9490]335            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j,k) read in eddy_viscosity_3D.nc file'
[5836]336            CALL iom_open( 'eddy_viscosity_3D.nc', inum )
337            CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt )
338            CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf )
339            CALL iom_close( inum )
340            !
341         CASE(  30  )       !==  fixed 3D shape  ==!
[9490]342            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth )'
343            IF(lwp) WRITE(numout,*) '           using a fixed viscous velocity = ', rn_Uv  ,' m/s   and   Ld = Max(e1,e2)'
344            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)'
345            CALL ldf_c2d( 'DYN', zUfac      , inn        , ahmt, ahmf )         ! surface value proportional to scale factor^inn
346            CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )  ! reduction with depth
[5836]347            !
348         CASE(  31  )       !==  time varying 3D field  ==!
[9490]349            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth , time )'
350            IF(lwp) WRITE(numout,*) '           proportional to the local velocity : 1/2 |u|e (lap) or 1/12 |u|e^3 (blp)'
[5836]351            !
352            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
353            !
[7646]354         CASE(  32  )       !==  time varying 3D field  ==!
[9490]355            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth , time )'
356            IF(lwp) WRITE(numout,*) '           proportional to the local deformation rate and gridscale (Smagorinsky)'
[7646]357            !
358            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
359            !
[9490]360            !                          ! allocate arrays used in ldf_dyn.
[10784]361            ALLOCATE( dtensq(jpi,jpj,jpk) , dshesq(jpi,jpj,jpk) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr )
[7646]362            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays')
363            !
[10784]364            DO jj = 1, jpj             ! Set local gridscale values
365               DO ji = 1, jpi
[11653]366                  esqt(ji,jj) = ( 2._wp * e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 
367                  esqf(ji,jj) = ( 2._wp * e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 
[7646]368               END DO
369            END DO
370            !
[5836]371         CASE DEFAULT
372            CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
373         END SELECT
[12535]374!CW      Add separate structure options for bi-Laplacian GM, to allow it to be used in combination with other types of diffusion, above.
375
376         IF( ln_dynldf_bgm ) THEN
377            SELECT CASE (nn_bhm_ijk_t)
378               !CW Define bi-Laplacian GM diffusivity and test the related computational stability criterion   
379
380            CASE(  11  )      !==  fixed profile: constant except for zero at top and bottom, so that eddy-induced velocity, w*=0  ==!
381               IF(lwp) WRITE(numout,*) '   ==>>>   bi-Laplacian GM eddy viscosity = constant in interior, zero at top and bottom'
382               IF(lwp) WRITE(numout,*) '           interior viscous coef. = constant = ', bhm0, cl_Units
383
384               !     First set to constant in interior, all levels
385               DO jj = 2, jpjm1
386                  DO ji = 2, jpim1
[12688]387                     bhm(ji,jj,:) = bhm0
[12535]388                  ENDDO
389               ENDDO
390
391               !     Surface BC : set to zero
[12688]392               bhm(:,:,1)   = 0._wp
[12535]393               !     Flat bottom BC : set to zero
[12688]394               bhm(:,:,jpk) = 0._wp 
395               !     Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm
[12535]396
397               !CW Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping
398               ! \frac{{\kappa \Delta t}{(\Delta z)^2 (\Delta x)^2}< 1/32
399               ! test at T-point - may need to revise this choice!!!!
400
401               ! Find domain maximum value of LHS, then test inequality.
402
403               ! initialise
404               lhscritmax=0._wp
405
406               DO jj = 2, jpjm1
407                  DO ji = 2, jpim1
[12688]408                     ikw = mbkt(ji,jj) !bottom last wet T-level
409                     DO jk = 2, ikw
[12535]410                        lhscrit=bhm0*rdt/(e3t_n(ji,jj,jk)**2*e1t(ji,jj)**2)
411                        IF(lhscrit .gt. lhscritmax) lhscritmax=lhscrit
412                     ENDDO
413                  ENDDO
414               ENDDO
415
416               IF ( lhscrit .ge. 1._wp/32._wp) THEN
417                  IF(lwp) THEN
418                     WRITE(numout,*)
419                     WRITE(numout,*) ' WARNING : Bi-Laplacian GM eddy viscosity is not',              &
420                          &            ' consistent with the model time step and grid setup: ', &
421                          &            ' LHS=',lhscrit, & 
422                          &            'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)'
423                     WRITE(numout,*) ' ========='
424                     WRITE(numout,*)
425                  ENDIF
426                  ! CW: warn, don't stop, as we may wish to violate this theoretical limit.
427                  ! CALL ctl_stop( 'ldf_dyn_init', 'Inconsistent Bi-Lap. GM diffusivity')
428               ELSE
429                  IF(lwp) THEN
430                     WRITE(numout,*)
431                     WRITE(numout,*) ' INFORMATION : Bi-Laplacian GM eddy viscosity is ',              &
432                          &            ' consistent with the model time step and grid setup: ', &
433                          &            ' LHS=',lhscrit, & 
434                          &            'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)'
435                     WRITE(numout,*) ' ========='
436                     WRITE(numout,*)
437                  ENDIF
438               ENDIF
439               !-----------------------------
440            CASE(  12  )      !==  fixed profile: constant in interior, zero at bottom and linearly-tapering to zero at top, over depth shallower than rn_bgmzcrit  ==!
441               IF(lwp) WRITE(numout,*) '   ==>>>   bi-Laplacian GM eddy viscosity = constant in interior, zero at bottom and'
442               IF(lwp) WRITE(numout,*) '           linearly-tapering to zero at top, over depth shallower than ', rn_bgmzcrit, ' metres.'
443               IF(lwp) WRITE(numout,*) '           Interior viscous coef. = constant = ', bhm0, cl_Units
444
445               !     Impose linear taper from interior value to zero at zero depth, for depths < rn_bgmzcrit
446               !     N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet.
447
448               DO jk = 1, jpk
449
450                  IF (gdept_1d(jk) .lt. rn_bgmzcrit) THEN
451
452                     DO jj = 2, jpjm1
453                        DO ji = 2, jpim1
[12688]454                           bhm(ji,jj,jk) = bhm0 * gdepw_1d(jk)/rn_bgmzcrit
[12535]455                        ENDDO
456                     ENDDO
457
458                  ELSE
459
460                     ! constant at depth rn_bgmzcrit and larger
461                     DO jj = 2, jpjm1
462                        DO ji = 2, jpim1
[12688]463                           bhm(ji,jj,jk) = bhm0
[12535]464                        ENDDO
465                     ENDDO
466
467                  ENDIF
468
469               ENDDO
470
471               !     Surface BC : set to zero
[12688]472               bhm(:,:,1)   = 0._wp 
[12535]473               !     Flat bottom BC : set to zero
[12688]474               bhm(:,:,jpk) = 0._wp 
475               !     Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm
[12535]476
477               !--------------------------------
478            CASE(  13  )      !==  fixed profile: steady profile of the form bhm0 X dx^2 X dz^2 in interior, zero at bottom and top  ==!
479
480               cl_Units = ' 1/s'
481
482               IF(lwp) WRITE(numout,*) '   ==>>>   bi-Laplacian GM eddy viscosity : steady profile of the form'
483               IF(lwp) WRITE(numout,*) '           bhm0 X dx^2 X dz^2 in interior, zero boundary conds. at bottom and top.'
484               IF(lwp) WRITE(numout,*) '           In this case, bhm0 is the constant of proportionality,'
485               IF(lwp) WRITE(numout,*) '           dimensioned as [diffusivity]/L^4, i.e. bhm0=', bhm0, cl_Units
486
487               cl_Units = ' m4/s'
488
489               !     N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet.
490
491               DO jk = 1, jpk
492                  ! constant at depth rn_bgmzcrit and larger
493                  DO jj = 2, jpjm1
494                     DO ji = 2, jpim1
[12688]495                        bhm(ji,jj,jk) = bhm0*e1t(ji,jj)**2*e3w_n(ji,jj,jk)**2
[12535]496                     ENDDO
497                  ENDDO
498               ENDDO
499
500               !     Surface BC : set to zero
[12688]501               bhm(:,:,1)   = 0._wp 
[12535]502               !     Flat bottom BC : set to zero
[12688]503               bhm(:,:,jpk) = 0._wp 
504               !     Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm
[12535]505
506               !--------------------------------
507
508            CASE DEFAULT
509               CALL ctl_stop('ldf_dyn_init: wrong choice for nn_bhm_ijk_t, the type of space-time variation of bhm')
510            END SELECT
511         ENDIF ! ln_dynldf_bgm
[5836]512         !
[9490]513         IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation
514            IF(     ln_dynldf_lap ) THEN                 !   laplacian operator (mask only)
515               ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1)
516               ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1)
517            ELSEIF( ln_dynldf_blp ) THEN                 ! bilaplacian operator (square root + mask)
518               ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1)
519               ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1)
[12535]520               !CW
521            ELSEIF( ln_dynldf_bgm ) THEN                 !bi-Laplacian GM operator (mask only)
[12688]522               bhm(:,:,1:jpkm1) =       bhm(:,:,1:jpkm1)   * wmask(:,:,1:jpkm1)
[12535]523               !
[9490]524            ENDIF
[5836]525         ENDIF
526         !
[71]527      ENDIF
[1601]528      !
[12535]529    END SUBROUTINE ldf_dyn_init
[71]530
531
[12535]532    SUBROUTINE ldf_dyn( kt )
[3]533      !!----------------------------------------------------------------------
[5836]534      !!                  ***  ROUTINE ldf_dyn  ***
535      !!
536      !! ** Purpose :   update at kt the momentum lateral mixing coeff. (ahmt and ahmf)
[3]537      !!
[5836]538      !! ** Method  :   time varying eddy viscosity coefficients:
[3]539      !!
[5836]540      !!    nn_ahm_ijk_t = 31    ahmt, ahmf = F(i,j,k,t) = F(local velocity)
541      !!                         ( |u|e /12  or  |u|e^3/12 for laplacian or bilaplacian operator )
[3]542      !!
[7646]543      !!    nn_ahm_ijk_t = 32    ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 
544      !!                         ( L^2|D|    or  L^4|D|/8  for laplacian or bilaplacian operator )
545      !!
546      !! ** note    :    in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian
547      !! ** action  :    ahmt, ahmf   updated at each time step
[3]548      !!----------------------------------------------------------------------
[5836]549      INTEGER, INTENT(in) ::   kt   ! time step index
[1601]550      !
[5836]551      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[10784]552      REAL(wp) ::   zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zemax   ! local scalar (option 31)
553      REAL(wp) ::   zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb     ! local scalar (option 32)
[5836]554      !!----------------------------------------------------------------------
555      !
[9019]556      IF( ln_timing )   CALL timing_start('ldf_dyn')
[5836]557      !
558      SELECT CASE(  nn_ahm_ijk_t  )       !== Eddy vicosity coefficients ==!
[12535]559         !
[5836]560      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity )
561         !
[7646]562         IF( ln_dynldf_lap   ) THEN        ! laplacian operator : |u| e /12 = |u/144| e
[5836]563            DO jk = 1, jpkm1
564               DO jj = 2, jpjm1
565                  DO ji = fs_2, fs_jpim1
566                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
567                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
[10784]568                     zemax = MAX( e1t(ji,jj) , e2t(ji,jj) )
569                     ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax * tmask(ji,jj,jk)      ! 288= 12*12 * 2
[5836]570                  END DO
571               END DO
[10784]572               DO jj = 1, jpjm1
573                  DO ji = 1, fs_jpim1
574                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
575                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
576                     zemax = MAX( e1f(ji,jj) , e2f(ji,jj) )
577                     ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax * fmask(ji,jj,jk)      ! 288= 12*12 * 2
578                  END DO
579               END DO
[5836]580            END DO
581         ELSEIF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e
582            DO jk = 1, jpkm1
583               DO jj = 2, jpjm1
584                  DO ji = fs_2, fs_jpim1
585                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
586                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
[10784]587                     zemax = MAX( e1t(ji,jj) , e2t(ji,jj) )
588                     ahmt(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax  ) * zemax * tmask(ji,jj,jk)
[5836]589                  END DO
590               END DO
[10784]591               DO jj = 1, jpjm1
592                  DO ji = 1, fs_jpim1
593                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
594                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
595                     zemax = MAX( e1f(ji,jj) , e2f(ji,jj) )
596                     ahmf(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax  ) * zemax * fmask(ji,jj,jk)
597                  END DO
598               END DO
[5836]599            END DO
600         ENDIF
601         !
[10425]602         CALL lbc_lnk_multi( 'ldfdyn', ahmt, 'T', 1.,  ahmf, 'F', 1. )
[5836]603         !
[7646]604         !
605      CASE(  32  )       !==  time varying 3D field  ==!   = F( local deformation rate and gridscale ) (Smagorinsky)
606         !
607         IF( ln_dynldf_lap .OR. ln_dynldf_blp  ) THEN        ! laplacian operator : (C_smag/pi)^2 L^2 |D|
608            !
[10784]609            zcmsmag   = (rn_csmc/rpi)**2                                            ! (C_smag/pi)^2
610            zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag )         ! lower limit stability factor scaling
611            zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt )               ! upper limit stability factor scaling
[7646]612            IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo          ! provide |U|L^3/12 lower limit instead
613            !                                                                       ! of |U|L^3/16 in blp case
614            DO jk = 1, jpkm1
615               !
[10784]616               DO jj = 2, jpjm1
617                  DO ji = 2, jpim1
618                     zdb =   ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) * r1_e1t(ji,jj) * e2t(ji,jj)  &
[12535]619                          &  - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) * r1_e2t(ji,jj) * e1t(ji,jj)
[10784]620                     dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk)
[7646]621                  END DO
622               END DO
623               !
624               DO jj = 1, jpjm1
625                  DO ji = 1, jpim1
[10784]626                     zdb =   ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) * r1_e2f(ji,jj) * e1f(ji,jj)  &
[12535]627                          &  + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) * r1_e1f(ji,jj) * e2f(ji,jj)
[10784]628                     dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk)
[7646]629                  END DO
630               END DO
631               !
[10784]632            END DO
633            !
634            CALL lbc_lnk_multi( 'ldfdyn', dtensq, 'T', 1. )  ! lbc_lnk on dshesq not needed
635            !
636            DO jk = 1, jpkm1
[12535]637               !
[10784]638               DO jj = 2, jpjm1                                ! T-point value
[7646]639                  DO ji = fs_2, fs_jpim1
640                     !
641                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
642                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
[10784]643                     !
[7646]644                     zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2
[10784]645                     ahmt(ji,jj,jk) = zdelta * SQRT(          dtensq(ji  ,jj,jk) +                         &
[12535]646                          &                            r1_4 * ( dshesq(ji  ,jj,jk) + dshesq(ji  ,jj-1,jk) +  &
647                          &                                     dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) )
[10784]648                     ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2
649                     ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt)
650                     !
651                  END DO
652               END DO
653               !
654               DO jj = 1, jpjm1                                ! F-point value
655                  DO ji = 1, fs_jpim1
656                     !
657                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
658                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
659                     !
[7646]660                     zdelta         = zcmsmag * esqf(ji,jj)                                        ! L^2 * (C_smag/pi)^2
[10784]661                     ahmf(ji,jj,jk) = zdelta * SQRT(          dshesq(ji  ,jj,jk) +                         &
[12535]662                          &                            r1_4 * ( dtensq(ji  ,jj,jk) + dtensq(ji  ,jj+1,jk) +  &
663                          &                                     dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) )
[10784]664                     ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2
665                     ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt)
[7646]666                     !
667                  END DO
668               END DO
[10784]669               !
[7646]670            END DO
671            !
672         ENDIF
673         !
674         IF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8)
[9490]675            !                          !                      = sqrt( A_lap_smag L^2/8 )
676            !                          ! stability limits already applied to laplacian values
677            !                          ! effective default limits are 1/12 |U|L^3 < B_hm < 1//(32*2dt) L^4
[7646]678            DO jk = 1, jpkm1
679               DO jj = 2, jpjm1
680                  DO ji = fs_2, fs_jpim1
[9490]681                     ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) )
[10784]682                  END DO
683               END DO
684               DO jj = 1, jpjm1
685                  DO ji = 1, fs_jpim1
[9490]686                     ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) )
[7646]687                  END DO
688               END DO
689            END DO
690            !
691         ENDIF
692         !
[10425]693         CALL lbc_lnk_multi( 'ldfdyn', ahmt, 'T', 1. , ahmf, 'F', 1. )
[7646]694         !
[5836]695      END SELECT
696      !
697      CALL iom_put( "ahmt_2d", ahmt(:,:,1) )   ! surface u-eddy diffusivity coeff.
698      CALL iom_put( "ahmf_2d", ahmf(:,:,1) )   ! surface v-eddy diffusivity coeff.
699      CALL iom_put( "ahmt_3d", ahmt(:,:,:) )   ! 3D      u-eddy diffusivity coeff.
700      CALL iom_put( "ahmf_3d", ahmf(:,:,:) )   ! 3D      v-eddy diffusivity coeff.
701      !
[9019]702      IF( ln_timing )   CALL timing_stop('ldf_dyn')
[5836]703      !
704   END SUBROUTINE ldf_dyn
[3]705
706   !!======================================================================
707END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.