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/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/LDF – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/LDF/ldfdyn.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
  • Property svn:keywords set to Id
File size: 30.6 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 "vectopt_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 jj = 1, jpj             ! Set local gridscale values
314               DO ji = 1, jpi
315                  esqt(ji,jj) = ( e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 
316                  esqf(ji,jj) = ( e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 
317               END DO
318            END DO
319            !
320         CASE DEFAULT
321            CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
322         END SELECT
323         !
324         IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation
325            IF(     ln_dynldf_lap ) THEN                 !   laplacian operator (mask only)
326               ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1)
327               ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1)
328            ELSEIF( ln_dynldf_blp ) THEN                 ! bilaplacian operator (square root + mask)
329               ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1)
330               ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1)
331            ENDIF
332         ENDIF
333         !
334      ENDIF
335      !
336   END SUBROUTINE ldf_dyn_init
337
338
339   SUBROUTINE ldf_dyn( kt )
340      !!----------------------------------------------------------------------
341      !!                  ***  ROUTINE ldf_dyn  ***
342      !!
343      !! ** Purpose :   update at kt the momentum lateral mixing coeff. (ahmt and ahmf)
344      !!
345      !! ** Method  :   time varying eddy viscosity coefficients:
346      !!
347      !!    nn_ahm_ijk_t = 31    ahmt, ahmf = F(i,j,k,t) = F(local velocity)
348      !!                         ( |u|e /12  or  |u|e^3/12 for laplacian or bilaplacian operator )
349      !!
350      !!    nn_ahm_ijk_t = 32    ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 
351      !!                         ( L^2|D|    or  L^4|D|/8  for laplacian or bilaplacian operator )
352      !!
353      !! ** note    :    in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian
354      !! ** action  :    ahmt, ahmf   updated at each time step
355      !!----------------------------------------------------------------------
356      INTEGER, INTENT(in) ::   kt   ! time step index
357      !
358      INTEGER  ::   ji, jj, jk   ! dummy loop indices
359      REAL(wp) ::   zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zemax   ! local scalar (option 31)
360      REAL(wp) ::   zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb     ! local scalar (option 32)
361      !!----------------------------------------------------------------------
362      !
363      IF( ln_timing )   CALL timing_start('ldf_dyn')
364      !
365      SELECT CASE(  nn_ahm_ijk_t  )       !== Eddy vicosity coefficients ==!
366      !
367      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity )
368         !
369         IF( ln_dynldf_lap   ) THEN        ! laplacian operator : |u| e /12 = |u/144| e
370            DO jk = 1, jpkm1
371               DO jj = 2, jpjm1
372                  DO ji = fs_2, fs_jpim1
373                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
374                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
375                     zemax = MAX( e1t(ji,jj) , e2t(ji,jj) )
376                     ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax * tmask(ji,jj,jk)      ! 288= 12*12 * 2
377                  END DO
378               END DO
379               DO jj = 1, jpjm1
380                  DO ji = 1, fs_jpim1
381                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
382                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
383                     zemax = MAX( e1f(ji,jj) , e2f(ji,jj) )
384                     ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax * fmask(ji,jj,jk)      ! 288= 12*12 * 2
385                  END DO
386               END DO
387            END DO
388         ELSEIF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e
389            DO jk = 1, jpkm1
390               DO jj = 2, jpjm1
391                  DO ji = fs_2, fs_jpim1
392                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
393                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
394                     zemax = MAX( e1t(ji,jj) , e2t(ji,jj) )
395                     ahmt(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax  ) * zemax * tmask(ji,jj,jk)
396                  END DO
397               END DO
398               DO jj = 1, jpjm1
399                  DO ji = 1, fs_jpim1
400                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
401                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
402                     zemax = MAX( e1f(ji,jj) , e2f(ji,jj) )
403                     ahmf(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax  ) * zemax * fmask(ji,jj,jk)
404                  END DO
405               END DO
406            END DO
407         ENDIF
408         !
409         CALL lbc_lnk_multi( 'ldfdyn', ahmt, 'T', 1.,  ahmf, 'F', 1. )
410         !
411         !
412      CASE(  32  )       !==  time varying 3D field  ==!   = F( local deformation rate and gridscale ) (Smagorinsky)
413         !
414         IF( ln_dynldf_lap .OR. ln_dynldf_blp  ) THEN        ! laplacian operator : (C_smag/pi)^2 L^2 |D|
415            !
416            zcmsmag   = (rn_csmc/rpi)**2                                            ! (C_smag/pi)^2
417            zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag )         ! lower limit stability factor scaling
418            zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt )               ! upper limit stability factor scaling
419            IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo          ! provide |U|L^3/12 lower limit instead
420            !                                                                       ! of |U|L^3/16 in blp case
421            DO jk = 1, jpkm1
422               !
423               DO jj = 2, jpjm1
424                  DO ji = 2, jpim1
425                     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)  &
426                        &  - ( 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)
427                     dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk)
428                  END DO
429               END DO
430               !
431               DO jj = 1, jpjm1
432                  DO ji = 1, jpim1
433                     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)  &
434                        &  + ( 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)
435                     dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk)
436                  END DO
437               END DO
438               !
439            END DO
440            !
441            CALL lbc_lnk_multi( 'ldfdyn', dtensq, 'T', 1. )  ! lbc_lnk on dshesq not needed
442            !
443            DO jk = 1, jpkm1
444              !
445               DO jj = 2, jpjm1                                ! T-point value
446                  DO ji = fs_2, fs_jpim1
447                     !
448                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
449                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
450                     !
451                     zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2
452                     ahmt(ji,jj,jk) = zdelta * SQRT(          dtensq(ji  ,jj,jk) +                         &
453                        &                            r1_4 * ( dshesq(ji  ,jj,jk) + dshesq(ji  ,jj-1,jk) +  &
454                        &                                     dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) )
455                     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
456                     ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt)
457                     !
458                  END DO
459               END DO
460               !
461               DO jj = 1, jpjm1                                ! F-point value
462                  DO ji = 1, fs_jpim1
463                     !
464                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
465                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
466                     !
467                     zdelta         = zcmsmag * esqf(ji,jj)                                        ! L^2 * (C_smag/pi)^2
468                     ahmf(ji,jj,jk) = zdelta * SQRT(          dshesq(ji  ,jj,jk) +                         &
469                        &                            r1_4 * ( dtensq(ji  ,jj,jk) + dtensq(ji  ,jj+1,jk) +  &
470                        &                                     dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) )
471                     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
472                     ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt)
473                     !
474                  END DO
475               END DO
476               !
477            END DO
478            !
479         ENDIF
480         !
481         IF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8)
482            !                          !                      = sqrt( A_lap_smag L^2/8 )
483            !                          ! stability limits already applied to laplacian values
484            !                          ! effective default limits are 1/12 |U|L^3 < B_hm < 1//(32*2dt) L^4
485            DO jk = 1, jpkm1
486               DO jj = 2, jpjm1
487                  DO ji = fs_2, fs_jpim1
488                     ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) )
489                  END DO
490               END DO
491               DO jj = 1, jpjm1
492                  DO ji = 1, fs_jpim1
493                     ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) )
494                  END DO
495               END DO
496            END DO
497            !
498         ENDIF
499         !
500         CALL lbc_lnk_multi( 'ldfdyn', ahmt, 'T', 1. , ahmf, 'F', 1. )
501         !
502      END SELECT
503      !
504      CALL iom_put( "ahmt_2d", ahmt(:,:,1) )   ! surface u-eddy diffusivity coeff.
505      CALL iom_put( "ahmf_2d", ahmf(:,:,1) )   ! surface v-eddy diffusivity coeff.
506      CALL iom_put( "ahmt_3d", ahmt(:,:,:) )   ! 3D      u-eddy diffusivity coeff.
507      CALL iom_put( "ahmf_3d", ahmf(:,:,:) )   ! 3D      v-eddy diffusivity coeff.
508      !
509      IF( ln_timing )   CALL timing_stop('ldf_dyn')
510      !
511   END SUBROUTINE ldf_dyn
512
513   !!======================================================================
514END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.