Context Navigation

calfis.F @ 5425

Last change on this file since 5425 was 1279, checked in by Laurent Fairhead, 15 years ago
Merged LMDZ4-dev branch changes r1241:1278 into the trunk Running trunk and LMDZ4-dev in LMDZOR configuration on local machine (sequential) and SX8 (4-proc) yields identical results (restart and restartphy are identical binarily) Log history from r1241 to r1278 is available by switching to source:LMDZ4/branches/LMDZ4-dev-20091210
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 1279 2009-12-10 09:02:56Z jyg $
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