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/trunk/src/OCE/LDF – NEMO

source: NEMO/trunk/src/OCE/LDF/ldfdyn.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 30.2 KB
Line 
1MODULE ldfdyn
2   !!======================================================================
3   !!                       ***  MODULE  ldfdyn  ***
4   !! Ocean physics:  lateral viscosity coefficient
5   !!=====================================================================
6   !! History :  OPA  ! 1997-07  (G. Madec)  multi dimensional coefficients
7   !!   NEMO     1.0  ! 2002-09  (G. Madec)  F90: Free form and module
8   !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
9   !!                 !                                  add velocity dependent coefficient and optional read in file
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   ldf_dyn_init  : initialization, namelist read, and parameters control
14   !!   ldf_dyn       : update lateral eddy viscosity coefficients at each time step
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers   
17   USE dom_oce         ! ocean space and time domain
18   USE phycst          ! physical constants
19   USE ldfslp          ! lateral diffusion: slopes of mixing orientation
20   USE ldfc1d_c2d      ! lateral diffusion: 1D and 2D cases
21   !
22   USE in_out_manager  ! I/O manager
23   USE iom             ! I/O module for ehanced bottom friction file
24   USE timing          ! Timing
25   USE lib_mpp         ! distribued memory computing library
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   ldf_dyn_init   ! called by nemogcm.F90
32   PUBLIC   ldf_dyn        ! called by step.F90
33
34   !                                    !!* Namelist namdyn_ldf : lateral mixing on momentum *
35   LOGICAL , PUBLIC ::   ln_dynldf_OFF   !: No operator (i.e. no explicit diffusion)
36   LOGICAL , PUBLIC ::   ln_dynldf_lap   !: laplacian operator
37   LOGICAL , PUBLIC ::   ln_dynldf_blp   !: bilaplacian operator
38   LOGICAL , PUBLIC ::   ln_dynldf_lev   !: iso-level direction
39   LOGICAL , PUBLIC ::   ln_dynldf_hor   !: horizontal (geopotential) direction
40!  LOGICAL , PUBLIC ::   ln_dynldf_iso   !: iso-neutral direction                        (see ldfslp)
41   INTEGER , PUBLIC ::   nn_ahm_ijk_t    !: choice of time & space variations of the lateral eddy viscosity coef.
42   !                                        !  time invariant coefficients:  aht = 1/2  Ud*Ld   (lap case)
43   !                                           !                             bht = 1/12 Ud*Ld^3 (blp case)
44   REAL(wp), PUBLIC ::   rn_Uv                 !: lateral viscous velocity  [m/s]
45   REAL(wp), PUBLIC ::   rn_Lv                 !: lateral viscous length    [m]
46   !                                        ! Smagorinsky viscosity  (nn_ahm_ijk_t = 32)
47   REAL(wp), PUBLIC ::   rn_csmc               !: Smagorinsky constant of proportionality
48   REAL(wp), PUBLIC ::   rn_minfac             !: Multiplicative factor of theorectical minimum Smagorinsky viscosity
49   REAL(wp), PUBLIC ::   rn_maxfac             !: Multiplicative factor of theorectical maximum Smagorinsky viscosity
50   !                                        ! iso-neutral laplacian (ln_dynldf_lap=ln_dynldf_iso=T)
51   REAL(wp), PUBLIC ::   rn_ahm_b              !: lateral laplacian background eddy viscosity  [m2/s]
52
53   !                                    !!* Parameter to control the type of lateral viscous operator
54   INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10                       !: error in setting the operator
55   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00                       !: without operator (i.e. no lateral viscous trend)
56   !                          !!      laplacian     !    bilaplacian    !
57   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  !: iso-level operator
58   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       !: iso-neutral or geopotential operator
59   !
60   INTEGER           , PUBLIC ::   nldf_dyn         !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals)
61   LOGICAL           , PUBLIC ::   l_ldfdyn_time    !: flag for time variation of the lateral eddy viscosity coef.
62
63   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy viscosity coef. at T- and F-points [m2/s or m4/s]
64   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dtensq       !: horizontal tension squared         (Smagorinsky only)
65   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dshesq       !: horizontal shearing strain squared (Smagorinsky only)
66   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   esqt, esqf   !: Square of the local gridscale (e1e2/(e1+e2))**2           
67
68   REAL(wp) ::   r1_2    = 0.5_wp            ! =1/2
69   REAL(wp) ::   r1_4    = 0.25_wp           ! =1/4
70   REAL(wp) ::   r1_8    = 0.125_wp          ! =1/8
71   REAL(wp) ::   r1_12   = 1._wp / 12._wp    ! =1/12
72   REAL(wp) ::   r1_288  = 1._wp / 288._wp   ! =1/( 12^2 * 2 )
73
74   !! * Substitutions
75#  include "do_loop_substitute.h90"
76   !!----------------------------------------------------------------------
77   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
78   !! $Id$
79   !! Software governed by the CeCILL license (see ./LICENSE)
80   !!----------------------------------------------------------------------
81CONTAINS
82
83   SUBROUTINE ldf_dyn_init
84      !!----------------------------------------------------------------------
85      !!                  ***  ROUTINE ldf_dyn_init  ***
86      !!                   
87      !! ** Purpose :   set the horizontal ocean dynamics physics
88      !!
89      !! ** Method  :   the eddy viscosity coef. specification depends on:
90      !!              - the operator:
91      !!             ln_dynldf_lap = T     laplacian operator
92      !!             ln_dynldf_blp = T   bilaplacian operator
93      !!              - the parameter nn_ahm_ijk_t:
94      !!    nn_ahm_ijk_t  =  0 => = constant
95      !!                  = 10 => = F(z) :     = constant with a reduction of 1/4 with depth
96      !!                  =-20 => = F(i,j)     = shape read in 'eddy_viscosity.nc' file
97      !!                  = 20    = F(i,j)     = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
98      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_viscosity.nc'  file
99      !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10)
100      !!                  = 31    = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator
101      !!                                                           or |u|e^3/12 bilaplacian operator )
102      !!                  = 32    = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 
103      !!                                                           (   L^2|D|      laplacian operator
104      !!                                                           or  L^4|D|/8  bilaplacian operator )
105      !!----------------------------------------------------------------------
106      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
107      INTEGER  ::   ioptio, ierr, inum, ios, inn   ! local integer
108      REAL(wp) ::   zah0, zah_max, zUfac           ! local scalar
109      CHARACTER(len=5) ::   cl_Units               ! units (m2/s or m4/s)
110      !!
111      NAMELIST/namdyn_ldf/ ln_dynldf_OFF, ln_dynldf_lap, ln_dynldf_blp,   &   ! type of operator
112         &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso,   &   ! acting direction of the operator
113         &                 nn_ahm_ijk_t , rn_Uv    , rn_Lv,   rn_ahm_b,   &   ! lateral eddy coefficient
114         &                 rn_csmc      , rn_minfac    , rn_maxfac            ! Smagorinsky settings
115      !!----------------------------------------------------------------------
116      !
117      READ  ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901)
118901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist' )
119
120      READ  ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 )
121902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist' )
122      IF(lwm) WRITE ( numond, namdyn_ldf )
123
124      IF(lwp) THEN                      ! Parameter print
125         WRITE(numout,*)
126         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
127         WRITE(numout,*) '~~~~~~~'
128         WRITE(numout,*) '   Namelist namdyn_ldf : set lateral mixing parameters'
129         !
130         WRITE(numout,*) '      type :'
131         WRITE(numout,*) '         no explicit diffusion                ln_dynldf_OFF = ', ln_dynldf_OFF
132         WRITE(numout,*) '         laplacian operator                   ln_dynldf_lap = ', ln_dynldf_lap
133         WRITE(numout,*) '         bilaplacian operator                 ln_dynldf_blp = ', ln_dynldf_blp
134         !
135         WRITE(numout,*) '      direction of action :'
136         WRITE(numout,*) '         iso-level                            ln_dynldf_lev = ', ln_dynldf_lev
137         WRITE(numout,*) '         horizontal (geopotential)            ln_dynldf_hor = ', ln_dynldf_hor
138         WRITE(numout,*) '         iso-neutral                          ln_dynldf_iso = ', ln_dynldf_iso
139         !
140         WRITE(numout,*) '      coefficients :'
141         WRITE(numout,*) '         type of time-space variation         nn_ahm_ijk_t  = ', nn_ahm_ijk_t
142         WRITE(numout,*) '         lateral viscous velocity  (if cst)      rn_Uv      = ', rn_Uv, ' m/s'
143         WRITE(numout,*) '         lateral viscous length    (if cst)      rn_Lv      = ', rn_Lv, ' m'
144         WRITE(numout,*) '         background viscosity (iso-lap case)     rn_ahm_b   = ', rn_ahm_b, ' m2/s'
145         !
146         WRITE(numout,*) '      Smagorinsky settings (nn_ahm_ijk_t  = 32) :'
147         WRITE(numout,*) '         Smagorinsky coefficient              rn_csmc       = ', rn_csmc
148         WRITE(numout,*) '         factor multiplier for eddy visc.'
149         WRITE(numout,*) '            lower limit (default 1.0)         rn_minfac    = ', rn_minfac
150         WRITE(numout,*) '            upper limit (default 1.0)         rn_maxfac    = ', rn_maxfac
151      ENDIF
152
153      !
154      !           !==  type of lateral operator used  ==!   (set nldf_dyn)
155      !           !=====================================!
156      !
157      nldf_dyn = np_ERROR
158      ioptio = 0
159      IF( ln_dynldf_OFF ) THEN   ;   nldf_dyn = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF
160      IF( ln_dynldf_lap ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF
161      IF( ln_dynldf_blp ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF
162      IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' )
163      !
164      IF(.NOT.ln_dynldf_OFF ) THEN     !==  direction ==>> type of operator  ==!
165         ioptio = 0
166         IF( ln_dynldf_lev )   ioptio = ioptio + 1
167         IF( ln_dynldf_hor )   ioptio = ioptio + 1
168         IF( ln_dynldf_iso )   ioptio = ioptio + 1
169         IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' )
170         !
171         !                             ! Set nldf_dyn, the type of lateral diffusion, from ln_dynldf_... logicals
172         ierr = 0
173         IF( ln_dynldf_lap ) THEN         ! laplacian operator
174            IF( ln_zco ) THEN                   ! z-coordinate
175               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
176               IF ( ln_dynldf_hor )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
177               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
178            ENDIF
179            IF( ln_zps ) THEN                   ! z-coordinate with partial step
180               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level              (no rotation)
181               IF ( ln_dynldf_hor )   nldf_dyn = np_lap     ! iso-level              (no rotation)
182               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
183            ENDIF
184            IF( ln_sco ) THEN                   ! s-coordinate
185               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
186               IF ( ln_dynldf_hor )   nldf_dyn = np_lap_i   ! horizontal             (   rotation)
187               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
188            ENDIF
189         ENDIF
190         !
191         IF( ln_dynldf_blp ) THEN         ! bilaplacian operator
192            IF( ln_zco ) THEN                   ! z-coordinate
193               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level = horizontal (no rotation)
194               IF( ln_dynldf_hor )   nldf_dyn = np_blp   ! iso-level = horizontal (no rotation)
195               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
196            ENDIF
197            IF( ln_zps ) THEN                   ! z-coordinate with partial step
198               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level              (no rotation)
199               IF( ln_dynldf_hor )   nldf_dyn = np_blp   ! iso-level              (no rotation)
200               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
201            ENDIF
202            IF( ln_sco ) THEN                   ! s-coordinate
203               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level              (no rotation)
204               IF( ln_dynldf_hor )   ierr = 2            ! horizontal             (   rotation)
205               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
206            ENDIF
207         ENDIF
208         !
209         IF( ierr == 2 )   CALL ctl_stop( 'rotated bi-laplacian operator does not exist' )
210         !
211         IF( nldf_dyn == np_lap_i )   l_ldfslp = .TRUE.  ! rotation require the computation of the slopes
212         !
213      ENDIF
214      !
215      IF(lwp) THEN
216         WRITE(numout,*)
217         SELECT CASE( nldf_dyn )
218         CASE( np_no_ldf )   ;   WRITE(numout,*) '   ==>>>   NO lateral viscosity'
219         CASE( np_lap    )   ;   WRITE(numout,*) '   ==>>>   iso-level laplacian operator'
220         CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background'
221         CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator'
222         END SELECT
223         WRITE(numout,*)
224      ENDIF
225     
226      !
227      !           !==  Space/time variation of eddy coefficients  ==!
228      !           !=================================================!
229      !
230      l_ldfdyn_time = .FALSE.                ! no time variation except in case defined below
231      !
232      IF( ln_dynldf_OFF ) THEN
233         IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt and ahmf are not allocated'
234         RETURN
235         !
236      ELSE                                   !==  a lateral diffusion operator is used  ==!
237         !
238         !                                         ! allocate the ahm arrays
239         ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr )
240         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays')
241         !
242         ahmt(:,:,:) = 0._wp                       ! init to 0 needed
243         ahmf(:,:,:) = 0._wp
244         !
245         !                                         ! value of lap/blp eddy mixing coef.
246         IF(     ln_dynldf_lap ) THEN   ;   zUfac = r1_2 *rn_Uv   ;   inn = 1   ;   cl_Units = ' m2/s'   !   laplacian
247         ELSEIF( ln_dynldf_blp ) THEN   ;   zUfac = r1_12*rn_Uv   ;   inn = 3   ;   cl_Units = ' m4/s'   ! bilaplacian
248         ENDIF
249         zah0    = zUfac *    rn_Lv**inn              ! mixing coefficient
250         zah_max = zUfac * (ra*rad)**inn              ! maximum reachable coefficient (value at the Equator)
251         !
252         SELECT CASE(  nn_ahm_ijk_t  )             !* Specification of space-time variations of ahmt, ahmf
253         !
254         CASE(   0  )      !==  constant  ==!
255            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity. = constant = ', zah0, cl_Units
256            ahmt(:,:,1:jpkm1) = zah0
257            ahmf(:,:,1:jpkm1) = zah0
258            !
259         CASE(  10  )      !==  fixed profile  ==!
260            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( depth )'
261            IF(lwp) WRITE(numout,*) '           surface viscous coef. = constant = ', zah0, cl_Units
262            ahmt(:,:,1) = zah0                        ! constant surface value
263            ahmf(:,:,1) = zah0
264            CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
265            !
266         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
267            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j) read in eddy_viscosity.nc file'
268            CALL iom_open( 'eddy_viscosity_2D.nc', inum )
269            CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) )
270            CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) )
271            CALL iom_close( inum )
272            DO jk = 2, jpkm1
273               ahmt(:,:,jk) = ahmt(:,:,1)
274               ahmf(:,:,jk) = ahmf(:,:,1)
275            END DO
276            !
277         CASE(  20  )      !== fixed horizontal shape  ==!
278            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)'
279            IF(lwp) WRITE(numout,*) '           using a fixed viscous velocity = ', rn_Uv  ,' m/s   and   Lv = Max(e1,e2)'
280            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)'
281            CALL ldf_c2d( 'DYN', zUfac      , inn        , ahmt, ahmf )         ! surface value proportional to scale factor^inn
282            !
283         CASE( -30  )      !== fixed 3D shape read in file  ==!
284            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j,k) read in eddy_viscosity_3D.nc file'
285            CALL iom_open( 'eddy_viscosity_3D.nc', inum )
286            CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt )
287            CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf )
288            CALL iom_close( inum )
289            !
290         CASE(  30  )       !==  fixed 3D shape  ==!
291            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth )'
292            IF(lwp) WRITE(numout,*) '           using a fixed viscous velocity = ', rn_Uv  ,' m/s   and   Ld = Max(e1,e2)'
293            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)'
294            CALL ldf_c2d( 'DYN', zUfac      , inn        , ahmt, ahmf )         ! surface value proportional to scale factor^inn
295            CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )  ! reduction with depth
296            !
297         CASE(  31  )       !==  time varying 3D field  ==!
298            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth , time )'
299            IF(lwp) WRITE(numout,*) '           proportional to the local velocity : 1/2 |u|e (lap) or 1/12 |u|e^3 (blp)'
300            !
301            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
302            !
303         CASE(  32  )       !==  time varying 3D field  ==!
304            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth , time )'
305            IF(lwp) WRITE(numout,*) '           proportional to the local deformation rate and gridscale (Smagorinsky)'
306            !
307            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
308            !
309            !                          ! allocate arrays used in ldf_dyn.
310            ALLOCATE( dtensq(jpi,jpj,jpk) , dshesq(jpi,jpj,jpk) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr )
311            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays')
312            !
313            DO_2D_11_11
314               esqt(ji,jj) = ( 2._wp * e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 
315               esqf(ji,jj) = ( 2._wp * e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 
316            END_2D
317            !
318         CASE DEFAULT
319            CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
320         END SELECT
321         !
322         IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation
323            IF(     ln_dynldf_lap ) THEN                 !   laplacian operator (mask only)
324               ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1)
325               ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1)
326            ELSEIF( ln_dynldf_blp ) THEN                 ! bilaplacian operator (square root + mask)
327               ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1)
328               ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1)
329            ENDIF
330         ENDIF
331         !
332      ENDIF
333      !
334   END SUBROUTINE ldf_dyn_init
335
336
337   SUBROUTINE ldf_dyn( kt, Kbb )
338      !!----------------------------------------------------------------------
339      !!                  ***  ROUTINE ldf_dyn  ***
340      !!
341      !! ** Purpose :   update at kt the momentum lateral mixing coeff. (ahmt and ahmf)
342      !!
343      !! ** Method  :   time varying eddy viscosity coefficients:
344      !!
345      !!    nn_ahm_ijk_t = 31    ahmt, ahmf = F(i,j,k,t) = F(local velocity)
346      !!                         ( |u|e /12  or  |u|e^3/12 for laplacian or bilaplacian operator )
347      !!
348      !!    nn_ahm_ijk_t = 32    ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 
349      !!                         ( L^2|D|    or  L^4|D|/8  for laplacian or bilaplacian operator )
350      !!
351      !! ** note    :    in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian
352      !! ** action  :    ahmt, ahmf   updated at each time step
353      !!----------------------------------------------------------------------
354      INTEGER, INTENT(in) ::   kt   ! time step index
355      INTEGER, INTENT(in) ::   Kbb  ! ocean time level indices
356      !
357      INTEGER  ::   ji, jj, jk   ! dummy loop indices
358      REAL(wp) ::   zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zemax   ! local scalar (option 31)
359      REAL(wp) ::   zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb     ! local scalar (option 32)
360      !!----------------------------------------------------------------------
361      !
362      IF( ln_timing )   CALL timing_start('ldf_dyn')
363      !
364      SELECT CASE(  nn_ahm_ijk_t  )       !== Eddy vicosity coefficients ==!
365      !
366      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity )
367         !
368         IF( ln_dynldf_lap   ) THEN        ! laplacian operator : |u| e /12 = |u/144| e
369            DO jk = 1, jpkm1
370               DO_2D_00_00
371                  zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
372                  zu2pv2_ij_m1 = uu(ji-1,jj  ,jk,Kbb) * uu(ji-1,jj  ,jk,Kbb) + vv(ji  ,jj-1,jk,Kbb) * vv(ji  ,jj-1,jk,Kbb)
373                  zemax = MAX( e1t(ji,jj) , e2t(ji,jj) )
374                  ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax * tmask(ji,jj,jk)      ! 288= 12*12 * 2
375               END_2D
376               DO_2D_10_10
377                  zu2pv2_ij_p1 = uu(ji  ,jj+1,jk, Kbb) * uu(ji  ,jj+1,jk, Kbb) + vv(ji+1,jj  ,jk, Kbb) * vv(ji+1,jj  ,jk, Kbb)
378                  zu2pv2_ij    = uu(ji  ,jj  ,jk, Kbb) * uu(ji  ,jj  ,jk, Kbb) + vv(ji  ,jj  ,jk, Kbb) * vv(ji  ,jj  ,jk, Kbb)
379                  zemax = MAX( e1f(ji,jj) , e2f(ji,jj) )
380                  ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax * fmask(ji,jj,jk)      ! 288= 12*12 * 2
381               END_2D
382            END DO
383         ELSEIF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e
384            DO jk = 1, jpkm1
385               DO_2D_00_00
386                  zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
387                  zu2pv2_ij_m1 = uu(ji-1,jj  ,jk,Kbb) * uu(ji-1,jj  ,jk,Kbb) + vv(ji  ,jj-1,jk,Kbb) * vv(ji  ,jj-1,jk,Kbb)
388                  zemax = MAX( e1t(ji,jj) , e2t(ji,jj) )
389                  ahmt(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax  ) * zemax * tmask(ji,jj,jk)
390               END_2D
391               DO_2D_10_10
392                  zu2pv2_ij_p1 = uu(ji  ,jj+1,jk, Kbb) * uu(ji  ,jj+1,jk, Kbb) + vv(ji+1,jj  ,jk, Kbb) * vv(ji+1,jj  ,jk, Kbb)
393                  zu2pv2_ij    = uu(ji  ,jj  ,jk, Kbb) * uu(ji  ,jj  ,jk, Kbb) + vv(ji  ,jj  ,jk, Kbb) * vv(ji  ,jj  ,jk, Kbb)
394                  zemax = MAX( e1f(ji,jj) , e2f(ji,jj) )
395                  ahmf(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax  ) * zemax * fmask(ji,jj,jk)
396               END_2D
397            END DO
398         ENDIF
399         !
400         CALL lbc_lnk_multi( 'ldfdyn', ahmt, 'T', 1.,  ahmf, 'F', 1. )
401         !
402         !
403      CASE(  32  )       !==  time varying 3D field  ==!   = F( local deformation rate and gridscale ) (Smagorinsky)
404         !
405         IF( ln_dynldf_lap .OR. ln_dynldf_blp  ) THEN        ! laplacian operator : (C_smag/pi)^2 L^2 |D|
406            !
407            zcmsmag   = (rn_csmc/rpi)**2                                            ! (C_smag/pi)^2
408            zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 12._wp * 12._wp * zcmsmag ) ! lower limit stability factor scaling
409            zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt )               ! upper limit stability factor scaling
410            IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo          ! provide |U|L^3/12 lower limit instead
411            !                                                                       ! of |U|L^3/16 in blp case
412            DO jk = 1, jpkm1
413               !
414               DO_2D_00_00
415                  zdb =    ( uu(ji,jj,jk,Kbb) * r1_e2u(ji,jj) -  uu(ji-1,jj,jk,Kbb) * r1_e2u(ji-1,jj) )  &
416                       &                      * r1_e1t(ji,jj) * e2t(ji,jj)                           &
417                       & - ( vv(ji,jj,jk,Kbb) * r1_e1v(ji,jj) -  vv(ji,jj-1,jk,Kbb) * r1_e1v(ji,jj-1) )  &
418                       &                      * r1_e2t(ji,jj) * e1t(ji,jj)
419                  dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk)
420               END_2D
421               !
422               DO_2D_10_10
423                  zdb =   (  uu(ji,jj+1,jk,Kbb) * r1_e1u(ji,jj+1) -  uu(ji,jj,jk,Kbb) * r1_e1u(ji,jj) )  &
424                       &                        * r1_e2f(ji,jj)   * e1f(ji,jj)                       &
425                       & + ( vv(ji+1,jj,jk,Kbb) * r1_e2v(ji+1,jj) -  vv(ji,jj,jk,Kbb) * r1_e2v(ji,jj) )  &
426                       &                        * r1_e1f(ji,jj)   * e2f(ji,jj)
427                  dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk)
428               END_2D
429               !
430            END DO
431            !
432            CALL lbc_lnk_multi( 'ldfdyn', dtensq, 'T', 1. )  ! lbc_lnk on dshesq not needed
433            !
434            DO jk = 1, jpkm1
435              !
436               DO_2D_00_00
437                  !
438                  zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
439                  zu2pv2_ij_m1 = uu(ji-1,jj  ,jk,Kbb) * uu(ji-1,jj  ,jk,Kbb) + vv(ji  ,jj-1,jk,Kbb) * vv(ji  ,jj-1,jk,Kbb)
440                  !
441                  zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2
442                  ahmt(ji,jj,jk) = zdelta * SQRT(          dtensq(ji  ,jj,jk) +                         &
443                     &                            r1_4 * ( dshesq(ji  ,jj,jk) + dshesq(ji  ,jj-1,jk) +  &
444                     &                                     dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) )
445                  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
446                  ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt)
447                  !
448               END_2D
449               !
450               DO_2D_10_10
451                  !
452                  zu2pv2_ij_p1 = uu(ji  ,jj+1,jk, kbb) * uu(ji  ,jj+1,jk, kbb) + vv(ji+1,jj  ,jk, kbb) * vv(ji+1,jj  ,jk, kbb)
453                  zu2pv2_ij    = uu(ji  ,jj  ,jk, kbb) * uu(ji  ,jj  ,jk, kbb) + vv(ji  ,jj  ,jk, kbb) * vv(ji  ,jj  ,jk, kbb)
454                  !
455                  zdelta         = zcmsmag * esqf(ji,jj)                                        ! L^2 * (C_smag/pi)^2
456                  ahmf(ji,jj,jk) = zdelta * SQRT(          dshesq(ji  ,jj,jk) +                         &
457                     &                            r1_4 * ( dtensq(ji  ,jj,jk) + dtensq(ji  ,jj+1,jk) +  &
458                     &                                     dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) )
459                  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
460                  ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt)
461                  !
462               END_2D
463               !
464            END DO
465            !
466         ENDIF
467         !
468         IF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8)
469            !                          !                      = sqrt( A_lap_smag L^2/8 )
470            !                          ! stability limits already applied to laplacian values
471            !                          ! effective default limits are 1/12 |U|L^3 < B_hm < 1//(32*2dt) L^4
472            DO jk = 1, jpkm1
473               DO_2D_00_00
474                  ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) )
475               END_2D
476               DO_2D_10_10
477                  ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) )
478               END_2D
479            END DO
480            !
481         ENDIF
482         !
483         CALL lbc_lnk_multi( 'ldfdyn', ahmt, 'T', 1. , ahmf, 'F', 1. )
484         !
485      END SELECT
486      !
487      CALL iom_put( "ahmt_2d", ahmt(:,:,1) )   ! surface u-eddy diffusivity coeff.
488      CALL iom_put( "ahmf_2d", ahmf(:,:,1) )   ! surface v-eddy diffusivity coeff.
489      CALL iom_put( "ahmt_3d", ahmt(:,:,:) )   ! 3D      u-eddy diffusivity coeff.
490      CALL iom_put( "ahmf_3d", ahmf(:,:,:) )   ! 3D      v-eddy diffusivity coeff.
491      !
492      IF( ln_timing )   CALL timing_stop('ldf_dyn')
493      !
494   END SUBROUTINE ldf_dyn
495
496   !!======================================================================
497END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.