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_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 5972

Last change on this file since 5972 was 5972, checked in by timgraham, 8 years ago

Upgraded to head of trunk (r5936)

  • 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.