Context Navigation

calfis.F @ 5306

Last change on this file since 5306 was 344, checked in by lmdz, 23 years ago
Inclusion des modifs de D. Hauglustaine pour la version 1 de INCA LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 16.6 KB

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

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

Download in other formats:

Original Format