Context Navigation

vlspltqs_loc.F @ 3863

Last change on this file since 3863 was 3852, checked in by dcugnet, 4 years ago

Extension of the tracers management.

The tracers files can be:

1) "traceur.def": old format, with:

the number of tracers on the first line
one line for each tracer: <tracer name> <hadv> <vadv> [<parent name>]

2) "tracer.def": new format with one section each model component.
3) "tracer_<name>.def": new format with a single section.

The formats 2 and 3 reading is driven by the "type_trac" key, which can be a

coma-separated list of components.

Format 2: read the sections from the "tracer.def" file.
format 3: read one section each "tracer_<section name>.def" file.
the first line of a section is "&<section name>
the other lines start with a tracer name followed by <key>=<val> pairs.
the "default" tracer name is reserved ; the other tracers of the section inherit its <key>=<val>, except for the keys that are redefined locally.

This format helps keeping the tracers files compact, thanks to the "default"
special tracer and the three levels of factorization:

on the tracers names: a tracer name can be a coma-separated list of tracers => all the tracers of the list have the same <key>=<val> properties
on the parents names: the value of the "parent" property can be a coma-separated list of tracers => only possible for geographic tagging tracers
on the phases: the property "phases" is [g](l][s] (gas/liquid/solid)

Read information is stored in the vector "tracers(:)", of derived type "tra".

"isotopes_params.def" is a similar file, with one section each isotopes family.
It contains a database of isotopes properties ; if there are second generation
tracers (isotopes), the corresponding sections are read.

Read information is stored in the vector "isotopes(:)", of derived type "iso".

The "getKey" function helps to get the values of the parameters stored in
"tracers" or "isotopes".

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.7 KB

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

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

Download in other formats:

Original Format