Context Navigation

advx.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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

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

Download in other formats:

Original Format