Context Navigation

vlspltqs_loc.F @ 4284

Last change on this file since 4284 was 4143, checked in by dcugnet, 3 years ago
Some variables are renamed or replaced by direct equivalents: iso_indnum -> tracers(:)%iso_iName niso_possibles -> niso iqiso -> iqIsoPha ; index_trac -> itZonIso ok_iso_verif -> isoCheck ntraceurs_zone -> nzone ; ntraciso -> ntiso qperemin -> min_qparent ; masseqmin -> min_qmass ; ratiomin -> min_ratio Some renamed variables are only aliased with the older name (using USE <module>, ONLY: <oldName> => <newName>) in routines where they are repeated many times. Few hard-coded indexes are now computed (examples: ilic, iso, ivap, irneb, iq_vap, iq_liq, iso_H2O, iso_HDO, iso_HTO, iso_O17, iso_O18). The IF(isoCheck) test is now embedded in the check_isotopes_seq and check_isotopes_loc routines (lighter calling).
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 4143 2022-05-09 10:35:40Z musat $
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,nqChilds(iq)=',
340	! & iq,ijb_x,tracers(iq)%nqChilds
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)%nqChilds
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,nqChilds(iq)=',iq,tracers(iq)%nqChilds
732
733	ijb=ij_begin-2*iip1
734	ije=ij_end+2*iip1
735	ijbm=ij_begin-iip1
736	ijem=ij_end+iip1
737	if (pole_nord) ijb=ij_begin
738	if (pole_sud) ije=ij_end
739	if (pole_nord) ijbm=ij_begin
740	if (pole_sud) ijem=ij_end
741
742	!write(lunout,*) 'vlspltqs 737: iq,ijb,ije=',iq,ijb,ije
743	!write(lunout,*) 'ij_begin,ij_end=',ij_begin,ij_end
744	!write(lunout,*) 'pole_nord,pole_sud=',pole_nord,pole_sud
745	do ifils=1,tracers(iq)%nqDescen
746	iq2=tracers(iq)%iqDescen(ifils)
747	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
748	DO l=1,llm
749	! modif des bornes: CRisi 16 nov 2020
750	! d'abord masse avec bornes corrigees
751	DO ij=ijbm,ijem
752	!MVals: veiller a ce qu'on n'ait pas de denominateur nul
753	masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
754	enddo !DO ij=ijbm,ijem
755
756	! ensuite Ratio avec anciennes bornes
757	DO ij=ijb,ije
758	!MVals: veiller a ce qu'on n'ait pas de denominateur nul
759	!write(lunout,*) 'ij,l,q(ij,l,iq)=',ij,l,q(ij,l,iq)
760	if (q(ij,l,iq).gt.min_qParent) then ! modif 13 nov 2020
761	Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
762	else
763	Ratio(ij,l,iq2)=min_ratio
764	endif
765	enddo !DO ij=ijbm,ijem
766	enddo !DO l=1,llm
767	c$OMP END DO NOWAIT
768	enddo
769	do ifils=1,tracers(iq)%nqChilds
770	iq2=tracers(iq)%iqDescen(ifils)
771	!write(lunout,*) 'vly: appel recursiv vly iq2=',iq2
772	call vly_loc(Ratio,pente_max,masse,qbyv,iq2)
773	enddo
774
775
776	! end CRisi
777
778	ijb=ij_begin
779	ije=ij_end
780	if (pole_nord) ijb=ij_begin+iip1
781	if (pole_sud) ije=ij_end-iip1
782
783	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
784	DO l=1,llm
785	DO ij=ijb,ije
786	newmasse=masse(ij,l,iq)
787	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
788	q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l,iq)
789	& -qbyv(ij-iip1,l,iq))/newmasse
790	masse(ij,l,iq)=newmasse
791	ENDDO
792	c.-. ancienne version
793
794	IF (pole_nord) THEN
795
796	convpn=SSUM(iim,qbyv(1,l,iq),1)/apoln
797	convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
798	DO ij = 1,iip1
799	newmasse=masse(ij,l,iq)+convmpn*aire(ij)
800	q(ij,l,iq)=(q(ij,l,iq)masse(ij,l,iq)+convpnaire(ij))/
801	& newmasse
802	masse(ij,l,iq)=newmasse
803	ENDDO
804
805	ENDIF
806
807	IF (pole_sud) THEN
808
809	convps = -SSUM(iim,qbyv(ip1jm-iim,l,iq),iq,1)/apols
810	convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
811	DO ij = ip1jm+1,ip1jmp1
812	newmasse=masse(ij,l,iq)+convmps*aire(ij)
813	q(ij,l,iq)=(q(ij,l,iq)masse(ij,l,iq)+convpsaire(ij))/
814	& newmasse
815	masse(ij,l,iq)=newmasse
816	ENDDO
817
818	ENDIF
819	c.-. fin ancienne version
820
821	c._. nouvelle version
822	c convpn=SSUM(iim,qbyv(1,l,iq),1)
823	c convmpn=ssum(iim,masse_adv_v(1,l),1)
824	c oldmasse=ssum(iim,masse(1,l,iq),1)
825	c newmasse=oldmasse+convmpn
826	c newq=(q(1,l,iq)*oldmasse+convpn)/newmasse
827	c newmasse=newmasse/apoln
828	c DO ij = 1,iip1
829	c q(ij,l,iq)=newq
830	c masse(ij,l,iq)=newmasse*aire(ij)
831	c ENDDO
832	c convps=-SSUM(iim,qbyv(ip1jm-iim,l,iq),1)
833	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
834	c oldmasse=ssum(iim,masse(ip1jm-iim,l,iq),1)
835	c newmasse=oldmasse+convmps
836	c newq=(q(ip1jmp1,l,iq)*oldmasse+convps)/newmasse
837	c newmasse=newmasse/apols
838	c DO ij = ip1jm+1,ip1jmp1
839	c q(ij,l,iq)=newq
840	c masse(ij,l,iq)=newmasse*aire(ij)
841	c ENDDO
842	c._. fin nouvelle version
843	ENDDO
844	c$OMP END DO NOWAIT
845
846	! retablir les fils en rapport de melange par rapport a l'air:
847	ijb=ij_begin
848	ije=ij_end
849	! if (pole_nord) ijb=ij_begin+iip1
850	! if (pole_sud) ije=ij_end-iip1
851
852	do ifils=1,tracers(iq)%nqDescen
853	iq2=tracers(iq)%iqDescen(ifils)
854	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
855	DO l=1,llm
856	DO ij=ijb,ije
857	q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
858	enddo
859	enddo
860	c$OMP END DO NOWAIT
861	enddo
862
863
864	RETURN
865	END

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

Download in other formats:

Original Format