Context Navigation

calfis.F @ 2087

Last change on this file since 2087 was 1676, checked in by aslmd, 12 years ago

Repare le probleme "un jour sans fin" des versions precedentes pour le modele generique. Pas d'impact sur modele terrestre.

NB: changements egalement faits sur dyn3dmem en prevision de son utilisation

Fix the bug "Groundhog Day" which was a problem for previous versions using generic physics. No impact whatsoever on terrestrial model.

NB: changes also done on dyn3dmem which is planned to be used soon

Author : AS equipe planeto

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 19.4 KB

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

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

Download in other formats:

Original Format