Context Navigation

calfis.F @ 1972

Last change on this file since 1972 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
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 1907 2013-11-26 13:10:46Z 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