Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

vlspltqs_loc.F90 @ 5248

Last change on this file since 5248 was 5246, checked in by abarral, 21 hours ago
Convert fixed-form to free-form sources .F -> .{f,F}90 (WIP: some .F remain, will be handled in subsequent commits)
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: 22.4 KB

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

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

Download in other formats:

Original Format