Context Navigation

calfis.F @ 5

Last change on this file since 5 was 5, checked in by slebonnois, 14 years ago

Ajout des fichiers .def Venus et Titan (tels qu'ils sont utilisés

actuellement) dans les deftanks.

Ajout d'une doc sur Cp(T).
Modifications dans dyn3d concernant Cp(T), cf le log (v5) dans chantiers
Premières modifs de l'appel à la physique dans dyn3d/calfis, cf log (v5)
Elimination des cpdet.* dans phytitan et phyvenus (remplacés par cpdet.F

dans dyn3d).

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

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

Download in other formats:

Original Format