Context Navigation

calfis.F @ 86

Last change on this file since 86 was 37, checked in by emillour, 14 years ago

Remise en route chantier compilation -- Ehouarn

Modifs et corrections pour pouvoir compiler le gcm (en séentiel, avec

makelmdz_fcm pour l'instant):

ajout de fichiers 'arch' pour linux-64 (pour Bellonzi, avec ioipsl et en r8)
modification de makelmdz_fcm, ajout de la cléPP_PHYS si on compile avec une physique
correction de quelques typos/bugs réléàa compilation:
infotrac.F90 : supression des appels àlnblnk' (remplacépar len_trim)
bilan_dyn.F : déaration des variables znom3,znom3l,zunites3, planet_type
cpdet.F : "use control_mod, ONLY: planet_type" mis aux bons endroits

(idem sur cpdet.F dans dyn3dpar)

leapfrog.F : declaration de ztetaec(), dtec, cpdet , itau_w, duspg()
diagedyn.F : correction typo; attention dans diagedyn.F il y a du

include "../phylmd/YOMCST.h"
Ca va poser problè qd on change de physique....

Avec ces modifs, la compilation marche sans physque, avec ou sans ioipsl et avec la physique terrestre phylmd.

File size: 20.2 KB

Line
1	!
2	! $Id: calfis.F 1407 2010-07-07 10:31:52Z fairhead $
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 (K/s)
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	#include "iniprint.h"
101
102	c Arguments :
103	c -----------
104	LOGICAL lafin
105
106
107	REAL pvcov(iip1,jjm,llm)
108	REAL pucov(iip1,jjp1,llm)
109	REAL pteta(iip1,jjp1,llm)
110	REAL pmasse(iip1,jjp1,llm)
111	REAL pq(iip1,jjp1,llm,nqtot)
112	REAL pphis(iip1,jjp1)
113	REAL pphi(iip1,jjp1,llm)
114	c
115	REAL pdvcov(iip1,jjm,llm)
116	REAL pducov(iip1,jjp1,llm)
117	REAL pdteta(iip1,jjp1,llm)
118	! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
119	! qui lui meme ne sert a rien dans la routine telle qu'elle est
120	! ecrite, et que j'ai donc commente....
121	REAL pdq(iip1,jjp1,llm,nqtot)
122	c
123	REAL pps(iip1,jjp1)
124	REAL pp(iip1,jjp1,llmp1)
125	REAL ppk(iip1,jjp1,llm)
126	c
127	c TENDENCIES in */s
128	REAL pdvfi(iip1,jjm,llm)
129	REAL pdufi(iip1,jjp1,llm)
130	REAL pdhfi(iip1,jjp1,llm)
131	REAL pdqfi(iip1,jjp1,llm,nqtot)
132	REAL pdpsfi(iip1,jjp1)
133
134	INTEGER longcles
135	PARAMETER ( longcles = 20 )
136	REAL clesphy0( longcles )
137
138
139	c Local variables :
140	c -----------------
141
142	INTEGER i,j,l,ig0,ig,iq,iiq
143	REAL zpsrf(ngridmx)
144	REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
145	REAL zphi(ngridmx,llm),zphis(ngridmx)
146	c
147	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
148	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
149	! ADAPTATION GCM POUR CP(T)
150	REAL zteta(ngridmx,llm)
151	REAL zpk(ngridmx,llm)
152	c
153	! RQ SL 13/10/10:
154	! Ces calculs ne servent pas.
155	! Si necessaire, decommenter ces variables et les calculs...
156	! REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
157	! REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
158	c
159	REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
160	REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
161	REAL zdpsrf(ngridmx)
162	c
163	REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
164	REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
165	REAL jH_cur_split,zdt_split
166	LOGICAL debut_split,lafin_split
167	INTEGER isplit
168
169	REAL zsin(iim),zcos(iim),z1(iim)
170	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
171	REAL unskap, pksurcp
172	c
173	cIM diagnostique PVteta, Amip2
174	INTEGER ntetaSTD
175	PARAMETER(ntetaSTD=3)
176	REAL rtetaSTD(ntetaSTD)
177	DATA rtetaSTD/350., 380., 405./
178	REAL PVteta(ngridmx,ntetaSTD)
179	c
180	REAL flxw(iip1,jjp1,llm) ! Flux de masse verticale sur la grille dynamique
181	REAL flxwfi(ngridmx,llm) ! Flux de masse verticale sur la grille physiq
182	c
183
184	REAL SSUM
185
186	LOGICAL firstcal, debut
187	DATA firstcal/.true./
188	SAVE firstcal,debut
189	! REAL rdayvrai
190	REAL, intent(in):: jD_cur, jH_cur
191	c
192	c-----------------------------------------------------------------------
193	c
194	c 1. Initialisations :
195	c --------------------
196	c
197	c
198	IF ( firstcal ) THEN
199	debut = .TRUE.
200	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
201	write(lunout,*) 'STOP dans calfis'
202	write(lunout,*)
203	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
204	write(lunout,*) ' ngridmx jjm iim '
205	write(lunout,*) ngridmx,jjm,iim
206	STOP
207	ENDIF
208	ELSE
209	debut = .FALSE.
210	ENDIF ! of IF (firstcal)
211
212	c
213	c
214	c-----------------------------------------------------------------------
215	c 40. transformation des variables dynamiques en variables physiques:
216	c ---------------------------------------------------------------
217
218	c 41. pressions au sol (en Pascals)
219	c ----------------------------------
220
221
222	zpsrf(1) = pps(1,1)
223
224	ig0 = 2
225	DO j = 2,jjm
226	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
227	ig0 = ig0+iim
228	ENDDO
229
230	zpsrf(ngridmx) = pps(1,jjp1)
231
232
233	c 42. pression intercouches et fonction d'Exner:
234	c
235	c -----------------------------------------------------------------
236	c .... zplev definis aux (llm +1) interfaces des couches ....
237	c .... zplay definis aux ( llm ) milieux des couches ....
238	c -----------------------------------------------------------------
239
240	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
241	c
242	unskap = 1./ kappa
243	c
244	! ADAPTATION GCM POUR CP(T)
245	DO l = 1, llm
246	zpk( 1,l ) = ppk(1,1,l)
247	zteta( 1,l ) = pteta(1,1,l)
248	zplev( 1,l ) = pp(1,1,l)
249	ig0 = 2
250	DO j = 2, jjm
251	DO i =1, iim
252	zpk( ig0,l ) = ppk(i,j,l)
253	zteta( ig0,l ) = pteta(i,j,l)
254	zplev( ig0,l ) = pp(i,j,l)
255	ig0 = ig0 +1
256	ENDDO
257	ENDDO
258	zpk( ngridmx,l ) = ppk(1,jjp1,l)
259	zteta( ngridmx,l ) = pteta(1,jjp1,l)
260	zplev( ngridmx,l ) = pp(1,jjp1,l)
261	ENDDO
262	zplev( 1,llmp1 ) = pp(1,1,llmp1)
263	ig0 = 2
264	DO j = 2, jjm
265	DO i =1, iim
266	zplev( ig0,llmp1 ) = pp(i,j,llmp1)
267	ig0 = ig0 +1
268	ENDDO
269	ENDDO
270	zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
271	c
272	c
273
274	c 43. temperature naturelle (en K) et pressions milieux couches .
275	c ---------------------------------------------------------------
276
277	! ADAPTATION GCM POUR CP(T)
278	call tpot2t(ngridmx*llm,zteta,ztfi,zpk)
279
280	DO l=1,llm
281
282	pksurcp = ppk(1,1,l) / cpp
283	zplay(1,l) = preff * pksurcp ** unskap
284	! pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
285	ig0 = 2
286
287	DO j = 2, jjm
288	DO i = 1, iim
289	pksurcp = ppk(i,j,l) / cpp
290	zplay(ig0,l) = preff * pksurcp ** unskap
291	! pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
292	ig0 = ig0 + 1
293	ENDDO
294	ENDDO
295
296	pksurcp = ppk(1,jjp1,l) / cpp
297	zplay(ig0,l) = preff * pksurcp ** unskap
298	! pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
299
300	ENDDO
301
302	c 43.bis traceurs (tous intensifs)
303	c ---------------
304	c
305	DO iq=1,nqtot
306	iiq=niadv(iq)
307	DO l=1,llm
308	zqfi(1,l,iq) = pq(1,1,l,iiq)
309	ig0 = 2
310	DO j=2,jjm
311	DO i = 1, iim
312	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
313	ig0 = ig0 + 1
314	ENDDO
315	ENDDO
316	zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
317	ENDDO
318	ENDDO ! boucle sur traceurs
319
320	!-----------------
321	! RQ SL 13/10/10:
322	! Ces calculs ne servent pas.
323	! Si necessaire, decommenter ces variables et les calculs...
324	!
325	! convergence dynamique pour les traceurs "EAU"
326	! Earth-specific treatment of first 2 tracers (water)
327	! if (planet_type=="earth") then
328	! DO iq=1,2
329	! DO l=1,llm
330	! pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
331	! ig0 = 2
332	! DO j=2,jjm
333	! DO i = 1, iim
334	! pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
335	! ig0 = ig0 + 1
336	! ENDDO
337	! ENDDO
338	! pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
339	! ENDDO
340	! ENDDO
341	! endif ! of if (planet_type=="earth")
342	!----------------
343
344	c Geopotentiel calcule par rapport a la surface locale:
345	c -----------------------------------------------------
346
347	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
348	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
349	DO l=1,llm
350	DO ig=1,ngridmx
351	zphi(ig,l)=zphi(ig,l)-zphis(ig)
352	ENDDO
353	ENDDO
354
355	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
356	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
357	c de masse est calclue dans advtrac.F
358	c DO l=1,llm
359	c pvervel(1,l)=pw(1,1,l) * g /apoln
360	c ig0=2
361	c DO j=2,jjm
362	c DO i = 1, iim
363	c pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
364	c ig0 = ig0 + 1
365	c ENDDO
366	c ENDDO
367	c pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
368	c ENDDO
369
370	c
371	c 45. champ u:
372	c ------------
373
374	DO 50 l=1,llm
375
376	DO 25 j=2,jjm
377	ig0 = 1+(j-2)*iim
378	zufi(ig0+1,l)= 0.5 *
379	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
380	! pcvgu(ig0+1,l)= 0.5 *
381	! $ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
382	DO 10 i=2,iim
383	zufi(ig0+i,l)= 0.5 *
384	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
385	! pcvgu(ig0+i,l)= 0.5 *
386	! $ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
387	10 CONTINUE
388	25 CONTINUE
389
390	50 CONTINUE
391
392
393	c 46.champ v:
394	c -----------
395
396	DO l=1,llm
397	DO j=2,jjm
398	ig0=1+(j-2)*iim
399	DO i=1,iim
400	zvfi(ig0+i,l)= 0.5 *
401	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
402	c pcvgv(ig0+i,l)= 0.5 *
403	c $ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
404	ENDDO
405	ENDDO
406	ENDDO
407
408
409	c 47. champs de vents aux pole nord
410	c ------------------------------
411	c U = 1 / pi * integrale [ v * cos(long) * d long ]
412	c V = 1 / pi * integrale [ v * sin(long) * d long ]
413
414	DO l=1,llm
415
416	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
417	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
418	DO i=2,iim
419	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
420	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
421	ENDDO
422
423	DO i=1,iim
424	zcos(i) = COS(rlonv(i))*z1(i)
425	zcosbis(i)= COS(rlonv(i))*z1bis(i)
426	zsin(i) = SIN(rlonv(i))*z1(i)
427	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
428	ENDDO
429
430	zufi(1,l) = SSUM(iim,zcos,1)/pi
431	! pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
432	zvfi(1,l) = SSUM(iim,zsin,1)/pi
433	! pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
434
435	ENDDO
436
437
438	c 48. champs de vents aux pole sud:
439	c ---------------------------------
440	c U = 1 / pi * integrale [ v * cos(long) * d long ]
441	c V = 1 / pi * integrale [ v * sin(long) * d long ]
442
443	DO l=1,llm
444
445	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
446	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
447	DO i=2,iim
448	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
449	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
450	ENDDO
451
452	DO i=1,iim
453	zcos(i) = COS(rlonv(i))*z1(i)
454	zcosbis(i) = COS(rlonv(i))*z1bis(i)
455	zsin(i) = SIN(rlonv(i))*z1(i)
456	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
457	ENDDO
458
459	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
460	! pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
461	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
462	! pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
463
464	ENDDO
465	c
466	if (planet_type=="earth") then
467	#ifdef CPP_EARTH
468	cIM calcul PV a teta=350, 380, 405K
469	CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
470	$ ztfi,zplay,zplev,
471	$ ntetaSTD,rtetaSTD,PVteta)
472	#endif
473	endif
474	c
475	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
476	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
477
478	c-----------------------------------------------------------------------
479	c Appel de la physique:
480	c ---------------------
481
482	! Appel de la physique: pose probleme quand on tourne
483	! SANS physique, car physiq.F est dans le repertoire phy[]...
484	! Il faut une cle CPP_PHYS
485
486	! Le fait que les arguments de physiq soient differents selon les planetes
487	! ne pose pas de probleme a priori.
488
489	#ifdef CPP_PHYS
490
491	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
492	zdt_split=dtphys/nsplit_phys
493	zdufic(:,:)=0.
494	zdvfic(:,:)=0.
495	zdtfic(:,:)=0.
496	zdqfic(:,:,:)=0.
497
498	do isplit=1,nsplit_phys
499
500	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
501	debut_split=debut.and.isplit==1
502	lafin_split=lafin.and.isplit==nsplit_phys
503
504	if (planet_type.eq."earth") then
505	CALL physiq (ngridmx,
506	. llm,
507	. debut_split,
508	. lafin_split,
509	. jD_cur,
510	. jH_cur_split,
511	. zdt_split,
512	. zplev,
513	. zplay,
514	. zphi,
515	. zphis,
516	. presnivs,
517	. clesphy0,
518	. zufi,
519	. zvfi,
520	. ztfi,
521	. zqfi,
522	. flxwfi,
523	. zdufi,
524	. zdvfi,
525	. zdtfi,
526	. zdqfi,
527	. zdpsrf,
528	cIM diagnostique PVteta, Amip2
529	. pducov,
530	. PVteta)
531
532	else ! a moduler pour Mars !!
533
534	CALL physiq (ngridmx,
535	. llm,
536	. nqtot,
537	. debut_split,
538	. lafin_split,
539	. jD_cur,
540	. jH_cur_split,
541	. zdt_split,
542	. zplev,
543	. zplay,
544	. zpk,
545	. zphi,
546	. zphis,
547	. presnivs,
548	. clesphy0,
549	. zufi,
550	. zvfi,
551	. ztfi,
552	. zqfi,
553	. flxwfi,
554	. zdufi,
555	. zdvfi,
556	. zdtfi,
557	. zdqfi,
558	. zdpsrf)
559
560	endif ! planet_type
561
562	zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
563	zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
564	ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
565	zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
566
567	zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
568	zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
569	zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
570	zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
571
572	enddo
573	zdufi(:,:)=zdufic(:,:)/nsplit_phys
574	zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
575	zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
576	zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
577
578	#endif
579	! #endif of #ifdef CPP_PHYS
580
581	500 CONTINUE
582
583	c-----------------------------------------------------------------------
584	c transformation des tendances physiques en tendances dynamiques:
585	c ---------------------------------------------------------------
586
587	c tendance sur la pression :
588	c -----------------------------------
589
590	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
591	c
592	c 62. enthalpie potentielle
593	c ---------------------
594
595	! ADAPTATION GCM POUR CP(T)
596	ztfi=ztfi+zdtfi*dtphys
597	call t2tpot(ngridmx*llm,ztfi,zteta,zpk)
598
599	DO i=1,iip1
600	pdhfi(i,1,:) = (zteta(1,:) - pteta(i,1,:))/dtphys
601	pdhfi(i,jjp1,:) = (zteta(ngridmx,:) - pteta(i,jjp1,:))/dtphys
602	ENDDO
603
604	DO j=2,jjm
605	ig0=1+(j-2)*iim
606	DO i=1,iim
607	pdhfi(i,j,:) = (zteta(ig0+i,:) - pteta(i,j,:))/dtphys
608	ENDDO
609	pdhfi(iip1,j,:) = pdhfi(1,j,:)
610	ENDDO
611
612
613	c 62. humidite specifique
614	c ---------------------
615	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
616	! DO iq=1,nqtot
617	! DO l=1,llm
618	! DO i=1,iip1
619	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
620	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
621	! ENDDO
622	! DO j=2,jjm
623	! ig0=1+(j-2)*iim
624	! DO i=1,iim
625	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
626	! ENDDO
627	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
628	! ENDDO
629	! ENDDO
630	! ENDDO
631
632	c 63. traceurs (tous en intensifs)
633	c ------------
634	C initialisation des tendances
635	pdqfi(:,:,:,:)=0.
636	C
637	DO iq=1,nqtot
638	iiq=niadv(iq)
639	DO l=1,llm
640	DO i=1,iip1
641	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
642	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
643	ENDDO
644	DO j=2,jjm
645	ig0=1+(j-2)*iim
646	DO i=1,iim
647	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
648	ENDDO
649	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
650	ENDDO
651	ENDDO
652	ENDDO
653
654	c 65. champ u:
655	c ------------
656
657	DO l=1,llm
658
659	DO i=1,iip1
660	pdufi(i,1,l) = 0.
661	pdufi(i,jjp1,l) = 0.
662	ENDDO
663
664	DO j=2,jjm
665	ig0=1+(j-2)*iim
666	DO i=1,iim-1
667	pdufi(i,j,l)=
668	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
669	ENDDO
670	pdufi(iim,j,l)=
671	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
672	pdufi(iip1,j,l)=pdufi(1,j,l)
673	ENDDO
674
675	ENDDO
676
677
678	c 67. champ v:
679	c ------------
680
681	DO l=1,llm
682
683	DO j=2,jjm-1
684	ig0=1+(j-2)*iim
685	DO i=1,iim
686	pdvfi(i,j,l)=
687	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
688	ENDDO
689	pdvfi(iip1,j,l) = pdvfi(1,j,l)
690	ENDDO
691	ENDDO
692
693
694	c 68. champ v pres des poles:
695	c ---------------------------
696	c v = U * cos(long) + V * SIN(long)
697
698	DO l=1,llm
699
700	DO i=1,iim
701	pdvfi(i,1,l)=
702	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
703	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
704	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
705	pdvfi(i,1,l)=
706	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
707	pdvfi(i,jjm,l)=
708	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
709	ENDDO
710
711	pdvfi(iip1,1,l) = pdvfi(1,1,l)
712	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
713
714	ENDDO
715
716	c-----------------------------------------------------------------------
717
718	700 CONTINUE
719
720	firstcal = .FALSE.
721
722	RETURN
723	END

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

Download in other formats:

Original Format