Context Navigation

calfis.F @ 3603

Last change on this file since 3603 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

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 1299 2010-01-20 14:27:21Z musat $
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	USE control_mod
34
35
36	IMPLICIT NONE
37	c=======================================================================
38	c
39	c 1. rearrangement des tableaux et transformation
40	c variables dynamiques > variables physiques
41	c 2. calcul des termes physiques
42	c 3. retransformation des tendances physiques en tendances dynamiques
43	c
44	c remarques:
45	c ----------
46	c
47	c - les vents sont donnes dans la physique par leurs composantes
48	c naturelles.
49	c - la variable thermodynamique de la physique est une variable
50	c intensive : T
51	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
52	c - les deux seules variables dependant de la geometrie necessaires
53	c pour la physique sont la latitude pour le rayonnement et
54	c l'aire de la maille quand on veut integrer une grandeur
55	c horizontalement.
56	c - les points de la physique sont les points scalaires de la
57	c la dynamique; numerotation:
58	c 1 pour le pole nord
59	c (jjm-1)*iim pour l'interieur du domaine
60	c ngridmx pour le pole sud
61	c ---> ngridmx=2+(jjm-1)*iim
62	c
63	c Input :
64	c -------
65	c pucov covariant zonal velocity
66	c pvcov covariant meridional velocity
67	c pteta potential temperature
68	c pps surface pressure
69	c pmasse masse d'air dans chaque maille
70	c pts surface temperature (K)
71	c callrad clef d'appel au rayonnement
72	c
73	c Output :
74	c --------
75	c pdufi tendency for the natural zonal velocity (ms-1)
76	c pdvfi tendency for the natural meridional velocity
77	c pdhfi tendency for the potential temperature
78	c pdtsfi tendency for the surface temperature
79	c
80	c pdtrad radiative tendencies \ both input
81	c pfluxrad radiative fluxes / and output
82	c
83	c=======================================================================
84	c
85	c-----------------------------------------------------------------------
86	c
87	c 0. Declarations :
88	c ------------------
89
90	#include "dimensions.h"
91	#include "paramet.h"
92	#include "temps.h"
93
94	INTEGER ngridmx
95	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
96
97	#include "comconst.h"
98	#include "comvert.h"
99	#include "comgeom2.h"
100
101	c Arguments :
102	c -----------
103	LOGICAL lafin
104
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	REAL, intent(in):: jD_cur, jH_cur
174	c
175	c-----------------------------------------------------------------------
176	c
177	c 1. Initialisations :
178	c --------------------
179	c
180	c
181	IF ( firstcal ) THEN
182	debut = .TRUE.
183	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
184	PRINT*,'STOP dans calfis'
185	PRINT,'La dimension ngridmx doit etre egale a 2 + (jjm-1)iim'
186	PRINT*,' ngridmx jjm iim '
187	PRINT*,ngridmx,jjm,iim
188	STOP
189	ENDIF
190	ELSE
191	debut = .FALSE.
192	ENDIF ! of IF (firstcal)
193
194	c
195	c
196	c-----------------------------------------------------------------------
197	c 40. transformation des variables dynamiques en variables physiques:
198	c ---------------------------------------------------------------
199
200	c 41. pressions au sol (en Pascals)
201	c ----------------------------------
202
203
204	zpsrf(1) = pps(1,1)
205
206	ig0 = 2
207	DO j = 2,jjm
208	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
209	ig0 = ig0+iim
210	ENDDO
211
212	zpsrf(ngridmx) = pps(1,jjp1)
213
214
215	c 42. pression intercouches :
216	c
217	c -----------------------------------------------------------------
218	c .... zplev definis aux (llm +1) interfaces des couches ....
219	c .... zplay definis aux ( llm ) milieux des couches ....
220	c -----------------------------------------------------------------
221
222	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
223	c
224	unskap = 1./ kappa
225	c
226	DO l = 1, llmp1
227	zplev( 1,l ) = pp(1,1,l)
228	ig0 = 2
229	DO j = 2, jjm
230	DO i =1, iim
231	zplev( ig0,l ) = pp(i,j,l)
232	ig0 = ig0 +1
233	ENDDO
234	ENDDO
235	zplev( ngridmx,l ) = pp(1,jjp1,l)
236	ENDDO
237	c
238	c
239
240	c 43. temperature naturelle (en K) et pressions milieux couches .
241	c ---------------------------------------------------------------
242
243	DO l=1,llm
244
245	pksurcp = ppk(1,1,l) / cpp
246	zplay(1,l) = preff * pksurcp ** unskap
247	ztfi(1,l) = pteta(1,1,l) * pksurcp
248	pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
249	ig0 = 2
250
251	DO j = 2, jjm
252	DO i = 1, iim
253	pksurcp = ppk(i,j,l) / cpp
254	zplay(ig0,l) = preff * pksurcp ** unskap
255	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
256	pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
257	ig0 = ig0 + 1
258	ENDDO
259	ENDDO
260
261	pksurcp = ppk(1,jjp1,l) / cpp
262	zplay(ig0,l) = preff * pksurcp ** unskap
263	ztfi (ig0,l) = pteta(1,jjp1,l) * pksurcp
264	pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
265
266	ENDDO
267
268	c 43.bis traceurs
269	c ---------------
270	c
271	DO iq=1,nqtot
272	iiq=niadv(iq)
273	DO l=1,llm
274	zqfi(1,l,iq) = pq(1,1,l,iiq)
275	ig0 = 2
276	DO j=2,jjm
277	DO i = 1, iim
278	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
279	ig0 = ig0 + 1
280	ENDDO
281	ENDDO
282	zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
283	ENDDO
284	ENDDO
285
286	c convergence dynamique pour les traceurs "EAU"
287	! Earth-specific treatment of first 2 tracers (water)
288	if (planet_type=="earth") then
289	DO iq=1,2
290	DO l=1,llm
291	pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
292	ig0 = 2
293	DO j=2,jjm
294	DO i = 1, iim
295	pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
296	ig0 = ig0 + 1
297	ENDDO
298	ENDDO
299	pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
300	ENDDO
301	ENDDO
302	endif ! of if (planet_type=="earth")
303
304
305	c Geopotentiel calcule par rapport a la surface locale:
306	c -----------------------------------------------------
307
308	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
309	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
310	DO l=1,llm
311	DO ig=1,ngridmx
312	zphi(ig,l)=zphi(ig,l)-zphis(ig)
313	ENDDO
314	ENDDO
315
316	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
317	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
318	c de masse est calclue dans advtrac.F
319	c DO l=1,llm
320	c pvervel(1,l)=pw(1,1,l) * g /apoln
321	c ig0=2
322	c DO j=2,jjm
323	c DO i = 1, iim
324	c pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
325	c ig0 = ig0 + 1
326	c ENDDO
327	c ENDDO
328	c pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
329	c ENDDO
330
331	c
332	c 45. champ u:
333	c ------------
334
335	DO 50 l=1,llm
336
337	DO 25 j=2,jjm
338	ig0 = 1+(j-2)*iim
339	zufi(ig0+1,l)= 0.5 *
340	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
341	pcvgu(ig0+1,l)= 0.5 *
342	$ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
343	DO 10 i=2,iim
344	zufi(ig0+i,l)= 0.5 *
345	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
346	pcvgu(ig0+i,l)= 0.5 *
347	$ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
348	10 CONTINUE
349	25 CONTINUE
350
351	50 CONTINUE
352
353
354	c 46.champ v:
355	c -----------
356
357	DO l=1,llm
358	DO j=2,jjm
359	ig0=1+(j-2)*iim
360	DO i=1,iim
361	zvfi(ig0+i,l)= 0.5 *
362	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
363	pcvgv(ig0+i,l)= 0.5 *
364	$ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
365	ENDDO
366	ENDDO
367	ENDDO
368
369
370	c 47. champs de vents aux pole nord
371	c ------------------------------
372	c U = 1 / pi * integrale [ v * cos(long) * d long ]
373	c V = 1 / pi * integrale [ v * sin(long) * d long ]
374
375	DO l=1,llm
376
377	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
378	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
379	DO i=2,iim
380	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
381	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
382	ENDDO
383
384	DO i=1,iim
385	zcos(i) = COS(rlonv(i))*z1(i)
386	zcosbis(i)= COS(rlonv(i))*z1bis(i)
387	zsin(i) = SIN(rlonv(i))*z1(i)
388	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
389	ENDDO
390
391	zufi(1,l) = SSUM(iim,zcos,1)/pi
392	pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
393	zvfi(1,l) = SSUM(iim,zsin,1)/pi
394	pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
395
396	ENDDO
397
398
399	c 48. champs de vents aux pole sud:
400	c ---------------------------------
401	c U = 1 / pi * integrale [ v * cos(long) * d long ]
402	c V = 1 / pi * integrale [ v * sin(long) * d long ]
403
404	DO l=1,llm
405
406	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
407	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
408	DO i=2,iim
409	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
410	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
411	ENDDO
412
413	DO i=1,iim
414	zcos(i) = COS(rlonv(i))*z1(i)
415	zcosbis(i) = COS(rlonv(i))*z1bis(i)
416	zsin(i) = SIN(rlonv(i))*z1(i)
417	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
418	ENDDO
419
420	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
421	pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
422	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
423	pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
424
425	ENDDO
426	c
427	if (planet_type=="earth") then
428	#ifdef CPP_EARTH
429	cIM calcul PV a teta=350, 380, 405K
430	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
431	$ ztfi,zplay,zplev,
432	$ ntetaSTD,rtetaSTD,PVteta)
433	#endif
434	endif
435	c
436	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
437	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
438
439	c-----------------------------------------------------------------------
440	c Appel de la physique:
441	c ---------------------
442
443
444	if (planet_type=="earth") then
445	#ifdef CPP_EARTH
446	CALL physiq (ngridmx,
447	. llm,
448	. debut,
449	. lafin,
450	. jD_cur,
451	. jH_cur,
452	. dtphys,
453	. zplev,
454	. zplay,
455	. zphi,
456	. zphis,
457	. presnivs,
458	. clesphy0,
459	. zufi,
460	. zvfi,
461	. ztfi,
462	. zqfi,
463	. flxwfi,
464	. zdufi,
465	. zdvfi,
466	. zdtfi,
467	. zdqfi,
468	. zdpsrf,
469	cIM diagnostique PVteta, Amip2
470	. pducov,
471	. PVteta)
472	#endif
473	endif !of if (planet_type=="earth")
474
475	500 CONTINUE
476
477	c-----------------------------------------------------------------------
478	c transformation des tendances physiques en tendances dynamiques:
479	c ---------------------------------------------------------------
480
481	c tendance sur la pression :
482	c -----------------------------------
483
484	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
485	c
486	c 62. enthalpie potentielle
487	c ---------------------
488
489	DO l=1,llm
490
491	DO i=1,iip1
492	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
493	pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
494	ENDDO
495
496	DO j=2,jjm
497	ig0=1+(j-2)*iim
498	DO i=1,iim
499	pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
500	ENDDO
501	pdhfi(iip1,j,l) = pdhfi(1,j,l)
502	ENDDO
503
504	ENDDO
505
506
507	c 62. humidite specifique
508	c ---------------------
509	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
510	! DO iq=1,nqtot
511	! DO l=1,llm
512	! DO i=1,iip1
513	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
514	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
515	! ENDDO
516	! DO j=2,jjm
517	! ig0=1+(j-2)*iim
518	! DO i=1,iim
519	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
520	! ENDDO
521	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
522	! ENDDO
523	! ENDDO
524	! ENDDO
525
526	c 63. traceurs
527	c ------------
528	C initialisation des tendances
529	pdqfi(:,:,:,:)=0.
530	C
531	DO iq=1,nqtot
532	iiq=niadv(iq)
533	DO l=1,llm
534	DO i=1,iip1
535	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
536	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
537	ENDDO
538	DO j=2,jjm
539	ig0=1+(j-2)*iim
540	DO i=1,iim
541	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
542	ENDDO
543	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
544	ENDDO
545	ENDDO
546	ENDDO
547
548	c 65. champ u:
549	c ------------
550
551	DO l=1,llm
552
553	DO i=1,iip1
554	pdufi(i,1,l) = 0.
555	pdufi(i,jjp1,l) = 0.
556	ENDDO
557
558	DO j=2,jjm
559	ig0=1+(j-2)*iim
560	DO i=1,iim-1
561	pdufi(i,j,l)=
562	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
563	ENDDO
564	pdufi(iim,j,l)=
565	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
566	pdufi(iip1,j,l)=pdufi(1,j,l)
567	ENDDO
568
569	ENDDO
570
571
572	c 67. champ v:
573	c ------------
574
575	DO l=1,llm
576
577	DO j=2,jjm-1
578	ig0=1+(j-2)*iim
579	DO i=1,iim
580	pdvfi(i,j,l)=
581	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
582	ENDDO
583	pdvfi(iip1,j,l) = pdvfi(1,j,l)
584	ENDDO
585	ENDDO
586
587
588	c 68. champ v pres des poles:
589	c ---------------------------
590	c v = U * cos(long) + V * SIN(long)
591
592	DO l=1,llm
593
594	DO i=1,iim
595	pdvfi(i,1,l)=
596	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
597	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
598	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
599	pdvfi(i,1,l)=
600	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
601	pdvfi(i,jjm,l)=
602	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
603	ENDDO
604
605	pdvfi(iip1,1,l) = pdvfi(1,1,l)
606	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
607
608	ENDDO
609
610	c-----------------------------------------------------------------------
611
612	700 CONTINUE
613
614	firstcal = .FALSE.
615
616	RETURN
617	END

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

Download in other formats:

Original Format