Context Navigation

calfis.F @ 2242

Last change on this file since 2242 was 2239, checked in by Ehouarn Millour, 9 years ago

Reorganizing physics/dynamics interface:

what is related to dynamics-physics interface is now in a seperate directory: dynlmdz_phy* for physics in phy*
1d model and related dependencies (including a couple from "dynamics", set up as symbolic links) is now in subdirectory "dyn1d" of phy*.
"bibio" directory is now "misc" and should only contain autonomous utilities.
"cosp" is now a subdirectory of phylmd.

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.8 KB

Line
1	!
2	! $Id: calfis.F 2239 2015-03-23 07:27:30Z emillour $
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	$ pdufi,
24	$ pdvfi,
25	$ pdhfi,
26	$ pdqfi,
27	$ pdpsfi)
28	c
29	c Auteur : P. Le Van, F. Hourdin
30	c .........
31	USE infotrac, ONLY: nqtot, niadv, tname
32	USE control_mod, ONLY: planet_type, nsplit_phys
33
34
35	IMPLICIT NONE
36	c=======================================================================
37	c
38	c 1. rearrangement des tableaux et transformation
39	c variables dynamiques > variables physiques
40	c 2. calcul des termes physiques
41	c 3. retransformation des tendances physiques en tendances dynamiques
42	c
43	c remarques:
44	c ----------
45	c
46	c - les vents sont donnes dans la physique par leurs composantes
47	c naturelles.
48	c - la variable thermodynamique de la physique est une variable
49	c intensive : T
50	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
51	c - les deux seules variables dependant de la geometrie necessaires
52	c pour la physique sont la latitude pour le rayonnement et
53	c l'aire de la maille quand on veut integrer une grandeur
54	c horizontalement.
55	c - les points de la physique sont les points scalaires de la
56	c la dynamique; numerotation:
57	c 1 pour le pole nord
58	c (jjm-1)*iim pour l'interieur du domaine
59	c ngridmx pour le pole sud
60	c ---> ngridmx=2+(jjm-1)*iim
61	c
62	c Input :
63	c -------
64	c pucov covariant zonal velocity
65	c pvcov covariant meridional velocity
66	c pteta potential temperature
67	c pps surface pressure
68	c pmasse masse d'air dans chaque maille
69	c pts surface temperature (K)
70	c callrad clef d'appel au rayonnement
71	c
72	c Output :
73	c --------
74	c pdufi tendency for the natural zonal velocity (ms-1)
75	c pdvfi tendency for the natural meridional velocity
76	c pdhfi tendency for the potential temperature
77	c pdtsfi tendency for the surface temperature
78	c
79	c pdtrad radiative tendencies \ both input
80	c pfluxrad radiative fluxes / and output
81	c
82	c=======================================================================
83	c
84	c-----------------------------------------------------------------------
85	c
86	c 0. Declarations :
87	c ------------------
88
89	#include "dimensions.h"
90	#include "paramet.h"
91	#include "temps.h"
92
93	INTEGER ngridmx
94	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
95
96	#include "comconst.h"
97	#include "comvert.h"
98	#include "comgeom2.h"
99	#include "iniprint.h"
100
101	c Arguments :
102	c -----------
103	LOGICAL,INTENT(IN) :: lafin ! .true. for the very last call to physics
104	REAL,INTENT(IN):: jD_cur, jH_cur
105	REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
106	REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
107	REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
108	REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
109	REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
110	REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
111	REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
112
113	REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
114	REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
115	REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
116	! NB: pdteta is used only to compute pcvgt which is in fact not used...
117	REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
118	! NB: pdq is only used to compute pcvgq which is in fact not used...
119
120	REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
121	REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
122	REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
123	REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on dynamics grid
124
125	! tendencies (in */s) from the physics
126	REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
127	REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
128	REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
129	REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
130	REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
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 zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
152	REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
153	REAL jH_cur_split,zdt_split
154	LOGICAL debut_split,lafin_split
155	INTEGER isplit
156
157	REAL zsin(iim),zcos(iim),z1(iim)
158	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
159	REAL unskap, pksurcp
160	c
161	REAL flxwfi(ngridmx,llm) ! Flux de masse verticale sur la grille physiq
162	c
163
164	REAL SSUM
165
166	LOGICAL,SAVE :: firstcal=.true., debut=.true.
167	! REAL rdayvrai
168
169	LOGICAL tracerdyn
170
171	c
172	c-----------------------------------------------------------------------
173	c
174	c 1. Initialisations :
175	c --------------------
176	c
177	c
178	IF ( firstcal ) THEN
179	debut = .TRUE.
180	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
181	write(lunout,*) 'STOP dans calfis'
182	write(lunout,*)
183	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
184	write(lunout,*) ' ngridmx jjm iim '
185	write(lunout,*) ngridmx,jjm,iim
186	STOP
187	ENDIF
188	ELSE
189	debut = .FALSE.
190	ENDIF ! of IF (firstcal)
191
192	c
193	c
194	c-----------------------------------------------------------------------
195	c 40. transformation des variables dynamiques en variables physiques:
196	c ---------------------------------------------------------------
197
198	c 41. pressions au sol (en Pascals)
199	c ----------------------------------
200
201
202	zpsrf(1) = pps(1,1)
203
204	ig0 = 2
205	DO j = 2,jjm
206	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
207	ig0 = ig0+iim
208	ENDDO
209
210	zpsrf(ngridmx) = pps(1,jjp1)
211
212
213	c 42. pression intercouches :
214	c
215	c -----------------------------------------------------------------
216	c .... zplev definis aux (llm +1) interfaces des couches ....
217	c .... zplay definis aux ( llm ) milieux des couches ....
218	c -----------------------------------------------------------------
219
220	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
221	c
222	unskap = 1./ kappa
223	c
224	DO l = 1, llmp1
225	zplev( 1,l ) = pp(1,1,l)
226	ig0 = 2
227	DO j = 2, jjm
228	DO i =1, iim
229	zplev( ig0,l ) = pp(i,j,l)
230	ig0 = ig0 +1
231	ENDDO
232	ENDDO
233	zplev( ngridmx,l ) = pp(1,jjp1,l)
234	ENDDO
235	c
236	c
237
238	c 43. temperature naturelle (en K) et pressions milieux couches .
239	c ---------------------------------------------------------------
240
241	DO l=1,llm
242
243	pksurcp = ppk(1,1,l) / cpp
244	zplay(1,l) = preff * pksurcp ** unskap
245	ztfi(1,l) = pteta(1,1,l) * pksurcp
246	pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
247	ig0 = 2
248
249	DO j = 2, jjm
250	DO i = 1, iim
251	pksurcp = ppk(i,j,l) / cpp
252	zplay(ig0,l) = preff * pksurcp ** unskap
253	ztfi(ig0,l) = pteta(i,j,l) * pksurcp
254	pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
255	ig0 = ig0 + 1
256	ENDDO
257	ENDDO
258
259	pksurcp = ppk(1,jjp1,l) / cpp
260	zplay(ig0,l) = preff * pksurcp ** unskap
261	ztfi (ig0,l) = pteta(1,jjp1,l) * pksurcp
262	pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
263
264	ENDDO
265
266	c 43.bis traceurs
267	c ---------------
268	c
269	DO iq=1,nqtot
270	iiq=niadv(iq)
271	DO l=1,llm
272	zqfi(1,l,iq) = pq(1,1,l,iiq)
273	ig0 = 2
274	DO j=2,jjm
275	DO i = 1, iim
276	zqfi(ig0,l,iq) = pq(i,j,l,iiq)
277	ig0 = ig0 + 1
278	ENDDO
279	ENDDO
280	zqfi(ig0,l,iq) = pq(1,jjp1,l,iiq)
281	ENDDO
282	ENDDO
283
284	c convergence dynamique pour les traceurs "EAU"
285	! Earth-specific treatment of first 2 tracers (water)
286	if (planet_type=="earth") then
287	DO iq=1,2
288	DO l=1,llm
289	pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
290	ig0 = 2
291	DO j=2,jjm
292	DO i = 1, iim
293	pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
294	ig0 = ig0 + 1
295	ENDDO
296	ENDDO
297	pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
298	ENDDO
299	ENDDO
300	endif ! of if (planet_type=="earth")
301
302
303	c Geopotentiel calcule par rapport a la surface locale:
304	c -----------------------------------------------------
305
306	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
307	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
308	DO l=1,llm
309	DO ig=1,ngridmx
310	zphi(ig,l)=zphi(ig,l)-zphis(ig)
311	ENDDO
312	ENDDO
313
314	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
315	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
316	c de masse est calclue dans advtrac.F
317	c DO l=1,llm
318	c pvervel(1,l)=pw(1,1,l) * g /apoln
319	c ig0=2
320	c DO j=2,jjm
321	c DO i = 1, iim
322	c pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
323	c ig0 = ig0 + 1
324	c ENDDO
325	c ENDDO
326	c pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
327	c ENDDO
328
329	c
330	c 45. champ u:
331	c ------------
332
333	DO 50 l=1,llm
334
335	DO 25 j=2,jjm
336	ig0 = 1+(j-2)*iim
337	zufi(ig0+1,l)= 0.5 *
338	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
339	pcvgu(ig0+1,l)= 0.5 *
340	$ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
341	DO 10 i=2,iim
342	zufi(ig0+i,l)= 0.5 *
343	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
344	pcvgu(ig0+i,l)= 0.5 *
345	$ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
346	10 CONTINUE
347	25 CONTINUE
348
349	50 CONTINUE
350
351
352	c 46.champ v:
353	c -----------
354
355	DO l=1,llm
356	DO j=2,jjm
357	ig0=1+(j-2)*iim
358	DO i=1,iim
359	zvfi(ig0+i,l)= 0.5 *
360	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
361	pcvgv(ig0+i,l)= 0.5 *
362	$ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
363	ENDDO
364	ENDDO
365	ENDDO
366
367
368	c 47. champs de vents aux pole nord
369	c ------------------------------
370	c U = 1 / pi * integrale [ v * cos(long) * d long ]
371	c V = 1 / pi * integrale [ v * sin(long) * d long ]
372
373	DO l=1,llm
374
375	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
376	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
377	DO i=2,iim
378	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
379	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
380	ENDDO
381
382	DO i=1,iim
383	zcos(i) = COS(rlonv(i))*z1(i)
384	zcosbis(i)= COS(rlonv(i))*z1bis(i)
385	zsin(i) = SIN(rlonv(i))*z1(i)
386	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
387	ENDDO
388
389	zufi(1,l) = SSUM(iim,zcos,1)/pi
390	pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
391	zvfi(1,l) = SSUM(iim,zsin,1)/pi
392	pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
393
394	ENDDO
395
396
397	c 48. champs de vents aux pole sud:
398	c ---------------------------------
399	c U = 1 / pi * integrale [ v * cos(long) * d long ]
400	c V = 1 / pi * integrale [ v * sin(long) * d long ]
401
402	DO l=1,llm
403
404	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
405	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
406	DO i=2,iim
407	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
408	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
409	ENDDO
410
411	DO i=1,iim
412	zcos(i) = COS(rlonv(i))*z1(i)
413	zcosbis(i) = COS(rlonv(i))*z1bis(i)
414	zsin(i) = SIN(rlonv(i))*z1(i)
415	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
416	ENDDO
417
418	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
419	pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
420	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
421	pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
422
423	ENDDO
424	c
425	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
426	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
427
428	c-----------------------------------------------------------------------
429	c Appel de la physique:
430	c ---------------------
431
432
433
434	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
435	zdt_split=dtphys/nsplit_phys
436	zdufic(:,:)=0.
437	zdvfic(:,:)=0.
438	zdtfic(:,:)=0.
439	zdqfic(:,:,:)=0.
440
441	#ifdef CPP_PHYS
442
443	do isplit=1,nsplit_phys
444
445	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
446	debut_split=debut.and.isplit==1
447	lafin_split=lafin.and.isplit==nsplit_phys
448
449	if (planet_type=="earth") then
450
451	CALL physiq (ngridmx,
452	. llm,
453	. debut_split,
454	. lafin_split,
455	. jD_cur,
456	. jH_cur_split,
457	. zdt_split,
458	. zplev,
459	. zplay,
460	. zphi,
461	. zphis,
462	. presnivs,
463	. zufi,
464	. zvfi,
465	. ztfi,
466	. zqfi,
467	. flxwfi,
468	. zdufi,
469	. zdvfi,
470	. zdtfi,
471	. zdqfi,
472	. zdpsrf,
473	. pducov)
474
475	else if ( planet_type=="generic" ) then
476
477	CALL physiq (ngridmx, !! ngrid
478	. llm, !! nlayer
479	. nqtot, !! nq
480	. tname, !! tracer names from dynamical core (given in infotrac)
481	. debut_split, !! firstcall
482	. lafin_split, !! lastcall
483	. jD_cur, !! pday. see leapfrog
484	. jH_cur_split, !! ptime "fraction of day"
485	. zdt_split, !! ptimestep
486	. zplev, !! pplev
487	. zplay, !! pplay
488	. zphi, !! pphi
489	. zufi, !! pu
490	. zvfi, !! pv
491	. ztfi, !! pt
492	. zqfi, !! pq
493	. flxwfi, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
494	. zdufi, !! pdu
495	. zdvfi, !! pdv
496	. zdtfi, !! pdt
497	. zdqfi, !! pdq
498	. zdpsrf, !! pdpsrf
499	. tracerdyn) !! tracerdyn <-- utilite ???
500
501	endif ! of if (planet_type=="earth")
502
503	zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
504	zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
505	ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
506	zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
507
508	zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
509	zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
510	zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
511	zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
512
513	enddo ! of do isplit=1,nsplit_phys
514
515	#endif
516	! of #ifdef CPP_PHYS
517
518	zdufi(:,:)=zdufic(:,:)/nsplit_phys
519	zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
520	zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
521	zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
522
523
524	500 CONTINUE
525
526	c-----------------------------------------------------------------------
527	c transformation des tendances physiques en tendances dynamiques:
528	c ---------------------------------------------------------------
529
530	c tendance sur la pression :
531	c -----------------------------------
532
533	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
534	c
535	c 62. enthalpie potentielle
536	c ---------------------
537
538	DO l=1,llm
539
540	DO i=1,iip1
541	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
542	pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
543	ENDDO
544
545	DO j=2,jjm
546	ig0=1+(j-2)*iim
547	DO i=1,iim
548	pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
549	ENDDO
550	pdhfi(iip1,j,l) = pdhfi(1,j,l)
551	ENDDO
552
553	ENDDO
554
555
556	c 62. humidite specifique
557	c ---------------------
558	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
559	! DO iq=1,nqtot
560	! DO l=1,llm
561	! DO i=1,iip1
562	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
563	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
564	! ENDDO
565	! DO j=2,jjm
566	! ig0=1+(j-2)*iim
567	! DO i=1,iim
568	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
569	! ENDDO
570	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
571	! ENDDO
572	! ENDDO
573	! ENDDO
574
575	c 63. traceurs
576	c ------------
577	C initialisation des tendances
578	pdqfi(:,:,:,:)=0.
579	C
580	DO iq=1,nqtot
581	iiq=niadv(iq)
582	DO l=1,llm
583	DO i=1,iip1
584	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
585	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
586	ENDDO
587	DO j=2,jjm
588	ig0=1+(j-2)*iim
589	DO i=1,iim
590	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
591	ENDDO
592	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
593	ENDDO
594	ENDDO
595	ENDDO
596
597	c 65. champ u:
598	c ------------
599
600	DO l=1,llm
601
602	DO i=1,iip1
603	pdufi(i,1,l) = 0.
604	pdufi(i,jjp1,l) = 0.
605	ENDDO
606
607	DO j=2,jjm
608	ig0=1+(j-2)*iim
609	DO i=1,iim-1
610	pdufi(i,j,l)=
611	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
612	ENDDO
613	pdufi(iim,j,l)=
614	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
615	pdufi(iip1,j,l)=pdufi(1,j,l)
616	ENDDO
617
618	ENDDO
619
620
621	c 67. champ v:
622	c ------------
623
624	DO l=1,llm
625
626	DO j=2,jjm-1
627	ig0=1+(j-2)*iim
628	DO i=1,iim
629	pdvfi(i,j,l)=
630	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
631	ENDDO
632	pdvfi(iip1,j,l) = pdvfi(1,j,l)
633	ENDDO
634	ENDDO
635
636
637	c 68. champ v pres des poles:
638	c ---------------------------
639	c v = U * cos(long) + V * SIN(long)
640
641	DO l=1,llm
642
643	DO i=1,iim
644	pdvfi(i,1,l)=
645	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
646	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
647	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
648	pdvfi(i,1,l)=
649	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
650	pdvfi(i,jjm,l)=
651	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
652	ENDDO
653
654	pdvfi(iip1,1,l) = pdvfi(1,1,l)
655	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
656
657	ENDDO
658
659	c-----------------------------------------------------------------------
660
661	700 CONTINUE
662
663	firstcal = .FALSE.
664
665	RETURN
666	END

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

Download in other formats:

Original Format