Context Navigation

calfis.F @ 1160

Last change on this file since 1160 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
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(lafin,
7	$ rdayvrai,
8	$ heure,
9	$ pucov,
10	$ pvcov,
11	$ pteta,
12	$ pq,
13	$ pmasse,
14	$ pps,
15	$ pp,
16	$ ppk,
17	$ pphis,
18	$ pphi,
19	$ pducov,
20	$ pdvcov,
21	$ pdteta,
22	$ pdq,
23	$ flxw,
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	USE infotrac
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
93	INTEGER ngridmx
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,nqtot)
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,nqtot)
118	c
119	REAL pps(iip1,jjp1)
120	REAL pp(iip1,jjp1,llmp1)
121	REAL ppk(iip1,jjp1,llm)
122	c
123	REAL pdvfi(iip1,jjm,llm)
124	REAL pdufi(iip1,jjp1,llm)
125	REAL pdhfi(iip1,jjp1,llm)
126	REAL pdqfi(iip1,jjp1,llm,nqtot)
127	REAL pdpsfi(iip1,jjp1)
128
129	INTEGER longcles
130	PARAMETER ( longcles = 20 )
131	REAL clesphy0( longcles )
132
133
134	c Local variables :
135	c -----------------
136
137	INTEGER i,j,l,ig0,ig,iq,iiq
138	REAL zpsrf(ngridmx)
139	REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
140	REAL zphi(ngridmx,llm),zphis(ngridmx)
141	c
142	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
143	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
144	c
145	REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
146	REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
147	c
148	REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
149	REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
150	REAL zdpsrf(ngridmx)
151	c
152	REAL zsin(iim),zcos(iim),z1(iim)
153	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
154	REAL unskap, pksurcp
155	c
156	cIM diagnostique PVteta, Amip2
157	INTEGER ntetaSTD
158	PARAMETER(ntetaSTD=3)
159	REAL rtetaSTD(ntetaSTD)
160	DATA rtetaSTD/350., 380., 405./
161	REAL PVteta(ngridmx,ntetaSTD)
162	c
163	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
164	REAL flxwfi(ngridmx,llm) ! Flux de masse verticale sur la grille physiq
165	c
166
167	REAL SSUM
168
169	LOGICAL firstcal, debut
170	DATA firstcal/.true./
171	SAVE firstcal,debut
172	REAL rdayvrai
173	c
174	c-----------------------------------------------------------------------
175	c
176	c 1. Initialisations :
177	c --------------------
178	c
179
180	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
181	PRINT*,'STOP dans calfis'
182	PRINT,'La dimension ngridmx doit etre egale a 2 + (jjm-1)iim'
183	PRINT*,' ngridmx jjm iim '
184	PRINT*,ngridmx,jjm,iim
185	STOP
186	ENDIF
187
188	c-----------------------------------------------------------------------
189	c latitude, longitude et aires des mailles pour la physique:
190	c ----------------------------------------------------------
191
192	c
193	IF ( firstcal ) THEN
194	debut = .TRUE.
195	ELSE
196	debut = .FALSE.
197	ENDIF
198
199	c
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 traceurs
274	c ---------------
275	c
276	DO iq=1,nqtot
277	iiq=niadv(iq)
278	DO l=1,llm
279	zqfi(1,l,iq) = pq(1,1,l,iiq)
280	ig0 = 2
281	DO j=2,jjm
282	DO i = 1, iim
283	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
284	ig0 = ig0 + 1
285	ENDDO
286	ENDDO
287	zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
288	ENDDO
289	ENDDO
290
291	c convergence dynamique pour les traceurs "EAU"
292
293	DO iq=1,2
294	DO l=1,llm
295	pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
296	ig0 = 2
297	DO j=2,jjm
298	DO i = 1, iim
299	pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
300	ig0 = ig0 + 1
301	ENDDO
302	ENDDO
303	pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
304	ENDDO
305	ENDDO
306
307
308	c Geopotentiel calcule par rapport a la surface locale:
309	c -----------------------------------------------------
310
311	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
312	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
313	DO l=1,llm
314	DO ig=1,ngridmx
315	zphi(ig,l)=zphi(ig,l)-zphis(ig)
316	ENDDO
317	ENDDO
318
319	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
320	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
321	c de masse est calclue dans advtrac.F
322	c DO l=1,llm
323	c pvervel(1,l)=pw(1,1,l) * g /apoln
324	c ig0=2
325	c DO j=2,jjm
326	c DO i = 1, iim
327	c pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
328	c ig0 = ig0 + 1
329	c ENDDO
330	c ENDDO
331	c pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
332	c ENDDO
333
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	c
430	cIM calcul PV a teta=350, 380, 405K
431	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
432	$ ztfi,zplay,zplev,
433	$ ntetaSTD,rtetaSTD,PVteta)
434	c
435	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
436	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
437
438	c-----------------------------------------------------------------------
439	c Appel de la physique:
440	c ---------------------
441
442
443	CALL physiq (ngridmx,
444	. llm,
445	. debut,
446	. lafin,
447	. rdayvrai,
448	. heure,
449	. dtphys,
450	. zplev,
451	. zplay,
452	. zphi,
453	. zphis,
454	. presnivs,
455	. clesphy0,
456	. zufi,
457	. zvfi,
458	. ztfi,
459	. zqfi,
460	. flxwfi,
461	. zdufi,
462	. zdvfi,
463	. zdtfi,
464	. zdqfi,
465	. zdpsrf,
466	cIM diagnostique PVteta, Amip2
467	. pducov,
468	. PVteta)
469
470	500 CONTINUE
471
472	c-----------------------------------------------------------------------
473	c transformation des tendances physiques en tendances dynamiques:
474	c ---------------------------------------------------------------
475
476	c tendance sur la pression :
477	c -----------------------------------
478
479	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
480	c
481	c 62. enthalpie potentielle
482	c ---------------------
483
484	DO l=1,llm
485
486	DO i=1,iip1
487	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
488	pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
489	ENDDO
490
491	DO j=2,jjm
492	ig0=1+(j-2)*iim
493	DO i=1,iim
494	pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
495	ENDDO
496	pdhfi(iip1,j,l) = pdhfi(1,j,l)
497	ENDDO
498
499	ENDDO
500
501
502	c 62. humidite specifique
503	c ---------------------
504
505	DO iq=1,nqtot
506	DO l=1,llm
507	DO i=1,iip1
508	pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
509	pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
510	ENDDO
511	DO j=2,jjm
512	ig0=1+(j-2)*iim
513	DO i=1,iim
514	pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
515	ENDDO
516	pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
517	ENDDO
518	ENDDO
519	ENDDO
520
521	c 63. traceurs
522	c ------------
523	C initialisation des tendances
524	pdqfi=0.
525	C
526	DO iq=1,nqtot
527	iiq=niadv(iq)
528	DO l=1,llm
529	DO i=1,iip1
530	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
531	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
532	ENDDO
533	DO j=2,jjm
534	ig0=1+(j-2)*iim
535	DO i=1,iim
536	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
537	ENDDO
538	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
539	ENDDO
540	ENDDO
541	ENDDO
542
543	c 65. champ u:
544	c ------------
545
546	DO l=1,llm
547
548	DO i=1,iip1
549	pdufi(i,1,l) = 0.
550	pdufi(i,jjp1,l) = 0.
551	ENDDO
552
553	DO j=2,jjm
554	ig0=1+(j-2)*iim
555	DO i=1,iim-1
556	pdufi(i,j,l)=
557	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
558	ENDDO
559	pdufi(iim,j,l)=
560	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
561	pdufi(iip1,j,l)=pdufi(1,j,l)
562	ENDDO
563
564	ENDDO
565
566
567	c 67. champ v:
568	c ------------
569
570	DO l=1,llm
571
572	DO j=2,jjm-1
573	ig0=1+(j-2)*iim
574	DO i=1,iim
575	pdvfi(i,j,l)=
576	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
577	ENDDO
578	pdvfi(iip1,j,l) = pdvfi(1,j,l)
579	ENDDO
580	ENDDO
581
582
583	c 68. champ v pres des poles:
584	c ---------------------------
585	c v = U * cos(long) + V * SIN(long)
586
587	DO l=1,llm
588
589	DO i=1,iim
590	pdvfi(i,1,l)=
591	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
592	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
593	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
594	pdvfi(i,1,l)=
595	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
596	pdvfi(i,jjm,l)=
597	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
598	ENDDO
599
600	pdvfi(iip1,1,l) = pdvfi(1,1,l)
601	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
602
603	ENDDO
604
605	c-----------------------------------------------------------------------
606
607	700 CONTINUE
608
609	firstcal = .FALSE.
610
611	RETURN
612	END

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

Download in other formats:

Original Format