Context Navigation

calfis.F @ 2343

Last change on this file since 2343 was 2333, checked in by lguez, 9 years ago

New parameterization of gravity wave drag due to front/jet systems, by

de la Camara and F. Lott. The new Camara-Lott parameterization

replaces the Hines parameterization so it is activated if not ok_hines
and ok_gwd_rando.

Also changed distribution of phase speeds in FLOTT_GWD_rando, from
uniform to Gaussian. Bug fix in sugwd_strato. Bug fix in the arguments
of the call to add_phys_tend for methane oxydation.

For the new Camara-Lott parameterization, we need to compute relative
vorticity in calfis and pass it as a new argument "rot" to
physiq. Interpolation of relative vorticity to the physics grid is not
optimal for now: it is not weighted by cell areas.

Alvaro de la Camara, Fran\c{}cois Lott

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

Line
1	!
2	! $Id: calfis.F 2333 2015-07-22 14:14:51Z 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 zrot(iip1,jjm,llm) ! AdlC May 2014
142	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm), zrfi(ngridmx,llm)
143	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
144	c
145	REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
146	REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
147	c
148	REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
149	REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
150	REAL zdpsrf(ngridmx)
151	c
152	REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
153	REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
154	REAL jH_cur_split,zdt_split
155	LOGICAL debut_split,lafin_split
156	INTEGER isplit
157
158	REAL zsin(iim),zcos(iim),z1(iim)
159	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
160	REAL unskap, pksurcp
161	c
162	REAL flxwfi(ngridmx,llm) ! Flux de masse verticale sur la grille physiq
163	c
164
165	REAL SSUM
166
167	LOGICAL,SAVE :: firstcal=.true., debut=.true.
168	! REAL rdayvrai
169
170	LOGICAL tracerdyn
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	if (planet_type=="earth") then
471
472	CALL physiq (ngridmx,
473	. llm,
474	. debut_split,
475	. lafin_split,
476	. jD_cur,
477	. jH_cur_split,
478	. zdt_split,
479	. zplev,
480	. zplay,
481	. zphi,
482	. zphis,
483	. presnivs,
484	. zufi,
485	. zvfi, zrfi,
486	. ztfi,
487	. zqfi,
488	. flxwfi,
489	. zdufi,
490	. zdvfi,
491	. zdtfi,
492	. zdqfi,
493	. zdpsrf,
494	. pducov)
495
496	else if ( planet_type=="generic" ) then
497
498	CALL physiq (ngridmx, !! ngrid
499	. llm, !! nlayer
500	. nqtot, !! nq
501	. tname, !! tracer names from dynamical core (given in infotrac)
502	. debut_split, !! firstcall
503	. lafin_split, !! lastcall
504	. jD_cur, !! pday. see leapfrog
505	. jH_cur_split, !! ptime "fraction of day"
506	. zdt_split, !! ptimestep
507	. zplev, !! pplev
508	. zplay, !! pplay
509	. zphi, !! pphi
510	. zufi, !! pu
511	. zvfi, !! pv
512	. ztfi, !! pt
513	. zqfi, !! pq
514	. flxwfi, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
515	. zdufi, !! pdu
516	. zdvfi, !! pdv
517	. zdtfi, !! pdt
518	. zdqfi, !! pdq
519	. zdpsrf, !! pdpsrf
520	. tracerdyn) !! tracerdyn <-- utilite ???
521
522	endif ! of if (planet_type=="earth")
523
524	zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
525	zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
526	ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
527	zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
528
529	zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
530	zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
531	zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
532	zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
533
534	enddo ! of do isplit=1,nsplit_phys
535
536	#endif
537	! of #ifdef CPP_PHYS
538
539	zdufi(:,:)=zdufic(:,:)/nsplit_phys
540	zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
541	zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
542	zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
543
544
545	500 CONTINUE
546
547	c-----------------------------------------------------------------------
548	c transformation des tendances physiques en tendances dynamiques:
549	c ---------------------------------------------------------------
550
551	c tendance sur la pression :
552	c -----------------------------------
553
554	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
555	c
556	c 62. enthalpie potentielle
557	c ---------------------
558
559	DO l=1,llm
560
561	DO i=1,iip1
562	pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l)
563	pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l)
564	ENDDO
565
566	DO j=2,jjm
567	ig0=1+(j-2)*iim
568	DO i=1,iim
569	pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l)
570	ENDDO
571	pdhfi(iip1,j,l) = pdhfi(1,j,l)
572	ENDDO
573
574	ENDDO
575
576
577	c 62. humidite specifique
578	c ---------------------
579	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
580	! DO iq=1,nqtot
581	! DO l=1,llm
582	! DO i=1,iip1
583	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
584	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
585	! ENDDO
586	! DO j=2,jjm
587	! ig0=1+(j-2)*iim
588	! DO i=1,iim
589	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
590	! ENDDO
591	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
592	! ENDDO
593	! ENDDO
594	! ENDDO
595
596	c 63. traceurs
597	c ------------
598	C initialisation des tendances
599	pdqfi(:,:,:,:)=0.
600	C
601	DO iq=1,nqtot
602	iiq=niadv(iq)
603	DO l=1,llm
604	DO i=1,iip1
605	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
606	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
607	ENDDO
608	DO j=2,jjm
609	ig0=1+(j-2)*iim
610	DO i=1,iim
611	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
612	ENDDO
613	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
614	ENDDO
615	ENDDO
616	ENDDO
617
618	c 65. champ u:
619	c ------------
620
621	DO l=1,llm
622
623	DO i=1,iip1
624	pdufi(i,1,l) = 0.
625	pdufi(i,jjp1,l) = 0.
626	ENDDO
627
628	DO j=2,jjm
629	ig0=1+(j-2)*iim
630	DO i=1,iim-1
631	pdufi(i,j,l)=
632	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
633	ENDDO
634	pdufi(iim,j,l)=
635	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
636	pdufi(iip1,j,l)=pdufi(1,j,l)
637	ENDDO
638
639	ENDDO
640
641
642	c 67. champ v:
643	c ------------
644
645	DO l=1,llm
646
647	DO j=2,jjm-1
648	ig0=1+(j-2)*iim
649	DO i=1,iim
650	pdvfi(i,j,l)=
651	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
652	ENDDO
653	pdvfi(iip1,j,l) = pdvfi(1,j,l)
654	ENDDO
655	ENDDO
656
657
658	c 68. champ v pres des poles:
659	c ---------------------------
660	c v = U * cos(long) + V * SIN(long)
661
662	DO l=1,llm
663
664	DO i=1,iim
665	pdvfi(i,1,l)=
666	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
667	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
668	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
669	pdvfi(i,1,l)=
670	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
671	pdvfi(i,jjm,l)=
672	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
673	ENDDO
674
675	pdvfi(iip1,1,l) = pdvfi(1,1,l)
676	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
677
678	ENDDO
679
680	c-----------------------------------------------------------------------
681
682	700 CONTINUE
683
684	firstcal = .FALSE.
685
686	RETURN
687	END

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

Download in other formats:

Original Format