Context Navigation

ldfdyn.F90 @ 8882

Last change on this file since 8882 was 8882, checked in by flavoni, 6 years ago
dev_CNRS_2017 branch: merged dev_r7881_ENHANCE09_RK3 with trunk r8864
Property svn:keywords set to `Id`
File size: 25.2 KB

Line
1	MODULE 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
27	IMPLICIT NONE
28	PRIVATE
29
30	PUBLIC ldf_dyn_init ! called by nemogcm.F90
31	PUBLIC ldf_dyn ! called by step.F90
32
33	! !!* Namelist namdyn_ldf : lateral mixing on momentum *
34	LOGICAL , PUBLIC :: ln_dynldf_NONE !: No operator (i.e. no explicit diffusion)
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	!! If nn_ahm_ijk_t = 32 a time and space varying Smagorinsky viscosity
45	!! will be computed.
46	REAL(wp), PUBLIC :: rn_csmc !: Smagorinsky constant of proportionality
47	REAL(wp), PUBLIC :: rn_minfac !: Multiplicative factor of theorectical minimum Smagorinsky viscosity
48	REAL(wp), PUBLIC :: rn_maxfac !: Multiplicative factor of theorectical maximum Smagorinsky viscosity
49
50	LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef.
51
52	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy diffusivity coef. at U- and V-points [m2/s or m4/s]
53	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dtensq !: horizontal tension squared (Smagorinsky only)
54	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dshesq !: horizontal shearing strain squared (Smagorinsky only)
55	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: esqt, esqf !: Square of the local gridscale (e1e2/(e1+e2))**2
56
57	REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12
58	REAL(wp) :: r1_4 = 0.25_wp ! =1/4
59	REAL(wp) :: r1_8 = 0.125_wp ! =1/8
60	REAL(wp) :: r1_288 = 1._wp / 288._wp ! =1/( 12^2 * 2 )
61
62	!! * Substitutions
63	# include "vectopt_loop_substitute.h90"
64	!!----------------------------------------------------------------------
65	!! NEMO/OPA 3.7 , NEMO Consortium (2014)
66	!! $Id$
67	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
68	!!----------------------------------------------------------------------
69	CONTAINS
70
71	SUBROUTINE ldf_dyn_init
72	!!----------------------------------------------------------------------
73	!! * ROUTINE ldf_dyn_init *
74	!!
75	!! ** Purpose : set the horizontal ocean dynamics physics
76	!!
77	!! ** Method : the eddy viscosity coef. specification depends on:
78	!! - the operator:
79	!! ln_dynldf_lap = T laplacian operator
80	!! ln_dynldf_blp = T bilaplacian operator
81	!! - the parameter nn_ahm_ijk_t:
82	!! nn_ahm_ijk_t = 0 => = constant
83	!! = 10 => = F(z) : = constant with a reduction of 1/4 with depth
84	!! =-20 => = F(i,j) = shape read in 'eddy_viscosity.nc' file
85	!! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
86	!! =-30 => = F(i,j,k) = shape read in 'eddy_viscosity.nc' file
87	!! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10)
88	!! = 31 = F(i,j,k,t) = F(local velocity) ( \|u\|e /12 laplacian operator
89	!! or \|u\|e^3/12 bilaplacian operator )
90	!! = 32 = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky)
91	!! ( L^2\|D\| laplacian operator
92	!! or L^4\|D\|/8 bilaplacian operator )
93	!!----------------------------------------------------------------------
94	INTEGER :: ji, jj, jk ! dummy loop indices
95	INTEGER :: ierr, inum, ios ! local integer
96	REAL(wp) :: zah0 ! local scalar
97	!
98	NAMELIST/namdyn_ldf/ ln_dynldf_NONE, ln_dynldf_lap, ln_dynldf_blp, & ! type of operator
99	& ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso , & ! acting direction of the operator
100	& nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0 , & ! lateral eddy coefficient
101	& rn_csmc , rn_minfac, rn_maxfac ! Smagorinsky settings
102	!!----------------------------------------------------------------------
103	!
104	REWIND( numnam_ref ) ! Namelist namdyn_ldf in reference namelist : Lateral physics
105	READ ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901)
106	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist', lwp )
107
108	REWIND( numnam_cfg ) ! Namelist namdyn_ldf in configuration namelist : Lateral physics
109	READ ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 )
110	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist', lwp )
111	IF(lwm) WRITE ( numond, namdyn_ldf )
112
113	IF(lwp) THEN ! Parameter print
114	WRITE(numout,*)
115	WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
116	WRITE(numout,*) '~~~~~~~'
117	WRITE(numout,*) ' Namelist namdyn_ldf : set lateral mixing parameters'
118	!
119	WRITE(numout,*) ' type :'
120	WRITE(numout,*) ' no explicit diffusion ln_dynldf_NONE= ', ln_dynldf_NONE
121	WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap
122	WRITE(numout,*) ' bilaplacian operator ln_dynldf_blp = ', ln_dynldf_blp
123	!
124	WRITE(numout,*) ' direction of action :'
125	WRITE(numout,*) ' iso-level ln_dynldf_lev = ', ln_dynldf_lev
126	WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor = ', ln_dynldf_hor
127	WRITE(numout,*) ' iso-neutral ln_dynldf_iso = ', ln_dynldf_iso
128	!
129	WRITE(numout,*) ' coefficients :'
130	WRITE(numout,*) ' type of time-space variation nn_ahm_ijk_t = ', nn_ahm_ijk_t
131	WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0 = ', rn_ahm_0, ' m2/s'
132	WRITE(numout,*) ' background viscosity (iso case) rn_ahm_b = ', rn_ahm_b, ' m2/s'
133	WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_bhm_0 = ', rn_bhm_0, ' m4/s'
134	WRITE(numout,*) ' Smagorinsky settings (nn_ahm_ijk_t = 32) :'
135	WRITE(numout,*) ' Smagorinsky coefficient rn_csmc = ', rn_csmc
136	WRITE(numout,*) ' factor multiplier for theorectical lower limit for '
137	WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_minfac = ', rn_minfac
138	WRITE(numout,*) ' factor multiplier for theorectical lower upper for '
139	WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_maxfac = ', rn_maxfac
140	ENDIF
141
142	! ! Parameter control
143	IF( ln_dynldf_NONE ) THEN
144	IF(lwp) WRITE(numout,*) ' No viscous operator selected. ahmt and ahmf are not allocated'
145	l_ldfdyn_time = .FALSE.
146	RETURN
147	ENDIF
148	!
149	IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN ! iso-neutral bilaplacian not implemented
150	CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' )
151	ENDIF
152
153	! ... Space/Time variation of eddy coefficients
154	! ! allocate the ahm arrays
155	ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr )
156	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays')
157	!
158	ahmt(:,:,jpk) = 0._wp ! last level always 0
159	ahmf(:,:,jpk) = 0._wp
160	!
161	! ! value of eddy mixing coef.
162	IF ( ln_dynldf_lap ) THEN ; zah0 = rn_ahm_0 ! laplacian operator
163	ELSEIF( ln_dynldf_blp ) THEN ; zah0 = ABS( rn_bhm_0 ) ! bilaplacian operator
164	ELSE ! NO viscous operator
165	CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' )
166	ENDIF
167	!
168	l_ldfdyn_time = .FALSE. ! no time variation except in case defined below
169	!
170	IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! only if a lateral diffusion operator is used
171	!
172	SELECT CASE( nn_ahm_ijk_t ) ! Specification of space time variations of ahmt, ahmf
173	!
174	CASE( 0 ) !== constant ==!
175	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = constant '
176	ahmt(:,:,:) = zah0 * tmask(:,:,:)
177	ahmf(:,:,:) = zah0 * fmask(:,:,:)
178	!
179	CASE( 10 ) !== fixed profile ==!
180	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( depth )'
181	ahmt(:,:,1) = zah0 * tmask(:,:,1) ! constant surface value
182	ahmf(:,:,1) = zah0 * fmask(:,:,1)
183	CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
184	!
185	CASE ( -20 ) !== fixed horizontal shape read in file ==!
186	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F(i,j) read in eddy_viscosity.nc file'
187	CALL iom_open( 'eddy_viscosity_2D.nc', inum )
188	CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) )
189	CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) )
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	!! better: check that the max is <=1 i.e. it is a shape from 0 to 1, not a coef that has physical dimension
194	DO jk = 2, jpkm1
195	ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk)
196	ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk)
197	END DO
198	!
199	CASE( 20 ) !== fixed horizontal shape ==!
200	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)'
201	IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor
202	IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor^3
203	!
204	CASE( -30 ) !== fixed 3D shape read in file ==!
205	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'
206	CALL iom_open( 'eddy_viscosity_3D.nc', inum )
207	CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt )
208	CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf )
209	CALL iom_close( inum )
210	!!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ????
211	!! do we introduce a scaling by the max value of the array, and then multiply by zah0 ????
212	DO jk = 1, jpkm1
213	ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk)
214	ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk)
215	END DO
216	!
217	CASE( 30 ) !== fixed 3D shape ==!
218	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth )'
219	IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor
220	IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor
221	! ! reduction with depth
222	CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
223	!
224	CASE( 31 ) !== time varying 3D field ==!
225	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth , time )'
226	IF(lwp) WRITE(numout,*) ' proportional to the velocity : \|u\|e/12 or \|u\|e^3/12'
227	!
228	l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90
229	!
230	CASE( 32 ) !== time varying 3D field ==!
231	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth , time )'
232	IF(lwp) WRITE(numout,*) ' proportional to the local deformation rate and gridscale (Smagorinsky)'
233	IF(lwp) WRITE(numout,*) ' : L^2\|D\| or L^4\|D\|/8'
234	!
235	l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90
236	!
237	! allocate arrays used in ldf_dyn.
238	ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr )
239	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays')
240	!
241	! Set local gridscale values
242	DO jj = 2, jpjm1
243	DO ji = fs_2, fs_jpim1
244	esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2
245	esqf(ji,jj) = ( e1e2f(ji,jj) /( e1f(ji,jj) + e2f(ji,jj) ) )**2
246	END DO
247	END DO
248	!
249	CASE DEFAULT
250	CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
251	END SELECT
252	!
253	IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN ! bilapcian and no time variation:
254	ahmt(:,:,:) = SQRT( ahmt(:,:,:) ) ! take the square root of the coefficient
255	ahmf(:,:,:) = SQRT( ahmf(:,:,:) )
256	ENDIF
257	!
258	ENDIF
259	!
260	END SUBROUTINE ldf_dyn_init
261
262
263	SUBROUTINE ldf_dyn( kt )
264	!!----------------------------------------------------------------------
265	!! * ROUTINE ldf_dyn *
266	!!
267	!! ** Purpose : update at kt the momentum lateral mixing coeff. (ahmt and ahmf)
268	!!
269	!! ** Method : time varying eddy viscosity coefficients:
270	!!
271	!! nn_ahm_ijk_t = 31 ahmt, ahmf = F(i,j,k,t) = F(local velocity)
272	!! ( \|u\|e /12 or \|u\|e^3/12 for laplacian or bilaplacian operator )
273	!!
274	!! nn_ahm_ijk_t = 32 ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky)
275	!! ( L^2\|D\| or L^4\|D\|/8 for laplacian or bilaplacian operator )
276	!!
277	!! ** note : in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian
278	!! ** action : ahmt, ahmf updated at each time step
279	!!----------------------------------------------------------------------
280	INTEGER, INTENT(in) :: kt ! time step index
281	!
282	INTEGER :: ji, jj, jk ! dummy loop indices
283	REAL(wp) :: zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax ! local scalar
284	REAL(wp) :: zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb ! local scalar
285	!!----------------------------------------------------------------------
286	!
287	IF( ln_timing ) CALL timing_start('ldf_dyn')
288	!
289	SELECT CASE( nn_ahm_ijk_t ) !== Eddy vicosity coefficients ==!
290	!
291	CASE( 31 ) !== time varying 3D field ==! = F( local velocity )
292	!
293	IF( ln_dynldf_lap ) THEN ! laplacian operator : \|u\| e /12 = \|u/144\| e
294	DO jk = 1, jpkm1
295	DO jj = 2, jpjm1
296	DO ji = fs_2, fs_jpim1
297	zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk)
298	zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk)
299	zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk)
300	zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
301	zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
302	ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax * tmask(ji,jj,jk) ! 288= 1212 2
303	ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax * fmask(ji,jj,jk)
304	END DO
305	END DO
306	END DO
307	ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( \|u\| e^3 /12 ) = sqrt( \|u/144\| e ) * e
308	DO jk = 1, jpkm1
309	DO jj = 2, jpjm1
310	DO ji = fs_2, fs_jpim1
311	zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk)
312	zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk)
313	zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk)
314	zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
315	zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
316	ahmt(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax ) * zetmax * tmask(ji,jj,jk)
317	ahmf(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax ) * zefmax * fmask(ji,jj,jk)
318	END DO
319	END DO
320	END DO
321	ENDIF
322	!
323	CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. )
324	!
325	!
326	CASE( 32 ) !== time varying 3D field ==! = F( local deformation rate and gridscale ) (Smagorinsky)
327	!
328	IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! laplacian operator : (C_smag/pi)^2 L^2 \|D\|
329	!
330	zcmsmag = (rn_csmc/rpi)**2 ! (C_smag/pi)^2
331	zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag ) ! lower limit stability factor scaling
332	zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt ) ! upper limit stability factor scaling
333	IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo ! provide \|U\|L^3/12 lower limit instead
334	! ! of \|U\|L^3/16 in blp case
335	DO jk = 1, jpkm1
336	!
337	DO jj = 2, jpj
338	DO ji = 2, jpi
339	zdb = ( ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) &
340	& * r1_e1t(ji,jj) * e2t(ji,jj) &
341	& - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) &
342	& * r1_e2t(ji,jj) * e1t(ji,jj) ) * tmask(ji,jj,jk)
343	dtensq(ji,jj) = zdb*zdb
344	END DO
345	END DO
346	!
347	DO jj = 1, jpjm1
348	DO ji = 1, jpim1
349	zdb = ( ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) &
350	& * r1_e2f(ji,jj) * e1f(ji,jj) &
351	& + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) &
352	& * r1_e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,jk)
353	dshesq(ji,jj) = zdb*zdb
354	END DO
355	END DO
356	!
357	DO jj = 2, jpjm1
358	DO ji = fs_2, fs_jpim1
359	!
360	zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk)
361	zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk)
362	zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk)
363	! T-point value
364	zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2
365	ahmt(ji,jj,jk) = zdelta * sqrt( dtensq(ji,jj) + &
366	& r1_4 * ( dshesq(ji,jj) + dshesq(ji,jj-1) + &
367	& dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) )
368	ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), &
369	& SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * \|U\|L/2
370	ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt)
371	! F-point value
372	zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2
373	ahmf(ji,jj,jk) = zdelta * sqrt( dshesq(ji,jj) + &
374	& r1_4 * ( dtensq(ji,jj) + dtensq(ji,jj+1) + &
375	& dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) )
376	ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), &
377	& SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * \|U\|L/2
378	ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt)
379	!
380	END DO
381	END DO
382	END DO
383	!
384	ENDIF
385	!
386	IF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 \|D\|/8)
387	! = sqrt( A_lap_smag L^2/8 )
388	! stability limits already applied to laplacian values
389	! effective default limits are \|U\|L^3/12 < B_hm < L^4/(32*2dt)
390	!
391	DO jk = 1, jpkm1
392	!
393	DO jj = 2, jpjm1
394	DO ji = fs_2, fs_jpim1
395	!
396	ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) )
397	ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) )
398	!
399	END DO
400	END DO
401	END DO
402	!
403	ENDIF
404	!
405	CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. )
406	!
407	END SELECT
408	!
409	CALL iom_put( "ahmt_2d", ahmt(:,:,1) ) ! surface u-eddy diffusivity coeff.
410	CALL iom_put( "ahmf_2d", ahmf(:,:,1) ) ! surface v-eddy diffusivity coeff.
411	CALL iom_put( "ahmt_3d", ahmt(:,:,:) ) ! 3D u-eddy diffusivity coeff.
412	CALL iom_put( "ahmf_3d", ahmf(:,:,:) ) ! 3D v-eddy diffusivity coeff.
413	!
414	IF( ln_timing ) CALL timing_stop('ldf_dyn')
415	!
416	END SUBROUTINE ldf_dyn
417
418	!!======================================================================
419	END MODULE ldfdyn

Note: See TracBrowser for help on using the repository browser.

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

Context Navigation

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 8882

Download in other formats: