Context Navigation

dynldf_iso.F90 @ 15540

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

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

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

Context Navigation

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DYN/dynldf_iso.F90 @ 15540

Download in other formats: