Context Navigation

vlspltqs_loc.F @ 4412

Last change on this file since 4412 was 4325, checked in by dcugnet, 2 years ago
simplify the parser usage: the getKey_init routine is now embedded in the readTracersFile routine. the initIsotopes routine is now embedded in the readIsotopesFile routine. the database is now unique, but can be changed using the get/setKeysDBase. the derived types descriptions, originally located in trac_types_mod, are moved to readTracFiles_mod. few checkings moved from infotrac to the routine testIsotopes, contained in the readIsotopesFile function from readTracFiles_mod. the readTracersFiles and readIsotopesFile routines no longer use a tracers/isotopes argument. remove tnat and alpha_ideal from infotrac ; use instead getKey to get them where they are used (check_isotopes, dynetat0, iniacademic) the trac_type field %Childs is renamed %Children move the isoSelect routine and the corresponding variables routine from infotrac and infotrac_phy to readTracFiles_mod infotrac_phy routine is now fully independant of the (very similar) routine infotrac (init_infotrac_phy has no arguments left). all the explicit keys of the trac_type are now included in the embedded keys database, accessible using the getKey function. the getKey/addKey routines are expanded to handle vectors of integers, reals, logicals or strings. few subroutines converted into functions with error return value. corrections for isotopic tagging tracers mode (to be continued).
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:keywords set to `Id`
File size: 24.4 KB

Line
1	!
2	! $Id: vlspltqs_loc.F 4325 2022-11-07 02:09:43Z evignon $
3	!
4	SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x,iq)
5	c
6	c Auteurs: P.Le Van, F.Hourdin, F.Forget
7	c
8	c ********************************************************************
9	c Shema d''advection " pseudo amont " .
10	c ********************************************************************
11	c
12	c --------------------------------------------------------------------
13	USE parallel_lmdz
14	USE infotrac, ONLY : nqtot,tracers, ! CRisi &
15	& min_qParent,min_qMass,min_ratio ! MVals et CRisi
16	IMPLICIT NONE
17	c
18	include "dimensions.h"
19	include "paramet.h"
20	c
21	c
22	c Arguments:
23	c ----------
24	REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
25	REAL u_m( ijb_u:ije_u,llm )
26	REAL q(ijb_u:ije_u,llm,nqtot)
27	REAL qsat(ijb_u:ije_u,llm)
28	INTEGER iq ! CRisi
29	c
30	c Local
31	c ---------
32	c
33	INTEGER ij,l,j,i,iju,ijq,indu(ijnb_u),niju
34	INTEGER n0,iadvplus(ijb_u:ije_u,llm),nl(llm)
35	c
36	REAL new_m,zu_m,zdum(ijb_u:ije_u,llm)
37	REAL dxq(ijb_u:ije_u,llm),dxqu(ijb_u:ije_u)
38	REAL zz(ijb_u:ije_u)
39	REAL adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm)
40	REAL u_mq(ijb_u:ije_u,llm)
41	REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
42	INTEGER ifils,iq2 ! CRisi
43
44
45	REAL SSUM
46
47
48	INTEGER ijb,ije,ijb_x,ije_x
49
50	!write(,) 'vlspltqs 58: entree vlxqs_loc, iq,ijb_x=',
51	! & iq,ijb_x
52
53	c calcul de la pente a droite et a gauche de la maille
54
55	c ijb=ij_begin
56	c ije=ij_end
57
58	ijb=ijb_x
59	ije=ije_x
60
61	if (pole_nord.and.ijb==1) ijb=ijb+iip1
62	if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1
63
64	IF (pente_max.gt.-1.e-5) THEN
65	c IF (pente_max.gt.10) THEN
66
67	c calcul des pentes avec limitation, Van Leer scheme I:
68	c -----------------------------------------------------
69
70	c calcul de la pente aux points u
71
72	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
73	DO l = 1, llm
74	DO ij=ijb,ije-1
75	dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
76	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
77	c sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
78	ENDDO
79	DO ij=ijb+iip1-1,ije,iip1
80	dxqu(ij)=dxqu(ij-iim)
81	c sigu(ij)=sigu(ij-iim)
82	ENDDO
83
84	DO ij=ijb,ije
85	adxqu(ij)=abs(dxqu(ij))
86	ENDDO
87
88	c calcul de la pente maximum dans la maille en valeur absolue
89
90	DO ij=ijb+1,ije
91	dxqmax(ij,l)=pente_max*
92	, min(adxqu(ij-1),adxqu(ij))
93	c limitation subtile
94	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
95
96
97	ENDDO
98
99	DO ij=ijb+iip1-1,ije,iip1
100	dxqmax(ij-iim,l)=dxqmax(ij,l)
101	ENDDO
102
103	DO ij=ijb+1,ije
104	#ifdef CRAY
105	dxq(ij,l)=
106	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
107	#else
108	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
109	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
110	ELSE
111	c extremum local
112	dxq(ij,l)=0.
113	ENDIF
114	#endif
115	dxq(ij,l)=0.5*dxq(ij,l)
116	dxq(ij,l)=
117	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
118	ENDDO
119
120	ENDDO ! l=1,llm
121	c$OMP END DO NOWAIT
122
123	ELSE ! (pente_max.lt.-1.e-5)
124
125	c Pentes produits:
126	c ----------------
127	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
128	DO l = 1, llm
129	DO ij=ijb,ije-1
130	dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
131	ENDDO
132	DO ij=ijb+iip1-1,ije,iip1
133	dxqu(ij)=dxqu(ij-iim)
134	ENDDO
135
136	DO ij=ijb+1,ije
137	zz(ij)=dxqu(ij-1)*dxqu(ij)
138	zz(ij)=zz(ij)+zz(ij)
139	IF(zz(ij).gt.0) THEN
140	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
141	ELSE
142	c extremum local
143	dxq(ij,l)=0.
144	ENDIF
145	ENDDO
146
147	ENDDO
148	c$OMP END DO NOWAIT
149	ENDIF ! (pente_max.lt.-1.e-5)
150
151	c bouclage de la pente en iip1:
152	c -----------------------------
153	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
154	DO l=1,llm
155	DO ij=ijb+iip1-1,ije,iip1
156	dxq(ij-iim,l)=dxq(ij,l)
157	ENDDO
158
159	DO ij=ijb,ije
160	iadvplus(ij,l)=0
161	ENDDO
162
163	ENDDO
164	c$OMP END DO NOWAIT
165
166	if (pole_nord) THEN
167	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
168	DO l=1,llm
169	iadvplus(1:iip1,l)=0
170	ENDDO
171	c$OMP END DO NOWAIT
172	endif
173
174	if (pole_sud) THEN
175	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
176	DO l=1,llm
177	iadvplus(ip1jm+1:ip1jmp1,l)=0
178	ENDDO
179	c$OMP END DO NOWAIT
180	endif
181
182	c calcul des flux a gauche et a droite
183
184	#ifdef CRAY
185	c--pas encore modification sur Qsat
186	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
187	DO l=1,llm
188	DO ij=ijb,ije-1
189	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
190	, 1.+u_m(ij,l)/masse(ij+1,l,iq),
191	, u_m(ij,l))
192	zdum(ij,l)=0.5*zdum(ij,l)
193	u_mq(ij,l)=cvmgp(
194	, q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
195	, q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
196	, u_m(ij,l))
197	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
198	ENDDO
199	ENDDO
200	c$OMP END DO NOWAIT
201
202	#else
203	c on cumule le flux correspondant a toutes les mailles dont la masse
204	c au travers de la paroi pENDant le pas de temps.
205	c le rapport de melange de l''air advecte est min(q_vanleer, Qsat_downwind)
206	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
207	DO l=1,llm
208	DO ij=ijb,ije-1
209	IF (u_m(ij,l).gt.0.) THEN
210	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
211	u_mq(ij,l)=u_m(ij,l)*
212	$ min(q(ij,l,iq)+0.5zdum(ij,l)dxq(ij,l),qsat(ij+1,l))
213	ELSE
214	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
215	u_mq(ij,l)=u_m(ij,l)*
216	$ min(q(ij+1,l,iq)-0.5zdum(ij,l)dxq(ij+1,l),qsat(ij,l))
217	ENDIF
218	ENDDO
219	ENDDO
220	c$OMP END DO NOWAIT
221	#endif
222
223
224	c detection des points ou on advecte plus que la masse de la
225	c maille
226	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
227	DO l=1,llm
228	DO ij=ijb,ije-1
229	IF(zdum(ij,l).lt.0) THEN
230	iadvplus(ij,l)=1
231	u_mq(ij,l)=0.
232	ENDIF
233	ENDDO
234	ENDDO
235	c$OMP END DO NOWAIT
236
237	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
238	DO l=1,llm
239	DO ij=ijb+iip1-1,ije,iip1
240	iadvplus(ij,l)=iadvplus(ij-iim,l)
241	ENDDO
242	ENDDO
243	c$OMP END DO NOWAIT
244
245
246
247	c traitement special pour le cas ou on advecte en longitude plus que le
248	c contenu de la maille.
249	c cette partie est mal vectorisee.
250
251	c pas d'influence de la pression saturante (pour l'instant)
252
253	c calcul du nombre de maille sur lequel on advecte plus que la maille.
254
255	n0=0
256	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
257	DO l=1,llm
258	nl(l)=0
259	DO ij=ijb,ije
260	nl(l)=nl(l)+iadvplus(ij,l)
261	ENDDO
262	n0=n0+nl(l)
263	ENDDO
264	c$OMP END DO NOWAIT
265
266	cym ATTENTION ICI en OpenMP reduction pas forcement necessaire
267	cym IF(n0.gt.1) THEN
268	cym IF(n0.gt.0) THEN
269	ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
270	ccc & ,'contenu de la maille : ',n0
271	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
272	DO l=1,llm
273	IF(nl(l).gt.0) THEN
274	iju=0
275	c indicage des mailles concernees par le traitement special
276	DO ij=ijb,ije
277	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
278	iju=iju+1
279	indu(iju)=ij
280	ENDIF
281	ENDDO
282	niju=iju
283	!PRINT*,'vlxqs 280: niju,nl',niju,nl(l)
284
285	c traitement des mailles
286	DO iju=1,niju
287	ij=indu(iju)
288	j=(ij-1)/iip1+1
289	zu_m=u_m(ij,l)
290	u_mq(ij,l)=0.
291	IF(zu_m.gt.0.) THEN
292	ijq=ij
293	i=ijq-(j-1)*iip1
294	c accumulation pour les mailles completements advectees
295	do while(zu_m.gt.masse(ijq,l,iq))
296	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
297	& *masse(ijq,l,iq)
298	zu_m=zu_m-masse(ijq,l,iq)
299	i=mod(i-2+iim,iim)+1
300	ijq=(j-1)*iip1+i
301	ENDDO
302	c ajout de la maille non completement advectee
303	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)
304	& +0.5(1.-zu_m/masse(ijq,l,iq))dxq(ijq,l))
305	ELSE
306	ijq=ij+1
307	i=ijq-(j-1)*iip1
308	c accumulation pour les mailles completements advectees
309	do while(-zu_m.gt.masse(ijq,l,iq))
310	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
311	& *masse(ijq,l,iq)
312	zu_m=zu_m+masse(ijq,l,iq)
313	i=mod(i,iim)+1
314	ijq=(j-1)*iip1+i
315	ENDDO
316	c ajout de la maille non completement advectee
317	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
318	& 0.5(1.+zu_m/masse(ijq,l,iq))dxq(ijq,l))
319	ENDIF
320	ENDDO
321	ENDIF
322	ENDDO
323	c$OMP END DO NOWAIT
324	cym ENDIF ! n0.gt.0
325
326
327
328	c bouclage en latitude
329	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
330	DO l=1,llm
331	DO ij=ijb+iip1-1,ije,iip1
332	u_mq(ij,l)=u_mq(ij-iim,l)
333	ENDDO
334	ENDDO
335	c$OMP END DO NOWAIT
336
337	! CRisi: appel recursif de l'advection sur les fils.
338	! Il faut faire ca avant d'avoir mis a jour q et masse
339	!write(,) 'vlspltqs 336: iq,ijb_x,nqChildren(iq)=',
340	! & iq,ijb_x,tracers(iq)%nqChildren
341
342	do ifils=1,tracers(iq)%nqDescen
343	iq2=tracers(iq)%iqDescen(ifils)
344	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
345	DO l=1,llm
346	DO ij=ijb,ije
347	!MVals: veiller a ce qu'on n'ait pas de denominateur nul
348	masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
349	if (q(ij,l,iq).gt.min_qParent) then ! modif 13 nov 2020
350	Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
351	else
352	Ratio(ij,l,iq2)=min_ratio
353	endif
354	enddo
355	enddo
356	c$OMP END DO NOWAIT
357	enddo
358	do ifils=1,tracers(iq)%nqChildren
359	iq2=tracers(iq)%iqDescen(ifils)
360	!write(,) 'vlxqs 349: on appelle vlx pour iq2=',iq2
361	call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2)
362	enddo
363	! end CRisi
364
365	!write(,) 'vlspltqs 360: iq,ijb_x=',iq,ijb_x
366
367	c calcul des tendances
368	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
369	DO l=1,llm
370	DO ij=ijb+1,ije
371	!MVals: veiller a ce qu'on n'ait pas de denominateur nul
372	new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),min_qMass)
373	q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
374	& u_mq(ij-1,l)-u_mq(ij,l))
375	& /new_m
376	masse(ij,l,iq)=new_m
377	ENDDO
378	c Modif Fred 22 03 96 correction d''un bug (les scopy ci-dessous)
379	DO ij=ijb+iip1-1,ije,iip1
380	q(ij-iim,l,iq)=q(ij,l,iq)
381	masse(ij-iim,l,iq)=masse(ij,l,iq)
382	ENDDO
383	ENDDO
384	c$OMP END DO NOWAIT
385
386	!write(,) 'vlspltqs 380: iq,ijb_x=',iq,ijb_x
387
388	! retablir les fils en rapport de melange par rapport a l'air:
389	do ifils=1,tracers(iq)%nqDescen
390	iq2=tracers(iq)%iqDescen(ifils)
391	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
392	DO l=1,llm
393	DO ij=ijb+1,ije
394	q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
395	enddo
396	DO ij=ijb+iip1-1,ije,iip1
397	q(ij-iim,l,iq2)=q(ij,l,iq2)
398	enddo ! DO ij=ijb+iip1-1,ije,iip1
399	enddo
400	c$OMP END DO NOWAIT
401	enddo
402
403	!write(,) 'vlspltqs 399: iq,ijb_x=',iq,ijb_x
404
405	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
406	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1,iq),iip1,masse(iip2,1,iq),iip1)
407
408
409	RETURN
410	END
411	SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat,iq)
412	c
413	c Auteurs: P.Le Van, F.Hourdin, F.Forget
414	c
415	c ********************************************************************
416	c Shema d'advection " pseudo amont " .
417	c ********************************************************************
418	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
419	c qsat est un argument de sortie pour le s-pg ....
420	c
421	c
422	c --------------------------------------------------------------------
423	USE parallel_lmdz
424	USE infotrac, ONLY : nqtot,tracers, ! CRisi &
425	& min_qParent,min_qMass,min_ratio ! MVals et CRisi
426	USE comconst_mod, ONLY: pi
427	IMPLICIT NONE
428	c
429	include "dimensions.h"
430	include "paramet.h"
431	include "comgeom.h"
432	include "iniprint.h"
433	c
434	c
435	c Arguments:
436	c ----------
437	REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
438	REAL masse_adv_v( ijb_v:ije_v,llm)
439	REAL q(ijb_u:ije_u,llm,nqtot)
440	REAL qsat(ijb_u:ije_u,llm)
441	INTEGER iq ! CRisi
442	c
443	c Local
444	c ---------
445	c
446	INTEGER i,ij,l
447	c
448	REAL airej2,airejjm,airescb(iim),airesch(iim)
449	REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v)
450	REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u)
451	REAL qbyv(ijb_v:ije_v,llm,nqtot)
452
453	REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
454	c REAL newq,oldmasse
455	Logical first
456	SAVE first
457	c$OMP THREADPRIVATE(first)
458	REAL convpn,convps,convmpn,convmps
459	REAL sinlon(iip1),sinlondlon(iip1)
460	REAL coslon(iip1),coslondlon(iip1)
461	SAVE sinlon,coslon,sinlondlon,coslondlon
462	SAVE airej2,airejjm
463	c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
464	c$OMP THREADPRIVATE(airej2,airejjm)
465	c
466	c
467	REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
468	INTEGER ifils,iq2 ! CRisi
469
470	REAL SSUM
471
472	DATA first/.true./
473	INTEGER ijb,ije
474	INTEGER ijbm,ijem
475
476	ijb=ij_begin-2*iip1
477	ije=ij_end+2*iip1
478	if (pole_nord) ijb=ij_begin
479	if (pole_sud) ije=ij_end
480	ij=3525
481	l=3
482	if ((ij.ge.ijb).and.(ij.le.ije)) then
483	!write(,) 'vlyqs 480: ij,l,iq,ijb,q(ij,l,:)=',
484	! & ij,l,iq,ijb,q(ij,l,:)
485	endif
486
487	IF(first) THEN
488	PRINT*,'Shema Amont nouveau appele dans Vanleer '
489	PRINT*,'vlyqs_loc, iq=',iq
490	first=.false.
491	do i=2,iip1
492	coslon(i)=cos(rlonv(i))
493	sinlon(i)=sin(rlonv(i))
494	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
495	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
496	ENDDO
497	coslon(1)=coslon(iip1)
498	coslondlon(1)=coslondlon(iip1)
499	sinlon(1)=sinlon(iip1)
500	sinlondlon(1)=sinlondlon(iip1)
501	airej2 = SSUM( iim, aire(iip2), 1 )
502	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
503	ENDIF
504
505	c
506
507	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
508	DO l = 1, llm
509	c
510	c --------------------------------
511	c CALCUL EN LATITUDE
512	c --------------------------------
513
514	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
515	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
516	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
517
518	if (pole_nord) then
519	DO i = 1, iim
520	airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
521	ENDDO
522	qpns = SSUM( iim, airescb ,1 ) / airej2
523	endif
524
525	if (pole_sud) then
526	DO i = 1, iim
527	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
528	ENDDO
529	qpsn = SSUM( iim, airesch ,1 ) / airejjm
530	endif
531
532
533	c calcul des pentes aux points v
534
535	ijb=ij_begin-2*iip1
536	ije=ij_end+iip1
537	if (pole_nord) ijb=ij_begin
538	if (pole_sud) ije=ij_end-iip1
539
540	DO ij=ijb,ije
541	dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
542	adyqv(ij)=abs(dyqv(ij))
543	ENDDO
544
545
546	c calcul des pentes aux points scalaires
547
548	ijb=ij_begin-iip1
549	ije=ij_end+iip1
550	if (pole_nord) ijb=ij_begin+iip1
551	if (pole_sud) ije=ij_end-iip1
552
553	DO ij=ijb,ije
554	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
555	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
556	dyqmax(ij)=pente_max*dyqmax(ij)
557	ENDDO
558
559	IF (pole_nord) THEN
560
561	c calcul des pentes aux poles
562	DO ij=1,iip1
563	dyq(ij,l)=qpns-q(ij+iip1,l,iq)
564	ENDDO
565
566	c filtrage de la derivee
567	dyn1=0.
568	dyn2=0.
569	DO ij=1,iim
570	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
571	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
572	ENDDO
573	DO ij=1,iip1
574	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
575	ENDDO
576
577	c calcul des pentes limites aux poles
578	fn=1.
579	DO ij=1,iim
580	IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
581	fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
582	ENDIF
583	ENDDO
584
585	DO ij=1,iip1
586	dyq(ij,l)=fn*dyq(ij,l)
587	ENDDO
588
589	ENDIF
590
591	IF (pole_sud) THEN
592
593	DO ij=1,iip1
594	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
595	ENDDO
596
597	dys1=0.
598	dys2=0.
599
600	DO ij=1,iim
601	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
602	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
603	ENDDO
604
605	DO ij=1,iip1
606	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
607	ENDDO
608
609	c calcul des pentes limites aux poles
610	fs=1.
611	DO ij=1,iim
612	IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
613	fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
614	ENDIF
615	ENDDO
616
617	DO ij=1,iip1
618	dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
619	ENDDO
620
621	ENDIF
622
623
624	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
625	C En memoire de dIFferents tests sur la
626	C limitation des pentes aux poles.
627	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
628	C PRINT*,dyq(1)
629	C PRINT*,dyqv(iip1+1)
630	C appn=abs(dyq(1)/dyqv(iip1+1))
631	C PRINT*,dyq(ip1jm+1)
632	C PRINT*,dyqv(ip1jm-iip1+1)
633	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
634	C DO ij=2,iim
635	C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
636	C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
637	C ENDDO
638	C appn=min(pente_max/appn,1.)
639	C apps=min(pente_max/apps,1.)
640	C
641	C
642	C cas ou on a un extremum au pole
643	C
644	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
645	C & appn=0.
646	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
647	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
648	C & apps=0.
649	C
650	C limitation des pentes aux poles
651	C DO ij=1,iip1
652	C dyq(ij)=appn*dyq(ij)
653	C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
654	C ENDDO
655	C
656	C test
657	C DO ij=1,iip1
658	C dyq(iip1+ij)=0.
659	C dyq(ip1jm+ij-iip1)=0.
660	C ENDDO
661	C DO ij=1,ip1jmp1
662	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
663	C ENDDO
664	C
665	C changement 10 07 96
666	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
667	C & THEN
668	C DO ij=1,iip1
669	C dyqmax(ij)=0.
670	C ENDDO
671	C ELSE
672	C DO ij=1,iip1
673	C dyqmax(ij)=pente_max*abs(dyqv(ij))
674	C ENDDO
675	C ENDIF
676	C
677	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
678	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
679	C &THEN
680	C DO ij=ip1jm+1,ip1jmp1
681	C dyqmax(ij)=0.
682	C ENDDO
683	C ELSE
684	C DO ij=ip1jm+1,ip1jmp1
685	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
686	C ENDDO
687	C ENDIF
688	C fin changement 10 07 96
689	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
690
691	c calcul des pentes limitees
692	ijb=ij_begin-iip1
693	ije=ij_end+iip1
694	if (pole_nord) ijb=ij_begin+iip1
695	if (pole_sud) ije=ij_end-iip1
696
697	DO ij=ijb,ije
698	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
699	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
700	ELSE
701	dyq(ij,l)=0.
702	ENDIF
703	ENDDO
704
705	ENDDO
706	c$OMP END DO NOWAIT
707
708	ijb=ij_begin-iip1
709	ije=ij_end
710	if (pole_nord) ijb=ij_begin
711	if (pole_sud) ije=ij_end-iip1
712
713	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
714	DO l=1,llm
715	DO ij=ijb,ije
716	IF( masse_adv_v(ij,l).GT.0. ) THEN
717	qbyv(ij,l,iq)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq ) +
718	, dyq(ij+iip1,l)0.5(1.-masse_adv_v(ij,l)
719	, /masse(ij+iip1,l,iq)))
720	ELSE
721	qbyv(ij,l,iq)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
722	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
723	ENDIF
724	qbyv(ij,l,iq) = masse_adv_v(ij,l)*qbyv(ij,l,iq)
725	ENDDO
726	ENDDO
727	c$OMP END DO NOWAIT
728
729	! CRisi: appel recursif de l'advection sur les fils.
730	! Il faut faire ca avant d'avoir mis a jour q et masse
731	! write(,)'vlyqs 689: iq,nqChildren(iq)=',iq,
732	! & tracers(iq)%nqChildren
733
734	ijb=ij_begin-2*iip1
735	ije=ij_end+2*iip1
736	ijbm=ij_begin-iip1
737	ijem=ij_end+iip1
738	if (pole_nord) ijb=ij_begin
739	if (pole_sud) ije=ij_end
740	if (pole_nord) ijbm=ij_begin
741	if (pole_sud) ijem=ij_end
742
743	!write(lunout,*) 'vlspltqs 737: iq,ijb,ije=',iq,ijb,ije
744	!write(lunout,*) 'ij_begin,ij_end=',ij_begin,ij_end
745	!write(lunout,*) 'pole_nord,pole_sud=',pole_nord,pole_sud
746	do ifils=1,tracers(iq)%nqDescen
747	iq2=tracers(iq)%iqDescen(ifils)
748	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
749	DO l=1,llm
750	! modif des bornes: CRisi 16 nov 2020
751	! d'abord masse avec bornes corrigees
752	DO ij=ijbm,ijem
753	!MVals: veiller a ce qu'on n'ait pas de denominateur nul
754	masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
755	enddo !DO ij=ijbm,ijem
756
757	! ensuite Ratio avec anciennes bornes
758	DO ij=ijb,ije
759	!MVals: veiller a ce qu'on n'ait pas de denominateur nul
760	!write(lunout,*) 'ij,l,q(ij,l,iq)=',ij,l,q(ij,l,iq)
761	if (q(ij,l,iq).gt.min_qParent) then ! modif 13 nov 2020
762	Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
763	else
764	Ratio(ij,l,iq2)=min_ratio
765	endif
766	enddo !DO ij=ijbm,ijem
767	enddo !DO l=1,llm
768	c$OMP END DO NOWAIT
769	enddo
770	do ifils=1,tracers(iq)%nqChildren
771	iq2=tracers(iq)%iqDescen(ifils)
772	!write(lunout,*) 'vly: appel recursiv vly iq2=',iq2
773	call vly_loc(Ratio,pente_max,masse,qbyv,iq2)
774	enddo
775
776
777	! end CRisi
778
779	ijb=ij_begin
780	ije=ij_end
781	if (pole_nord) ijb=ij_begin+iip1
782	if (pole_sud) ije=ij_end-iip1
783
784	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
785	DO l=1,llm
786	DO ij=ijb,ije
787	newmasse=masse(ij,l,iq)
788	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
789	q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l,iq)
790	& -qbyv(ij-iip1,l,iq))/newmasse
791	masse(ij,l,iq)=newmasse
792	ENDDO
793	c.-. ancienne version
794
795	IF (pole_nord) THEN
796
797	convpn=SSUM(iim,qbyv(1,l,iq),1)/apoln
798	convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
799	DO ij = 1,iip1
800	newmasse=masse(ij,l,iq)+convmpn*aire(ij)
801	q(ij,l,iq)=(q(ij,l,iq)masse(ij,l,iq)+convpnaire(ij))/
802	& newmasse
803	masse(ij,l,iq)=newmasse
804	ENDDO
805
806	ENDIF
807
808	IF (pole_sud) THEN
809
810	convps = -SSUM(iim,qbyv(ip1jm-iim,l,iq),iq,1)/apols
811	convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
812	DO ij = ip1jm+1,ip1jmp1
813	newmasse=masse(ij,l,iq)+convmps*aire(ij)
814	q(ij,l,iq)=(q(ij,l,iq)masse(ij,l,iq)+convpsaire(ij))/
815	& newmasse
816	masse(ij,l,iq)=newmasse
817	ENDDO
818
819	ENDIF
820	c.-. fin ancienne version
821
822	c._. nouvelle version
823	c convpn=SSUM(iim,qbyv(1,l,iq),1)
824	c convmpn=ssum(iim,masse_adv_v(1,l),1)
825	c oldmasse=ssum(iim,masse(1,l,iq),1)
826	c newmasse=oldmasse+convmpn
827	c newq=(q(1,l,iq)*oldmasse+convpn)/newmasse
828	c newmasse=newmasse/apoln
829	c DO ij = 1,iip1
830	c q(ij,l,iq)=newq
831	c masse(ij,l,iq)=newmasse*aire(ij)
832	c ENDDO
833	c convps=-SSUM(iim,qbyv(ip1jm-iim,l,iq),1)
834	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
835	c oldmasse=ssum(iim,masse(ip1jm-iim,l,iq),1)
836	c newmasse=oldmasse+convmps
837	c newq=(q(ip1jmp1,l,iq)*oldmasse+convps)/newmasse
838	c newmasse=newmasse/apols
839	c DO ij = ip1jm+1,ip1jmp1
840	c q(ij,l,iq)=newq
841	c masse(ij,l,iq)=newmasse*aire(ij)
842	c ENDDO
843	c._. fin nouvelle version
844	ENDDO
845	c$OMP END DO NOWAIT
846
847	! retablir les fils en rapport de melange par rapport a l'air:
848	ijb=ij_begin
849	ije=ij_end
850	! if (pole_nord) ijb=ij_begin+iip1
851	! if (pole_sud) ije=ij_end-iip1
852
853	do ifils=1,tracers(iq)%nqDescen
854	iq2=tracers(iq)%iqDescen(ifils)
855	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
856	DO l=1,llm
857	DO ij=ijb,ije
858	q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
859	enddo
860	enddo
861	c$OMP END DO NOWAIT
862	enddo
863
864
865	RETURN
866	END

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

Download in other formats:

Original Format