Context Navigation

calfis.F @ 1625

Last change on this file since 1625 was 1617, checked in by Ehouarn Millour, 13 years ago
Minor correction: with previous commits, compilation without physics was no longer possible (because PVtheta calls routine tetalevel, which is in the physics directory). Note also that the levels (hard coded in calfis and calfis_p) given as inputs to PVtheta are Earth-specific. This needs to be cleaned and improved... EM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 18.2 KB

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

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

Download in other formats:

Original Format