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 branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 9499

Last change on this file since 9499 was 9499, checked in by davestorkey, 6 years ago

branches/UKMO/dev_merge_2017_CICE_interface : clear SVN keywords.

File size: 30.1 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_NONE  !: 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 "vectopt_loop_substitute.h90"
76   !!----------------------------------------------------------------------
77   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
78   !! $Id$
79   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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_NONE, 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      REWIND( numnam_ref )              ! Namelist namdyn_ldf in reference namelist : Lateral physics
118      READ  ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901)
119901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist', lwp )
120
121      REWIND( numnam_cfg )              ! Namelist namdyn_ldf in configuration namelist : Lateral physics
122      READ  ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 )
123902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist', lwp )
124      IF(lwm) WRITE ( numond, namdyn_ldf )
125
126      IF(lwp) THEN                      ! Parameter print
127         WRITE(numout,*)
128         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
129         WRITE(numout,*) '~~~~~~~'
130         WRITE(numout,*) '   Namelist namdyn_ldf : set lateral mixing parameters'
131         !
132         WRITE(numout,*) '      type :'
133         WRITE(numout,*) '         no explicit diffusion                ln_dynldf_NONE= ', ln_dynldf_NONE
134         WRITE(numout,*) '         laplacian operator                   ln_dynldf_lap = ', ln_dynldf_lap
135         WRITE(numout,*) '         bilaplacian operator                 ln_dynldf_blp = ', ln_dynldf_blp
136         !
137         WRITE(numout,*) '      direction of action :'
138         WRITE(numout,*) '         iso-level                            ln_dynldf_lev = ', ln_dynldf_lev
139         WRITE(numout,*) '         horizontal (geopotential)            ln_dynldf_hor = ', ln_dynldf_hor
140         WRITE(numout,*) '         iso-neutral                          ln_dynldf_iso = ', ln_dynldf_iso
141         !
142         WRITE(numout,*) '      coefficients :'
143         WRITE(numout,*) '         type of time-space variation         nn_ahm_ijk_t  = ', nn_ahm_ijk_t
144         WRITE(numout,*) '         lateral viscous velocity  (if cst)      rn_Uv      = ', rn_Uv, ' m/s'
145         WRITE(numout,*) '         lateral viscous length    (if cst)      rn_Lv      = ', rn_Lv, ' m'
146         WRITE(numout,*) '         background viscosity (iso-lap case)     rn_ahm_b   = ', rn_ahm_b, ' m2/s'
147         !
148         WRITE(numout,*) '      Smagorinsky settings (nn_ahm_ijk_t  = 32) :'
149         WRITE(numout,*) '         Smagorinsky coefficient              rn_csmc       = ', rn_csmc
150         WRITE(numout,*) '         factor multiplier for eddy visc.'
151         WRITE(numout,*) '            lower limit (default 1.0)         rn_minfac    = ', rn_minfac
152         WRITE(numout,*) '            upper limit (default 1.0)         rn_maxfac    = ', rn_maxfac
153      ENDIF
154
155      !
156      !           !==  type of lateral operator used  ==!   (set nldf_dyn)
157      !           !=====================================!
158      !
159      nldf_dyn = np_ERROR
160      ioptio = 0
161      IF( ln_dynldf_NONE ) THEN   ;   nldf_dyn = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF
162      IF( ln_dynldf_lap  ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF
163      IF( ln_dynldf_blp  ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF
164      IF( ioptio /= 1    )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' )
165      !
166      IF(.NOT.ln_dynldf_NONE ) THEN    !==  direction ==>> type of operator  ==!
167         ioptio = 0
168         IF( ln_dynldf_lev )   ioptio = ioptio + 1
169         IF( ln_dynldf_hor )   ioptio = ioptio + 1
170         IF( ln_dynldf_iso )   ioptio = ioptio + 1
171         IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' )
172         !
173         !                             ! Set nldf_dyn, the type of lateral diffusion, from ln_dynldf_... logicals
174         ierr = 0
175         IF( ln_dynldf_lap ) THEN         ! laplacian operator
176            IF( ln_zco ) THEN                   ! z-coordinate
177               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
178               IF ( ln_dynldf_hor )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
179               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
180            ENDIF
181            IF( ln_zps ) THEN                   ! z-coordinate with partial step
182               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level              (no rotation)
183               IF ( ln_dynldf_hor )   nldf_dyn = np_lap     ! iso-level              (no rotation)
184               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
185            ENDIF
186            IF( ln_sco ) THEN                   ! s-coordinate
187               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
188               IF ( ln_dynldf_hor )   nldf_dyn = np_lap_i   ! horizontal             (   rotation)
189               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
190            ENDIF
191         ENDIF
192         !
193         IF( ln_dynldf_blp ) THEN         ! bilaplacian operator
194            IF( ln_zco ) THEN                   ! z-coordinate
195               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level = horizontal (no rotation)
196               IF( ln_dynldf_hor )   nldf_dyn = np_blp   ! iso-level = horizontal (no rotation)
197               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
198            ENDIF
199            IF( ln_zps ) THEN                   ! z-coordinate with partial step
200               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level              (no rotation)
201               IF( ln_dynldf_hor )   nldf_dyn = np_blp   ! iso-level              (no rotation)
202               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
203            ENDIF
204            IF( ln_sco ) THEN                   ! s-coordinate
205               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level              (no rotation)
206               IF( ln_dynldf_hor )   ierr = 2            ! horizontal             (   rotation)
207               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
208            ENDIF
209         ENDIF
210         !
211         IF( ierr == 2 )   CALL ctl_stop( 'rotated bi-laplacian operator does not exist' )
212         !
213         IF( nldf_dyn == np_lap_i )   l_ldfslp = .TRUE.  ! rotation require the computation of the slopes
214         !
215      ENDIF
216      !
217      IF(lwp) THEN
218         WRITE(numout,*)
219         SELECT CASE( nldf_dyn )
220         CASE( np_no_ldf )   ;   WRITE(numout,*) '   ==>>>   NO lateral viscosity'
221         CASE( np_lap    )   ;   WRITE(numout,*) '   ==>>>   iso-level laplacian operator'
222         CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background'
223         CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator'
224         END SELECT
225         WRITE(numout,*)
226      ENDIF
227     
228      !
229      !           !==  Space/time variation of eddy coefficients  ==!
230      !           !=================================================!
231      !
232      l_ldfdyn_time = .FALSE.                ! no time variation except in case defined below
233      !
234      IF( ln_dynldf_NONE ) THEN
235         IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt and ahmf are not allocated'
236         RETURN
237         !
238      ELSE                                   !==  a lateral diffusion operator is used  ==!
239         !
240         !                                         ! allocate the ahm arrays
241         ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr )
242         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays')
243         !
244         ahmt(:,:,jpk) = 0._wp                     ! last level always 0 
245         ahmf(:,:,jpk) = 0._wp
246         !
247         !                                         ! value of lap/blp eddy mixing coef.
248         IF(     ln_dynldf_lap ) THEN   ;   zUfac = r1_2 *rn_Uv   ;   inn = 1   ;   cl_Units = ' m2/s'   !   laplacian
249         ELSEIF( ln_dynldf_blp ) THEN   ;   zUfac = r1_12*rn_Uv   ;   inn = 3   ;   cl_Units = ' m4/s'   ! bilaplacian
250         ENDIF
251         zah0    = zUfac *    rn_Lv**inn              ! mixing coefficient
252         zah_max = zUfac * (ra*rad)**inn              ! maximum reachable coefficient (value at the Equator)
253         !
254         SELECT CASE(  nn_ahm_ijk_t  )             !* Specification of space-time variations of ahmt, ahmf
255         !
256         CASE(   0  )      !==  constant  ==!
257            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity. = constant = ', zah0, cl_Units
258            ahmt(:,:,1:jpkm1) = zah0
259            ahmf(:,:,1:jpkm1) = zah0
260            !
261         CASE(  10  )      !==  fixed profile  ==!
262            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( depth )'
263            IF(lwp) WRITE(numout,*) '           surface viscous coef. = constant = ', zah0, cl_Units
264            ahmt(:,:,1) = zah0                        ! constant surface value
265            ahmf(:,:,1) = zah0
266            CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
267            !
268         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
269            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j) read in eddy_viscosity.nc file'
270            CALL iom_open( 'eddy_viscosity_2D.nc', inum )
271            CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) )
272            CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) )
273            CALL iom_close( inum )
274            DO jk = 2, jpkm1
275               ahmt(:,:,jk) = ahmt(:,:,1)
276               ahmf(:,:,jk) = ahmf(:,:,1)
277            END DO
278            !
279         CASE(  20  )      !== fixed horizontal shape  ==!
280            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)'
281            IF(lwp) WRITE(numout,*) '           using a fixed viscous velocity = ', rn_Uv  ,' m/s   and   Lv = Max(e1,e2)'
282            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)'
283            CALL ldf_c2d( 'DYN', zUfac      , inn        , ahmt, ahmf )         ! surface value proportional to scale factor^inn
284            !
285         CASE( -30  )      !== fixed 3D shape read in file  ==!
286            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j,k) read in eddy_viscosity_3D.nc file'
287            CALL iom_open( 'eddy_viscosity_3D.nc', inum )
288            CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt )
289            CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf )
290            CALL iom_close( inum )
291            !
292         CASE(  30  )       !==  fixed 3D shape  ==!
293            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth )'
294            IF(lwp) WRITE(numout,*) '           using a fixed viscous velocity = ', rn_Uv  ,' m/s   and   Ld = Max(e1,e2)'
295            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)'
296            CALL ldf_c2d( 'DYN', zUfac      , inn        , ahmt, ahmf )         ! surface value proportional to scale factor^inn
297            CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )  ! reduction with depth
298            !
299         CASE(  31  )       !==  time varying 3D field  ==!
300            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth , time )'
301            IF(lwp) WRITE(numout,*) '           proportional to the local velocity : 1/2 |u|e (lap) or 1/12 |u|e^3 (blp)'
302            !
303            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
304            !
305         CASE(  32  )       !==  time varying 3D field  ==!
306            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth , time )'
307            IF(lwp) WRITE(numout,*) '           proportional to the local deformation rate and gridscale (Smagorinsky)'
308            !
309            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
310            !
311            !                          ! allocate arrays used in ldf_dyn.
312            ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) ,  esqf(jpi,jpj) , STAT=ierr )
313            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays')
314            !
315            DO jj = 2, jpjm1           ! Set local gridscale values
316               DO ji = fs_2, fs_jpim1
317                  esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2 
318                  esqf(ji,jj) = ( e1e2f(ji,jj) /( e1f(ji,jj) + e2f(ji,jj) ) )**2 
319               END DO
320            END DO
321            !
322         CASE DEFAULT
323            CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
324         END SELECT
325         !
326         IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation
327            IF(     ln_dynldf_lap ) THEN                 !   laplacian operator (mask only)
328               ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1)
329               ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1)
330            ELSEIF( ln_dynldf_blp ) THEN                 ! bilaplacian operator (square root + mask)
331               ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1)
332               ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1)
333            ENDIF
334         ENDIF
335         !
336      ENDIF
337      !
338   END SUBROUTINE ldf_dyn_init
339
340
341   SUBROUTINE ldf_dyn( kt )
342      !!----------------------------------------------------------------------
343      !!                  ***  ROUTINE ldf_dyn  ***
344      !!
345      !! ** Purpose :   update at kt the momentum lateral mixing coeff. (ahmt and ahmf)
346      !!
347      !! ** Method  :   time varying eddy viscosity coefficients:
348      !!
349      !!    nn_ahm_ijk_t = 31    ahmt, ahmf = F(i,j,k,t) = F(local velocity)
350      !!                         ( |u|e /12  or  |u|e^3/12 for laplacian or bilaplacian operator )
351      !!
352      !!    nn_ahm_ijk_t = 32    ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 
353      !!                         ( L^2|D|    or  L^4|D|/8  for laplacian or bilaplacian operator )
354      !!
355      !! ** note    :    in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian
356      !! ** action  :    ahmt, ahmf   updated at each time step
357      !!----------------------------------------------------------------------
358      INTEGER, INTENT(in) ::   kt   ! time step index
359      !
360      INTEGER  ::   ji, jj, jk   ! dummy loop indices
361      REAL(wp) ::   zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax   ! local scalar
362      REAL(wp) ::   zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb              ! local scalar
363      !!----------------------------------------------------------------------
364      !
365      IF( ln_timing )   CALL timing_start('ldf_dyn')
366      !
367      SELECT CASE(  nn_ahm_ijk_t  )       !== Eddy vicosity coefficients ==!
368      !
369      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity )
370         !
371         IF( ln_dynldf_lap   ) THEN        ! laplacian operator : |u| e /12 = |u/144| e
372            DO jk = 1, jpkm1
373               DO jj = 2, jpjm1
374                  DO ji = fs_2, fs_jpim1
375                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
376                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
377                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
378                     zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
379                     zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
380                     ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax * tmask(ji,jj,jk)      ! 288= 12*12 * 2
381                     ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax * fmask(ji,jj,jk)
382                  END DO
383               END DO
384            END DO
385         ELSEIF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e
386            DO jk = 1, jpkm1
387               DO jj = 2, jpjm1
388                  DO ji = fs_2, fs_jpim1
389                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
390                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
391                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
392                     zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
393                     zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
394                     ahmt(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax  ) * zetmax * tmask(ji,jj,jk)
395                     ahmf(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax  ) * zefmax * fmask(ji,jj,jk)
396                  END DO
397               END DO
398            END DO
399         ENDIF
400         !
401         CALL lbc_lnk_multi( ahmt, 'T', 1.,  ahmf, 'F', 1. )
402         !
403         !
404      CASE(  32  )       !==  time varying 3D field  ==!   = F( local deformation rate and gridscale ) (Smagorinsky)
405         !
406         IF( ln_dynldf_lap .OR. ln_dynldf_blp  ) THEN        ! laplacian operator : (C_smag/pi)^2 L^2 |D|
407            !
408            zcmsmag = (rn_csmc/rpi)**2                                              ! (C_smag/pi)^2
409            zstabf_lo  = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag )        ! lower limit stability factor scaling
410            zstabf_up  = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt )              ! upper limit stability factor scaling
411            IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo          ! provide |U|L^3/12 lower limit instead
412            !                                                                       ! of |U|L^3/16 in blp case
413            DO jk = 1, jpkm1
414               !
415               DO jj = 2, jpj
416                  DO ji = 2, jpi
417                     zdb = ( (  ub(ji,jj,jk) * r1_e2u(ji,jj) -  ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) )  &
418                          &                  * r1_e1t(ji,jj) * e2t(ji,jj)                           &
419                          & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) -  vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) )  &
420                          &                  * r1_e2t(ji,jj) * e1t(ji,jj)    ) * tmask(ji,jj,jk)
421                     dtensq(ji,jj) = zdb * zdb
422                  END DO
423               END DO
424               !
425               DO jj = 1, jpjm1
426                  DO ji = 1, jpim1
427                     zdb = ( (  ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) -  ub(ji,jj,jk) * r1_e1u(ji,jj) )  &
428                          &                    * r1_e2f(ji,jj)   * e1f(ji,jj)                       &
429                          & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) -  vb(ji,jj,jk) * r1_e2v(ji,jj) )  &
430                          &                    * r1_e1f(ji,jj)   * e2f(ji,jj)  ) * fmask(ji,jj,jk)
431                     dshesq(ji,jj) = zdb * zdb
432                  END DO
433               END DO
434               !
435               DO jj = 2, jpjm1
436                  DO ji = fs_2, fs_jpim1
437                     !
438                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
439                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
440                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
441                                                     ! T-point value
442                     zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2
443                     ahmt(ji,jj,jk) = zdelta * sqrt(          dtensq(ji,jj)   +                        &
444                                     &               r1_4 * ( dshesq(ji,jj)   + dshesq(ji,jj-1)   +    &
445                                     &                        dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) )
446                     ahmt(ji,jj,jk) =   MAX( ahmt(ji,jj,jk),   &
447                                     &   SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2
448                     ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up )                    ! Impose upper limit == maxfac  * L^2/(4*2dt)
449                                                     ! F-point value
450                     zdelta         = zcmsmag * esqf(ji,jj)                                        ! L^2 * (C_smag/pi)^2
451                     ahmf(ji,jj,jk) = zdelta * sqrt(          dshesq(ji,jj)   +                        &
452                                     &               r1_4 * ( dtensq(ji,jj)   + dtensq(ji,jj+1)   +    &
453                                     &                        dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) )
454                     ahmf(ji,jj,jk) =   MAX( ahmf(ji,jj,jk),   &
455                                     &   SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2
456                     ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up )                    ! Impose upper limit == maxfac  * L^2/(4*2dt)
457                     !
458                  END DO
459               END DO
460            END DO
461            !
462         ENDIF
463         !
464         IF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8)
465            !                          !                      = sqrt( A_lap_smag L^2/8 )
466            !                          ! stability limits already applied to laplacian values
467            !                          ! effective default limits are 1/12 |U|L^3 < B_hm < 1//(32*2dt) L^4
468            DO jk = 1, jpkm1
469               DO jj = 2, jpjm1
470                  DO ji = fs_2, fs_jpim1
471                     ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) )
472                     ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) )
473                  END DO
474               END DO
475            END DO
476            !
477         ENDIF
478         !
479         CALL lbc_lnk_multi( ahmt, 'T', 1. , ahmf, 'F', 1. )
480         !
481      END SELECT
482      !
483      CALL iom_put( "ahmt_2d", ahmt(:,:,1) )   ! surface u-eddy diffusivity coeff.
484      CALL iom_put( "ahmf_2d", ahmf(:,:,1) )   ! surface v-eddy diffusivity coeff.
485      CALL iom_put( "ahmt_3d", ahmt(:,:,:) )   ! 3D      u-eddy diffusivity coeff.
486      CALL iom_put( "ahmf_3d", ahmf(:,:,:) )   ! 3D      v-eddy diffusivity coeff.
487      !
488      IF( ln_timing )   CALL timing_stop('ldf_dyn')
489      !
490   END SUBROUTINE ldf_dyn
491
492   !!======================================================================
493END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.