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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/LDF/ldfdyn.F90 @ 12340

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

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

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