Context Navigation

calfis.F @ 1017

Last change on this file since 1017 was 1017, checked in by emillour, 11 years ago

Common dynamics: (and collateral adaptations in Venus physics)
Improved cpdet routines in and additional sponge mode:

Additionnal sponge mode (trigered with "callsponge" flag), in line with the one used in the Generic and Martian GCM. This sponge is called whenever there is a dissipation step.
Improvement of the cpdet routines : created routines tpot2t_glo_p and t2tpot_glo_p which handle fields on the whole dynamics (scaler) grid, which are more efficient than calling tpot2t_p or t2tpot_p with slabs of data (generated use of intermediate copies of these chunks of data at every call)
Turned cpdet.F into a module cpdet_mod.F90 (and correspondingly adapted all routines in the Venus physics).

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

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

Download in other formats:

Original Format