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

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

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