Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

ldfdyn.F90 @ 12286

Last change on this file since 12286 was 12286, checked in by jcastill, 4 years ago
Remove svn keywords
File size: 25.0 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	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	!! 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_lap, ln_dynldf_blp, &
99	& ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso, &
100	& nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0, &
101	& rn_csmc , rn_minfac, rn_maxfac
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,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap
121	WRITE(numout,*) ' bilaplacian operator ln_dynldf_blp = ', ln_dynldf_blp
122	!
123	WRITE(numout,*) ' direction of action :'
124	WRITE(numout,*) ' iso-level ln_dynldf_lev = ', ln_dynldf_lev
125	WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor = ', ln_dynldf_hor
126	WRITE(numout,*) ' iso-neutral ln_dynldf_iso = ', ln_dynldf_iso
127	!
128	WRITE(numout,*) ' coefficients :'
129	WRITE(numout,*) ' type of time-space variation nn_ahm_ijk_t = ', nn_ahm_ijk_t
130	WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0 = ', rn_ahm_0, ' m2/s'
131	WRITE(numout,*) ' background viscosity (iso case) rn_ahm_b = ', rn_ahm_b, ' m2/s'
132	WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_bhm_0 = ', rn_bhm_0, ' m4/s'
133	WRITE(numout,*) ' smagorinsky settings (nn_ahm_ijk_t = 32) :'
134	WRITE(numout,*) ' Smagorinsky coefficient rn_csmc = ', rn_csmc
135	WRITE(numout,*) ' factor multiplier for theorectical lower limit for '
136	WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_minfac = ', rn_minfac
137	WRITE(numout,*) ' factor multiplier for theorectical lower upper for '
138	WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_maxfac = ', rn_maxfac
139	ENDIF
140
141	! ! Parameter control
142	IF( .NOT.ln_dynldf_lap .AND. .NOT.ln_dynldf_blp ) THEN
143	IF(lwp) WRITE(numout,*) ' No viscous operator selected. ahmt and ahmf are not allocated'
144	l_ldfdyn_time = .FALSE.
145	RETURN
146	ENDIF
147	!
148	IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN ! iso-neutral bilaplacian not implemented
149	CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' )
150	ENDIF
151
152	! ... Space/Time variation of eddy coefficients
153	! ! allocate the ahm arrays
154	ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr )
155	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays')
156	!
157	ahmt(:,:,jpk) = 0._wp ! last level always 0
158	ahmf(:,:,jpk) = 0._wp
159	!
160	! ! value of eddy mixing coef.
161	IF ( ln_dynldf_lap ) THEN ; zah0 = rn_ahm_0 ! laplacian operator
162	ELSEIF( ln_dynldf_blp ) THEN ; zah0 = ABS( rn_bhm_0 ) ! bilaplacian operator
163	ELSE ! NO viscous operator
164	CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' )
165	ENDIF
166	!
167	l_ldfdyn_time = .FALSE. ! no time variation except in case defined below
168	!
169	IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! only if a lateral diffusion operator is used
170	!
171	SELECT CASE( nn_ahm_ijk_t ) ! Specification of space time variations of ahmt, ahmf
172	!
173	CASE( 0 ) !== constant ==!
174	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = constant '
175	ahmt(:,:,:) = zah0 * tmask(:,:,:)
176	ahmf(:,:,:) = zah0 * fmask(:,:,:)
177	!
178	CASE( 10 ) !== fixed profile ==!
179	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( depth )'
180	ahmt(:,:,1) = zah0 * tmask(:,:,1) ! constant surface value
181	ahmf(:,:,1) = zah0 * fmask(:,:,1)
182	CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
183	!
184	CASE ( -20 ) !== fixed horizontal shape read in file ==!
185	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F(i,j) read in eddy_viscosity.nc file'
186	CALL iom_open( 'eddy_viscosity_2D.nc', inum )
187	CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) )
188	CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) )
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	!! better: check that the max is <=1 i.e. it is a shape from 0 to 1, not a coef that has physical dimension
193	DO jk = 2, jpkm1
194	ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk)
195	ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk)
196	END DO
197	!
198	CASE( 20 ) !== fixed horizontal shape ==!
199	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)'
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^3
202	!
203	CASE( -30 ) !== fixed 3D shape read in file ==!
204	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'
205	CALL iom_open( 'eddy_viscosity_3D.nc', inum )
206	CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt )
207	CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf )
208	CALL iom_close( inum )
209	!!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ????
210	!! do we introduce a scaling by the max value of the array, and then multiply by zah0 ????
211	DO jk = 1, jpkm1
212	ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk)
213	ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk)
214	END DO
215	!
216	CASE( 30 ) !== fixed 3D shape ==!
217	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth )'
218	IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor
219	IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor
220	! ! reduction with depth
221	CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
222	!
223	CASE( 31 ) !== time varying 3D field ==!
224	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth , time )'
225	IF(lwp) WRITE(numout,*) ' proportional to the velocity : \|u\|e/12 or \|u\|e^3/12'
226	!
227	l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90
228	!
229	CASE( 32 ) !== time varying 3D field ==!
230	IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth , time )'
231	IF(lwp) WRITE(numout,*) ' proportional to the local deformation rate and gridscale (Smagorinsky)'
232	IF(lwp) WRITE(numout,*) ' : L^2\|D\| or L^4\|D\|/8'
233	!
234	l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90
235	!
236	! allocate arrays used in ldf_dyn.
237	ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr )
238	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays')
239	!
240	! Set local gridscale values
241	DO jj = 2, jpjm1
242	DO ji = fs_2, fs_jpim1
243	esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2
244	esqf(ji,jj) = ( e1e2f(ji,jj) /( e1f(ji,jj) + e2f(ji,jj) ) )**2
245	END DO
246	END DO
247	!
248	CASE DEFAULT
249	CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
250	END SELECT
251	!
252	IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN ! bilapcian and no time variation:
253	ahmt(:,:,:) = SQRT( ahmt(:,:,:) ) ! take the square root of the coefficient
254	ahmf(:,:,:) = SQRT( ahmf(:,:,:) )
255	ENDIF
256	!
257	ENDIF
258	!
259	END SUBROUTINE ldf_dyn_init
260
261
262	SUBROUTINE ldf_dyn( kt )
263	!!----------------------------------------------------------------------
264	!! * ROUTINE ldf_dyn *
265	!!
266	!! ** Purpose : update at kt the momentum lateral mixing coeff. (ahmt and ahmf)
267	!!
268	!! ** Method : time varying eddy viscosity coefficients:
269	!!
270	!! nn_ahm_ijk_t = 31 ahmt, ahmf = F(i,j,k,t) = F(local velocity)
271	!! ( \|u\|e /12 or \|u\|e^3/12 for laplacian or bilaplacian operator )
272	!!
273	!! nn_ahm_ijk_t = 32 ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky)
274	!! ( L^2\|D\| or L^4\|D\|/8 for laplacian or bilaplacian operator )
275	!!
276	!! ** note : in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian
277	!! ** action : ahmt, ahmf updated at each time step
278	!!----------------------------------------------------------------------
279	INTEGER, INTENT(in) :: kt ! time step index
280	!
281	INTEGER :: ji, jj, jk ! dummy loop indices
282	REAL(wp) :: zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax ! local scalar
283	REAL(wp) :: zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb ! local scalar
284	!!----------------------------------------------------------------------
285	!
286	IF( nn_timing == 1 ) CALL timing_start('ldf_dyn')
287	!
288	SELECT CASE( nn_ahm_ijk_t ) !== Eddy vicosity coefficients ==!
289	!
290	CASE( 31 ) !== time varying 3D field ==! = F( local velocity )
291	!
292	IF( ln_dynldf_lap ) THEN ! laplacian operator : \|u\| e /12 = \|u/144\| e
293	DO jk = 1, jpkm1
294	DO jj = 2, jpjm1
295	DO ji = fs_2, fs_jpim1
296	zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk)
297	zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk)
298	zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk)
299	zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
300	zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
301	ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax * tmask(ji,jj,jk) ! 288= 1212 2
302	ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax * fmask(ji,jj,jk)
303	END DO
304	END DO
305	END DO
306	ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( \|u\| e^3 /12 ) = sqrt( \|u/144\| e ) * e
307	DO jk = 1, jpkm1
308	DO jj = 2, jpjm1
309	DO ji = fs_2, fs_jpim1
310	zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk)
311	zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk)
312	zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk)
313	zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) )
314	zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) )
315	ahmt(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax ) * zetmax * tmask(ji,jj,jk)
316	ahmf(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax ) * zefmax * fmask(ji,jj,jk)
317	END DO
318	END DO
319	END DO
320	ENDIF
321	!
322	CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. )
323	!
324	!
325	CASE( 32 ) !== time varying 3D field ==! = F( local deformation rate and gridscale ) (Smagorinsky)
326	!
327	IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! laplacian operator : (C_smag/pi)^2 L^2 \|D\|
328	!
329	zcmsmag = (rn_csmc/rpi)**2 ! (C_smag/pi)^2
330	zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag ) ! lower limit stability factor scaling
331	zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt ) ! upper limit stability factor scaling
332	IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo ! provide \|U\|L^3/12 lower limit instead
333	! ! of \|U\|L^3/16 in blp case
334	DO jk = 1, jpkm1
335	!
336	DO jj = 2, jpj
337	DO ji = 2, jpi
338	zdb = ( ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) &
339	& * r1_e1t(ji,jj) * e2t(ji,jj) &
340	& - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) &
341	& * r1_e2t(ji,jj) * e1t(ji,jj) ) * tmask(ji,jj,jk)
342	dtensq(ji,jj) = zdb*zdb
343	END DO
344	END DO
345	!
346	DO jj = 1, jpjm1
347	DO ji = 1, jpim1
348	zdb = ( ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) &
349	& * r1_e2f(ji,jj) * e1f(ji,jj) &
350	& + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) &
351	& * r1_e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,jk)
352	dshesq(ji,jj) = zdb*zdb
353	END DO
354	END DO
355	!
356	DO jj = 2, jpjm1
357	DO ji = fs_2, fs_jpim1
358	!
359	zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk)
360	zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk)
361	zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk)
362	! T-point value
363	zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2
364	ahmt(ji,jj,jk) = zdelta * sqrt( dtensq(ji,jj) + &
365	& r1_4 * ( dshesq(ji,jj) + dshesq(ji,jj-1) + &
366	& dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) )
367	ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), &
368	& SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * \|U\|L/2
369	ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt)
370	! F-point value
371	zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2
372	ahmf(ji,jj,jk) = zdelta * sqrt( dshesq(ji,jj) + &
373	& r1_4 * ( dtensq(ji,jj) + dtensq(ji,jj+1) + &
374	& dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) )
375	ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), &
376	& SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * \|U\|L/2
377	ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt)
378	!
379	END DO
380	END DO
381	END DO
382	!
383	ENDIF
384	!
385	IF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 \|D\|/8)
386	! = sqrt( A_lap_smag L^2/8 )
387	! stability limits already applied to laplacian values
388	! effective default limits are \|U\|L^3/12 < B_hm < L^4/(32*2dt)
389	!
390	DO jk = 1, jpkm1
391	!
392	DO jj = 2, jpjm1
393	DO ji = fs_2, fs_jpim1
394	!
395	ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) )
396	ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) )
397	!
398	END DO
399	END DO
400	END DO
401	!
402	ENDIF
403	!
404	CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. )
405	!
406	END SELECT
407	!
408	CALL iom_put( "ahmt_2d", ahmt(:,:,1) ) ! surface u-eddy diffusivity coeff.
409	CALL iom_put( "ahmf_2d", ahmf(:,:,1) ) ! surface v-eddy diffusivity coeff.
410	CALL iom_put( "ahmt_3d", ahmt(:,:,:) ) ! 3D u-eddy diffusivity coeff.
411	CALL iom_put( "ahmf_3d", ahmf(:,:,:) ) ! 3D v-eddy diffusivity coeff.
412	!
413	IF( nn_timing == 1 ) CALL timing_stop('ldf_dyn')
414	!
415	END SUBROUTINE ldf_dyn
416
417	!!======================================================================
418	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/UKMO/r8395_cpl_tauwav/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 @ 12286

Download in other formats: