Context Navigation

calfis.F @ 1006

Last change on this file since 1006 was 960, checked in by lsce, 16 years ago

Ajoute du parametre config_inca dans conf_gcm.F config_inca='none'(sans INCA, par defaut) config_inca='chem'(avec INCA config chemie) config_inca='aero'(avec INCA config aerosol)
Menage parmis les cles CPP INCA
Enleve le calcul d'omega dans calfis.F et active le calcul correspondant dans advtrac.F(avant uniquement pour INCA).

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

File size: 16.7 KB

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

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

Download in other formats:

Original Format