Context Navigation

calfis.F @ 201

Last change on this file since 201 was 177, checked in by lmdzadmin, 24 years ago

Lots of stuff, plus particulierement:

appel a ORCHIDEE en etat de marche (pb de grille subsiste)
modifs de Pascale sur soil dans le cas ou ok_veget=false -

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

File size: 15.7 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	REAL zcufi(ngridmx),zcvfi(ngridmx)
131	c
132	REAL zsin(iim),zcos(iim),z1(iim)
133	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
134	REAL unskap, pksurcp
135	c
136
137	EXTERNAL gr_dyn_fi,gr_fi_dyn
138	EXTERNAL physiq,multipl
139	REAL SSUM
140	EXTERNAL SSUM
141
142	REAL latfi(ngridmx),lonfi(ngridmx)
143	REAL airefi(ngridmx)
144	SAVE latfi, lonfi, airefi
145
146	LOGICAL firstcal, debut
147	DATA firstcal/.true./
148	SAVE firstcal,debut
149	REAL rdayvrai,rday_ecri
150	c
151	c-----------------------------------------------------------------------
152	c
153	c 1. Initialisations :
154	c --------------------
155	c
156
157	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
158	PRINT*,'STOP dans calfis'
159	PRINT,'La dimension ngridmx doit etre egale a 2 + (jjm-1)iim'
160	PRINT*,' ngridmx jjm iim '
161	PRINT*,ngridmx,jjm,iim
162	STOP
163	ENDIF
164
165	c-----------------------------------------------------------------------
166	c latitude, longitude et aires des mailles pour la physique:
167	c ----------------------------------------------------------
168
169	c
170	IF ( firstcal ) THEN
171	debut = .TRUE.
172	ELSE
173	debut = .FALSE.
174	ENDIF
175
176	c
177	IF (firstcal) THEN
178	latfi(1)=rlatu(1)
179	lonfi(1)=0.
180	zcufi(1) = cu(1,1)
181	zcvfi(1) = cv(1,1)
182	DO j=2,jjm
183	DO i=1,iim
184	latfi((j-2)*iim+1+i)= rlatu(j)
185	lonfi((j-2)*iim+1+i)= rlonv(i)
186	zcufi((j-2)*iim+1+i) = cu(i,j)
187	zcvfi((j-2)*iim+1+i) = cv(i,j)
188	ENDDO
189	ENDDO
190	latfi(ngridmx)= rlatu(jjp1)
191	lonfi(ngridmx)= 0.
192	zcufi(ngridmx) = cu(1,jjp1)
193	zcvfi(ngridmx) = cv(1,jjm)
194	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
195	PRINT*,'WARNING!!! vitesse verticale nulle dans la physique'
196	CALL inifis(ngridmx,llm,daysec,day_ini,dtphys ,
197	, latfi,lonfi,airefi,rad,g,r,cpp )
198	ENDIF
199
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 humidite specifique (en kg/kg)
274	c -------------------------------------
275
276	DO iq=1,nqmx
277	DO l=1,llm
278	zqfi(1,l,iq) = pq(1,1,l,iq)
279	ig0 = 2
280	DO j=2,jjm
281	DO i = 1, iim
282	zqfi(ig0,l,iq) = pq(i,j,l,iq)
283	ig0 = ig0 + 1
284	ENDDO
285	ENDDO
286	zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
287	ENDDO
288	ENDDO
289
290	c convergence dynamique pour les traceurs "EAU"
291
292	DO iq=1,2
293	DO l=1,llm
294	pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
295	ig0 = 2
296	DO j=2,jjm
297	DO i = 1, iim
298	pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
299	ig0 = ig0 + 1
300	ENDDO
301	ENDDO
302	pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
303	ENDDO
304	ENDDO
305
306
307	c Geopotentiel calcule par rapport a la surface locale:
308	c -----------------------------------------------------
309
310	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
311	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
312	DO l=1,llm
313	DO ig=1,ngridmx
314	zphi(ig,l)=zphi(ig,l)-zphis(ig)
315	ENDDO
316	ENDDO
317
318	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
319	c
320	DO l=1,llm
321	pvervel(1,l)=pw(1,1,l)/apoln
322	ig0=2
323	DO j=2,jjm
324	DO i = 1, iim
325	pvervel(ig0,l) = pw(i,j,l) * unsaire(i,j)
326	ig0 = ig0 + 1
327	ENDDO
328	ENDDO
329	pvervel(ig0,l)=pw(1,jjp1,l)/apols
330	ENDDO
331
332	c CALL initial0( ngridmx*llm,pvervel )
333	c
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
430	c-----------------------------------------------------------------------
431	c Appel de la physique:
432	c ---------------------
433
434
435	CALL physiq (ngridmx,llm,nq,debut,lafin,
436	, rdayvrai,rday_ecri,heure,dtphys,zplev,zplay,zphi,zphis,airefi,
437	, presnivs,clesphy0, zufi, zvfi,ztfi, zqfi,
438	ccc , pcvgu, pcvgv, pcvgt, pcvgq,
439	, pvervel, zcufi, zcvfi,
440	C - sorties
441	s zdufi, zdvfi, zdtfi, zdqfi,zdpsrf )
442
443	500 CONTINUE
444
445	c-----------------------------------------------------------------------
446	c transformation des tendances physiques en tendances dynamiques:
447	c ---------------------------------------------------------------
448
449	c tendance sur la pression :
450	c -----------------------------------
451
452	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
453	c
454	ccc CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi)
455
456	c 62. enthalpie potentielle
457	c ---------------------
458
459	DO l=1,llm
460
461	DO i=1,iip1
462	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
463	pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
464	ENDDO
465
466	DO j=2,jjm
467	ig0=1+(j-2)*iim
468	DO i=1,iim
469	pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
470	ENDDO
471	pdhfi(iip1,j,l) = pdhfi(1,j,l)
472	ENDDO
473
474	ENDDO
475
476
477	c 62. humidite specifique
478	c ---------------------
479
480	DO iq=1,nqmx
481	DO l=1,llm
482	DO i=1,iip1
483	pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
484	pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
485	ENDDO
486	DO j=2,jjm
487	ig0=1+(j-2)*iim
488	DO i=1,iim
489	pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
490	ENDDO
491	pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
492	ENDDO
493	ENDDO
494	ENDDO
495
496	c 65. champ u:
497	c ------------
498
499	DO l=1,llm
500
501	DO i=1,iip1
502	pdufi(i,1,l) = 0.
503	pdufi(i,jjp1,l) = 0.
504	ENDDO
505
506	DO j=2,jjm
507	ig0=1+(j-2)*iim
508	DO i=1,iim-1
509	pdufi(i,j,l)=
510	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
511	ENDDO
512	pdufi(iim,j,l)=
513	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
514	pdufi(iip1,j,l)=pdufi(1,j,l)
515	ENDDO
516
517	ENDDO
518
519
520	c 67. champ v:
521	c ------------
522
523	DO l=1,llm
524
525	DO j=2,jjm-1
526	ig0=1+(j-2)*iim
527	DO i=1,iim
528	pdvfi(i,j,l)=
529	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
530	ENDDO
531	pdvfi(iip1,j,l) = pdvfi(1,j,l)
532	ENDDO
533	ENDDO
534
535
536	c 68. champ v pres des poles:
537	c ---------------------------
538	c v = U * cos(long) + V * SIN(long)
539
540	DO l=1,llm
541
542	DO i=1,iim
543	pdvfi(i,1,l)=
544	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
545	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
546	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
547	pdvfi(i,1,l)=
548	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
549	pdvfi(i,jjm,l)=
550	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
551	ENDDO
552
553	pdvfi(iip1,1,l) = pdvfi(1,1,l)
554	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
555
556	ENDDO
557
558	c-----------------------------------------------------------------------
559
560	700 CONTINUE
561
562	firstcal = .FALSE.
563
564	RETURN
565	END

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

Download in other formats:

Original Format