Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

calfis.F @ 5440

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

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

Download in other formats:

Original Format