Context Navigation

advx.F @ 5082

Last change on this file since 5082 was 5077, checked in by abarral, 5 months ago
Fix incorrect alignment of variables in COMMON block comdissip.h Replace obsolete DO with shared termination in advx.F, advxp.F (minor) replace obsolete bool operators
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: 11.9 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE advx(limit,dtx,pbaru,sm,s0,
5	$ sx,sy,sz,lati,latf)
6	IMPLICIT NONE
7
8	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
9	C C
10	C first-order moments (FOM) advection of tracer in X direction C
11	C C
12	C Source : Pascal Simon (Meteo,CNRM) C
13	C Adaptation : A.Armengaud (LGGE) juin 94 C
14	C C
15	C limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz C
16	C sont des arguments d'entree pour le s-pg... C
17	C C
18	C sm,s0,sx,sy,sz C
19	C sont les arguments de sortie pour le s-pg C
20	C C
21	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
22	C
23	C parametres principaux du modele
24	C
25	include "dimensions.h"
26	include "paramet.h"
27
28	C Arguments :
29	C -----------
30	C dtx : frequence fictive d'appel du transport
31	C pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1
32
33	INTEGER ntra
34	PARAMETER (ntra = 1)
35
36	C ATTENTION partout ou on trouve ntra, insertion de boucle
37	C possible dans l'avenir.
38
39	REAL dtx
40	REAL pbaru ( iip1,jjp1,llm )
41
42	C moments: SM total mass in each grid box
43	C S0 mass of tracer in each grid box
44	C Si 1rst order moment in i direction
45	C
46	REAL SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra)
47	REAL sx(iip1,jjp1,llm,ntra)
48	$ ,sy(iip1,jjp1,llm,ntra)
49	REAL sz(iip1,jjp1,llm,ntra)
50
51	C Local :
52	C -------
53
54	C mass fluxes across the boundaries (UGRI,VGRI,WGRI)
55	C mass fluxes in kg
56	C declaration :
57
58	REAL UGRI(iip1,jjp1,llm)
59
60	C Rem : VGRI et WGRI ne sont pas utilises dans
61	C cette subroutine ( advection en x uniquement )
62	C
63	C Ti are the moments for the current latitude and level
64	C
65	REAL TM(iim)
66	REAL T0(iim,ntra),TX(iim,ntra)
67	REAL TY(iim,ntra),TZ(iim,ntra)
68	REAL TEMPTM ! just a temporary variable
69	C
70	C the moments F are similarly defined and used as temporary
71	C storage for portions of the grid boxes in transit
72	C
73	REAL FM(iim)
74	REAL F0(iim,ntra),FX(iim,ntra)
75	REAL FY(iim,ntra),FZ(iim,ntra)
76	C
77	C work arrays
78	C
79	REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
80	C
81	REAL SMNEW(iim),UEXT(iim)
82	C
83	REAL sqi,sqf
84
85	LOGICAL LIMIT
86	INTEGER NUM(jjp1),LONK,NUMK
87	INTEGER lon,lati,latf,niv
88	INTEGER i,i2,i3,j,jv,l,k,itrac
89
90	lon = iim
91	niv = llm
92
93	C * Test de passage d'arguments ****
94
95
96	C -------------------------------------
97	DO j = 1,jjp1
98	NUM(j) = 1
99	END DO
100	sqi = 0.
101	sqf = 0.
102
103	DO l = 1,llm
104	DO j = 1,jjp1
105	DO i = 1,iim
106	cIM 240305 sqi = sqi + S0(i,j,l,9)
107	sqi = sqi + S0(i,j,l,ntra)
108	ENDDO
109	ENDDO
110	ENDDO
111	PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
112	PRINT*,'sqi=',sqi
113
114
115	C Interface : adaptation nouveau modele
116	C -------------------------------------
117	C
118	C ---------------------------------------------------------
119	C Conversion des flux de masses en kg/s
120	C pbaru est en N/s d'ou :
121	C ugri est en kg/s
122
123	DO l = 1,llm
124	DO j = 1,jjm+1
125	DO i = 1,iip1
126	C ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
127	ugri (i,j,llm+1-l) = pbaru (i,j,l)
128	END DO
129	END DO
130	END DO
131
132
133	C ---------------------------------------------------------
134	C ---------------------------------------------------------
135	C ---------------------------------------------------------
136
137	C start here
138	C
139	C boucle principale sur les niveaux et les latitudes
140	C
141	DO L=1,NIV
142	DO K=lati,latf
143	C
144	C initialisation
145	C
146	C program assumes periodic boundaries in X
147	C
148	DO I=2,LON
149	SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
150	END DO
151	SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
152	C
153	C modifications for extended polar zones
154	C
155	NUMK=NUM(K)
156	LONK=LON/NUMK
157	C
158	IF(NUMK>1) THEN
159	C
160	DO I=1,LON
161	TM(I)=0.
162	END DO
163	DO JV=1,NTRA
164	DO I=1,LON
165	T0(I,JV)=0.
166	TX(I,JV)=0.
167	TY(I,JV)=0.
168	TZ(I,JV)=0.
169	END DO
170	END DO
171	C
172	DO I2=1,NUMK
173	C
174	DO I=1,LONK
175	I3=(I-1)*NUMK+I2
176	TM(I)=TM(I)+SM(I3,K,L)
177	ALF(I)=SM(I3,K,L)/TM(I)
178	ALF1(I)=1.-ALF(I)
179	END DO
180	C
181	DO JV=1,NTRA
182	DO I=1,LONK
183	I3=(I-1)*NUMK+I2
184	TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)
185	$ *S0(I3,K,L,JV)
186	T0(I,JV)=T0(I,JV)+S0(I3,K,L,JV)
187	TX(I,JV)=ALF(I) *sx(I3,K,L,JV)+
188	$ ALF1(I)TX(I,JV) +3.TEMPTM
189	TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV)
190	TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV)
191	ENDDO
192	ENDDO
193	C
194	END DO
195	C
196	ELSE
197	C
198	DO I=1,LON
199	TM(I)=SM(I,K,L)
200	END DO
201	DO JV=1,NTRA
202	DO I=1,LON
203	T0(I,JV)=S0(I,K,L,JV)
204	TX(I,JV)=sx(I,K,L,JV)
205	TY(I,JV)=sy(I,K,L,JV)
206	TZ(I,JV)=sz(I,K,L,JV)
207	END DO
208	END DO
209	C
210	ENDIF
211	C
212	DO I=1,LONK
213	UEXT(I)=UGRI(I*NUMK,K,L)
214	END DO
215	C
216	C place limits on appropriate moments before transport
217	C (if flux-limiting is to be applied)
218	C
219	IF(.NOT.LIMIT) GO TO 13
220	C
221	DO JV=1,NTRA
222	DO I=1,LONK
223	TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
224	END DO
225	END DO
226	C
227	13 CONTINUE
228	C
229	C calculate flux and moments between adjacent boxes
230	C 1- create temporary moments/masses for partial boxes in transit
231	C 2- reajusts moments remaining in the box
232	C
233	C flux from IP to I if U(I).lt.0
234	C
235	DO I=1,LONK-1
236	IF(UEXT(I)<0.) THEN
237	FM(I)=-UEXT(I)*DTX
238	ALF(I)=FM(I)/TM(I+1)
239	TM(I+1)=TM(I+1)-FM(I)
240	ENDIF
241	END DO
242	C
243	I=LONK
244	IF(UEXT(I)<0.) THEN
245	FM(I)=-UEXT(I)*DTX
246	ALF(I)=FM(I)/TM(1)
247	TM(1)=TM(1)-FM(I)
248	ENDIF
249	C
250	C flux from I to IP if U(I).gt.0
251	C
252	DO I=1,LONK
253	IF(UEXT(I)>=0.) THEN
254	FM(I)=UEXT(I)*DTX
255	ALF(I)=FM(I)/TM(I)
256	TM(I)=TM(I)-FM(I)
257	ENDIF
258	END DO
259	C
260	DO I=1,LONK
261	ALFQ(I)=ALF(I)*ALF(I)
262	ALF1(I)=1.-ALF(I)
263	ALF1Q(I)=ALF1(I)*ALF1(I)
264	END DO
265	C
266	DO JV=1,NTRA
267	DO I=1,LONK-1
268	C
269	IF(UEXT(I)<0.) THEN
270	C
271	F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
272	FX(I,JV)=ALFQ(I)*TX(I+1,JV)
273	FY(I,JV)=ALF (I)*TY(I+1,JV)
274	FZ(I,JV)=ALF (I)*TZ(I+1,JV)
275	C
276	T0(I+1,JV)=T0(I+1,JV)-F0(I,JV)
277	TX(I+1,JV)=ALF1Q(I)*TX(I+1,JV)
278	TY(I+1,JV)=TY(I+1,JV)-FY(I,JV)
279	TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,JV)
280	C
281	ENDIF
282	C
283	END DO
284	END DO
285	C
286	I=LONK
287	IF(UEXT(I)<0.) THEN
288	C
289	DO JV=1,NTRA
290	C
291	F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
292	FX (I,JV)=ALFQ(I)*TX(1,JV)
293	FY (I,JV)=ALF (I)*TY(1,JV)
294	FZ (I,JV)=ALF (I)*TZ(1,JV)
295	C
296	T0(1,JV)=T0(1,JV)-F0(I,JV)
297	TX(1,JV)=ALF1Q(I)*TX(1,JV)
298	TY(1,JV)=TY(1,JV)-FY(I,JV)
299	TZ(1,JV)=TZ(1,JV)-FZ(I,JV)
300	C
301	END DO
302	C
303	ENDIF
304	C
305	DO JV=1,NTRA
306	DO I=1,LONK
307	C
308	IF(UEXT(I)>=0.) THEN
309	C
310	F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
311	FX(I,JV)=ALFQ(I)*TX(I,JV)
312	FY(I,JV)=ALF (I)*TY(I,JV)
313	FZ(I,JV)=ALF (I)*TZ(I,JV)
314	C
315	T0(I,JV)=T0(I,JV)-F0(I,JV)
316	TX(I,JV)=ALF1Q(I)*TX(I,JV)
317	TY(I,JV)=TY(I,JV)-FY(I,JV)
318	TZ(I,JV)=TZ(I,JV)-FZ(I,JV)
319	C
320	ENDIF
321	C
322	END DO
323	END DO
324	C
325	C puts the temporary moments Fi into appropriate neighboring boxes
326	C
327	DO I=1,LONK
328	IF(UEXT(I)<0.) THEN
329	TM(I)=TM(I)+FM(I)
330	ALF(I)=FM(I)/TM(I)
331	ENDIF
332	END DO
333	C
334	DO I=1,LONK-1
335	IF(UEXT(I)>=0.) THEN
336	TM(I+1)=TM(I+1)+FM(I)
337	ALF(I)=FM(I)/TM(I+1)
338	ENDIF
339	END DO
340	C
341	I=LONK
342	IF(UEXT(I)>=0.) THEN
343	TM(1)=TM(1)+FM(I)
344	ALF(I)=FM(I)/TM(1)
345	ENDIF
346	C
347	DO I=1,LONK
348	ALF1(I)=1.-ALF(I)
349	END DO
350	C
351	DO JV=1,NTRA
352	DO I=1,LONK
353	C
354	IF(UEXT(I)<0.) THEN
355	C
356	TEMPTM=-ALF(I)T0(I,JV)+ALF1(I)F0(I,JV)
357	T0(I,JV)=T0(I,JV)+F0(I,JV)
358	TX(I,JV)=ALF(I)FX(I,JV)+ALF1(I)TX(I,JV)+3.*TEMPTM
359	TY(I,JV)=TY(I,JV)+FY(I,JV)
360	TZ(I,JV)=TZ(I,JV)+FZ(I,JV)
361	C
362	ENDIF
363	C
364	END DO
365	END DO
366	C
367	DO JV=1,NTRA
368	DO I=1,LONK-1
369	C
370	IF(UEXT(I)>=0.) THEN
371	C
372	TEMPTM=ALF(I)T0(I+1,JV)-ALF1(I)F0(I,JV)
373	T0(I+1,JV)=T0(I+1,JV)+F0(I,JV)
374	TX(I+1,JV)=ALF(I)FX(I,JV)+ALF1(I)TX(I+1,JV)+3.*TEMPTM
375	TY(I+1,JV)=TY(I+1,JV)+FY(I,JV)
376	TZ(I+1,JV)=TZ(I+1,JV)+FZ(I,JV)
377	C
378	ENDIF
379	C
380	END DO
381	END DO
382	C
383	I=LONK
384	IF(UEXT(I)>=0.) THEN
385	DO JV=1,NTRA
386	TEMPTM=ALF(I)T0(1,JV)-ALF1(I)F0(I,JV)
387	T0(1,JV)=T0(1,JV)+F0(I,JV)
388	TX(1,JV)=ALF(I)FX(I,JV)+ALF1(I)TX(1,JV)+3.*TEMPTM
389	TY(1,JV)=TY(1,JV)+FY(I,JV)
390	TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
391	END DO
392	ENDIF
393	C
394	C retour aux mailles d'origine (passage des Tij aux Sij)
395	C
396	IF(NUMK>1) THEN
397	C
398	DO I2=1,NUMK
399	C
400	DO I=1,LONK
401	C
402	I3=I2+(I-1)*NUMK
403	SM(I3,K,L)=SMNEW(I3)
404	ALF(I)=SMNEW(I3)/TM(I)
405	TM(I)=TM(I)-SMNEW(I3)
406	C
407	ALFQ(I)=ALF(I)*ALF(I)
408	ALF1(I)=1.-ALF(I)
409	ALF1Q(I)=ALF1(I)*ALF1(I)
410	C
411	END DO
412	END DO
413	C
414	DO JV=1,NTRA
415	DO I=1,LONK
416	C
417	I3=I2+(I-1)*NUMK
418	S0(I3,K,L,JV)=ALF (I)
419	$ * (T0(I,JV)-ALF1(I)*TX(I,JV))
420	sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV)
421	sy(I3,K,L,JV)=ALF (I)*TY(I,JV)
422	sz(I3,K,L,JV)=ALF (I)*TZ(I,JV)
423	C
424	C reajusts moments remaining in the box
425	C
426	T0(I,JV)=T0(I,JV)-S0(I3,K,L,JV)
427	TX(I,JV)=ALF1Q(I)*TX(I,JV)
428	TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV)
429	TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV)
430	ENDDO
431	ENDDO
432	C
433	C
434	ELSE
435	C
436	DO I=1,LON
437	SM(I,K,L)=TM(I)
438	END DO
439	DO JV=1,NTRA
440	DO I=1,LON
441	S0(I,K,L,JV)=T0(I,JV)
442	sx(I,K,L,JV)=TX(I,JV)
443	sy(I,K,L,JV)=TY(I,JV)
444	sz(I,K,L,JV)=TZ(I,JV)
445	END DO
446	END DO
447	C
448	ENDIF
449	C
450	END DO
451	END DO
452	C
453	C ----------- AA Test en fin de ADVX ------ Controle des S*
454	c OK
455	c DO 9998 l = 1, llm
456	c DO 9998 j = 1, jjp1
457	c DO 9998 i = 1, iip1
458	c IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
459	c PRINT*, '-------------------'
460	c PRINT*, 'En fin de ADVX'
461	c PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
462	c PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
463	c print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
464	c print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
465	c print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
466	c WRITE (,) 'On arrete !! - pbl en fin de ADVX1'
467	cc STOP
468	c ENDIF
469	c 9998 CONTINUE
470	c
471	C ---------- bouclage cyclique
472	DO itrac=1,ntra
473	DO l = 1,llm
474	DO j = lati,latf
475	SM(iip1,j,l) = SM(1,j,l)
476	S0(iip1,j,l,itrac) = S0(1,j,l,itrac)
477	sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
478	sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
479	sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
480	END DO
481	END DO
482	ENDDO
483
484	c ----------- qqtite totale de traceur dans tte l'atmosphere
485	DO l = 1, llm
486	DO j = 1, jjp1
487	DO i = 1, iim
488	cIM 240405 sqf = sqf + S0(i,j,l,9)
489	sqf = sqf + S0(i,j,l,ntra)
490	END DO
491	END DO
492	END DO
493	c
494	PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
495	PRINT*,'sqf=',sqf
496	c-------------
497
498	RETURN
499	END
500	C_________________________________________________________________
501	C_________________________________________________________________

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

Download in other formats:

Original Format