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

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 4596

Last change on this file since 4596 was 4596, checked in by gm, 10 years ago

#1260: LDF simplification + bilap iso-neutral for TRA and GYRE
  • Property svn:keywords set to Id
File size: 15.9 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_c3d   : 3D eddy viscosity coefficient initialization
15   !!   ldf_dyn_c2d   : 2D eddy viscosity coefficient initialization
16   !!   ldf_dyn_c1d   : 1D eddy viscosity coefficient initialization
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers   
19   USE dom_oce         ! ocean space and time domain
20   USE phycst          ! physical constants
21   USE ldfc1d          ! lateral diffusion: 1D case
22   USE ldfc2d          ! lateral diffusion: 2D case
23!   USE ldfc3d          ! lateral diffusion: 3D case
24   !
25   USE in_out_manager  ! I/O manager
26   USE iom             ! I/O module for ehanced bottom friction file
27   USE timing          ! Timing
28   USE lib_mpp         ! distribued memory computing library
29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
30   USE wrk_nemo        ! Memory Allocation
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   ldf_dyn_init   ! called by opa.F90
36
37   !                                                !!* Namelist namdyn_ldf : lateral mixing on momentum *
38   LOGICAL , PUBLIC ::   ln_dynldf_lap = .TRUE.      !: laplacian operator
39   LOGICAL , PUBLIC ::   ln_dynldf_blp = .FALSE.     !: bilaplacian operator
40   LOGICAL , PUBLIC ::   ln_dynldf_lev = .FALSE.     !: iso-level direction
41   LOGICAL , PUBLIC ::   ln_dynldf_hor = .TRUE.      !: horizontal (geopotential) direction
42   LOGICAL , PUBLIC ::   ln_dynldf_iso = .FALSE.     !: iso-neutral direction
43   INTEGER , PUBLIC ::   nn_ahm_ijk_t  = 0           !:   ??????  !!gm
44   REAL(wp), PUBLIC ::   rn_ahm_0      = 40000._wp   !: lateral laplacian eddy viscosity            [m2/s]
45   REAL(wp), PUBLIC ::   rn_ahm_b      =     0._wp   !: lateral laplacian background eddy viscosity [m2/s]
46   REAL(wp), PUBLIC ::   rn_bhm_0      = 5.0e+11_wp  !: lateral bilaplacian eddy viscosity          [m4/s]
47
48   LOGICAL , PUBLIC ::   l_ldfdyn_time = .FALSE.     !: flag for time variation of the lateral eddy viscosity coef.
49
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy diffusivity coef. at U- and V-points   [m2/s or m4/s]
51
52   REAL(wp) ::   r1_12   = 1._wp / 12._wp    ! =1/12
53   REAL(wp) ::   r1_4    = 0.25_wp           ! =1/4
54   REAL(wp) ::   r1_288  = 1._wp / 288._wp   ! =1/( 12^2 * 2 )
55
56   !! * Substitutions
57#  include "domzgr_substitute.h90"
58#  include "vectopt_loop_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE ldf_dyn_init
67      !!----------------------------------------------------------------------
68      !!                  ***  ROUTINE ldf_dyn_init  ***
69      !!                   
70      !! ** Purpose :   set the horizontal ocean dynamics physics
71      !!
72      !! ** Method  :   the eddy viscosity coef. specification depends on:
73      !!
74      !!    ln_dynldf_lap = T     laplacian operator
75      !!    ln_dynldf_blp = T   bilaplacian operator
76      !!
77      !!    nn_ahm_ijk_t  =  0 => = constant
78      !!                  !
79      !!                  = 10 => = F(z) : constant with a reduction of 1/4 with depth
80      !!                  !
81      !!                  =-20 => = F(i,j)   = shape read in 'eddy_viscosity.nc' file
82      !!                  = 20    = F(i,j)   = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
83      !!                  !
84      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_viscosity.nc'  file
85      !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10)
86      !!                  = 31    = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator
87      !!                                                           or |u|e^3/12 bilaplacian operator )
88      !!
89      !!----------------------------------------------------------------------
90      INTEGER  ::   jk                ! dummy loop indices
91      INTEGER  ::   ierr, inum, ios   ! local integer
92      REAL(wp) ::   zah0              ! local scalar
93      !
94      NAMELIST/namdyn_ldf/ ln_dynldf_lap, ln_dynldf_blp,                  &
95         &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso,   &
96         &                 nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0
97      !!----------------------------------------------------------------------
98!!      NAMELIST/namdyn_ldf/ ln_dynldf_lap  , ln_dynldf_bilap,                  &
99!!         &                 ln_dynldf_level, ln_dynldf_hor  , ln_dynldf_iso,   &
100!!         &                 rn_ahm_0_lap   , rn_ahmb_0      , rn_ahm_0_blp ,   &
101!
102!!         &                 rn_cmsmag_1    , rn_cmsmag_2    , rn_cmsh,         &
103!
104!!         &                 rn_ahm_m_lap   , rn_ahm_m_blp
105   !!----------------------------------------------------------------------
106
107      REWIND( numnam_ref )              ! Namelist namdyn_ldf in reference namelist : Lateral physics
108      READ  ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901)
109901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist', lwp )
110
111      REWIND( numnam_cfg )              ! Namelist namdyn_ldf in configuration namelist : Lateral physics
112      READ  ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 )
113902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist', lwp )
114      WRITE ( numond, namdyn_ldf )
115
116      IF(lwp) THEN                      ! Parameter print
117         WRITE(numout,*)
118         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
119         WRITE(numout,*) '~~~~~~~'
120         WRITE(numout,*) '   Namelist nam_dynldf : set lateral mixing parameters'
121         !
122         WRITE(numout,*) '      type :'
123         WRITE(numout,*) '         laplacian operator                   ln_dynldf_lap = ', ln_dynldf_lap
124         WRITE(numout,*) '         bilaplacian operator                 ln_dynldf_blp = ', ln_dynldf_blp
125         !
126         WRITE(numout,*) '      direction of action :'
127         WRITE(numout,*) '         iso-level                            ln_dynldf_lev = ', ln_dynldf_lev
128         WRITE(numout,*) '         horizontal (geopotential)            ln_dynldf_hor = ', ln_dynldf_hor
129         WRITE(numout,*) '         iso-neutral                          ln_dynldf_iso = ', ln_dynldf_iso
130         !
131         WRITE(numout,*) '      coefficients :'
132         WRITE(numout,*) '         type of time-space variation         nn_ahm_ijk_t  = ', nn_ahm_ijk_t
133         WRITE(numout,*) '         lateral laplacian eddy viscosity     rn_ahm_0_lap  = ', rn_ahm_0, ' m2/s'
134         WRITE(numout,*) '         background viscosity (iso case)      rn_ahm_b      = ', rn_ahm_b, ' m2/s'
135         WRITE(numout,*) '         lateral bilaplacian eddy viscosity   rn_ahm_0_blp  = ', rn_bhm_0, ' m4/s'
136      ENDIF
137
138      !                                ! Parameter control
139      IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN     ! iso-neutral bilaplacian not implemented
140         CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' ) 
141      ENDIF
142
143      ! ... Space/Time variation of eddy coefficients
144      !                                               ! allocate the ahm arrays
145      ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr )
146      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays')
147      !
148      ahmt(:,:,jpk) = 0._wp                           ! last level always 0 
149      ahmf(:,:,jpk) = 0._wp
150      !
151      !                                               ! value of eddy mixing coef.
152      IF    ( ln_dynldf_lap ) THEN   ;   zah0 =            rn_ahm_0         !   laplacian operator
153      ELSEIF( ln_dynldf_blp ) THEN   ;   zah0 = SQRT( ABS( rn_bhm_0 ) )     ! bilaplacian operator
154      ELSE                                                                  ! NO viscous  operator
155         CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' )
156      ENDIF
157      !
158      l_ldfdyn_time = .FALSE.                         ! no time variation except in case defined below
159      !
160      IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN     ! only if a lateral diffusion operator is used
161         !
162         SELECT CASE(  nn_ahm_ijk_t  )                ! Specification of space time variations of ahmt, ahmf
163         !
164         CASE(   0  )      !==  constant  ==!
165            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = constant '
166            ahmt(:,:,:) = zah0 * tmask(:,:,:)
167            ahmf(:,:,:) = zah0 * fmask(:,:,:)
168            !
169         CASE(  10  )      !==  fixed profile  ==!
170            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( depth )'
171            ahmt(:,:,1) = zah0 * tmask(:,:,1)                      ! constant surface value
172            ahmf(:,:,1) = zah0 * fmask(:,:,1)
173            CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
174            !
175         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
176            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F(i,j) read in eddy_viscosity.nc file'
177            CALL iom_open( 'eddy_viscosity.nc', inum )
178            CALL iom_get ( inum, jpdom_data, 'ahmt_2D', ahmt(:,:,1) )
179            CALL iom_get ( inum, jpdom_data, 'ahmf_2D', ahmf(:,:,1) )
180            CALL iom_close( inum )
181            DO jk = 2, jpkm1
182               ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk)
183               ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk)
184            END DO
185            !
186         CASE(  20  )      !== fixed horizontal shape  ==!
187            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)'
188            IF( ln_dynldf_lap )   CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor
189            IF( ln_dynldf_blp )   CALL ldf_c2d( 'DYN', 'BiL', zah0, ahmt, ahmf )    ! surface value proportional to scale factor^3
190            !
191         CASE( -30  )      !== fixed 3D shape read in file  ==!
192            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'
193            CALL iom_open( 'eddy_diffusivity_3D.nc', inum )
194            CALL iom_get ( inum, jpdom_data, 'ahmt', ahmt )
195            CALL iom_get ( inum, jpdom_data, 'ahmf', ahmf )
196            CALL iom_close( inum )
197            DO jk = 1, jpkm1
198               ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk)
199               ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk)
200            END DO
201            !
202         CASE(  30  )       !==  fixed 3D shape  ==!
203            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth )'
204            IF( ln_dynldf_lap )   CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf )    ! surface value proportional to scale factor
205            IF( ln_dynldf_blp )   CALL ldf_c2d( 'DYN', 'BiL', zah0, ahmt, ahmf )    ! surface value proportional to scale factor
206            !                                                    ! reduction with depth
207            CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
208            !
209         CASE(  31  )       !==  time varying 3D field  ==!
210            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth , time )'
211            IF(lwp) WRITE(numout,*) '                                proportional to the velocity : |u|e/12 or |u|e^3/12'
212            !
213            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
214            !
215            CALL ctl_stop( 'STOP', 'ldf_dyn_init: ahm=F(velocity) not yet implemented')
216            !
217         CASE DEFAULT
218            CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
219         END SELECT
220         !
221      ENDIF
222      !
223   END SUBROUTINE ldf_dyn_init
224
225
226   SUBROUTINE ldf_dyn( kt )
227      !!----------------------------------------------------------------------
228      !!                  ***  ROUTINE ldf_dyn_init  ***
229      !!
230      !! ** Purpose :   update at kt the momentum lateral mixing coeff. (ahmt and ahmf)
231      !!
232      !! ** Method  :   time varying eddy viscosity coefficients:
233      !!
234      !!    nn_ahm_ijk_t = 31    ahmt, ahmf = F(i,j,k,t) = F(local velocity)
235      !!                         ( |u|e /12  or  |u|e^3/12 for laplacian or bilaplacian operator )
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_i_j_p1, zu2pv2_i_j, zu2pv2_i_j_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
252            DO jk = 1, jpkm1
253               DO jj = 2, jpjm1
254                  DO ji = fs_2, fs_jpim1
255!!gm should probably be defined as an average of " e1u*u + e2v*v " not the
256                     zu2pv2_i_j_p1 = ub(ji+1,jj,jk) * ub(ji+1,jj,jk) + vb(ji,jj+1,jk) * vb(ji,jj+1,jk)
257                     zu2pv2_i_j    = ub(ji  ,jj,jk) * ub(ji  ,jj,jk) + vb(ji,jj  ,jk) * vb(ji,jj  ,jk)
258                     zu2pv2_i_j_m1 = ub(ji-1,jj,jk) * ub(ji-1,jj,jk) + vb(ji,jj-1,jk) * vb(ji,jj-1,jk)
259                     zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
260                     zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
261                     ahmt(ji,jj,jk) = SQRT( zu2pv2_i_j + zu2pv2_i_j_m1 * r1_288 ) * zetmax * tmask(ji,jj,jk)      ! 288= 12*12 * 2
262                     ahmf(ji,jj,jk) = SQRT( zu2pv2_i_j + zu2pv2_i_j_p1 * r1_288 ) * zefmax * fmask(ji,jj,jk)
263                  END DO
264               END DO
265            END DO
266         ELSEIF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( |u| e^3 /12 )
267            DO jk = 1, jpkm1
268               DO jj = 2, jpjm1
269                  DO ji = fs_2, fs_jpim1
270                     zu2pv2_i_j_p1 = ub(ji+1,jj,jk) * ub(ji+1,jj,jk) + vb(ji,jj+1,jk) * vb(ji,jj+1,jk)
271                     zu2pv2_i_j    = ub(ji  ,jj,jk) * ub(ji  ,jj,jk) + vb(ji,jj  ,jk) * vb(ji,jj  ,jk)
272                     zu2pv2_i_j_m1 = ub(ji-1,jj,jk) * ub(ji-1,jj,jk) + vb(ji,jj-1,jk) * vb(ji,jj-1,jk)
273                     zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
274                     zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
275                     ahmt(ji,jj,jk) = SQRT(  SQRT( zu2pv2_i_j + zu2pv2_i_j_m1 * r1_288 ) * zetmax  ) * zetmax * tmask(ji,jj,jk)
276                     ahmf(ji,jj,jk) = SQRT(  SQRT( zu2pv2_i_j + zu2pv2_i_j_p1 * r1_288 ) * zefmax  ) * zefmax * fmask(ji,jj,jk)
277                  END DO
278               END DO
279            END DO
280         ENDIF
281         !
282         CALL lbc_lnk( ahmt, 'U', 1. )   ;   CALL lbc_lnk( ahmf, 'V', 1. )
283         !
284      END SELECT
285      !
286      IF( nn_timing == 1 )  CALL timing_stop('ldf_dyn')
287     !
288   END SUBROUTINE ldf_dyn
289
290   !!======================================================================
291END MODULE ldfdyn
Note: See TracBrowser for help on using the repository browser.