Context Navigation

calfis.F @ 1206

Last change on this file since 1206 was 1201, checked in by Laurent Fairhead, 15 years ago
Modifications nécessaires a l'inclusion d'un calendrier réaliste. La date courante est calculée dans leapfrog.F et exprimée en Jour Julien (modifié). On en a profité pour faire un peu de ménage dans la gestion des dates du modèle. Dans la physique, on utilise les routines de passages entre calendrier Julien et Gregorien incluses dans IOIPSL pour calculer le nombre de jours écoulés depuis le 1er janvier (pour les conditions aux limites) ou l'equinoxe (pour le calcul de la longitude solaire). Le calcul de l'orbite reprend celui du gcm planétaire (codé par FH) On décide du calendrier à utiliser à l'aide du paramètre calend du run.def. Par défaut celui-ci est à earth_360d 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	! $Id: calfis.F 1201 2009-07-07 14:01:00Z emillour $
3	!
4	C
5	C
6	SUBROUTINE calfis(lafin,
7	$ jD_cur, jH_cur,
8	$ pucov,
9	$ pvcov,
10	$ pteta,
11	$ pq,
12	$ pmasse,
13	$ pps,
14	$ pp,
15	$ ppk,
16	$ pphis,
17	$ pphi,
18	$ pducov,
19	$ pdvcov,
20	$ pdteta,
21	$ pdq,
22	$ flxw,
23	$ clesphy0,
24	$ pdufi,
25	$ pdvfi,
26	$ pdhfi,
27	$ pdqfi,
28	$ pdpsfi)
29	c
30	c Auteur : P. Le Van, F. Hourdin
31	c .........
32	USE infotrac
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 pucov covariant zonal velocity
64	c pvcov covariant meridional velocity
65	c pteta potential temperature
66	c pps surface pressure
67	c pmasse masse d'air dans chaque maille
68	c pts surface temperature (K)
69	c callrad clef d'appel au rayonnement
70	c
71	c Output :
72	c --------
73	c pdufi tendency for the natural zonal velocity (ms-1)
74	c pdvfi tendency for the natural meridional velocity
75	c pdhfi tendency for the potential temperature
76	c pdtsfi tendency for the surface temperature
77	c
78	c pdtrad radiative tendencies \ both input
79	c pfluxrad radiative fluxes / and output
80	c
81	c=======================================================================
82	c
83	c-----------------------------------------------------------------------
84	c
85	c 0. Declarations :
86	c ------------------
87
88	#include "dimensions.h"
89	#include "paramet.h"
90	#include "temps.h"
91
92	INTEGER ngridmx
93	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
94
95	#include "comconst.h"
96	#include "comvert.h"
97	#include "comgeom2.h"
98	#include "control.h"
99
100	c Arguments :
101	c -----------
102	LOGICAL lafin
103
104
105	REAL pvcov(iip1,jjm,llm)
106	REAL pucov(iip1,jjp1,llm)
107	REAL pteta(iip1,jjp1,llm)
108	REAL pmasse(iip1,jjp1,llm)
109	REAL pq(iip1,jjp1,llm,nqtot)
110	REAL pphis(iip1,jjp1)
111	REAL pphi(iip1,jjp1,llm)
112	c
113	REAL pdvcov(iip1,jjm,llm)
114	REAL pducov(iip1,jjp1,llm)
115	REAL pdteta(iip1,jjp1,llm)
116	REAL pdq(iip1,jjp1,llm,nqtot)
117	c
118	REAL pps(iip1,jjp1)
119	REAL pp(iip1,jjp1,llmp1)
120	REAL ppk(iip1,jjp1,llm)
121	c
122	REAL pdvfi(iip1,jjm,llm)
123	REAL pdufi(iip1,jjp1,llm)
124	REAL pdhfi(iip1,jjp1,llm)
125	REAL pdqfi(iip1,jjp1,llm,nqtot)
126	REAL pdpsfi(iip1,jjp1)
127
128	INTEGER longcles
129	PARAMETER ( longcles = 20 )
130	REAL clesphy0( longcles )
131
132
133	c Local variables :
134	c -----------------
135
136	INTEGER i,j,l,ig0,ig,iq,iiq
137	REAL zpsrf(ngridmx)
138	REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
139	REAL zphi(ngridmx,llm),zphis(ngridmx)
140	c
141	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
142	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
143	c
144	REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
145	REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
146	c
147	REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
148	REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
149	REAL zdpsrf(ngridmx)
150	c
151	REAL zsin(iim),zcos(iim),z1(iim)
152	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
153	REAL unskap, pksurcp
154	c
155	cIM diagnostique PVteta, Amip2
156	INTEGER ntetaSTD
157	PARAMETER(ntetaSTD=3)
158	REAL rtetaSTD(ntetaSTD)
159	DATA rtetaSTD/350., 380., 405./
160	REAL PVteta(ngridmx,ntetaSTD)
161	c
162	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
163	REAL flxwfi(ngridmx,llm) ! Flux de masse verticale sur la grille physiq
164	c
165
166	REAL SSUM
167
168	LOGICAL firstcal, debut
169	DATA firstcal/.true./
170	SAVE firstcal,debut
171	! REAL rdayvrai
172	REAL :: jD_cur, jH_cur
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	. jD_cur,
448	. jH_cur,
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