Context Navigation

calfis.F @ 932

Last change on this file since 932 was 644, checked in by Laurent Fairhead, 19 years ago
Synchronisation avec tous les diagnostiques de Ionela IM Inclusion du slab ocean IM 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(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