Context Navigation

calfis.F @ 958

Last change on this file since 958 was 956, checked in by lmdzadmin, 17 years ago

Nettoyage du controle des parametres physiques. FH

Les parametres cycle_diurne, soil_model, new_oliq, ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad et iflag_con
sont maintenant geres par la physique uniquement.
ecritphy est elimine.
dimphy.F90 et clesphys.h ne sont plus utilises par le code dynamique.
Le test academique obtenu en compilant avec
makegcm -p nophys gcm
fonctionne. FH
IM

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 pucov covariant zonal velocity
68	c pvcov covariant meridional velocity
69	c pteta potential temperature
70	c pps surface pressure
71	c pmasse masse d'air dans chaque maille
72	c pts surface temperature (K)
73	c callrad clef d'appel au rayonnement
74	c
75	c Output :
76	c --------
77	c pdufi tendency for the natural zonal velocity (ms-1)
78	c pdvfi tendency for the natural meridional velocity
79	c pdhfi tendency for the potential temperature
80	c pdtsfi tendency for the surface temperature
81	c
82	c pdtrad radiative tendencies \ both input
83	c pfluxrad radiative fluxes / and output
84	c
85	c=======================================================================
86	c
87	c-----------------------------------------------------------------------
88	c
89	c 0. Declarations :
90	c ------------------
91
92	#include "dimensions.h"
93	#include "paramet.h"
94	#include "temps.h"
95	#include "advtrac.h"
96
97	INTEGER ngridmx,nq
98	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
99
100	#include "comconst.h"
101	#include "comvert.h"
102	#include "comgeom2.h"
103	#include "control.h"
104
105	c Arguments :
106	c -----------
107	LOGICAL lafin
108	REAL heure
109
110	REAL pvcov(iip1,jjm,llm)
111	REAL pucov(iip1,jjp1,llm)
112	REAL pteta(iip1,jjp1,llm)
113	REAL pmasse(iip1,jjp1,llm)
114	REAL pq(iip1,jjp1,llm,nqmx)
115	REAL pphis(iip1,jjp1)
116	REAL pphi(iip1,jjp1,llm)
117	c
118	REAL pdvcov(iip1,jjm,llm)
119	REAL pducov(iip1,jjp1,llm)
120	REAL pdteta(iip1,jjp1,llm)
121	REAL pdq(iip1,jjp1,llm,nqmx)
122	c
123	REAL pw(iip1,jjp1,llm)
124
125	REAL pps(iip1,jjp1)
126	REAL pp(iip1,jjp1,llmp1)
127	REAL ppk(iip1,jjp1,llm)
128	c
129	REAL pdvfi(iip1,jjm,llm)
130	REAL pdufi(iip1,jjp1,llm)
131	REAL pdhfi(iip1,jjp1,llm)
132	REAL pdqfi(iip1,jjp1,llm,nqmx)
133	REAL pdpsfi(iip1,jjp1)
134
135	INTEGER longcles
136	PARAMETER ( longcles = 20 )
137	REAL clesphy0( longcles )
138
139
140	c Local variables :
141	c -----------------
142
143	INTEGER i,j,l,ig0,ig,iq,iiq
144	REAL zpsrf(ngridmx)
145	REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
146	REAL zphi(ngridmx,llm),zphis(ngridmx)
147	c
148	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
149	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqmx)
150	c
151	REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
152	REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
153	c
154	REAL pvervel(ngridmx,llm)
155	c
156	REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
157	REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqmx)
158	REAL zdpsrf(ngridmx)
159	c
160	REAL zsin(iim),zcos(iim),z1(iim)
161	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
162	REAL unskap, pksurcp
163	c
164	cIM diagnostique PVteta, Amip2
165	INTEGER ntetaSTD
166	PARAMETER(ntetaSTD=3)
167	REAL rtetaSTD(ntetaSTD)
168	DATA rtetaSTD/350., 380., 405./
169	REAL PVteta(ngridmx,ntetaSTD)
170	c
171	#ifdef INCA
172	REAL flxw(iip1,jjp1,llm)
173	REAL flxwfi(ngridmx,llm)
174	#endif
175	c
176
177	REAL SSUM
178
179	LOGICAL firstcal, debut
180	DATA firstcal/.true./
181	SAVE firstcal,debut
182	REAL rdayvrai
183	c
184	c-----------------------------------------------------------------------
185	c
186	c 1. Initialisations :
187	c --------------------
188	c
189
190	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
191	PRINT*,'STOP dans calfis'
192	PRINT,'La dimension ngridmx doit etre egale a 2 + (jjm-1)iim'
193	PRINT*,' ngridmx jjm iim '
194	PRINT*,ngridmx,jjm,iim
195	STOP
196	ENDIF
197
198	c-----------------------------------------------------------------------
199	c latitude, longitude et aires des mailles pour la physique:
200	c ----------------------------------------------------------
201
202	c
203	IF ( firstcal ) THEN
204	debut = .TRUE.
205	ELSE
206	debut = .FALSE.
207	ENDIF
208
209	c
210	c
211	c-----------------------------------------------------------------------
212	c 40. transformation des variables dynamiques en variables physiques:
213	c ---------------------------------------------------------------
214
215	c 41. pressions au sol (en Pascals)
216	c ----------------------------------
217
218
219	zpsrf(1) = pps(1,1)
220
221	ig0 = 2
222	DO j = 2,jjm
223	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
224	ig0 = ig0+iim
225	ENDDO
226
227	zpsrf(ngridmx) = pps(1,jjp1)
228
229
230	c 42. pression intercouches :
231	c
232	c -----------------------------------------------------------------
233	c .... zplev definis aux (llm +1) interfaces des couches ....
234	c .... zplay definis aux ( llm ) milieux des couches ....
235	c -----------------------------------------------------------------
236
237	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
238	c
239	unskap = 1./ kappa
240	c
241	DO l = 1, llmp1
242	zplev( 1,l ) = pp(1,1,l)
243	ig0 = 2
244	DO j = 2, jjm
245	DO i =1, iim
246	zplev( ig0,l ) = pp(i,j,l)
247	ig0 = ig0 +1
248	ENDDO
249	ENDDO
250	zplev( ngridmx,l ) = pp(1,jjp1,l)
251	ENDDO
252	c
253	c
254
255	c 43. temperature naturelle (en K) et pressions milieux couches .
256	c ---------------------------------------------------------------
257
258	DO l=1,llm
259
260	pksurcp = ppk(1,1,l) / cpp
261	zplay(1,l) = preff * pksurcp ** unskap
262	ztfi(1,l) = pteta(1,1,l) * pksurcp
263	pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
264	ig0 = 2
265
266	DO j = 2, jjm
267	DO i = 1, iim
268	pksurcp = ppk(i,j,l) / cpp
269	zplay(ig0,l) = preff * pksurcp ** unskap
270	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
271	pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
272	ig0 = ig0 + 1
273	ENDDO
274	ENDDO
275
276	pksurcp = ppk(1,jjp1,l) / cpp
277	zplay(ig0,l) = preff * pksurcp ** unskap
278	ztfi (ig0,l) = pteta(1,jjp1,l) * pksurcp
279	pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
280
281	ENDDO
282
283	c 43.bis traceurs
284	c ---------------
285	c
286	DO iq=1,nq
287	iiq=niadv(iq)
288	DO l=1,llm
289	zqfi(1,l,iq) = pq(1,1,l,iiq)
290	ig0 = 2
291	DO j=2,jjm
292	DO i = 1, iim
293	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
294	ig0 = ig0 + 1
295	ENDDO
296	ENDDO
297	zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
298	ENDDO
299	ENDDO
300
301	c convergence dynamique pour les traceurs "EAU"
302
303	DO iq=1,2
304	DO l=1,llm
305	pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
306	ig0 = 2
307	DO j=2,jjm
308	DO i = 1, iim
309	pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
310	ig0 = ig0 + 1
311	ENDDO
312	ENDDO
313	pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
314	ENDDO
315	ENDDO
316
317
318	c Geopotentiel calcule par rapport a la surface locale:
319	c -----------------------------------------------------
320
321	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
322	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
323	DO l=1,llm
324	DO ig=1,ngridmx
325	zphi(ig,l)=zphi(ig,l)-zphis(ig)
326	ENDDO
327	ENDDO
328
329	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
330	c
331	DO l=1,llm
332	pvervel(1,l)=pw(1,1,l) * g /apoln
333	ig0=2
334	DO j=2,jjm
335	DO i = 1, iim
336	pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
337	ig0 = ig0 + 1
338	ENDDO
339	ENDDO
340	pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
341	ENDDO
342
343	c
344	c 45. champ u:
345	c ------------
346
347	DO 50 l=1,llm
348
349	DO 25 j=2,jjm
350	ig0 = 1+(j-2)*iim
351	zufi(ig0+1,l)= 0.5 *
352	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
353	pcvgu(ig0+1,l)= 0.5 *
354	$ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
355	DO 10 i=2,iim
356	zufi(ig0+i,l)= 0.5 *
357	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
358	pcvgu(ig0+i,l)= 0.5 *
359	$ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
360	10 CONTINUE
361	25 CONTINUE
362
363	50 CONTINUE
364
365
366	c 46.champ v:
367	c -----------
368
369	DO l=1,llm
370	DO j=2,jjm
371	ig0=1+(j-2)*iim
372	DO i=1,iim
373	zvfi(ig0+i,l)= 0.5 *
374	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
375	pcvgv(ig0+i,l)= 0.5 *
376	$ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
377	ENDDO
378	ENDDO
379	ENDDO
380
381
382	c 47. champs de vents aux pole nord
383	c ------------------------------
384	c U = 1 / pi * integrale [ v * cos(long) * d long ]
385	c V = 1 / pi * integrale [ v * sin(long) * d long ]
386
387	DO l=1,llm
388
389	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
390	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
391	DO i=2,iim
392	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
393	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
394	ENDDO
395
396	DO i=1,iim
397	zcos(i) = COS(rlonv(i))*z1(i)
398	zcosbis(i)= COS(rlonv(i))*z1bis(i)
399	zsin(i) = SIN(rlonv(i))*z1(i)
400	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
401	ENDDO
402
403	zufi(1,l) = SSUM(iim,zcos,1)/pi
404	pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
405	zvfi(1,l) = SSUM(iim,zsin,1)/pi
406	pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
407
408	ENDDO
409
410
411	c 48. champs de vents aux pole sud:
412	c ---------------------------------
413	c U = 1 / pi * integrale [ v * cos(long) * d long ]
414	c V = 1 / pi * integrale [ v * sin(long) * d long ]
415
416	DO l=1,llm
417
418	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
419	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
420	DO i=2,iim
421	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
422	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
423	ENDDO
424
425	DO i=1,iim
426	zcos(i) = COS(rlonv(i))*z1(i)
427	zcosbis(i) = COS(rlonv(i))*z1bis(i)
428	zsin(i) = SIN(rlonv(i))*z1(i)
429	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
430	ENDDO
431
432	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
433	pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
434	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
435	pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
436
437	ENDDO
438	c
439	cIM calcul PV a teta=350, 380, 405K
440	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
441	$ ztfi,zplay,zplev,
442	$ ntetaSTD,rtetaSTD,PVteta)
443	c
444	#ifdef INCA
445	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
446	#endif
447
448
449	c-----------------------------------------------------------------------
450	c Appel de la physique:
451	c ---------------------
452
453
454	CALL physiq (ngridmx,
455	. llm,
456	. nq,
457	. debut,
458	. lafin,
459	. rdayvrai,
460	. heure,
461	. dtphys,
462	. zplev,
463	. zplay,
464	. zphi,
465	. zphis,
466	. presnivs,
467	. clesphy0,
468	. zufi,
469	. zvfi,
470	. ztfi,
471	. zqfi,
472	. pvervel,
473	#ifdef INCA
474	. flxwfi,
475	#endif
476	. zdufi,
477	. zdvfi,
478	. zdtfi,
479	. zdqfi,
480	. zdpsrf,
481	cIM diagnostique PVteta, Amip2
482	. pducov,
483	. PVteta)
484
485	500 CONTINUE
486
487	c-----------------------------------------------------------------------
488	c transformation des tendances physiques en tendances dynamiques:
489	c ---------------------------------------------------------------
490
491	c tendance sur la pression :
492	c -----------------------------------
493
494	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
495	c
496	c 62. enthalpie potentielle
497	c ---------------------
498
499	DO l=1,llm
500
501	DO i=1,iip1
502	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
503	pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
504	ENDDO
505
506	DO j=2,jjm
507	ig0=1+(j-2)*iim
508	DO i=1,iim
509	pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
510	ENDDO
511	pdhfi(iip1,j,l) = pdhfi(1,j,l)
512	ENDDO
513
514	ENDDO
515
516
517	c 62. humidite specifique
518	c ---------------------
519
520	DO iq=1,nqmx
521	DO l=1,llm
522	DO i=1,iip1
523	pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
524	pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
525	ENDDO
526	DO j=2,jjm
527	ig0=1+(j-2)*iim
528	DO i=1,iim
529	pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
530	ENDDO
531	pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
532	ENDDO
533	ENDDO
534	ENDDO
535
536	c 63. traceurs
537	c ------------
538	C initialisation des tendances
539	pdqfi=0.
540	C
541	DO iq=1,nq
542	iiq=niadv(iq)
543	DO l=1,llm
544	DO i=1,iip1
545	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
546	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
547	ENDDO
548	DO j=2,jjm
549	ig0=1+(j-2)*iim
550	DO i=1,iim
551	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
552	ENDDO
553	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
554	ENDDO
555	ENDDO
556	ENDDO
557
558	c 65. champ u:
559	c ------------
560
561	DO l=1,llm
562
563	DO i=1,iip1
564	pdufi(i,1,l) = 0.
565	pdufi(i,jjp1,l) = 0.
566	ENDDO
567
568	DO j=2,jjm
569	ig0=1+(j-2)*iim
570	DO i=1,iim-1
571	pdufi(i,j,l)=
572	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
573	ENDDO
574	pdufi(iim,j,l)=
575	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
576	pdufi(iip1,j,l)=pdufi(1,j,l)
577	ENDDO
578
579	ENDDO
580
581
582	c 67. champ v:
583	c ------------
584
585	DO l=1,llm
586
587	DO j=2,jjm-1
588	ig0=1+(j-2)*iim
589	DO i=1,iim
590	pdvfi(i,j,l)=
591	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
592	ENDDO
593	pdvfi(iip1,j,l) = pdvfi(1,j,l)
594	ENDDO
595	ENDDO
596
597
598	c 68. champ v pres des poles:
599	c ---------------------------
600	c v = U * cos(long) + V * SIN(long)
601
602	DO l=1,llm
603
604	DO i=1,iim
605	pdvfi(i,1,l)=
606	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
607	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
608	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
609	pdvfi(i,1,l)=
610	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
611	pdvfi(i,jjm,l)=
612	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
613	ENDDO
614
615	pdvfi(iip1,1,l) = pdvfi(1,1,l)
616	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
617
618	ENDDO
619
620	c-----------------------------------------------------------------------
621
622	700 CONTINUE
623
624	firstcal = .FALSE.
625
626	RETURN
627	END

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

Download in other formats:

Original Format