Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

dynldf_iso.F90

Last change on this file was 15094, checked in by clem, 3 years ago
repair isoneutral diffusion on momentum (ln_dynldf_iso=T)
Property svn:keywords set to `Id`
File size: 21.4 KB

Line
1	MODULE dynldf_iso
2	!!======================================================================
3	!! * MODULE dynldf_iso *
4	!! Ocean dynamics: lateral viscosity trend (rotated laplacian operator)
5	!!======================================================================
6	!! History : OPA ! 97-07 (G. Madec) Original code
7	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
8	!! - ! 2004-08 (C. Talandier) New trends organization
9	!! 2.0 ! 2005-11 (G. Madec) s-coordinate: horizontal diffusion
10	!! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification,
11	!! ! add velocity dependent coefficient and optional read in file
12	!!----------------------------------------------------------------------
13
14	!!----------------------------------------------------------------------
15	!! dyn_ldf_iso : update the momentum trend with the horizontal part
16	!! of the lateral diffusion using isopycnal or horizon-
17	!! tal s-coordinate laplacian operator.
18	!!----------------------------------------------------------------------
19	USE oce ! ocean dynamics and tracers
20	USE dom_oce ! ocean space and time domain
21	USE ldfdyn ! lateral diffusion: eddy viscosity coef.
22	USE ldftra ! lateral physics: eddy diffusivity
23	USE zdf_oce ! ocean vertical physics
24	USE ldfslp ! iso-neutral slopes
25	!
26	USE in_out_manager ! I/O manager
27	USE lib_mpp ! MPP library
28	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
29	USE prtctl ! Print control
30	#if defined key_loop_fusion
31	USE dynldf_iso_lf, ONLY: dyn_ldf_iso_lf ! lateral mixing - loop fusion version (dyn_ldf_iso routine )
32	#endif
33
34	IMPLICIT NONE
35	PRIVATE
36
37	PUBLIC dyn_ldf_iso ! called by step.F90
38	PUBLIC dyn_ldf_iso_alloc ! called by nemogcm.F90
39
40	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akzu, akzv !: vertical component of rotated lateral viscosity
41
42	!! * Substitutions
43	# include "do_loop_substitute.h90"
44	# include "domzgr_substitute.h90"
45	!!----------------------------------------------------------------------
46	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
47	!! $Id$
48	!! Software governed by the CeCILL license (see ./LICENSE)
49	!!----------------------------------------------------------------------
50	CONTAINS
51
52	INTEGER FUNCTION dyn_ldf_iso_alloc()
53	!!----------------------------------------------------------------------
54	!! * ROUTINE dyn_ldf_iso_alloc *
55	!!----------------------------------------------------------------------
56	dyn_ldf_iso_alloc = 0
57	IF( .NOT. ALLOCATED( akzu ) ) THEN
58	ALLOCATE( akzu(jpi,jpj,jpk), akzv(jpi,jpj,jpk), STAT=dyn_ldf_iso_alloc )
59	!
60	IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
61	ENDIF
62	END FUNCTION dyn_ldf_iso_alloc
63
64
65	SUBROUTINE dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs )
66	!!----------------------------------------------------------------------
67	!! * ROUTINE dyn_ldf_iso *
68	!!
69	!! ** Purpose : Compute the before trend of the rotated laplacian
70	!! operator of lateral momentum diffusion except the diagonal
71	!! vertical term that will be computed in dynzdf module. Add it
72	!! to the general trend of momentum equation.
73	!!
74	!! ** Method :
75	!! The momentum lateral diffusive trend is provided by a 2nd
76	!! order operator rotated along neutral or geopotential surfaces
77	!! (in s-coordinates).
78	!! It is computed using before fields (forward in time) and isopyc-
79	!! nal or geopotential slopes computed in routine ldfslp.
80	!! Here, u and v components are considered as 2 independent scalar
81	!! fields. Therefore, the property of splitting divergent and rota-
82	!! tional part of the flow of the standard, z-coordinate laplacian
83	!! momentum diffusion is lost.
84	!! horizontal fluxes associated with the rotated lateral mixing:
85	!! u-component:
86	!! ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t di[ uu ]
87	!! - ahmt e2t * mi-1(uslp) dk[ mi(mk(uu)) ]
88	!! zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f dj[ uu ]
89	!! - ahmf e1f * mi(vslp) dk[ mj(mk(uu)) ]
90	!! v-component:
91	!! zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t di[ vv ]
92	!! - ahmf e2t * mj(uslp) dk[ mi(mk(vv)) ]
93	!! zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f dj[ vv ]
94	!! - ahmt e1f * mj-1(vslp) dk[ mj(mk(vv)) ]
95	!! take the horizontal divergence of the fluxes:
96	!! diffu = 1/(e1ue2ue3u) { di [ ziut ] + dj-1[ zjuf ] }
97	!! diffv = 1/(e1ve2ve3v) { di-1[ zivf ] + dj [ zjvt ] }
98	!! Add this trend to the general trend (uu(rhs),vv(rhs)):
99	!! uu(rhs) = uu(rhs) + diffu
100	!! CAUTION: here the isopycnal part is with a coeff. of aht. This
101	!! should be modified for applications others than orca_r2 (!!bug)
102	!!
103	!! ** Action :
104	!! -(puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the before geopotential harmonic mixing trend
105	!! -(akzu,akzv) to accompt for the diagonal vertical component
106	!! of the rotated operator in dynzdf module
107	!!----------------------------------------------------------------------
108	INTEGER , INTENT( in ) :: kt ! ocean time-step index
109	INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices
110	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
111	!
112	INTEGER :: ji, jj, jk ! dummy loop indices
113	REAL(wp) :: zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj ! local scalars
114	REAL(wp) :: zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj ! - -
115	REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0 ! - -
116	REAL(wp), DIMENSION(A2D(nn_hls)) :: ziut, zivf, zdku, zdk1u ! 2D workspace
117	REAL(wp), DIMENSION(A2D(nn_hls)) :: zjuf, zjvt, zdkv, zdk1v ! - -
118	REAL(wp), DIMENSION(A1Di(nn_hls),jpk) :: zfuw, zdiu, zdju, zdj1u ! - -
119	REAL(wp), DIMENSION(A1Di(nn_hls),jpk) :: zfvw, zdiv, zdjv, zdj1v ! - -
120	!!----------------------------------------------------------------------
121	!
122	#if defined key_loop_fusion
123	CALL dyn_ldf_iso_lf( kt, Kbb, Kmm, puu, pvv, Krhs )
124	#else
125
126	IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
127	IF( kt == nit000 ) THEN
128	IF(lwp) WRITE(numout,*)
129	IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
130	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate horizontal diffusive operator'
131	! ! allocate dyn_ldf_iso arrays
132	IF( dyn_ldf_iso_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
133	!
134	DO_2D_OVR( 0, 0, 0, 0 )
135	akzu(ji,jj,1) = 0._wp
136	akzu(ji,jj,jpk) = 0._wp
137	akzv(ji,jj,1) = 0._wp
138	akzv(ji,jj,jpk) = 0._wp
139	END_2D
140	!
141	ENDIF
142	ENDIF
143
144	!!gm bug is dyn_ldf_iso called before tra_ldf_iso .... <<<<<===== TO BE CHECKED
145	! s-coordinate: Iso-level diffusion on momentum but not on tracer
146	IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
147	!
148	DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) ! set the slopes of iso-level
149	uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
150	vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
151	wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kbb) - gdepw(ji-1,jj,jk,Kbb) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5
152	wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kbb) - gdepw(ji,jj-1,jk,Kbb) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5
153	END_3D
154	! Lateral boundary conditions on the slopes
155	IF (nn_hls == 1) CALL lbc_lnk( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp )
156	!
157	ENDIF
158
159	zaht_0 = 0.5_wp * rn_Ud * rn_Ld ! aht_0 from namtra_ldf = zaht_max
160
161	! ! ===============
162	DO jk = 1, jpkm1 ! Horizontal slab
163	! ! ===============
164
165	! Vertical u- and v-shears at level jk and jk+1
166	! ---------------------------------------------
167	! surface boundary condition: zdku(jk=1)=zdku(jk=2)
168	! zdkv(jk=1)=zdkv(jk=2)
169
170	DO_2D( 1, 1, 1, 1 )
171	zdk1u(ji,jj) = ( puu(ji,jj,jk,Kbb) -puu(ji,jj,jk+1,Kbb) ) * umask(ji,jj,jk+1)
172	zdk1v(ji,jj) = ( pvv(ji,jj,jk,Kbb) -pvv(ji,jj,jk+1,Kbb) ) * vmask(ji,jj,jk+1)
173	END_2D
174
175	IF( jk == 1 ) THEN
176	zdku(:,:) = zdk1u(:,:)
177	zdkv(:,:) = zdk1v(:,:)
178	ELSE
179	DO_2D( 1, 1, 1, 1 )
180	zdku(ji,jj) = ( puu(ji,jj,jk-1,Kbb) - puu(ji,jj,jk,Kbb) ) * umask(ji,jj,jk)
181	zdkv(ji,jj) = ( pvv(ji,jj,jk-1,Kbb) - pvv(ji,jj,jk,Kbb) ) * vmask(ji,jj,jk)
182	END_2D
183	ENDIF
184
185	! -----f-----
186	! Horizontal fluxes on U \|
187	! --------------------=== t u t
188	! \|
189	! i-flux at t-point -----f-----
190
191	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
192	DO_2D( 0, 1, 0, 0 )
193	zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) &
194	& * MIN( e3u(ji ,jj,jk,Kmm), &
195	& e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
196
197	zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
198	& + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ) , 1._wp )
199
200	zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
201
202	ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) &
203	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
204	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
205	END_2D
206	ELSE ! other coordinate system (zco or sco) : e3t
207	DO_2D( 0, 1, 0, 0 )
208	zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) &
209	& * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
210
211	zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) &
212	& + umask(ji-1,jj,jk+1) + umask(ji,jj,jk ) , 1._wp )
213
214	zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
215
216	ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) &
217	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
218	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
219	END_2D
220	ENDIF
221
222	! j-flux at f-point
223	DO_2D( 1, 0, 1, 0 )
224	zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) &
225	& * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
226
227	zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
228	& + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ) , 1._wp )
229
230	zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
231
232	zjuf(ji,jj) = ( zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) ) &
233	& + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
234	& +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) * fmask(ji,jj,jk)
235	END_2D
236
237	! \| t \|
238	! Horizontal fluxes on V \| \|
239	! --------------------=== f---v---f
240	! \| \|
241	! i-flux at f-point \| t \|
242
243	DO_2D( 1, 0, 0, 0 )
244	zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) &
245	& * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
246
247	zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
248	& + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ) , 1._wp )
249
250	zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
251
252	zivf(ji,jj) = ( zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) ) &
253	& + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) &
254	& + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk)
255	END_2D
256
257	! j-flux at t-point
258	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
259	DO_2D( 1, 0, 0, 1 )
260	zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) &
261	& * MIN( e3v(ji,jj ,jk,Kmm), &
262	& e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj)
263
264	zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
265	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ) , 1._wp )
266
267	zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
268
269	zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) &
270	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
271	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
272	END_2D
273	ELSE ! other coordinate system (zco or sco) : e3t
274	DO_2D( 1, 0, 0, 1 )
275	zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) &
276	& * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
277
278	zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
279	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )
280
281	zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
282
283	zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) &
284	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
285	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
286	END_2D
287	ENDIF
288
289
290	! Second derivative (divergence) and add to the general trend
291	! -----------------------------------------------------------
292	DO_2D( 0, 0, 0, 0 ) !!gm Question vectop possible??? !!bug
293	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( ziut(ji+1,jj) - ziut(ji,jj ) &
294	& + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) &
295	& / e3u(ji,jj,jk,Kmm)
296	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zivf(ji,jj ) - zivf(ji-1,jj) &
297	& + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) &
298	& / e3v(ji,jj,jk,Kmm)
299	END_2D
300	! ! ===============
301	END DO ! End of slab
302	! ! ===============
303
304	! print sum trends (used for debugging)
305	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldfh - Ua: ', mask1=umask, &
306	& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
307
308
309	! ! ===============
310	DO jj = ntsj, ntej ! Vertical slab
311	! ! ===============
312
313
314	! I. vertical trends associated with the lateral mixing
315	! =====================================================
316	! (excluding the vertical flux proportional to dk[t]
317
318
319	! I.1 horizontal momentum gradient
320	! --------------------------------
321
322	DO jk = 1, jpk
323	DO ji = ntsi, ntei + nn_hls
324	! i-gradient of u at jj
325	zdiu (ji,jk) = tmask(ji,jj ,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji-1,jj ,jk,Kbb) )
326	! j-gradient of u and v at jj
327	zdju (ji,jk) = fmask(ji,jj ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji ,jj ,jk,Kbb) )
328	zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pvv(ji,jj ,jk,Kbb) - pvv(ji ,jj-1,jk,Kbb) )
329	! j-gradient of u and v at jj+1
330	zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji ,jj-1,jk,Kbb) )
331	zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) )
332	END DO
333	END DO
334	DO jk = 1, jpk
335	DO ji = ntsi - nn_hls, ntei
336	! i-gradient of v at jj
337	zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) )
338	END DO
339	END DO
340
341
342	! I.2 Vertical fluxes
343	! -------------------
344
345	! Surface and bottom vertical fluxes set to zero
346	DO ji = ntsi - nn_hls, ntei + nn_hls
347	zfuw(ji, 1 ) = 0.e0
348	zfvw(ji, 1 ) = 0.e0
349	zfuw(ji,jpk) = 0.e0
350	zfvw(ji,jpk) = 0.e0
351	END DO
352
353	! interior (2=<jk=<jpk-1) on U field
354	DO jk = 2, jpkm1
355	DO ji = ntsi, ntei
356	zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk)
357	!
358	zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
359	zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
360	!
361	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) &
362	+ tmask(ji,jj,jk )+tmask(ji+1,jj,jk ) , 1. )
363	zmkf = 1./MAX( fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1) &
364	+ fmask(ji,jj-1,jk ) + fmask(ji,jj,jk ) , 1. )
365
366	zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
367	zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
368	! vertical flux on u field
369	zfuw(ji,jk) = zcof3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &
370	& + zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &
371	& + zcof4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &
372	& + zdj1u(ji,jk ) + zdju (ji ,jk ) )
373	! vertical mixing coefficient (akzu)
374	! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
375	akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
376	END DO
377	END DO
378
379	! interior (2=<jk=<jpk-1) on V field
380	DO jk = 2, jpkm1
381	DO ji = ntsi, ntei
382	zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk)
383	!
384	zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
385	zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
386	!
387	zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) &
388	& + fmask(ji-1,jj,jk )+fmask(ji,jj,jk ) , 1. )
389	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) &
390	& + tmask(ji,jj,jk )+tmask(ji,jj+1,jk ) , 1. )
391
392	zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
393	zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
394	! vertical flux on v field
395	zfvw(ji,jk) = zcof3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &
396	& + zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &
397	& + zcof4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &
398	& + zdjv (ji,jk ) + zdj1v(ji ,jk ) )
399	! vertical mixing coefficient (akzv)
400	! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
401	akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
402	END DO
403	END DO
404
405
406	! I.3 Divergence of vertical fluxes added to the general tracer trend
407	! -------------------------------------------------------------------
408	DO jk = 1, jpkm1
409	DO ji = ntsi, ntei
410	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) &
411	& / e3u(ji,jj,jk,Kmm)
412	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) &
413	& / e3v(ji,jj,jk,Kmm)
414	END DO
415	END DO
416	! ! ===============
417	END DO ! End of slab
418	! ! ===============
419	#endif
420	END SUBROUTINE dyn_ldf_iso
421
422	!!======================================================================
423	END MODULE dynldf_iso

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

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

Context Navigation

source: NEMO/trunk/src/OCE/DYN/dynldf_iso.F90

Download in other formats: