New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
ldfdyn.F90 in branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 7508

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

changes on code duplication and workshare construct
  • Property svn:keywords set to Id
File size: 17.2 KB
Line 
1MODULE ldfdyn
2   !!======================================================================
3   !!                       ***  MODULE  ldfdyn  ***
4   !! Ocean physics:  lateral viscosity coefficient
5   !!=====================================================================
6   !! History :  OPA  ! 1997-07  (G. Madec)  multi dimensional coefficients
7   !!   NEMO     1.0  ! 2002-09  (G. Madec)  F90: Free form and module
8   !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
9   !!                 !                                  add velocity dependent coefficient and optional read in file
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   ldf_dyn_init  : initialization, namelist read, and parameters control
14   !!   ldf_dyn       : update lateral eddy viscosity coefficients at each time step
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers   
17   USE dom_oce         ! ocean space and time domain
18   USE phycst          ! physical constants
19   USE 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
45   LOGICAL , PUBLIC ::   l_ldfdyn_time   !: flag for time variation of the lateral eddy viscosity coef.
46
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy diffusivity coef. at U- and V-points   [m2/s or m4/s]
48
49   REAL(wp) ::   r1_12   = 1._wp / 12._wp    ! =1/12
50   REAL(wp) ::   r1_4    = 0.25_wp           ! =1/4
51   REAL(wp) ::   r1_288  = 1._wp / 288._wp   ! =1/( 12^2 * 2 )
52
53   !! * Substitutions
54#  include "vectopt_loop_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
57   !! $Id$
58   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE ldf_dyn_init
63      !!----------------------------------------------------------------------
64      !!                  ***  ROUTINE ldf_dyn_init  ***
65      !!                   
66      !! ** Purpose :   set the horizontal ocean dynamics physics
67      !!
68      !! ** Method  :   the eddy viscosity coef. specification depends on:
69      !!              - the operator:
70      !!             ln_dynldf_lap = T     laplacian operator
71      !!             ln_dynldf_blp = T   bilaplacian operator
72      !!              - the parameter nn_ahm_ijk_t:
73      !!    nn_ahm_ijk_t  =  0 => = constant
74      !!                  = 10 => = F(z) :     = constant with a reduction of 1/4 with depth
75      !!                  =-20 => = F(i,j)     = shape read in 'eddy_viscosity.nc' file
76      !!                  = 20    = F(i,j)     = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
77      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_viscosity.nc'  file
78      !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10)
79      !!                  = 31    = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator
80      !!                                                           or |u|e^3/12 bilaplacian operator )
81      !!----------------------------------------------------------------------
82      INTEGER  ::   jk, jj, ji        ! dummy loop indices
83      INTEGER  ::   ierr, inum, ios   ! local integer
84      REAL(wp) ::   zah0              ! local scalar
85      !
86      NAMELIST/namdyn_ldf/ ln_dynldf_lap, ln_dynldf_blp,                  &
87         &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso,   &
88         &                 nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0
89      !!----------------------------------------------------------------------
90      !
91      REWIND( numnam_ref )              ! Namelist namdyn_ldf in reference namelist : Lateral physics
92      READ  ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901)
93901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist', lwp )
94
95      REWIND( numnam_cfg )              ! Namelist namdyn_ldf in configuration namelist : Lateral physics
96      READ  ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 )
97902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist', lwp )
98      IF(lwm) WRITE ( numond, namdyn_ldf )
99
100      IF(lwp) THEN                      ! Parameter print
101         WRITE(numout,*)
102         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
103         WRITE(numout,*) '~~~~~~~'
104         WRITE(numout,*) '   Namelist namdyn_ldf : set lateral mixing parameters'
105         !
106         WRITE(numout,*) '      type :'
107         WRITE(numout,*) '         laplacian operator                   ln_dynldf_lap = ', ln_dynldf_lap
108         WRITE(numout,*) '         bilaplacian operator                 ln_dynldf_blp = ', ln_dynldf_blp
109         !
110         WRITE(numout,*) '      direction of action :'
111         WRITE(numout,*) '         iso-level                            ln_dynldf_lev = ', ln_dynldf_lev
112         WRITE(numout,*) '         horizontal (geopotential)            ln_dynldf_hor = ', ln_dynldf_hor
113         WRITE(numout,*) '         iso-neutral                          ln_dynldf_iso = ', ln_dynldf_iso
114         !
115         WRITE(numout,*) '      coefficients :'
116         WRITE(numout,*) '         type of time-space variation         nn_ahm_ijk_t  = ', nn_ahm_ijk_t
117         WRITE(numout,*) '         lateral laplacian eddy viscosity     rn_ahm_0_lap  = ', rn_ahm_0, ' m2/s'
118         WRITE(numout,*) '         background viscosity (iso case)      rn_ahm_b      = ', rn_ahm_b, ' m2/s'
119         WRITE(numout,*) '         lateral bilaplacian eddy viscosity   rn_ahm_0_blp  = ', rn_bhm_0, ' m4/s'
120      ENDIF
121
122      !                                ! Parameter control
123      IF( .NOT.ln_dynldf_lap .AND. .NOT.ln_dynldf_blp ) THEN
124         IF(lwp) WRITE(numout,*) '   No viscous operator selected. ahmt and ahmf are not allocated'
125         l_ldfdyn_time = .FALSE.
126         RETURN
127      ENDIF
128      !
129      IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN     ! iso-neutral bilaplacian not implemented
130         CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' ) 
131      ENDIF
132
133      ! ... Space/Time variation of eddy coefficients
134      !                                               ! allocate the ahm arrays
135      ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr )
136      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays')
137      !
138!$OMP PARALLEL DO schedule(static) private(jj, ji)
139      DO jj = 1, jpj
140         DO ji = 1, jpi
141            ahmt(ji,jj,jpk) = 0._wp                           ! last level always 0 
142            ahmf(ji,jj,jpk) = 0._wp
143         END DO
144      END DO
145      !
146      !                                               ! value of eddy mixing coef.
147      IF    ( ln_dynldf_lap ) THEN   ;   zah0 =      rn_ahm_0         !   laplacian operator
148      ELSEIF( ln_dynldf_blp ) THEN   ;   zah0 = ABS( rn_bhm_0 )       ! bilaplacian operator
149      ELSE                                                                  ! NO viscous  operator
150         CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' )
151      ENDIF
152      !
153      l_ldfdyn_time = .FALSE.                         ! no time variation except in case defined below
154      !
155      IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN     ! only if a lateral diffusion operator is used
156         !
157         SELECT CASE(  nn_ahm_ijk_t  )                ! Specification of space time variations of ahmt, ahmf
158         !
159         CASE(   0  )      !==  constant  ==!
160            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = constant '
161!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
162      DO jk = 1, jpk
163         DO jj = 1, jpj
164            DO ji = 1, jpi
165               ahmt(ji,jj,jk) = zah0 * tmask(ji,jj,jk)
166               ahmf(ji,jj,jk) = zah0 * fmask(ji,jj,jk)
167            END DO
168         END DO
169      END DO
170            !
171         CASE(  10  )      !==  fixed profile  ==!
172            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( depth )'
173!$OMP PARALLEL DO schedule(static) private(jj, ji)
174               DO jj = 1, jpj
175                  DO ji = 1, jpi
176                     ahmt(ji,jj,1) = zah0 * tmask(ji,jj,1)                      ! constant surface value
177                     ahmf(ji,jj,1) = zah0 * fmask(ji,jj,1)
178                  END DO
179               END DO
180            CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
181            !
182         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
183            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F(i,j) read in eddy_viscosity.nc file'
184            CALL iom_open( 'eddy_viscosity_2D.nc', inum )
185            CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) )
186            CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) )
187            CALL iom_close( inum )
188!!gm Question : info for LAP or BLP case  to take into account the SQRT in the bilaplacian case ???
189!!              do we introduce a scaling by the max value of the array, and then multiply by zah0 ????
190!!              better:  check that the max is <=1  i.e. it is a shape from 0 to 1, not a coef that has physical dimension
191!$OMP PARALLEL DO schedule(static) private(jk)
192            DO jk = 2, jpkm1
193               ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk)
194               ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk)
195            END DO
196            !
197         CASE(  20  )      !== fixed horizontal shape  ==!
198            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)'
199            IF( ln_dynldf_lap )   CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor
200            IF( ln_dynldf_blp )   CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor^3
201            !
202         CASE( -30  )      !== fixed 3D shape read in file  ==!
203            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'
204            CALL iom_open( 'eddy_viscosity_3D.nc', inum )
205            CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt )
206            CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf )
207            CALL iom_close( inum )
208!!gm Question : info for LAP or BLP case  to take into account the SQRT in the bilaplacian case ????
209!!              do we introduce a scaling by the max value of the array, and then multiply by zah0 ????
210!$OMP PARALLEL DO schedule(static) private(jk)
211            DO jk = 1, jpkm1
212               ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk)
213               ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk)
214            END DO
215            !
216         CASE(  30  )       !==  fixed 3D shape  ==!
217            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth )'
218            IF( ln_dynldf_lap )   CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor
219            IF( ln_dynldf_blp )   CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor
220            !                                                    ! reduction with depth
221            CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
222            !
223         CASE(  31  )       !==  time varying 3D field  ==!
224            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth , time )'
225            IF(lwp) WRITE(numout,*) '                                proportional to the velocity : |u|e/12 or |u|e^3/12'
226            !
227            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
228            !
229         CASE DEFAULT
230            CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
231         END SELECT
232         !
233         IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN       ! bilapcian and no time variation:
234            ahmt(:,:,:) = SQRT( ahmt(:,:,:) )                     ! take the square root of the coefficient
235            ahmf(:,:,:) = SQRT( ahmf(:,:,:) )
236         ENDIF
237         !
238      ENDIF
239      !
240   END SUBROUTINE ldf_dyn_init
241
242
243   SUBROUTINE ldf_dyn( kt )
244      !!----------------------------------------------------------------------
245      !!                  ***  ROUTINE ldf_dyn  ***
246      !!
247      !! ** Purpose :   update at kt the momentum lateral mixing coeff. (ahmt and ahmf)
248      !!
249      !! ** Method  :   time varying eddy viscosity coefficients:
250      !!
251      !!    nn_ahm_ijk_t = 31    ahmt, ahmf = F(i,j,k,t) = F(local velocity)
252      !!                         ( |u|e /12  or  |u|e^3/12 for laplacian or bilaplacian operator )
253      !!                BLP case : sqrt of the eddy coef, since bilaplacian is en re-entrant laplacian
254      !!
255      !! ** action  :    ahmt, ahmf   update at each time step
256      !!----------------------------------------------------------------------
257      INTEGER, INTENT(in) ::   kt   ! time step index
258      !
259      INTEGER  ::   ji, jj, jk   ! dummy loop indices
260      REAL(wp) ::   zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax   ! local scalar
261      !!----------------------------------------------------------------------
262      !
263      IF( nn_timing == 1 )  CALL timing_start('ldf_dyn')
264      !
265      SELECT CASE(  nn_ahm_ijk_t  )       !== Eddy vicosity coefficients ==!
266      !
267      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity )
268         !
269         IF( ln_dynldf_lap   ) THEN          !   laplacian operator : |u| e /12 = |u/144| e
270            DO jk = 1, jpkm1
271               DO jj = 2, jpjm1
272                  DO ji = fs_2, fs_jpim1
273                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
274                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
275                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
276                     zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
277                     zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
278                     ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax * tmask(ji,jj,jk)      ! 288= 12*12 * 2
279                     ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax * fmask(ji,jj,jk)
280                  END DO
281               END DO
282            END DO
283         ELSEIF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e
284            DO jk = 1, jpkm1
285               DO jj = 2, jpjm1
286                  DO ji = fs_2, fs_jpim1
287                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
288                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
289                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
290                     zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
291                     zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
292                     ahmt(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax  ) * zetmax * tmask(ji,jj,jk)
293                     ahmf(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax  ) * zefmax * fmask(ji,jj,jk)
294                  END DO
295               END DO
296            END DO
297         ENDIF
298         !
299         CALL lbc_lnk( ahmt, 'T', 1. )   ;   CALL lbc_lnk( ahmf, 'F', 1. )
300         !
301      END SELECT
302      !
303      CALL iom_put( "ahmt_2d", ahmt(:,:,1) )   ! surface u-eddy diffusivity coeff.
304      CALL iom_put( "ahmf_2d", ahmf(:,:,1) )   ! surface v-eddy diffusivity coeff.
305      CALL iom_put( "ahmt_3d", ahmt(:,:,:) )   ! 3D      u-eddy diffusivity coeff.
306      CALL iom_put( "ahmf_3d", ahmf(:,:,:) )   ! 3D      v-eddy diffusivity coeff.
307      !
308      IF( nn_timing == 1 )  CALL timing_stop('ldf_dyn')
309      !
310   END SUBROUTINE ldf_dyn
311
312   !!======================================================================
313END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.