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 trunk/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
  • Property svn:keywords set to Id
File size: 26.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 ldfc1d_c2d      ! lateral diffusion: 1D and 2D cases
20   !
21   USE in_out_manager  ! I/O manager
22   USE iom             ! I/O module for ehanced bottom friction file
23   USE timing          ! Timing
24   USE lib_mpp         ! distribued memory computing library
25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
26   USE wrk_nemo        ! Memory Allocation
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_lap   !: laplacian operator
36   LOGICAL , PUBLIC ::   ln_dynldf_blp   !: bilaplacian operator
37   LOGICAL , PUBLIC ::   ln_dynldf_lev   !: iso-level direction
38   LOGICAL , PUBLIC ::   ln_dynldf_hor   !: horizontal (geopotential) direction
39   LOGICAL , PUBLIC ::   ln_dynldf_iso   !: iso-neutral direction
40   INTEGER , PUBLIC ::   nn_ahm_ijk_t    !: choice of time & space variations of the lateral eddy viscosity coef.
41   REAL(wp), PUBLIC ::   rn_ahm_0        !: lateral laplacian eddy viscosity            [m2/s]
42   REAL(wp), PUBLIC ::   rn_ahm_b        !: lateral laplacian background eddy viscosity [m2/s]
43   REAL(wp), PUBLIC ::   rn_bhm_0        !: lateral bilaplacian eddy viscosity          [m4/s]
44                                         !! If nn_ahm_ijk_t = 32 a time and space varying Smagorinsky viscosity
45                                         !! will be computed.
46   REAL(wp), PUBLIC ::   rn_csmc         !: Smagorinsky constant of proportionality
47   REAL(wp), PUBLIC ::   rn_minfac       !: Multiplicative factor of theorectical minimum Smagorinsky viscosity
48   REAL(wp), PUBLIC ::   rn_maxfac       !: Multiplicative factor of theorectical maximum Smagorinsky viscosity
49
50   LOGICAL , PUBLIC ::   l_ldfdyn_time   !: flag for time variation of the lateral eddy viscosity coef.
51
52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy diffusivity coef. at U- and V-points   [m2/s or m4/s]
53   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dtensq       !: horizontal tension squared         (Smagorinsky only)
54   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dshesq       !: horizontal shearing strain squared (Smagorinsky only)
55   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   esqt, esqf   !: Square of the local gridscale (e1e2/(e1+e2))**2           
56
57   REAL(wp) ::   r1_12   = 1._wp / 12._wp    ! =1/12
58   REAL(wp) ::   r1_4    = 0.25_wp           ! =1/4
59   REAL(wp) ::   r1_8    = 0.125_wp          ! =1/8
60   REAL(wp) ::   r1_288  = 1._wp / 288._wp   ! =1/( 12^2 * 2 )
61
62   !! * Substitutions
63#  include "vectopt_loop_substitute.h90"
64   !!----------------------------------------------------------------------
65   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
66   !! $Id$
67   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
68   !!----------------------------------------------------------------------
69CONTAINS
70
71   SUBROUTINE ldf_dyn_init
72      !!----------------------------------------------------------------------
73      !!                  ***  ROUTINE ldf_dyn_init  ***
74      !!                   
75      !! ** Purpose :   set the horizontal ocean dynamics physics
76      !!
77      !! ** Method  :   the eddy viscosity coef. specification depends on:
78      !!              - the operator:
79      !!             ln_dynldf_lap = T     laplacian operator
80      !!             ln_dynldf_blp = T   bilaplacian operator
81      !!              - the parameter nn_ahm_ijk_t:
82      !!    nn_ahm_ijk_t  =  0 => = constant
83      !!                  = 10 => = F(z) :     = constant with a reduction of 1/4 with depth
84      !!                  =-20 => = F(i,j)     = shape read in 'eddy_viscosity.nc' file
85      !!                  = 20    = F(i,j)     = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
86      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_viscosity.nc'  file
87      !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10)
88      !!                  = 31    = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator
89      !!                                                           or |u|e^3/12 bilaplacian operator )
90      !!                  = 32    = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 
91      !!                                                           (   L^2|D|      laplacian operator
92      !!                                                           or  L^4|D|/8  bilaplacian operator )
93      !!----------------------------------------------------------------------
94      INTEGER  ::   ji, jj, jk        ! dummy loop indices
95      INTEGER  ::   ierr, inum, ios   ! local integer
96      REAL(wp) ::   zah0              ! local scalar
97      !
98      NAMELIST/namdyn_ldf/ ln_dynldf_lap, ln_dynldf_blp,                  &
99         &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso,   &
100         &                 nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0,   &
101         &                 rn_csmc      , rn_minfac, rn_maxfac
102      !!----------------------------------------------------------------------
103      !
104      REWIND( numnam_ref )              ! Namelist namdyn_ldf in reference namelist : Lateral physics
105      READ  ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901)
106901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist', lwp )
107
108      REWIND( numnam_cfg )              ! Namelist namdyn_ldf in configuration namelist : Lateral physics
109      READ  ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 )
110902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist', lwp )
111      IF(lwm) WRITE ( numond, namdyn_ldf )
112
113      IF(lwp) THEN                      ! Parameter print
114         WRITE(numout,*)
115         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
116         WRITE(numout,*) '~~~~~~~'
117         WRITE(numout,*) '   Namelist namdyn_ldf : set lateral mixing parameters'
118         !
119         WRITE(numout,*) '      type :'
120         WRITE(numout,*) '         laplacian operator                   ln_dynldf_lap = ', ln_dynldf_lap
121         WRITE(numout,*) '         bilaplacian operator                 ln_dynldf_blp = ', ln_dynldf_blp
122         !
123         WRITE(numout,*) '      direction of action :'
124         WRITE(numout,*) '         iso-level                            ln_dynldf_lev = ', ln_dynldf_lev
125         WRITE(numout,*) '         horizontal (geopotential)            ln_dynldf_hor = ', ln_dynldf_hor
126         WRITE(numout,*) '         iso-neutral                          ln_dynldf_iso = ', ln_dynldf_iso
127         !
128         WRITE(numout,*) '      coefficients :'
129         WRITE(numout,*) '         type of time-space variation         nn_ahm_ijk_t  = ', nn_ahm_ijk_t
130         WRITE(numout,*) '         lateral laplacian eddy viscosity     rn_ahm_0      = ', rn_ahm_0, ' m2/s'
131         WRITE(numout,*) '         background viscosity (iso case)      rn_ahm_b      = ', rn_ahm_b, ' m2/s'
132         WRITE(numout,*) '         lateral bilaplacian eddy viscosity   rn_bhm_0      = ', rn_bhm_0, ' m4/s'
133         WRITE(numout,*) '      smagorinsky settings (nn_ahm_ijk_t  = 32) :'
134         WRITE(numout,*) '         Smagorinsky coefficient              rn_csmc       = ', rn_csmc
135         WRITE(numout,*) '         factor multiplier for theorectical lower limit for '
136         WRITE(numout,*) '           Smagorinsky eddy visc (def. 1.0)   rn_minfac    = ', rn_minfac
137         WRITE(numout,*) '         factor multiplier for theorectical lower upper for '
138         WRITE(numout,*) '           Smagorinsky eddy visc (def. 1.0)   rn_maxfac    = ', rn_maxfac
139      ENDIF
140
141      !                                ! Parameter control
142      IF( .NOT.ln_dynldf_lap .AND. .NOT.ln_dynldf_blp ) THEN
143         IF(lwp) WRITE(numout,*) '   No viscous operator selected. ahmt and ahmf are not allocated'
144         l_ldfdyn_time = .FALSE.
145         RETURN
146      ENDIF
147      !
148      IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN     ! iso-neutral bilaplacian not implemented
149         CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' ) 
150      ENDIF
151
152      ! ... Space/Time variation of eddy coefficients
153      !                                               ! allocate the ahm arrays
154      ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr )
155      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays')
156      !
157!$OMP PARALLEL DO schedule(static) private(jj, ji)
158      DO jj = 1, jpj
159         DO ji = 1, jpi
160            ahmt(ji,jj,jpk) = 0._wp                           ! last level always 0 
161            ahmf(ji,jj,jpk) = 0._wp
162         END DO
163      END DO
164      !
165      !                                               ! value of eddy mixing coef.
166      IF    ( ln_dynldf_lap ) THEN   ;   zah0 =      rn_ahm_0         !   laplacian operator
167      ELSEIF( ln_dynldf_blp ) THEN   ;   zah0 = ABS( rn_bhm_0 )       ! bilaplacian operator
168      ELSE                                                                  ! NO viscous  operator
169         CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' )
170      ENDIF
171      !
172      l_ldfdyn_time = .FALSE.                         ! no time variation except in case defined below
173      !
174      IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN     ! only if a lateral diffusion operator is used
175         !
176         SELECT CASE(  nn_ahm_ijk_t  )                ! Specification of space time variations of ahmt, ahmf
177         !
178         CASE(   0  )      !==  constant  ==!
179            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = constant '
180!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
181            DO jk = 1, jpk
182               DO jj = 1, jpj
183                  DO ji = 1, jpi
184                     ahmt(ji,jj,jk) = zah0 * tmask(ji,jj,jk)
185                     ahmf(ji,jj,jk) = zah0 * fmask(ji,jj,jk)
186                  END DO
187               END DO
188            END DO
189            !
190         CASE(  10  )      !==  fixed profile  ==!
191            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( depth )'
192!$OMP PARALLEL DO schedule(static) private(jj, ji)
193            DO jj = 1, jpj
194               DO ji = 1, jpi
195                  ahmt(ji,jj,1) = zah0 * tmask(ji,jj,1)            ! constant surface value
196                  ahmf(ji,jj,1) = zah0 * fmask(ji,jj,1)
197               END DO
198            END DO
199            CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
200            !
201         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
202            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F(i,j) read in eddy_viscosity.nc file'
203            CALL iom_open( 'eddy_viscosity_2D.nc', inum )
204            CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) )
205            CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) )
206            CALL iom_close( inum )
207!!gm Question : info for LAP or BLP case  to take into account the SQRT in the bilaplacian case ???
208!!              do we introduce a scaling by the max value of the array, and then multiply by zah0 ????
209!!              better:  check that the max is <=1  i.e. it is a shape from 0 to 1, not a coef that has physical dimension
210!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
211            DO jk = 2, jpkm1
212               DO jj = 1, jpj
213                  DO ji = 1, jpi
214                     ahmt(ji,jj,jk) = ahmt(ji,jj,1) * tmask(ji,jj,jk)
215                     ahmf(ji,jj,jk) = ahmf(ji,jj,1) * fmask(ji,jj,jk)
216                  END DO
217               END DO
218            END DO
219            !
220         CASE(  20  )      !== fixed horizontal shape  ==!
221            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)'
222            IF( ln_dynldf_lap )   CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor
223            IF( ln_dynldf_blp )   CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor^3
224            !
225         CASE( -30  )      !== fixed 3D shape read in file  ==!
226            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'
227            CALL iom_open( 'eddy_viscosity_3D.nc', inum )
228            CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt )
229            CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf )
230            CALL iom_close( inum )
231!!gm Question : info for LAP or BLP case  to take into account the SQRT in the bilaplacian case ????
232!!              do we introduce a scaling by the max value of the array, and then multiply by zah0 ????
233!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
234            DO jk = 1, jpkm1
235               DO jj = 1, jpj
236                  DO ji = 1, jpi
237                     ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * tmask(ji,jj,jk)
238                     ahmf(ji,jj,jk) = ahmf(ji,jj,jk) * fmask(ji,jj,jk)
239                  END DO
240               END DO
241            END DO
242            !
243         CASE(  30  )       !==  fixed 3D shape  ==!
244            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth )'
245            IF( ln_dynldf_lap )   CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor
246            IF( ln_dynldf_blp )   CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor
247            !                                                    ! reduction with depth
248            CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
249            !
250         CASE(  31  )       !==  time varying 3D field  ==!
251            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth , time )'
252            IF(lwp) WRITE(numout,*) '                                proportional to the velocity : |u|e/12 or |u|e^3/12'
253            !
254            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
255            !
256         CASE(  32  )       !==  time varying 3D field  ==!
257            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth , time )'
258            IF(lwp) WRITE(numout,*) '             proportional to the local deformation rate and gridscale (Smagorinsky)'
259            IF(lwp) WRITE(numout,*) '                                                             : L^2|D| or L^4|D|/8'
260            !
261            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
262            !
263            ! allocate arrays used in ldf_dyn.
264            ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) ,  esqf(jpi,jpj) , STAT=ierr )
265            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays')
266            !
267            ! Set local gridscale values
268!$OMP PARALLEL DO schedule(static) private(jj,ji)
269            DO jj = 2, jpjm1
270               DO ji = fs_2, fs_jpim1
271                  esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2 
272                  esqf(ji,jj) = ( e1e2f(ji,jj) /( e1f(ji,jj) + e2f(ji,jj) ) )**2 
273               END DO
274            END DO
275            !
276         CASE DEFAULT
277            CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
278         END SELECT
279         !
280         IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN       ! bilapcian and no time variation:
281!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
282            DO jk = 1, jpk
283               DO jj = 1, jpj
284                  DO ji = 1, jpi
285                     ahmt(ji,jj,jk) = SQRT( ahmt(ji,jj,jk) )      ! take the square root of the coefficient
286                     ahmf(ji,jj,jk) = SQRT( ahmf(ji,jj,jk) )
287                  END DO
288               END DO
289            END DO
290         ENDIF
291         !
292      ENDIF
293      !
294   END SUBROUTINE ldf_dyn_init
295
296
297   SUBROUTINE ldf_dyn( kt )
298      !!----------------------------------------------------------------------
299      !!                  ***  ROUTINE ldf_dyn  ***
300      !!
301      !! ** Purpose :   update at kt the momentum lateral mixing coeff. (ahmt and ahmf)
302      !!
303      !! ** Method  :   time varying eddy viscosity coefficients:
304      !!
305      !!    nn_ahm_ijk_t = 31    ahmt, ahmf = F(i,j,k,t) = F(local velocity)
306      !!                         ( |u|e /12  or  |u|e^3/12 for laplacian or bilaplacian operator )
307      !!
308      !!    nn_ahm_ijk_t = 32    ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 
309      !!                         ( L^2|D|    or  L^4|D|/8  for laplacian or bilaplacian operator )
310      !!
311      !! ** note    :    in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian
312      !! ** action  :    ahmt, ahmf   updated at each time step
313      !!----------------------------------------------------------------------
314      INTEGER, INTENT(in) ::   kt   ! time step index
315      !
316      INTEGER  ::   ji, jj, jk   ! dummy loop indices
317      REAL(wp) ::   zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax   ! local scalar
318      REAL(wp) ::   zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb              ! local scalar
319      !!----------------------------------------------------------------------
320      !
321      IF( nn_timing == 1 )  CALL timing_start('ldf_dyn')
322      !
323      SELECT CASE(  nn_ahm_ijk_t  )       !== Eddy vicosity coefficients ==!
324      !
325      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity )
326         !
327         IF( ln_dynldf_lap   ) THEN        ! laplacian operator : |u| e /12 = |u/144| e
328            DO jk = 1, jpkm1
329               DO jj = 2, jpjm1
330                  DO ji = fs_2, fs_jpim1
331                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
332                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
333                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
334                     zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
335                     zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
336                     ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax * tmask(ji,jj,jk)      ! 288= 12*12 * 2
337                     ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax * fmask(ji,jj,jk)
338                  END DO
339               END DO
340            END DO
341         ELSEIF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e
342            DO jk = 1, jpkm1
343               DO jj = 2, jpjm1
344                  DO ji = fs_2, fs_jpim1
345                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
346                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
347                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
348                     zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
349                     zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
350                     ahmt(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax  ) * zetmax * tmask(ji,jj,jk)
351                     ahmf(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax  ) * zefmax * fmask(ji,jj,jk)
352                  END DO
353               END DO
354            END DO
355         ENDIF
356         !
357         CALL lbc_lnk( ahmt, 'T', 1. )   ;   CALL lbc_lnk( ahmf, 'F', 1. )
358         !
359         !
360      CASE(  32  )       !==  time varying 3D field  ==!   = F( local deformation rate and gridscale ) (Smagorinsky)
361         !
362         IF( ln_dynldf_lap .OR. ln_dynldf_blp  ) THEN        ! laplacian operator : (C_smag/pi)^2 L^2 |D|
363            !
364            zcmsmag = (rn_csmc/rpi)**2                                              ! (C_smag/pi)^2
365            zstabf_lo  = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag )        ! lower limit stability factor scaling
366            zstabf_up  = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt )              ! upper limit stability factor scaling
367            IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo          ! provide |U|L^3/12 lower limit instead
368            !                                                                       ! of |U|L^3/16 in blp case
369            DO jk = 1, jpkm1
370               !
371               DO jj = 2, jpj
372                  DO ji = 2, jpi
373                     zdb = ( (  ub(ji,jj,jk) * r1_e2u(ji,jj) -  ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) )  &
374                          &                  * r1_e1t(ji,jj) * e2t(ji,jj)                           &
375                          & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) -  vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) )  &
376                          &                  * r1_e2t(ji,jj) * e1t(ji,jj)    ) * tmask(ji,jj,jk)
377                     dtensq(ji,jj) = zdb*zdb
378                  END DO
379               END DO
380               !
381               DO jj = 1, jpjm1
382                  DO ji = 1, jpim1
383                     zdb = ( (  ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) -  ub(ji,jj,jk) * r1_e1u(ji,jj) )  &
384                          &                    * r1_e2f(ji,jj)   * e1f(ji,jj)                       &
385                          & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) -  vb(ji,jj,jk) * r1_e2v(ji,jj) )  &
386                          &                    * r1_e1f(ji,jj)   * e2f(ji,jj)  ) * fmask(ji,jj,jk)
387                     dshesq(ji,jj) = zdb*zdb
388                  END DO
389               END DO
390               !
391               DO jj = 2, jpjm1
392                  DO ji = fs_2, fs_jpim1
393                     !
394                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
395                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
396                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
397                                                     ! T-point value
398                     zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2
399                     ahmt(ji,jj,jk) = zdelta * sqrt(          dtensq(ji,jj)   +                        &
400                                     &               r1_4 * ( dshesq(ji,jj)   + dshesq(ji,jj-1)   +    &
401                                     &                        dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) )
402                     ahmt(ji,jj,jk) =   MAX( ahmt(ji,jj,jk),   &
403                                     &   SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2
404                     ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up )                    ! Impose upper limit == maxfac  * L^2/(4*2dt)
405                                                     ! F-point value
406                     zdelta         = zcmsmag * esqf(ji,jj)                                        ! L^2 * (C_smag/pi)^2
407                     ahmf(ji,jj,jk) = zdelta * sqrt(          dshesq(ji,jj)   +                        &
408                                     &               r1_4 * ( dtensq(ji,jj)   + dtensq(ji,jj+1)   +    &
409                                     &                        dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) )
410                     ahmf(ji,jj,jk) =   MAX( ahmf(ji,jj,jk),   &
411                                     &   SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2
412                     ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up )                    ! Impose upper limit == maxfac  * L^2/(4*2dt)
413                     !
414                  END DO
415               END DO
416            END DO
417            !
418         ENDIF
419         !
420         IF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8)
421                                       !                      = sqrt( A_lap_smag L^2/8 )
422                                       ! stability limits already applied to laplacian values
423                                       ! effective default limits are |U|L^3/12 < B_hm < L^4/(32*2dt)
424            !
425            DO jk = 1, jpkm1
426               !
427               DO jj = 2, jpjm1
428                  DO ji = fs_2, fs_jpim1
429                     !
430                     ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) )
431                     ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) )
432                     !
433                  END DO
434               END DO
435            END DO
436            !
437         ENDIF
438         !
439         CALL lbc_lnk( ahmt, 'T', 1. )   ;   CALL lbc_lnk( ahmf, 'F', 1. )
440         !
441      END SELECT
442      !
443      CALL iom_put( "ahmt_2d", ahmt(:,:,1) )   ! surface u-eddy diffusivity coeff.
444      CALL iom_put( "ahmf_2d", ahmf(:,:,1) )   ! surface v-eddy diffusivity coeff.
445      CALL iom_put( "ahmt_3d", ahmt(:,:,:) )   ! 3D      u-eddy diffusivity coeff.
446      CALL iom_put( "ahmf_3d", ahmf(:,:,:) )   ! 3D      v-eddy diffusivity coeff.
447      !
448      IF( nn_timing == 1 )  CALL timing_stop('ldf_dyn')
449      !
450   END SUBROUTINE ldf_dyn
451
452   !!======================================================================
453END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.