Context Navigation

calfis.F @ 2011

Last change on this file since 2011 was 1987, checked in by Ehouarn Millour, 11 years ago
Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency). Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid). EM
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: 20.5 KB

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

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

Download in other formats:

Original Format