Context Navigation

calfis.F @ 115

Last change on this file since 115 was 79, checked in by (none), 25 years ago
This commit was manufactured by cvs2svn to create branch 'rel-LF'.
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 15.4 KB

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

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

Download in other formats:

Original Format