Context Navigation

calfis.F @ 2597

Last change on this file since 2597 was 2597, checked in by Ehouarn Millour, 8 years ago
Cleanup in the dynamics: get rid of comconst.h, make it a module comconst_mod. 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: 21.3 KB

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

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

Download in other formats:

Original Format