Context Navigation

calfis.F @ 5065

Last change on this file since 5065 was 1225, checked in by lguez, 15 years ago

Added some "intent" attributes in declarations.

In "phyredem", "dtime" is not declared. It is not in any included
file. Probably accepted by compilers as an intrinsic non-standard
function. Removed this element of "tab_cntrl".

Added some "only" clauses in "use" statements.

If the ozone field is read from a file, it is now updated every
360th of the length of the current year, regardless of that length.

In "physiq", "omega" was output before it was defined. Moved the
output instruction after the definition.

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

File size: 16.8 KB

Line
1	!
2	! $Id: calfis.F 1225 2009-08-07 16:29:23Z abarral $
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, intent(in):: jD_cur, jH_cur
173	c
174	c-----------------------------------------------------------------------
175	c
176	c 1. Initialisations :
177	c --------------------
178	c
179	c
180	IF ( firstcal ) THEN
181	debut = .TRUE.
182	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
183	PRINT*,'STOP dans calfis'
184	PRINT,'La dimension ngridmx doit etre egale a 2 + (jjm-1)iim'
185	PRINT*,' ngridmx jjm iim '
186	PRINT*,ngridmx,jjm,iim
187	STOP
188	ENDIF
189	ELSE
190	debut = .FALSE.
191	ENDIF ! of IF (firstcal)
192
193	c
194	c
195	c-----------------------------------------------------------------------
196	c 40. transformation des variables dynamiques en variables physiques:
197	c ---------------------------------------------------------------
198
199	c 41. pressions au sol (en Pascals)
200	c ----------------------------------
201
202
203	zpsrf(1) = pps(1,1)
204
205	ig0 = 2
206	DO j = 2,jjm
207	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
208	ig0 = ig0+iim
209	ENDDO
210
211	zpsrf(ngridmx) = pps(1,jjp1)
212
213
214	c 42. pression intercouches :
215	c
216	c -----------------------------------------------------------------
217	c .... zplev definis aux (llm +1) interfaces des couches ....
218	c .... zplay definis aux ( llm ) milieux des couches ....
219	c -----------------------------------------------------------------
220
221	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
222	c
223	unskap = 1./ kappa
224	c
225	DO l = 1, llmp1
226	zplev( 1,l ) = pp(1,1,l)
227	ig0 = 2
228	DO j = 2, jjm
229	DO i =1, iim
230	zplev( ig0,l ) = pp(i,j,l)
231	ig0 = ig0 +1
232	ENDDO
233	ENDDO
234	zplev( ngridmx,l ) = pp(1,jjp1,l)
235	ENDDO
236	c
237	c
238
239	c 43. temperature naturelle (en K) et pressions milieux couches .
240	c ---------------------------------------------------------------
241
242	DO l=1,llm
243
244	pksurcp = ppk(1,1,l) / cpp
245	zplay(1,l) = preff * pksurcp ** unskap
246	ztfi(1,l) = pteta(1,1,l) * pksurcp
247	pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
248	ig0 = 2
249
250	DO j = 2, jjm
251	DO i = 1, iim
252	pksurcp = ppk(i,j,l) / cpp
253	zplay(ig0,l) = preff * pksurcp ** unskap
254	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
255	pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
256	ig0 = ig0 + 1
257	ENDDO
258	ENDDO
259
260	pksurcp = ppk(1,jjp1,l) / cpp
261	zplay(ig0,l) = preff * pksurcp ** unskap
262	ztfi (ig0,l) = pteta(1,jjp1,l) * pksurcp
263	pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
264
265	ENDDO
266
267	c 43.bis traceurs
268	c ---------------
269	c
270	DO iq=1,nqtot
271	iiq=niadv(iq)
272	DO l=1,llm
273	zqfi(1,l,iq) = pq(1,1,l,iiq)
274	ig0 = 2
275	DO j=2,jjm
276	DO i = 1, iim
277	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
278	ig0 = ig0 + 1
279	ENDDO
280	ENDDO
281	zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
282	ENDDO
283	ENDDO
284
285	c convergence dynamique pour les traceurs "EAU"
286	! Earth-specific treatment of first 2 tracers (water)
287	if (planet_type=="earth") then
288	DO iq=1,2
289	DO l=1,llm
290	pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
291	ig0 = 2
292	DO j=2,jjm
293	DO i = 1, iim
294	pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
295	ig0 = ig0 + 1
296	ENDDO
297	ENDDO
298	pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
299	ENDDO
300	ENDDO
301	endif ! of if (planet_type=="earth")
302
303
304	c Geopotentiel calcule par rapport a la surface locale:
305	c -----------------------------------------------------
306
307	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
308	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
309	DO l=1,llm
310	DO ig=1,ngridmx
311	zphi(ig,l)=zphi(ig,l)-zphis(ig)
312	ENDDO
313	ENDDO
314
315	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
316	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
317	c de masse est calclue dans advtrac.F
318	c DO l=1,llm
319	c pvervel(1,l)=pw(1,1,l) * g /apoln
320	c ig0=2
321	c DO j=2,jjm
322	c DO i = 1, iim
323	c pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
324	c ig0 = ig0 + 1
325	c ENDDO
326	c ENDDO
327	c pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
328	c ENDDO
329
330	c
331	c 45. champ u:
332	c ------------
333
334	DO 50 l=1,llm
335
336	DO 25 j=2,jjm
337	ig0 = 1+(j-2)*iim
338	zufi(ig0+1,l)= 0.5 *
339	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
340	pcvgu(ig0+1,l)= 0.5 *
341	$ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
342	DO 10 i=2,iim
343	zufi(ig0+i,l)= 0.5 *
344	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
345	pcvgu(ig0+i,l)= 0.5 *
346	$ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
347	10 CONTINUE
348	25 CONTINUE
349
350	50 CONTINUE
351
352
353	c 46.champ v:
354	c -----------
355
356	DO l=1,llm
357	DO j=2,jjm
358	ig0=1+(j-2)*iim
359	DO i=1,iim
360	zvfi(ig0+i,l)= 0.5 *
361	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
362	pcvgv(ig0+i,l)= 0.5 *
363	$ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
364	ENDDO
365	ENDDO
366	ENDDO
367
368
369	c 47. champs de vents aux pole nord
370	c ------------------------------
371	c U = 1 / pi * integrale [ v * cos(long) * d long ]
372	c V = 1 / pi * integrale [ v * sin(long) * d long ]
373
374	DO l=1,llm
375
376	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
377	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
378	DO i=2,iim
379	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
380	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
381	ENDDO
382
383	DO i=1,iim
384	zcos(i) = COS(rlonv(i))*z1(i)
385	zcosbis(i)= COS(rlonv(i))*z1bis(i)
386	zsin(i) = SIN(rlonv(i))*z1(i)
387	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
388	ENDDO
389
390	zufi(1,l) = SSUM(iim,zcos,1)/pi
391	pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
392	zvfi(1,l) = SSUM(iim,zsin,1)/pi
393	pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
394
395	ENDDO
396
397
398	c 48. champs de vents aux pole sud:
399	c ---------------------------------
400	c U = 1 / pi * integrale [ v * cos(long) * d long ]
401	c V = 1 / pi * integrale [ v * sin(long) * d long ]
402
403	DO l=1,llm
404
405	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
406	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
407	DO i=2,iim
408	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
409	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
410	ENDDO
411
412	DO i=1,iim
413	zcos(i) = COS(rlonv(i))*z1(i)
414	zcosbis(i) = COS(rlonv(i))*z1bis(i)
415	zsin(i) = SIN(rlonv(i))*z1(i)
416	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
417	ENDDO
418
419	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
420	pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
421	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
422	pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
423
424	ENDDO
425	c
426	if (planet_type=="earth") then
427	#ifdef CPP_EARTH
428	cIM calcul PV a teta=350, 380, 405K
429	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
430	$ ztfi,zplay,zplev,
431	$ ntetaSTD,rtetaSTD,PVteta)
432	#endif
433	endif
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	if (planet_type=="earth") then
444	#ifdef CPP_EARTH
445	CALL physiq (ngridmx,
446	. llm,
447	. debut,
448	. lafin,
449	. jD_cur,
450	. jH_cur,
451	. dtphys,
452	. zplev,
453	. zplay,
454	. zphi,
455	. zphis,
456	. presnivs,
457	. clesphy0,
458	. zufi,
459	. zvfi,
460	. ztfi,
461	. zqfi,
462	. flxwfi,
463	. zdufi,
464	. zdvfi,
465	. zdtfi,
466	. zdqfi,
467	. zdpsrf,
468	cIM diagnostique PVteta, Amip2
469	. pducov,
470	. PVteta)
471	#endif
472	endif !of if (planet_type=="earth")
473
474	500 CONTINUE
475
476	c-----------------------------------------------------------------------
477	c transformation des tendances physiques en tendances dynamiques:
478	c ---------------------------------------------------------------
479
480	c tendance sur la pression :
481	c -----------------------------------
482
483	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
484	c
485	c 62. enthalpie potentielle
486	c ---------------------
487
488	DO l=1,llm
489
490	DO i=1,iip1
491	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
492	pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
493	ENDDO
494
495	DO j=2,jjm
496	ig0=1+(j-2)*iim
497	DO i=1,iim
498	pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
499	ENDDO
500	pdhfi(iip1,j,l) = pdhfi(1,j,l)
501	ENDDO
502
503	ENDDO
504
505
506	c 62. humidite specifique
507	c ---------------------
508	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
509	! DO iq=1,nqtot
510	! DO l=1,llm
511	! DO i=1,iip1
512	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
513	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
514	! ENDDO
515	! DO j=2,jjm
516	! ig0=1+(j-2)*iim
517	! DO i=1,iim
518	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
519	! ENDDO
520	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
521	! ENDDO
522	! ENDDO
523	! ENDDO
524
525	c 63. traceurs
526	c ------------
527	C initialisation des tendances
528	pdqfi(:,:,:,:)=0.
529	C
530	DO iq=1,nqtot
531	iiq=niadv(iq)
532	DO l=1,llm
533	DO i=1,iip1
534	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
535	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
536	ENDDO
537	DO j=2,jjm
538	ig0=1+(j-2)*iim
539	DO i=1,iim
540	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
541	ENDDO
542	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
543	ENDDO
544	ENDDO
545	ENDDO
546
547	c 65. champ u:
548	c ------------
549
550	DO l=1,llm
551
552	DO i=1,iip1
553	pdufi(i,1,l) = 0.
554	pdufi(i,jjp1,l) = 0.
555	ENDDO
556
557	DO j=2,jjm
558	ig0=1+(j-2)*iim
559	DO i=1,iim-1
560	pdufi(i,j,l)=
561	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
562	ENDDO
563	pdufi(iim,j,l)=
564	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
565	pdufi(iip1,j,l)=pdufi(1,j,l)
566	ENDDO
567
568	ENDDO
569
570
571	c 67. champ v:
572	c ------------
573
574	DO l=1,llm
575
576	DO j=2,jjm-1
577	ig0=1+(j-2)*iim
578	DO i=1,iim
579	pdvfi(i,j,l)=
580	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
581	ENDDO
582	pdvfi(iip1,j,l) = pdvfi(1,j,l)
583	ENDDO
584	ENDDO
585
586
587	c 68. champ v pres des poles:
588	c ---------------------------
589	c v = U * cos(long) + V * SIN(long)
590
591	DO l=1,llm
592
593	DO i=1,iim
594	pdvfi(i,1,l)=
595	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
596	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
597	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
598	pdvfi(i,1,l)=
599	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
600	pdvfi(i,jjm,l)=
601	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
602	ENDDO
603
604	pdvfi(iip1,1,l) = pdvfi(1,1,l)
605	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
606
607	ENDDO
608
609	c-----------------------------------------------------------------------
610
611	700 CONTINUE
612
613	firstcal = .FALSE.
614
615	RETURN
616	END

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

Download in other formats:

Original Format