Context Navigation

advn.f90 @ 5134

Last change on this file since 5134 was 5134, checked in by abarral, 4 months ago
Replace academic.h, alpale.h, comdissip.h, comdissipn.h, comdissnew.h by modules Remove unused clesph0.h
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:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 25.1 KB

Line
1	! $Header$
2
3	SUBROUTINE advn(q, masse, w, pbaru, pbarv, pdt, mode)
4	!
5	! Auteur : F. Hourdin
6	!
7	! ********************************************************************
8	! Shema d'advection " pseudo amont " .
9	! ********************************************************************
10	! q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
11	!
12	! pbaru,pbarv,w flux de masse en u ,v ,w
13	! pdt pas de temps
14	!
15	! --------------------------------------------------------------------
16	USE lmdz_iniprint, ONLY: lunout, prt_level
17	USE lmdz_ssum_scopy, ONLY: ssum
18
19	IMPLICIT NONE
20	!
21	INCLUDE "dimensions.h"
22	INCLUDE "paramet.h"
23	INCLUDE "comgeom.h"
24
25	!
26	! Arguments:
27	! ----------
28	INTEGER :: mode
29	REAL :: masse(ip1jmp1, llm)
30	REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
31	REAL :: q(ip1jmp1, llm)
32	REAL :: w(ip1jmp1, llm), pdt
33	!
34	! Local
35	! ---------
36	!
37	INTEGER :: i, ij, l, j, ii
38	INTEGER :: ijlqmin, iqmin, jqmin, lqmin
39	INTEGER :: ismin
40	!
41	REAL :: zm(ip1jmp1, llm), newmasse
42	REAL :: mu(ip1jmp1, llm)
43	REAL :: mv(ip1jm, llm)
44	REAL :: mw(ip1jmp1, llm + 1)
45	REAL :: zq(ip1jmp1, llm), zz, qpn, qps
46	REAL :: zqg(ip1jmp1, llm), zqd(ip1jmp1, llm)
47	REAL :: zqs(ip1jmp1, llm), zqn(ip1jmp1, llm)
48	REAL :: zqh(ip1jmp1, llm), zqb(ip1jmp1, llm)
49	REAL :: temps0, temps1, temps2, temps3
50	REAL :: ztemps1, ztemps2
51	save temps1, temps2, temps3
52	REAL :: zzpbar, zzw
53
54	REAL :: qmin, qmax
55	data qmin, qmax/0., 1./
56	data temps1, temps2, temps3/0., 0., 0./
57
58	zzpbar = 0.5 * pdt
59	zzw = pdt
60
61	DO l = 1, llm
62	DO ij = iip2, ip1jm
63	mu(ij, l) = pbaru(ij, l) * zzpbar
64	ENDDO
65	DO ij = 1, ip1jm
66	mv(ij, l) = pbarv(ij, l) * zzpbar
67	ENDDO
68	DO ij = 1, ip1jmp1
69	mw(ij, l) = w(ij, l) * zzw
70	ENDDO
71	ENDDO
72
73	DO ij = 1, ip1jmp1
74	mw(ij, llm + 1) = 0.
75	ENDDO
76
77	do l = 1, llm
78	qpn = 0.
79	qps = 0.
80	do ij = 1, iim
81	qpn = qpn + q(ij, l) * masse(ij, l)
82	qps = qps + q(ip1jm + ij, l) * masse(ip1jm + ij, l)
83	enddo
84	qpn = qpn / ssum(iim, masse(1, l), 1)
85	qps = qps / ssum(iim, masse(ip1jm + 1, l), 1)
86	do ij = 1, iip1
87	q(ij, l) = qpn
88	q(ip1jm + ij, l) = qps
89	enddo
90	enddo
91
92	do ij = 1, ip1jmp1
93	mw(ij, llm + 1) = 0.
94	enddo
95	do l = 1, llm
96	do ij = 1, ip1jmp1
97	zq(ij, l) = q(ij, l)
98	zm(ij, l) = masse(ij, l)
99	enddo
100	enddo
101
102	! CALL minmaxq(zq,qmin,qmax,'avant vlx ')
103	CALL advnqx(zq, zqg, zqd)
104	CALL advnx(zq, zqg, zqd, zm, mu, mode)
105	CALL advnqy(zq, zqs, zqn)
106	CALL advny(zq, zqs, zqn, zm, mv)
107	CALL advnqz(zq, zqh, zqb)
108	CALL advnz(zq, zqh, zqb, zm, mw)
109	! CALL vlz(zq,0.,zm,mw)
110	CALL advnqy(zq, zqs, zqn)
111	CALL advny(zq, zqs, zqn, zm, mv)
112	CALL advnqx(zq, zqg, zqd)
113	CALL advnx(zq, zqg, zqd, zm, mu, mode)
114	! CALL minmaxq(zq,qmin,qmax,'apres vlx ')
115
116	do l = 1, llm
117	do ij = 1, ip1jmp1
118	q(ij, l) = zq(ij, l)
119	enddo
120	do ij = 1, ip1jm + 1, iip1
121	q(ij + iim, l) = q(ij, l)
122	enddo
123	enddo
124
125	RETURN
126	END SUBROUTINE advn
127
128	SUBROUTINE advnqx(q, qg, qd)
129	!
130	! Auteurs: Calcul des valeurs de q aux point u.
131	!
132	! --------------------------------------------------------------------
133	USE lmdz_iniprint, ONLY: lunout, prt_level
134	IMPLICIT NONE
135	!
136	INCLUDE "dimensions.h"
137	INCLUDE "paramet.h"
138	!
139	!
140	! Arguments:
141	! ----------
142	REAL :: q(ip1jmp1, llm), qg(ip1jmp1, llm), qd(ip1jmp1, llm)
143	!
144	! Local
145	! ---------
146	!
147	INTEGER :: ij, l
148	!
149	REAL :: dxqu(ip1jmp1), zqu(ip1jmp1)
150	REAL :: zqmax(ip1jmp1), zqmin(ip1jmp1)
151	LOGICAL :: extremum(ip1jmp1)
152
153	INTEGER :: mode
154	save mode
155	data mode/1/
156
157	! calcul des pentes en u:
158	! -----------------------
159	IF (mode==0) THEN
160	do l = 1, llm
161	do ij = 1, ip1jm
162	qd(ij, l) = q(ij, l)
163	qg(ij, l) = q(ij, l)
164	enddo
165	enddo
166	else
167	do l = 1, llm
168	do ij = iip2, ip1jm - 1
169	dxqu(ij) = q(ij + 1, l) - q(ij, l)
170	zqu(ij) = 0.5 * (q(ij + 1, l) + q(ij, l))
171	enddo
172	do ij = iip1 + iip1, ip1jm, iip1
173	dxqu(ij) = dxqu(ij - iim)
174	zqu(ij) = zqu(ij - iim)
175	enddo
176	do ij = iip2, ip1jm - 1
177	zqu(ij) = zqu(ij) - dxqu(ij + 1) / 12.
178	enddo
179	do ij = iip1 + iip1, ip1jm, iip1
180	zqu(ij) = zqu(ij - iim)
181	enddo
182	do ij = iip2 + 1, ip1jm
183	zqu(ij) = zqu(ij) + dxqu(ij - 1) / 12.
184	enddo
185	do ij = iip1 + iip1, ip1jm, iip1
186	zqu(ij - iim) = zqu(ij)
187	enddo
188
189	! calcul des valeurs max et min acceptees aux interfaces
190
191	do ij = iip2, ip1jm - 1
192	zqmax(ij) = max(q(ij + 1, l), q(ij, l))
193	zqmin(ij) = min(q(ij + 1, l), q(ij, l))
194	enddo
195	do ij = iip1 + iip1, ip1jm, iip1
196	zqmax(ij) = zqmax(ij - iim)
197	zqmin(ij) = zqmin(ij - iim)
198	enddo
199	do ij = iip2 + 1, ip1jm
200	extremum(ij) = dxqu(ij) * dxqu(ij - 1)<=0.
201	enddo
202	do ij = iip1 + iip1, ip1jm, iip1
203	extremum(ij - iim) = extremum(ij)
204	enddo
205	do ij = iip2, ip1jm
206	zqu(ij) = min(max(zqmin(ij), zqu(ij)), zqmax(ij))
207	enddo
208	do ij = iip2 + 1, ip1jm
209	IF(extremum(ij)) THEN
210	qg(ij, l) = q(ij, l)
211	qd(ij, l) = q(ij, l)
212	else
213	qd(ij, l) = zqu(ij)
214	qg(ij, l) = zqu(ij - 1)
215	endif
216	enddo
217	do ij = iip1 + iip1, ip1jm, iip1
218	qd(ij - iim, l) = qd(ij, l)
219	qg(ij - iim, l) = qg(ij, l)
220	enddo
221
222	goto 8888
223
224	do ij = iip2 + 1, ip1jm
225	IF(extremum(ij).and..not.extremum(ij - 1)) &
226	qd(ij - 1, l) = q(ij, l)
227	enddo
228
229	do ij = iip1 + iip1, ip1jm, iip1
230	qd(ij - iim, l) = qd(ij, l)
231	enddo
232	do ij = iip2, ip1jm - 1
233	IF (extremum(ij).and..not.extremum(ij + 1)) &
234	qg(ij + 1, l) = q(ij, l)
235	enddo
236
237	do ij = iip1 + iip1, ip1jm, iip1
238	qg(ij, l) = qg(ij - iim, l)
239	enddo
240	8888 continue
241	enddo
242	ENDIF
243	RETURN
244	END SUBROUTINE advnqx
245	SUBROUTINE advnqy(q, qs, qn)
246	!
247	! Auteurs: Calcul des valeurs de q aux point v.
248	!
249	! --------------------------------------------------------------------
250	USE lmdz_iniprint, ONLY: lunout, prt_level
251	IMPLICIT NONE
252	!
253	INCLUDE "dimensions.h"
254	INCLUDE "paramet.h"
255	!
256	!
257	! Arguments:
258	! ----------
259	REAL :: q(ip1jmp1, llm), qs(ip1jmp1, llm), qn(ip1jmp1, llm)
260	!
261	! Local
262	! ---------
263	!
264	INTEGER :: ij, l
265	!
266	REAL :: dyqv(ip1jm), zqv(ip1jm, llm)
267	REAL :: zqmax(ip1jm), zqmin(ip1jm)
268	LOGICAL :: extremum(ip1jmp1)
269
270	INTEGER :: mode
271	save mode
272	data mode/1/
273
274	IF (mode==0) THEN
275	do l = 1, llm
276	do ij = 1, ip1jmp1
277	qn(ij, l) = q(ij, l)
278	qs(ij, l) = q(ij, l)
279	enddo
280	enddo
281	else
282
283	! calcul des pentes en u:
284	! -----------------------
285	do l = 1, llm
286	do ij = 1, ip1jm
287	dyqv(ij) = q(ij, l) - q(ij + iip1, l)
288	enddo
289
290	do ij = iip2, ip1jm - iip1
291	zqv(ij, l) = 0.5 * (q(ij + iip1, l) + q(ij, l))
292	zqv(ij, l) = zqv(ij, l) + (dyqv(ij + iip1) - dyqv(ij - iip1)) / 12.
293	enddo
294
295	do ij = iip2, ip1jm
296	extremum(ij) = dyqv(ij) * dyqv(ij - iip1)<=0.
297	enddo
298
299	! Pas de pentes aux poles
300	do ij = 1, iip1
301	zqv(ij, l) = q(ij, l)
302	zqv(ip1jm - iip1 + ij, l) = q(ip1jm + ij, l)
303	extremum(ij) = .TRUE.
304	extremum(ip1jmp1 - iip1 + ij) = .TRUE.
305	enddo
306
307	! calcul des valeurs max et min acceptees aux interfaces
308	do ij = 1, ip1jm
309	zqmax(ij) = max(q(ij + iip1, l), q(ij, l))
310	zqmin(ij) = min(q(ij + iip1, l), q(ij, l))
311	enddo
312
313	do ij = 1, ip1jm
314	zqv(ij, l) = min(max(zqmin(ij), zqv(ij, l)), zqmax(ij))
315	enddo
316
317	do ij = iip2, ip1jm
318	IF(extremum(ij)) THEN
319	qs(ij, l) = q(ij, l)
320	qn(ij, l) = q(ij, l)
321	! if (.NOT.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
322	! if (.NOT.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
323	else
324	qs(ij, l) = zqv(ij, l)
325	qn(ij, l) = zqv(ij - iip1, l)
326	endif
327	enddo
328
329	do ij = 1, iip1
330	qs(ij, l) = q(ij, l)
331	qn(ij, l) = q(ij, l)
332	qs(ip1jm + ij, l) = q(ip1jm + ij, l)
333	qn(ip1jm + ij, l) = q(ip1jm + ij, l)
334	enddo
335
336	enddo
337	ENDIF
338	RETURN
339	END SUBROUTINE advnqy
340
341	SUBROUTINE advnqz(q, qh, qb)
342	!
343	! Auteurs: Calcul des valeurs de q aux point v.
344	!
345	! --------------------------------------------------------------------
346	USE lmdz_iniprint, ONLY: lunout, prt_level
347	IMPLICIT NONE
348	!
349	INCLUDE "dimensions.h"
350	INCLUDE "paramet.h"
351	!
352	!
353	! Arguments:
354	! ----------
355	REAL :: q(ip1jmp1, llm), qh(ip1jmp1, llm), qb(ip1jmp1, llm)
356	!
357	! Local
358	! ---------
359	!
360	INTEGER :: ij, l
361	!
362	REAL :: dzqw(ip1jmp1, llm + 1), zqw(ip1jmp1, llm + 1)
363	REAL :: zqmax(ip1jmp1, llm), zqmin(ip1jmp1, llm)
364	LOGICAL :: extremum(ip1jmp1, llm)
365
366	INTEGER :: mode
367	save mode
368
369	data mode/1/
370
371	! calcul des pentes en u:
372	! -----------------------
373
374	IF (mode==0) THEN
375	do l = 1, llm
376	do ij = 1, ip1jmp1
377	qb(ij, l) = q(ij, l)
378	qh(ij, l) = q(ij, l)
379	enddo
380	enddo
381	else
382	do l = 2, llm
383	do ij = 1, ip1jmp1
384	dzqw(ij, l) = q(ij, l - 1) - q(ij, l)
385	zqw(ij, l) = 0.5 * (q(ij, l - 1) + q(ij, l))
386	enddo
387	enddo
388	do ij = 1, ip1jmp1
389	dzqw(ij, 1) = 0.
390	dzqw(ij, llm + 1) = 0.
391	enddo
392	do l = 2, llm
393	do ij = 1, ip1jmp1
394	zqw(ij, l) = zqw(ij, l) + (dzqw(ij, l + 1) - dzqw(ij, l - 1)) / 12.
395	enddo
396	enddo
397	do l = 2, llm - 1
398	do ij = 1, ip1jmp1
399	extremum(ij, l) = dzqw(ij, l) * dzqw(ij, l + 1)<=0.
400	enddo
401	enddo
402
403	! Pas de pentes en bas et en haut
404	do ij = 1, ip1jmp1
405	zqw(ij, 2) = q(ij, 1)
406	zqw(ij, llm) = q(ij, llm)
407	extremum(ij, 1) = .TRUE.
408	extremum(ij, llm) = .TRUE.
409	enddo
410
411	! calcul des valeurs max et min acceptees aux interfaces
412	do l = 2, llm
413	do ij = 1, ip1jmp1
414	zqmax(ij, l) = max(q(ij, l - 1), q(ij, l))
415	zqmin(ij, l) = min(q(ij, l - 1), q(ij, l))
416	enddo
417	enddo
418
419	do l = 2, llm
420	do ij = 1, ip1jmp1
421	zqw(ij, l) = min(max(zqmin(ij, l), zqw(ij, l)), zqmax(ij, l))
422	enddo
423	enddo
424
425	do l = 2, llm - 1
426	do ij = 1, ip1jmp1
427	IF(extremum(ij, l)) THEN
428	qh(ij, l) = q(ij, l)
429	qb(ij, l) = q(ij, l)
430	else
431	qh(ij, l) = zqw(ij, l + 1)
432	qb(ij, l) = zqw(ij, l)
433	endif
434	enddo
435	enddo
436	! do l=2,llm-1
437	! do ij=1,ip1jmp1
438	! IF(extremum(ij,l)) THEN
439	! if (.NOT.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
440	! if (.NOT.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
441	! endif
442	! enddo
443	! enddo
444
445	do ij = 1, ip1jmp1
446	qb(ij, 1) = q(ij, 1)
447	qh(ij, 1) = q(ij, 1)
448	qb(ij, llm) = q(ij, llm)
449	qh(ij, llm) = q(ij, llm)
450	enddo
451
452	ENDIF
453
454	RETURN
455	END SUBROUTINE advnqz
456
457	SUBROUTINE advnx(q, qg, qd, masse, u_m, mode)
458	!
459	! Auteur : F. Hourdin
460	!
461	! ********************************************************************
462	! Shema d'advection " pseudo amont " .
463	! ********************************************************************
464	! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
465	!
466	!
467	! --------------------------------------------------------------------
468	USE lmdz_iniprint, ONLY: lunout, prt_level
469	IMPLICIT NONE
470	!
471	INCLUDE "dimensions.h"
472	INCLUDE "paramet.h"
473	!
474	!
475	! Arguments:
476	! ----------
477	INTEGER :: mode
478	REAL :: masse(ip1jmp1, llm)
479	REAL :: u_m(ip1jmp1, llm)
480	REAL :: q(ip1jmp1, llm), qd(ip1jmp1, llm), qg(ip1jmp1, llm)
481	!
482	! Local
483	! ---------
484	!
485	INTEGER :: i, j, ij, l, indu(ip1jmp1), niju, iju, ijq
486	INTEGER :: n0, nl(llm)
487	!
488	REAL :: new_m, zu_m, zdq, zz
489	REAL :: zsigg(ip1jmp1, llm), zsigd(ip1jmp1, llm), zsig
490	REAL :: u_mq(ip1jmp1, llm)
491
492	REAL :: zm, zq, zsigm, zsigp, zqm, zqp, zu
493
494	LOGICAL :: ladvplus(ip1jmp1, llm)
495
496	REAL :: prec
497	save prec
498
499	data prec/1.e-15/
500
501	do l = 1, llm
502	do ij = iip2, ip1jm
503	zdq = qd(ij, l) - qg(ij, l)
504	! if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) THEN
505	! PRINT*,'probleme au point ij=',ij,' l=',l
506	! PRINT*,qd(ij,l),q(ij,l),qg(ij,l)
507	! qd(ij,l)=q(ij,l)
508	! qg(ij,l)=q(ij,l)
509	! END IF
510	IF(abs(zdq)>prec) THEN
511	zsigd(ij, l) = (q(ij, l) - qg(ij, l)) / zdq
512	zsigg(ij, l) = 1. - zsigd(ij, l)
513	! IF(.NOT.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .AND.
514	! s zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) THEN
515	! PRINT*,'probleme au point ij=',ij,' l=',l
516	! PRINT*,'sigg=',zsigg(ij,l),' sigd=',zsigd(ij,l)
517	! PRINT*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
518	! stop
519	! END IF
520	else
521	zsigd(ij, l) = 0.5
522	zsigg(ij, l) = 0.5
523	qd(ij, l) = q(ij, l)
524	qg(ij, l) = q(ij, l)
525	endif
526	enddo
527	enddo
528
529	! calcul de la pente maximum dans la maille en valeur absolue
530
531	do l = 1, llm
532	do ij = iip2, ip1jm - 1
533	IF (u_m(ij, l)>=0.) THEN
534	zsigp = zsigd(ij, l)
535	zsigm = zsigg(ij, l)
536	zqp = qd(ij, l)
537	zqm = qg(ij, l)
538	zm = masse(ij, l)
539	zq = q(ij, l)
540	else
541	zsigm = zsigd(ij + 1, l)
542	zsigp = zsigg(ij + 1, l)
543	zqm = qd(ij + 1, l)
544	zqp = qg(ij + 1, l)
545	zm = masse(ij + 1, l)
546	zq = q(ij + 1, l)
547	endif
548	zu = abs(u_m(ij, l))
549	ladvplus(ij, l) = zu>zm
550	zsig = zu / zm
551	IF(zsig==0.) zsigp = 0.1
552	IF (mode==1) THEN
553	IF (zsig<=zsigp) THEN
554	u_mq(ij, l) = u_m(ij, l) * zqp
555	ELSE IF (mode==1) THEN
556	u_mq(ij, l) = &
557	sign(zm, u_m(ij, l)) * (zsigp * zqp + (zsig - zsigp) * zqm)
558	endif
559	else
560	IF (zsig<=zsigp) THEN
561	u_mq(ij, l) = u_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq))
562	else
563	zz = 0.5 * (zsig - zsigp) / zsigm
564	u_mq(ij, l) = sign(zm, u_m(ij, l)) * (0.5 * (zq + zqp) * zsigp &
565	+ (zsig - zsigp) * (zq + zz * (zqm - zq)))
566	endif
567	endif
568	! IF(zsig.lt.0.) THEN
569	! PRINT*,'au point ij=',ij,' l=',l,' sig=',zsig
570	! stop
571	! END IF
572	enddo
573	enddo
574
575	do l = 1, llm
576	do ij = iip1 + iip1, ip1jm, iip1
577	u_mq(ij, l) = u_mq(ij - iim, l)
578	ladvplus(ij, l) = ladvplus(ij - iim, l)
579	enddo
580	enddo
581
582	!=================================================================
583	! SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
584	!=================================================================
585	! tris des regions a traiter
586	n0 = 0
587	do l = 1, llm
588	nl(l) = 0
589	do ij = iip2, ip1jm
590	IF(ladvplus(ij, l)) THEN
591	nl(l) = nl(l) + 1
592	u_mq(ij, l) = 0.
593	endif
594	enddo
595	n0 = n0 + nl(l)
596	enddo
597
598	IF(n0>1) THEN
599	IF (prt_level > 9) WRITE(lunout, *) &
600	'Nombre de points pour lesquels on advect plus que le' &
601	, 'contenu de la maille : ', n0
602
603	do l = 1, llm
604	IF(nl(l)>0) THEN
605	iju = 0
606	! indicage des mailles concernees par le traitement special
607	do ij = iip2, ip1jm
608	IF(ladvplus(ij, l).AND.mod(ij, iip1)/=0) THEN
609	iju = iju + 1
610	indu(iju) = ij
611	endif
612	enddo
613	niju = iju
614	! PRINT*,'niju,nl',niju,nl(l)
615
616	! traitement des mailles
617	do iju = 1, niju
618	ij = indu(iju)
619	j = (ij - 1) / iip1 + 1
620	zu_m = u_m(ij, l)
621	u_mq(ij, l) = 0.
622	IF(zu_m>0.) THEN
623	ijq = ij
624	i = ijq - (j - 1) * iip1
625	! accumulation pour les mailles completements advectees
626	do while(zu_m>masse(ijq, l))
627	u_mq(ij, l) = u_mq(ij, l) + q(ijq, l) * masse(ijq, l)
628	zu_m = zu_m - masse(ijq, l)
629	i = mod(i - 2 + iim, iim) + 1
630	ijq = (j - 1) * iip1 + i
631	enddo
632	! MODIFS SPECIFIQUES DU SCHEMA
633	! ajout de la maille non completement advectee
634	zsig = zu_m / masse(ijq, l)
635	IF(zsig<=zsigd(ijq, l)) THEN
636	u_mq(ij, l) = u_mq(ij, l) + zu_m * (qd(ijq, l) &
637	- 0.5 * zsig / zsigd(ijq, l) * (qd(ijq, l) - q(ijq, l)))
638	else
639	! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
640	! goto 8888
641	zz = 0.5 * (zsig - zsigd(ijq, l)) / zsigg(ijq, l)
642	IF(.NOT.(zz>0..and.zz<=0.5)) THEN
643	WRITE(lunout, *)'probleme2 au point ij=', ij, &
644	' l=', l
645	WRITE(lunout, *)'zz=', zz
646	stop
647	endif
648	u_mq(ij, l) = u_mq(ij, l) + masse(ijq, l) * (&
649	0.5 * (q(ijq, l) + qd(ijq, l)) * zsigd(ijq, l) &
650	+ (zsig - zsigd(ijq, l)) * (q(ijq, l) + zz * (qg(ijq, l) - q(ijq, l))))
651	endif
652	else
653	ijq = ij + 1
654	i = ijq - (j - 1) * iip1
655	! accumulation pour les mailles completements advectees
656	do while(-zu_m>masse(ijq, l))
657	u_mq(ij, l) = u_mq(ij, l) - q(ijq, l) * masse(ijq, l)
658	zu_m = zu_m + masse(ijq, l)
659	i = mod(i, iim) + 1
660	ijq = (j - 1) * iip1 + i
661	enddo
662	! ajout de la maille non completement advectee
663	! 2eme MODIF SPECIFIQUE
664	zsig = -zu_m / masse(ij + 1, l)
665	IF(zsig<=zsigg(ijq, l)) THEN
666	u_mq(ij, l) = u_mq(ij, l) + zu_m * (qg(ijq, l) &
667	- 0.5 * zsig / zsigg(ijq, l) * (qg(ijq, l) - q(ijq, l)))
668	else
669	! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
670	! goto 9999
671	zz = 0.5 * (zsig - zsigg(ijq, l)) / zsigd(ijq, l)
672	IF(.NOT.(zz>0..and.zz<=0.5)) THEN
673	WRITE(lunout, *)'probleme22 au point ij=', ij &
674	, ' l=', l
675	WRITE(lunout, *)'zz=', zz
676	stop
677	endif
678	u_mq(ij, l) = u_mq(ij, l) - masse(ijq, l) * (&
679	0.5 * (q(ijq, l) + qg(ijq, l)) * zsigg(ijq, l) &
680	+ (zsig - zsigg(ijq, l)) * &
681	(q(ijq, l) + zz * (qd(ijq, l) - q(ijq, l))))
682	endif
683	! fin de la modif
684	endif
685	enddo
686	endif
687	enddo
688	ENDIF ! n0.gt.0
689
690	! bouclage en latitude
691	do l = 1, llm
692	do ij = iip1 + iip1, ip1jm, iip1
693	u_mq(ij, l) = u_mq(ij - iim, l)
694	enddo
695	enddo
696
697	!=================================================================
698	! CALCUL DE LA CONVERGENCE DES FLUX
699	!=================================================================
700
701	do l = 1, llm
702	do ij = iip2 + 1, ip1jm
703	new_m = masse(ij, l) + u_m(ij - 1, l) - u_m(ij, l)
704	q(ij, l) = (q(ij, l) * masse(ij, l) + &
705	u_mq(ij - 1, l) - u_mq(ij, l)) &
706	/ new_m
707	masse(ij, l) = new_m
708	enddo
709	! Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
710	do ij = iip1 + iip1, ip1jm, iip1
711	q(ij - iim, l) = q(ij, l)
712	masse(ij - iim, l) = masse(ij, l)
713	enddo
714	enddo
715
716	RETURN
717	END SUBROUTINE advnx
718	SUBROUTINE advny(q, qs, qn, masse, v_m)
719	!
720	! Auteur : F. Hourdin
721	!
722	! ********************************************************************
723	! Shema d'advection " pseudo amont " .
724	! ********************************************************************
725	! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
726	!
727	!
728	! --------------------------------------------------------------------
729	USE lmdz_iniprint, ONLY: lunout, prt_level
730	USE lmdz_ssum_scopy, ONLY: ssum
731
732	IMPLICIT NONE
733	!
734	INCLUDE "dimensions.h"
735	INCLUDE "paramet.h"
736	INCLUDE "comgeom.h"
737	!
738	!
739	! Arguments:
740	! ----------
741	REAL :: masse(ip1jmp1, llm)
742	REAL :: v_m(ip1jm, llm)
743	REAL :: q(ip1jmp1, llm), qn(ip1jmp1, llm), qs(ip1jmp1, llm)
744	!
745	! Local
746	! ---------
747	!
748	INTEGER :: ij, l
749	!
750	REAL :: new_m, zdq, zz
751	REAL :: zsigs(ip1jmp1), zsign(ip1jmp1), zsig
752	REAL :: v_mq(ip1jm, llm)
753	REAL :: convpn, convps, convmpn, convmps, massen, masses
754	REAL :: zm, zq, zsigm, zsigp, zqm, zqp
755	REAL :: prec
756	save prec
757
758	data prec/1.e-15/
759	do l = 1, llm
760	do ij = 1, ip1jmp1
761	zdq = qn(ij, l) - qs(ij, l)
762	! if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) THEN
763	! PRINT*,'probleme au point ij=',ij,' l=',l,' advnqx'
764	! PRINT*,qn(ij,l),q(ij,l),qs(ij,l)
765	! qn(ij,l)=q(ij,l)
766	! qs(ij,l)=q(ij,l)
767	! END IF
768	IF(abs(zdq)>prec) THEN
769	zsign(ij) = (q(ij, l) - qs(ij, l)) / zdq
770	zsigs(ij) = 1. - zsign(ij)
771	! IF(.NOT.(zsign(ij).ge.0..and.zsign(ij).le.1. .AND.
772	! s zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) THEN
773	! PRINT*,'probleme au point ij=',ij,' l=',l
774	! PRINT*,'sigs=',zsigs(ij),' sign=',zsign(ij)
775	! stop
776	! END IF
777	else
778	zsign(ij) = 0.5
779	zsigs(ij) = 0.5
780	endif
781	enddo
782
783	! calcul de la pente maximum dans la maille en valeur absolue
784
785	do ij = 1, ip1jm
786	IF (v_m(ij, l)>=0.) THEN
787	zsigp = zsign(ij + iip1)
788	zsigm = zsigs(ij + iip1)
789	zqp = qn(ij + iip1, l)
790	zqm = qs(ij + iip1, l)
791	zm = masse(ij + iip1, l)
792	zq = q(ij + iip1, l)
793	else
794	zsigm = zsign(ij)
795	zsigp = zsigs(ij)
796	zqm = qn(ij, l)
797	zqp = qs(ij, l)
798	zm = masse(ij, l)
799	zq = q(ij, l)
800	endif
801	zsig = abs(v_m(ij, l)) / zm
802	IF(zsig==0.) zsigp = 0.1
803	IF (zsig<=zsigp) THEN
804	v_mq(ij, l) = v_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq))
805	else
806	zz = 0.5 * (zsig - zsigp) / zsigm
807	v_mq(ij, l) = sign(zm, v_m(ij, l)) * (0.5 * (zq + zqp) * zsigp &
808	+ (zsig - zsigp) * (zq + zz * (zqm - zq)))
809	endif
810	enddo
811	enddo
812
813	do l = 1, llm
814	do ij = iip2, ip1jm
815	new_m = masse(ij, l) &
816	+ v_m(ij, l) - v_m(ij - iip1, l)
817	q(ij, l) = (q(ij, l) * masse(ij, l) + v_mq(ij, l) - v_mq(ij - iip1, l)) &
818	/ new_m
819	masse(ij, l) = new_m
820	enddo
821	!.-. ancienne version
822	convpn = SSUM(iim, v_mq(1, l), 1)
823	convmpn = ssum(iim, v_m(1, l), 1)
824	massen = ssum(iim, masse(1, l), 1)
825	new_m = massen + convmpn
826	q(1, l) = (q(1, l) * massen + convpn) / new_m
827	do ij = 1, iip1
828	q(ij, l) = q(1, l)
829	masse(ij, l) = new_m * aire(ij) / apoln
830	enddo
831
832	convps = -SSUM(iim, v_mq(ip1jm - iim, l), 1)
833	convmps = -ssum(iim, v_m(ip1jm - iim, l), 1)
834	masses = ssum(iim, masse(ip1jm + 1, l), 1)
835	new_m = masses + convmps
836	q(ip1jm + 1, l) = (q(ip1jm + 1, l) * masses + convps) / new_m
837	do ij = ip1jm + 1, ip1jmp1
838	q(ij, l) = q(ip1jm + 1, l)
839	masse(ij, l) = new_m * aire(ij) / apols
840	enddo
841	enddo
842
843	RETURN
844	END SUBROUTINE advny
845	SUBROUTINE advnz(q, qh, qb, masse, w_m)
846	!
847	! Auteurs: F.Hourdin
848	!
849	! ********************************************************************
850	! Shema d'advection " pseudo amont " .
851	! b designe le bas et h le haut
852	! il y a une correspondance entre le b en z et le d en x
853	! ********************************************************************
854	!
855	!
856	! --------------------------------------------------------------------
857	USE lmdz_iniprint, ONLY: lunout, prt_level
858	IMPLICIT NONE
859	!
860	INCLUDE "dimensions.h"
861	INCLUDE "paramet.h"
862	INCLUDE "comgeom.h"
863	!
864	!
865	! Arguments:
866	! ----------
867	REAL :: masse(ip1jmp1, llm)
868	REAL :: w_m(ip1jmp1, llm + 1)
869	REAL :: q(ip1jmp1, llm), qb(ip1jmp1, llm), qh(ip1jmp1, llm)
870
871	!
872	! Local
873	! ---------
874	!
875	INTEGER :: ij, l
876	!
877	REAL :: new_m, zdq, zz
878	REAL :: zsigh(ip1jmp1, llm), zsigb(ip1jmp1, llm), zsig
879	REAL :: w_mq(ip1jmp1, llm + 1)
880	REAL :: zm, zq, zsigm, zsigp, zqm, zqp
881	REAL :: prec
882	save prec
883
884	data prec/1.e-13/
885
886	do l = 1, llm
887	do ij = 1, ip1jmp1
888	zdq = qb(ij, l) - qh(ij, l)
889	! if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) THEN
890	! PRINT*,'probleme au point ij=',ij,' l=',l
891	! PRINT*,qh(ij,l),q(ij,l),qb(ij,l)
892	! qh(ij,l)=q(ij,l)
893	! qb(ij,l)=q(ij,l)
894	! END IF
895
896	IF(abs(zdq)>prec) THEN
897	zsigb(ij, l) = (q(ij, l) - qh(ij, l)) / zdq
898	zsigh(ij, l) = 1. - zsigb(ij, l)
899	zsigb(ij, l) = min(max(zsigb(ij, l), 0.), 1.)
900	else
901	zsigb(ij, l) = 0.5
902	zsigh(ij, l) = 0.5
903	endif
904	enddo
905	enddo
906
907	! PRINT*,'ok1'
908	! calcul de la pente maximum dans la maille en valeur absolue
909	do l = 2, llm
910	do ij = 1, ip1jmp1
911	IF (w_m(ij, l)>=0.) THEN
912	zsigp = zsigb(ij, l)
913	zsigm = zsigh(ij, l)
914	zqp = qb(ij, l)
915	zqm = qh(ij, l)
916	zm = masse(ij, l)
917	zq = q(ij, l)
918	else
919	zsigm = zsigb(ij, l - 1)
920	zsigp = zsigh(ij, l - 1)
921	zqm = qb(ij, l - 1)
922	zqp = qh(ij, l - 1)
923	zm = masse(ij, l - 1)
924	zq = q(ij, l - 1)
925	endif
926	zsig = abs(w_m(ij, l)) / zm
927	IF(zsig==0.) zsigp = 0.1
928	IF (zsig<=zsigp) THEN
929	w_mq(ij, l) = w_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq))
930	else
931	zz = 0.5 * (zsig - zsigp) / zsigm
932	w_mq(ij, l) = sign(zm, w_m(ij, l)) * (0.5 * (zq + zqp) * zsigp &
933	+ (zsig - zsigp) * (zq + zz * (zqm - zq)))
934	endif
935	enddo
936	enddo
937
938	do ij = 1, ip1jmp1
939	w_mq(ij, llm + 1) = 0.
940	w_mq(ij, 1) = 0.
941	enddo
942
943	do l = 1, llm
944	do ij = 1, ip1jmp1
945	new_m = masse(ij, l) + w_m(ij, l + 1) - w_m(ij, l)
946	q(ij, l) = (q(ij, l) * masse(ij, l) + w_mq(ij, l + 1) - w_mq(ij, l)) &
947	/ new_m
948	masse(ij, l) = new_m
949	enddo
950	enddo
951	! PRINT*,'ok3'
952	RETURN
953	END SUBROUTINE advnz

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

Download in other formats:

Original Format