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

source: branches/UKMO/dev_r6393_CO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 7019

Last change on this file since 7019 was 7019, checked in by deazer, 8 years ago

Cleared svn keywords
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 "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                ! 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      ahmt(:,:,jpk) = 0._wp                           ! last level always 0 
139      ahmf(:,:,jpk) = 0._wp
140      !
141      !                                               ! value of eddy mixing coef.
142      IF    ( ln_dynldf_lap ) THEN   ;   zah0 =      rn_ahm_0         !   laplacian operator
143      ELSEIF( ln_dynldf_blp ) THEN   ;   zah0 = ABS( rn_bhm_0 )       ! bilaplacian operator
144      ELSE                                                                  ! NO viscous  operator
145         CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' )
146      ENDIF
147      !
148      l_ldfdyn_time = .FALSE.                         ! no time variation except in case defined below
149      !
150      IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN     ! only if a lateral diffusion operator is used
151         !
152         SELECT CASE(  nn_ahm_ijk_t  )                ! Specification of space time variations of ahmt, ahmf
153         !
154         CASE(   0  )      !==  constant  ==!
155            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = constant '
156            ahmt(:,:,:) = zah0 * tmask(:,:,:)
157            ahmf(:,:,:) = zah0 * fmask(:,:,:)
158            !
159         CASE(  10  )      !==  fixed profile  ==!
160            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( depth )'
161            ahmt(:,:,1) = zah0 * tmask(:,:,1)                      ! constant surface value
162            ahmf(:,:,1) = zah0 * fmask(:,:,1)
163            CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
164            !
165         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
166            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F(i,j) read in eddy_viscosity.nc file'
167            CALL iom_open( 'eddy_viscosity_2D.nc', inum )
168            CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) )
169            CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) )
170            CALL iom_close( inum )
171!!gm Question : info for LAP or BLP case  to take into account the SQRT in the bilaplacian case ???
172!!              do we introduce a scaling by the max value of the array, and then multiply by zah0 ????
173!!              better:  check that the max is <=1  i.e. it is a shape from 0 to 1, not a coef that has physical dimension
174            DO jk = 2, jpkm1
175               ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk)
176               ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk)
177            END DO
178            !
179         CASE(  20  )      !== fixed horizontal shape  ==!
180            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)'
181            IF( ln_dynldf_lap )   CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor
182            IF( ln_dynldf_blp )   CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor^3
183            !
184         CASE( -30  )      !== fixed 3D shape read in file  ==!
185            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'
186            CALL iom_open( 'eddy_viscosity_3D.nc', inum )
187            CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt )
188            CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf )
189            CALL iom_close( inum )
190!!gm Question : info for LAP or BLP case  to take into account the SQRT in the bilaplacian case ????
191!!              do we introduce a scaling by the max value of the array, and then multiply by zah0 ????
192            DO jk = 1, jpkm1
193               ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk)
194               ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk)
195            END DO
196            !
197         CASE(  30  )       !==  fixed 3D shape  ==!
198            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth )'
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
201            !                                                    ! reduction with depth
202            CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
203            !
204         CASE(  31  )       !==  time varying 3D field  ==!
205            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth , time )'
206            IF(lwp) WRITE(numout,*) '                                proportional to the velocity : |u|e/12 or |u|e^3/12'
207            !
208            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
209            !
210         CASE DEFAULT
211            CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
212         END SELECT
213         !
214         IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN       ! bilapcian and no time variation:
215            ahmt(:,:,:) = SQRT( ahmt(:,:,:) )                     ! take the square root of the coefficient
216            ahmf(:,:,:) = SQRT( ahmf(:,:,:) )
217         ENDIF
218         !
219      ENDIF
220      !
221   END SUBROUTINE ldf_dyn_init
222
223
224   SUBROUTINE ldf_dyn( kt )
225      !!----------------------------------------------------------------------
226      !!                  ***  ROUTINE ldf_dyn  ***
227      !!
228      !! ** Purpose :   update at kt the momentum lateral mixing coeff. (ahmt and ahmf)
229      !!
230      !! ** Method  :   time varying eddy viscosity coefficients:
231      !!
232      !!    nn_ahm_ijk_t = 31    ahmt, ahmf = F(i,j,k,t) = F(local velocity)
233      !!                         ( |u|e /12  or  |u|e^3/12 for laplacian or bilaplacian operator )
234      !!                BLP case : sqrt of the eddy coef, since bilaplacian is en re-entrant laplacian
235      !!
236      !! ** action  :    ahmt, ahmf   update at each time step
237      !!----------------------------------------------------------------------
238      INTEGER, INTENT(in) ::   kt   ! time step index
239      !
240      INTEGER  ::   ji, jj, jk   ! dummy loop indices
241      REAL(wp) ::   zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax   ! local scalar
242      !!----------------------------------------------------------------------
243      !
244      IF( nn_timing == 1 )  CALL timing_start('ldf_dyn')
245      !
246      SELECT CASE(  nn_ahm_ijk_t  )       !== Eddy vicosity coefficients ==!
247      !
248      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity )
249         !
250         IF( ln_dynldf_lap   ) THEN          !   laplacian operator : |u| e /12 = |u/144| e
251            DO jk = 1, jpkm1
252               DO jj = 2, jpjm1
253                  DO ji = fs_2, fs_jpim1
254                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
255                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
256                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
257                     zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
258                     zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
259                     ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax * tmask(ji,jj,jk)      ! 288= 12*12 * 2
260                     ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax * fmask(ji,jj,jk)
261                  END DO
262               END DO
263            END DO
264         ELSEIF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e
265            DO jk = 1, jpkm1
266               DO jj = 2, jpjm1
267                  DO ji = fs_2, fs_jpim1
268                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk)
269                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk)
270                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk)
271                     zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
272                     zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
273                     ahmt(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax  ) * zetmax * tmask(ji,jj,jk)
274                     ahmf(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax  ) * zefmax * fmask(ji,jj,jk)
275                  END DO
276               END DO
277            END DO
278         ENDIF
279         !
280         CALL lbc_lnk( ahmt, 'T', 1. )   ;   CALL lbc_lnk( ahmf, 'F', 1. )
281         !
282      END SELECT
283      !
284      CALL iom_put( "ahmt_2d", ahmt(:,:,1) )   ! surface u-eddy diffusivity coeff.
285      CALL iom_put( "ahmf_2d", ahmf(:,:,1) )   ! surface v-eddy diffusivity coeff.
286      CALL iom_put( "ahmt_3d", ahmt(:,:,:) )   ! 3D      u-eddy diffusivity coeff.
287      CALL iom_put( "ahmf_3d", ahmf(:,:,:) )   ! 3D      v-eddy diffusivity coeff.
288      !
289      IF( nn_timing == 1 )  CALL timing_stop('ldf_dyn')
290      !
291   END SUBROUTINE ldf_dyn
292
293   !!======================================================================
294END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.