Context Navigation

calfis.F @ 1198

Last change on this file since 1198 was 1114, checked in by jghattas, 16 years ago

Creation du module infotrac:

contient les variables de advtrac.h
contient la subroutine iniadvtrac renommer en infotrac_init
le nombre des traceurs est lu dans tracer.def en dynamique (ou par default ou recu par INCA)
ce module est utilise dans la dynamique et la physique
contient aussi la variable nbtr qui avant etait stockee dans dimphy

Le fichier advtrac.h n'existe plus.
La compilation ne prend plus en compte le nombre de traceur.

/JG

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 16.6 KB

Line
1	!
2	! $Header$
3	!
4	C
5	C
6	SUBROUTINE calfis(lafin,
7	$ rdayvrai,
8	$ heure,
9	$ pucov,
10	$ pvcov,
11	$ pteta,
12	$ pq,
13	$ pmasse,
14	$ pps,
15	$ pp,
16	$ ppk,
17	$ pphis,
18	$ pphi,
19	$ pducov,
20	$ pdvcov,
21	$ pdteta,
22	$ pdq,
23	$ flxw,
24	$ clesphy0,
25	$ pdufi,
26	$ pdvfi,
27	$ pdhfi,
28	$ pdqfi,
29	$ pdpsfi)
30	c
31	c Auteur : P. Le Van, F. Hourdin
32	c .........
33	USE infotrac
34
35	IMPLICIT NONE
36	c=======================================================================
37	c
38	c 1. rearrangement des tableaux et transformation
39	c variables dynamiques > variables physiques
40	c 2. calcul des termes physiques
41	c 3. retransformation des tendances physiques en tendances dynamiques
42	c
43	c remarques:
44	c ----------
45	c
46	c - les vents sont donnes dans la physique par leurs composantes
47	c naturelles.
48	c - la variable thermodynamique de la physique est une variable
49	c intensive : T
50	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
51	c - les deux seules variables dependant de la geometrie necessaires
52	c pour la physique sont la latitude pour le rayonnement et
53	c l'aire de la maille quand on veut integrer une grandeur
54	c horizontalement.
55	c - les points de la physique sont les points scalaires de la
56	c la dynamique; numerotation:
57	c 1 pour le pole nord
58	c (jjm-1)*iim pour l'interieur du domaine
59	c ngridmx pour le pole sud
60	c ---> ngridmx=2+(jjm-1)*iim
61	c
62	c Input :
63	c -------
64	c pucov covariant zonal velocity
65	c pvcov covariant meridional velocity
66	c pteta potential temperature
67	c pps surface pressure
68	c pmasse masse d'air dans chaque maille
69	c pts surface temperature (K)
70	c callrad clef d'appel au rayonnement
71	c
72	c Output :
73	c --------
74	c pdufi tendency for the natural zonal velocity (ms-1)
75	c pdvfi tendency for the natural meridional velocity
76	c pdhfi tendency for the potential temperature
77	c pdtsfi tendency for the surface temperature
78	c
79	c pdtrad radiative tendencies \ both input
80	c pfluxrad radiative fluxes / and output
81	c
82	c=======================================================================
83	c
84	c-----------------------------------------------------------------------
85	c
86	c 0. Declarations :
87	c ------------------
88
89	#include "dimensions.h"
90	#include "paramet.h"
91	#include "temps.h"
92
93	INTEGER ngridmx
94	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
95
96	#include "comconst.h"
97	#include "comvert.h"
98	#include "comgeom2.h"
99	#include "control.h"
100
101	c Arguments :
102	c -----------
103	LOGICAL lafin
104	REAL heure
105
106	REAL pvcov(iip1,jjm,llm)
107	REAL pucov(iip1,jjp1,llm)
108	REAL pteta(iip1,jjp1,llm)
109	REAL pmasse(iip1,jjp1,llm)
110	REAL pq(iip1,jjp1,llm,nqtot)
111	REAL pphis(iip1,jjp1)
112	REAL pphi(iip1,jjp1,llm)
113	c
114	REAL pdvcov(iip1,jjm,llm)
115	REAL pducov(iip1,jjp1,llm)
116	REAL pdteta(iip1,jjp1,llm)
117	REAL pdq(iip1,jjp1,llm,nqtot)
118	c
119	REAL pps(iip1,jjp1)
120	REAL pp(iip1,jjp1,llmp1)
121	REAL ppk(iip1,jjp1,llm)
122	c
123	REAL pdvfi(iip1,jjm,llm)
124	REAL pdufi(iip1,jjp1,llm)
125	REAL pdhfi(iip1,jjp1,llm)
126	REAL pdqfi(iip1,jjp1,llm,nqtot)
127	REAL pdpsfi(iip1,jjp1)
128
129	INTEGER longcles
130	PARAMETER ( longcles = 20 )
131	REAL clesphy0( longcles )
132
133
134	c Local variables :
135	c -----------------
136
137	INTEGER i,j,l,ig0,ig,iq,iiq
138	REAL zpsrf(ngridmx)
139	REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
140	REAL zphi(ngridmx,llm),zphis(ngridmx)
141	c
142	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
143	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
144	c
145	REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
146	REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
147	c
148	REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
149	REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
150	REAL zdpsrf(ngridmx)
151	c
152	REAL zsin(iim),zcos(iim),z1(iim)
153	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
154	REAL unskap, pksurcp
155	c
156	cIM diagnostique PVteta, Amip2
157	INTEGER ntetaSTD
158	PARAMETER(ntetaSTD=3)
159	REAL rtetaSTD(ntetaSTD)
160	DATA rtetaSTD/350., 380., 405./
161	REAL PVteta(ngridmx,ntetaSTD)
162	c
163	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
164	REAL flxwfi(ngridmx,llm) ! Flux de masse verticale sur la grille physiq
165	c
166
167	REAL SSUM
168
169	LOGICAL firstcal, debut
170	DATA firstcal/.true./
171	SAVE firstcal,debut
172	REAL rdayvrai
173	c
174	c-----------------------------------------------------------------------
175	c
176	c 1. Initialisations :
177	c --------------------
178	c
179
180	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
181	PRINT*,'STOP dans calfis'
182	PRINT,'La dimension ngridmx doit etre egale a 2 + (jjm-1)iim'
183	PRINT*,' ngridmx jjm iim '
184	PRINT*,ngridmx,jjm,iim
185	STOP
186	ENDIF
187
188	c-----------------------------------------------------------------------
189	c latitude, longitude et aires des mailles pour la physique:
190	c ----------------------------------------------------------
191
192	c
193	IF ( firstcal ) THEN
194	debut = .TRUE.
195	ELSE
196	debut = .FALSE.
197	ENDIF
198
199	c
200	c
201	c-----------------------------------------------------------------------
202	c 40. transformation des variables dynamiques en variables physiques:
203	c ---------------------------------------------------------------
204
205	c 41. pressions au sol (en Pascals)
206	c ----------------------------------
207
208
209	zpsrf(1) = pps(1,1)
210
211	ig0 = 2
212	DO j = 2,jjm
213	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
214	ig0 = ig0+iim
215	ENDDO
216
217	zpsrf(ngridmx) = pps(1,jjp1)
218
219
220	c 42. pression intercouches :
221	c
222	c -----------------------------------------------------------------
223	c .... zplev definis aux (llm +1) interfaces des couches ....
224	c .... zplay definis aux ( llm ) milieux des couches ....
225	c -----------------------------------------------------------------
226
227	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
228	c
229	unskap = 1./ kappa
230	c
231	DO l = 1, llmp1
232	zplev( 1,l ) = pp(1,1,l)
233	ig0 = 2
234	DO j = 2, jjm
235	DO i =1, iim
236	zplev( ig0,l ) = pp(i,j,l)
237	ig0 = ig0 +1
238	ENDDO
239	ENDDO
240	zplev( ngridmx,l ) = pp(1,jjp1,l)
241	ENDDO
242	c
243	c
244
245	c 43. temperature naturelle (en K) et pressions milieux couches .
246	c ---------------------------------------------------------------
247
248	DO l=1,llm
249
250	pksurcp = ppk(1,1,l) / cpp
251	zplay(1,l) = preff * pksurcp ** unskap
252	ztfi(1,l) = pteta(1,1,l) * pksurcp
253	pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
254	ig0 = 2
255
256	DO j = 2, jjm
257	DO i = 1, iim
258	pksurcp = ppk(i,j,l) / cpp
259	zplay(ig0,l) = preff * pksurcp ** unskap
260	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
261	pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
262	ig0 = ig0 + 1
263	ENDDO
264	ENDDO
265
266	pksurcp = ppk(1,jjp1,l) / cpp
267	zplay(ig0,l) = preff * pksurcp ** unskap
268	ztfi (ig0,l) = pteta(1,jjp1,l) * pksurcp
269	pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
270
271	ENDDO
272
273	c 43.bis traceurs
274	c ---------------
275	c
276	DO iq=1,nqtot
277	iiq=niadv(iq)
278	DO l=1,llm
279	zqfi(1,l,iq) = pq(1,1,l,iiq)
280	ig0 = 2
281	DO j=2,jjm
282	DO i = 1, iim
283	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
284	ig0 = ig0 + 1
285	ENDDO
286	ENDDO
287	zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
288	ENDDO
289	ENDDO
290
291	c convergence dynamique pour les traceurs "EAU"
292
293	DO iq=1,2
294	DO l=1,llm
295	pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
296	ig0 = 2
297	DO j=2,jjm
298	DO i = 1, iim
299	pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
300	ig0 = ig0 + 1
301	ENDDO
302	ENDDO
303	pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
304	ENDDO
305	ENDDO
306
307
308	c Geopotentiel calcule par rapport a la surface locale:
309	c -----------------------------------------------------
310
311	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
312	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
313	DO l=1,llm
314	DO ig=1,ngridmx
315	zphi(ig,l)=zphi(ig,l)-zphis(ig)
316	ENDDO
317	ENDDO
318
319	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
320	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
321	c de masse est calclue dans advtrac.F
322	c DO l=1,llm
323	c pvervel(1,l)=pw(1,1,l) * g /apoln
324	c ig0=2
325	c DO j=2,jjm
326	c DO i = 1, iim
327	c pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
328	c ig0 = ig0 + 1
329	c ENDDO
330	c ENDDO
331	c pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
332	c ENDDO
333
334	c
335	c 45. champ u:
336	c ------------
337
338	DO 50 l=1,llm
339
340	DO 25 j=2,jjm
341	ig0 = 1+(j-2)*iim
342	zufi(ig0+1,l)= 0.5 *
343	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
344	pcvgu(ig0+1,l)= 0.5 *
345	$ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
346	DO 10 i=2,iim
347	zufi(ig0+i,l)= 0.5 *
348	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
349	pcvgu(ig0+i,l)= 0.5 *
350	$ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
351	10 CONTINUE
352	25 CONTINUE
353
354	50 CONTINUE
355
356
357	c 46.champ v:
358	c -----------
359
360	DO l=1,llm
361	DO j=2,jjm
362	ig0=1+(j-2)*iim
363	DO i=1,iim
364	zvfi(ig0+i,l)= 0.5 *
365	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
366	pcvgv(ig0+i,l)= 0.5 *
367	$ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
368	ENDDO
369	ENDDO
370	ENDDO
371
372
373	c 47. champs de vents aux pole nord
374	c ------------------------------
375	c U = 1 / pi * integrale [ v * cos(long) * d long ]
376	c V = 1 / pi * integrale [ v * sin(long) * d long ]
377
378	DO l=1,llm
379
380	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
381	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
382	DO i=2,iim
383	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
384	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
385	ENDDO
386
387	DO i=1,iim
388	zcos(i) = COS(rlonv(i))*z1(i)
389	zcosbis(i)= COS(rlonv(i))*z1bis(i)
390	zsin(i) = SIN(rlonv(i))*z1(i)
391	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
392	ENDDO
393
394	zufi(1,l) = SSUM(iim,zcos,1)/pi
395	pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
396	zvfi(1,l) = SSUM(iim,zsin,1)/pi
397	pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
398
399	ENDDO
400
401
402	c 48. champs de vents aux pole sud:
403	c ---------------------------------
404	c U = 1 / pi * integrale [ v * cos(long) * d long ]
405	c V = 1 / pi * integrale [ v * sin(long) * d long ]
406
407	DO l=1,llm
408
409	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
410	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
411	DO i=2,iim
412	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
413	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
414	ENDDO
415
416	DO i=1,iim
417	zcos(i) = COS(rlonv(i))*z1(i)
418	zcosbis(i) = COS(rlonv(i))*z1bis(i)
419	zsin(i) = SIN(rlonv(i))*z1(i)
420	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
421	ENDDO
422
423	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
424	pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
425	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
426	pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
427
428	ENDDO
429	c
430	cIM calcul PV a teta=350, 380, 405K
431	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
432	$ ztfi,zplay,zplev,
433	$ ntetaSTD,rtetaSTD,PVteta)
434	c
435	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
436	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
437
438	c-----------------------------------------------------------------------
439	c Appel de la physique:
440	c ---------------------
441
442
443	CALL physiq (ngridmx,
444	. llm,
445	. debut,
446	. lafin,
447	. rdayvrai,
448	. heure,
449	. dtphys,
450	. zplev,
451	. zplay,
452	. zphi,
453	. zphis,
454	. presnivs,
455	. clesphy0,
456	. zufi,
457	. zvfi,
458	. ztfi,
459	. zqfi,
460	. flxwfi,
461	. zdufi,
462	. zdvfi,
463	. zdtfi,
464	. zdqfi,
465	. zdpsrf,
466	cIM diagnostique PVteta, Amip2
467	. pducov,
468	. PVteta)
469
470	500 CONTINUE
471
472	c-----------------------------------------------------------------------
473	c transformation des tendances physiques en tendances dynamiques:
474	c ---------------------------------------------------------------
475
476	c tendance sur la pression :
477	c -----------------------------------
478
479	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
480	c
481	c 62. enthalpie potentielle
482	c ---------------------
483
484	DO l=1,llm
485
486	DO i=1,iip1
487	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
488	pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
489	ENDDO
490
491	DO j=2,jjm
492	ig0=1+(j-2)*iim
493	DO i=1,iim
494	pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
495	ENDDO
496	pdhfi(iip1,j,l) = pdhfi(1,j,l)
497	ENDDO
498
499	ENDDO
500
501
502	c 62. humidite specifique
503	c ---------------------
504
505	DO iq=1,nqtot
506	DO l=1,llm
507	DO i=1,iip1
508	pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
509	pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
510	ENDDO
511	DO j=2,jjm
512	ig0=1+(j-2)*iim
513	DO i=1,iim
514	pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
515	ENDDO
516	pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
517	ENDDO
518	ENDDO
519	ENDDO
520
521	c 63. traceurs
522	c ------------
523	C initialisation des tendances
524	pdqfi=0.
525	C
526	DO iq=1,nqtot
527	iiq=niadv(iq)
528	DO l=1,llm
529	DO i=1,iip1
530	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
531	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
532	ENDDO
533	DO j=2,jjm
534	ig0=1+(j-2)*iim
535	DO i=1,iim
536	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
537	ENDDO
538	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
539	ENDDO
540	ENDDO
541	ENDDO
542
543	c 65. champ u:
544	c ------------
545
546	DO l=1,llm
547
548	DO i=1,iip1
549	pdufi(i,1,l) = 0.
550	pdufi(i,jjp1,l) = 0.
551	ENDDO
552
553	DO j=2,jjm
554	ig0=1+(j-2)*iim
555	DO i=1,iim-1
556	pdufi(i,j,l)=
557	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
558	ENDDO
559	pdufi(iim,j,l)=
560	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
561	pdufi(iip1,j,l)=pdufi(1,j,l)
562	ENDDO
563
564	ENDDO
565
566
567	c 67. champ v:
568	c ------------
569
570	DO l=1,llm
571
572	DO j=2,jjm-1
573	ig0=1+(j-2)*iim
574	DO i=1,iim
575	pdvfi(i,j,l)=
576	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
577	ENDDO
578	pdvfi(iip1,j,l) = pdvfi(1,j,l)
579	ENDDO
580	ENDDO
581
582
583	c 68. champ v pres des poles:
584	c ---------------------------
585	c v = U * cos(long) + V * SIN(long)
586
587	DO l=1,llm
588
589	DO i=1,iim
590	pdvfi(i,1,l)=
591	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
592	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
593	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
594	pdvfi(i,1,l)=
595	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
596	pdvfi(i,jjm,l)=
597	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
598	ENDDO
599
600	pdvfi(iip1,1,l) = pdvfi(1,1,l)
601	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
602
603	ENDDO
604
605	c-----------------------------------------------------------------------
606
607	700 CONTINUE
608
609	firstcal = .FALSE.
610
611	RETURN
612	END

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

Download in other formats:

Original Format