Context Navigation

calfis.F @ 1459

Last change on this file since 1459 was 1459, checked in by emillour, 9 years ago

Common dynamics: A couple of bug fixes

calfis[_p].F array boundaries must be explicitely specified when underlying arrays are of different sizes.
advect_new_p.F : missing initializations of intermediate variables at topmost layer.

File size: 29.6 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, ONLY: nqtot, niadv, tname
32	USE control_mod, ONLY: planet_type, nsplit_phys
33	USE write_field
34	USE cpdet_mod, only: t2tpot,tpot2t
35
36	! used only for zonal averages
37	USE moyzon_mod
38	USE comvert_mod, ONLY: presnivs,preff
39	USE comconst_mod, ONLY: daysec,dtvr,dtphys,kappa,cpp,g,rad,pi
40	USE logic_mod, ONLY: moyzon_ch,moyzon_mu
41
42	IMPLICIT NONE
43	c=======================================================================
44	c
45	c 1. rearrangement des tableaux et transformation
46	c variables dynamiques > variables physiques
47	c 2. calcul des termes physiques
48	c 3. retransformation des tendances physiques en tendances dynamiques
49	c
50	c remarques:
51	c ----------
52	c
53	c - les vents sont donnes dans la physique par leurs composantes
54	c naturelles.
55	c - la variable thermodynamique de la physique est une variable
56	c intensive : T
57	c pour la dynamique on prend T * ( preff / p(l) ) **kappa
58	c - les deux seules variables dependant de la geometrie necessaires
59	c pour la physique sont la latitude pour le rayonnement et
60	c l'aire de la maille quand on veut integrer une grandeur
61	c horizontalement.
62	c - les points de la physique sont les points scalaires de la
63	c la dynamique; numerotation:
64	c 1 pour le pole nord
65	c (jjm-1)*iim pour l'interieur du domaine
66	c ngridmx pour le pole sud
67	c ---> ngridmx=2+(jjm-1)*iim
68	c
69	c Input :
70	c -------
71	c pucov covariant zonal velocity
72	c pvcov covariant meridional velocity
73	c pteta potential temperature
74	c pps surface pressure
75	c pmasse masse d'air dans chaque maille
76	c pts surface temperature (K)
77	c callrad clef d'appel au rayonnement
78	c
79	c Output :
80	c --------
81	c pdufi tendency for the natural zonal velocity (ms-1)
82	c pdvfi tendency for the natural meridional velocity
83	c pdhfi tendency for the potential temperature (K/s)
84	c pdtsfi tendency for the surface temperature
85	c
86	c pdtrad radiative tendencies \ both input
87	c pfluxrad radiative fluxes / and output
88	c
89	c=======================================================================
90	c
91	c-----------------------------------------------------------------------
92	c
93	c 0. Declarations :
94	c ------------------
95
96	include "dimensions.h"
97	include "paramet.h"
98
99	INTEGER ngridmx
100	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
101
102	include "comgeom2.h"
103	include "iniprint.h"
104
105	c Arguments :
106	c -----------
107	LOGICAL,INTENT(IN) :: lafin ! .true. for the very last call to physics
108	REAL,INTENT(IN) :: jD_cur, jH_cur
109	REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
110	REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
111	REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
112	REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
113	REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
114	REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
115	REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
116
117	REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
118	REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
119	REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
120	! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
121	! qui lui meme ne sert a rien dans la routine telle qu'elle est
122	! ecrite, et que j'ai donc commente....
123	REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
124	! NB: pdq is only used to compute pcvgq which is in fact not used...
125
126	REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
127	REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
128	REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
129	REAL,INTENT(IN) :: flxw(iip1,jjp1,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) (on llm because flxw(:,:,llm+1)=0)
130
131	! tendencies (in */s) from the physics
132	REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
133	REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
134	REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
135	REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
136	REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
137
138	#ifndef CPP_PARA
139	c Local variables :
140	c -----------------
141
142	INTEGER i,j,l,ig0,ig,iq,iiq
143	REAL zpsrf(ngridmx)
144	REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm)
145	REAL zphi(ngridmx,llm),zphis(ngridmx)
146
147	REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
148	REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
149	! ADAPTATION GCM POUR CP(T)
150	REAL zteta(ngridmx,llm)
151	REAL zpk(ngridmx,llm)
152
153	! RQ SL 13/10/10:
154	! Ces calculs ne servent pas.
155	! Si necessaire, decommenter ces variables et les calculs...
156	! REAL pcvgu(ngridmx,llm), pcvgv(ngridmx,llm)
157	! REAL pcvgt(ngridmx,llm), pcvgq(ngridmx,llm,2)
158
159	REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
160	REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
161	REAL zdpsrf(ngridmx)
162
163	REAL zdufic(ngridmx,llm),zdvfic(ngridmx,llm)
164	REAL zdtfic(ngridmx,llm),zdqfic(ngridmx,llm,nqtot)
165	REAL jH_cur_split,zdt_split
166	LOGICAL debut_split,lafin_split
167	INTEGER isplit
168
169	REAL zsin(iim),zcos(iim),z1(iim)
170	REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
171	REAL unskap, pksurcp
172	save unskap
173
174	REAL flxwfi(ngridmx,llm) ! Flux de masse verticale sur la grille physiq
175
176	REAL SSUM
177
178	LOGICAL,SAVE :: firstcal=.true., debut=.true.
179	! REAL rdayvrai
180
181	LOGICAL tracerdyn ! for generic/mars physics call ; possibly to get rid of
182
183	! For Titan only right now:
184	! to allow for 2D computation of microphys and chemistry
185	LOGICAL,save :: flag_moyzon
186	REAL,allocatable,save :: tmpvar(:,:)
187	REAL,allocatable,save :: tmpvarp1(:,:)
188	REAL,allocatable,save :: tmpvarbar(:)
189	REAL,allocatable,save :: tmpvarbarp1(:)
190	real :: zz1,zz2
191
192	c-----------------------------------------------------------------------
193
194	c 1. Initialisations :
195	c --------------------
196
197
198	IF ( firstcal ) THEN
199	debut = .TRUE.
200	IF (ngridmx.NE.2+(jjm-1)*iim) THEN
201	write(lunout,*) 'STOP dans calfis'
202	write(lunout,*)
203	& 'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
204	write(lunout,*) ' ngridmx jjm iim '
205	write(lunout,*) ngridmx,jjm,iim
206	STOP
207	ENDIF
208
209	unskap = 1./ kappa
210
211	flag_moyzon = .false.
212	if(moyzon_ch.or.moyzon_mu) then
213	flag_moyzon = .true.
214	allocate(tmpvar(iip1,llm))
215	allocate(tmpvarp1(iip1,llmp1))
216	allocate(tmpvarbar(llm))
217	allocate(tmpvarbarp1(llmp1))
218	endif
219
220	if (flag_moyzon) call moyzon_init(ngridmx,llm,nqtot)
221
222	c----------------------------------------------
223	c moyennes globales pour le profil de pression
224	if(planet_type.eq."titan".or.planet_type.eq."venus") then
225	ALLOCATE(plevmoy(llm+1))
226	ALLOCATE(playmoy(llm))
227	ALLOCATE(tmoy(llm))
228	ALLOCATE(tetamoy(llm))
229	ALLOCATE(pkmoy(llm))
230	ALLOCATE(phimoy(0:llm))
231	ALLOCATE(zlevmoy(llm+1))
232	ALLOCATE(zlaymoy(llm))
233	plevmoy=0.
234	do l=1,llmp1
235	do i=1,iip1
236	do j=1,jjp1
237	plevmoy(l)=plevmoy(l)+pp(i,j,l)/(iip1*jjp1)
238	enddo
239	enddo
240	enddo
241	tetamoy=0.
242	pkmoy=0.
243	phimoy=0.
244	do i=1,iip1
245	do j=1,jjp1
246	phimoy(0)=phimoy(0)+pphis(i,j)/(iip1*jjp1)
247	enddo
248	enddo
249	do l=1,llm
250	do i=1,iip1
251	do j=1,jjp1
252	tetamoy(l)=tetamoy(l)+pteta(i,j,l)/(iip1*jjp1)
253	pkmoy(l)=pkmoy(l)+ppk(i,j,l)/(iip1*jjp1)
254	phimoy(l)=phimoy(l)+pphi(i,j,l)/(iip1*jjp1)
255	enddo
256	enddo
257	enddo
258	playmoy(:) = preff * (pkmoy(:)/cpp) ** unskap
259	call tpot2t(llm,tetamoy,tmoy,pkmoy)
260	c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
261	zlaymoy(1:llm) = gradrad/(g*rad-phimoy(1:llm))-rad
262	zlevmoy(1) = phimoy(0)/g
263	DO l=2,llm
264	zz1=(playmoy(l-1)+plevmoy(l))/(playmoy(l-1)-plevmoy(l))
265	zz2=(plevmoy(l) +playmoy(l))/(plevmoy(l) -playmoy(l))
266	zlevmoy(l)=(zz1zlaymoy(l-1)+zz2zlaymoy(l))/(zz1+zz2)
267	ENDDO
268	zlevmoy(llmp1)=zlaymoy(llm)+(zlaymoy(llm)-zlevmoy(llm))
269	c-------------------
270	c + lat index
271	allocate(klat(ngridmx))
272	klat=0
273	klat(1) = 1
274	ig0 = 2
275	DO j = 2,jjm
276	do i=0,iim-1
277	klat(ig0+i) = j
278	enddo
279	ig0 = ig0+iim
280	ENDDO
281	klat(ngridmx) = jjp1
282	endif ! planet_type=titan
283	c----------------------------------------------
284	ELSE
285	debut = .FALSE.
286	ENDIF ! of IF (firstcal)
287
288
289	c-----------------------------------------------------------------------
290	c 40. transformation des variables dynamiques en variables physiques:
291	c ---------------------------------------------------------------
292
293	c 41. pressions au sol (en Pascals)
294	c ----------------------------------
295
296	zpsrf(1) = pps(1,1)
297
298	ig0 = 2
299	DO j = 2,jjm
300	CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 )
301	ig0 = ig0+iim
302	ENDDO
303
304	zpsrf(ngridmx) = pps(1,jjp1)
305
306	c 42. pression intercouches et fonction d'Exner:
307
308	c -----------------------------------------------------------------
309	c .... zplev definis aux (llm +1) interfaces des couches ....
310	c .... zplay definis aux ( llm ) milieux des couches ....
311	c -----------------------------------------------------------------
312
313	c ... Exner = cp * ( p(l) / preff ) ** kappa ....
314
315	! ADAPTATION GCM POUR CP(T)
316	DO l = 1, llm
317	zpk( 1,l ) = ppk(1,1,l)
318	zteta( 1,l ) = pteta(1,1,l)
319	zplev( 1,l ) = pp(1,1,l)
320	ig0 = 2
321	DO j = 2, jjm
322	DO i =1, iim
323	zpk( ig0,l ) = ppk(i,j,l)
324	zteta( ig0,l ) = pteta(i,j,l)
325	zplev( ig0,l ) = pp(i,j,l)
326	ig0 = ig0 +1
327	ENDDO
328	ENDDO
329	zpk( ngridmx,l ) = ppk(1,jjp1,l)
330	zteta( ngridmx,l ) = pteta(1,jjp1,l)
331	zplev( ngridmx,l ) = pp(1,jjp1,l)
332	ENDDO
333	zplev( 1,llmp1 ) = pp(1,1,llmp1)
334	ig0 = 2
335	DO j = 2, jjm
336	DO i =1, iim
337	zplev( ig0,llmp1 ) = pp(i,j,llmp1)
338	ig0 = ig0 +1
339	ENDDO
340	ENDDO
341	zplev( ngridmx,llmp1 ) = pp(1,jjp1,llmp1)
342
343	if (flag_moyzon) then
344	tmpvarp1(:,:) = pp(:,1,:)
345	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
346	zplevbar(1,:) = tmpvarbarp1
347	tmpvar(:,:) = ppk(:,1,:)
348	call moyzon(llm,tmpvar,tmpvarbar)
349	zpkbar(1,:) = tmpvarbar
350	tmpvar(:,:) = pteta(:,1,:)
351	call moyzon(llm,tmpvar,tmpvarbar)
352	ztetabar(1,:) = tmpvarbar
353	call tpot2t(llm,ztetabar(1,:),ztfibar(1,:),zpkbar(1,:))
354	ig0 = 2
355	do j = 2, jjm
356	tmpvarp1(:,:) = pp(:,j,:)
357	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
358	zplevbar(ig0,:) = tmpvarbarp1
359	tmpvar(:,:) = ppk(:,j,:)
360	call moyzon(llm,tmpvar,tmpvarbar)
361	zpkbar(ig0,:) = tmpvarbar
362	tmpvar(:,:) = pteta(:,j,:)
363	call moyzon(llm,tmpvar,tmpvarbar)
364	ztetabar(ig0,:) = tmpvarbar
365	call tpot2t(llm,ztetabar(ig0,:),ztfibar(ig0,:),zpkbar(ig0,:))
366	ig0 = ig0+1
367	do i=2,iim
368	zplevbar(ig0,:) = zplevbar(ig0-1,:)
369	zpkbar(ig0,:) = zpkbar(ig0-1,:)
370	ztetabar(ig0,:) = ztetabar(ig0-1,:)
371	ztfibar(ig0,:) = ztfibar(ig0-1,:)
372	ig0 = ig0+1
373	enddo
374	enddo
375	tmpvarp1(:,:) = pp(:,jjp1,:)
376	call moyzon(llmp1,tmpvarp1,tmpvarbarp1)
377	zplevbar(ngridmx,:) = tmpvarbarp1
378	tmpvar(:,:) = ppk(:,jjp1,:)
379	call moyzon(llm,tmpvar,tmpvarbar)
380	zpkbar(ngridmx,:) = tmpvarbar
381	tmpvar(:,:) = pteta(:,jjp1,:)
382	call moyzon(llm,tmpvar,tmpvarbar)
383	ztetabar(ngridmx,:) = tmpvarbar
384	call tpot2t(llm,ztetabar(ngridmx,:),
385	. ztfibar(ngridmx,:),zpkbar(ngridmx,:))
386	endif
387
388	c 43. temperature naturelle (en K) et pressions milieux couches .
389	c ---------------------------------------------------------------
390
391	! ADAPTATION GCM POUR CP(T)
392	call tpot2t(ngridmx*llm,zteta,ztfi,zpk)
393
394	DO l=1,llm
395
396	pksurcp = ppk(1,1,l) / cpp
397	zplay(1,l) = preff * pksurcp ** unskap
398	! pcvgt(1,l) = pdteta(1,1,l) * pksurcp / pmasse(1,1,l)
399	ig0 = 2
400
401	DO j = 2, jjm
402	DO i = 1, iim
403	pksurcp = ppk(i,j,l) / cpp
404	zplay(ig0,l) = preff * pksurcp ** unskap
405	! pcvgt(ig0,l) = pdteta(i,j,l) * pksurcp / pmasse(i,j,l)
406	ig0 = ig0 + 1
407	ENDDO
408	ENDDO
409
410	pksurcp = ppk(1,jjp1,l) / cpp
411	zplay(ig0,l) = preff * pksurcp ** unskap
412	! pcvgt(ig0,l) = pdteta(1,jjp1,l) * pksurcp/ pmasse(1,jjp1,l)
413
414	ENDDO
415
416	if (flag_moyzon) then
417	zplaybar(:,:) = preff * (zpkbar(:,:)/cpp)**unskap
418	endif
419
420	c 43.bis traceurs (tous intensifs)
421	c ---------------
422
423	DO iq=1,nqtot
424	DO l=1,llm
425	zqfi(1,l,iq) = pq(1,1,l,iq)
426	ig0 = 2
427	DO j=2,jjm
428	DO i = 1, iim
429	zqfi(ig0,l,iq) = pq(i,j,l,iq)
430	ig0 = ig0 + 1
431	ENDDO
432	ENDDO
433	zqfi(ig0,l,iq) = pq(1,jjp1,l,iq)
434	ENDDO
435	ENDDO ! boucle sur traceurs
436
437	if (flag_moyzon) then
438	DO iq=1,nqtot
439	! RQ: REVOIR A QUOI CA SERT... ET VERIFIER...
440	! iiq=niadv(iq)
441	! en fait, iiq=iq...
442	! FIN RQ
443	tmpvar(:,:) = pq(:,1,:,iq)
444	call moyzon(llm,tmpvar,tmpvarbar)
445	zqfibar(1,:,iq) = tmpvarbar
446	ig0 = 2
447	do j = 2, jjm
448	tmpvar(:,:) = pq(:,j,:,iq)
449	call moyzon(llm,tmpvar,tmpvarbar)
450	zqfibar(ig0,:,iq) = tmpvarbar
451	ig0 = ig0+1
452	do i=2,iim
453	zqfibar(ig0,:,iq) = zqfibar(ig0-1,:,iq)
454	ig0 = ig0+1
455	enddo
456	enddo
457	tmpvar(:,:) = pq(:,jjp1,:,iq)
458	call moyzon(llm,tmpvar,tmpvarbar)
459	zqfibar(ngridmx,:,iq) = tmpvarbar
460	ENDDO ! of DO iq=1,nqtot
461	endif
462
463	! DEBUG
464	! do ig0=1,ngridmx
465	! write(*,'(6(e13.5,1x))') zqfibar(ig0,1,10),zqfi(ig0,1,10),
466	! . zqfibar(ig0,llm/2,10),zqfi(ig0,llm/2,10),
467	! . zqfibar(ig0,llm,10),zqfi(ig0,llm,10)
468	! enddo
469	! stop
470
471	!-----------------
472	! RQ SL 13/10/10:
473	! Ces calculs ne servent pas.
474	! Si necessaire, decommenter ces variables et les calculs...
475	!
476	! convergence dynamique pour les traceurs "EAU"
477	! Earth-specific treatment of first 2 tracers (water)
478	! if (planet_type=="earth") then
479	! DO iq=1,2
480	! DO l=1,llm
481	! pcvgq(1,l,iq)= pdq(1,1,l,iq) / pmasse(1,1,l)
482	! ig0 = 2
483	! DO j=2,jjm
484	! DO i = 1, iim
485	! pcvgq(ig0,l,iq) = pdq(i,j,l,iq) / pmasse(i,j,l)
486	! ig0 = ig0 + 1
487	! ENDDO
488	! ENDDO
489	! pcvgq(ig0,l,iq)= pdq(1,jjp1,l,iq) / pmasse(1,jjp1,l)
490	! ENDDO
491	! ENDDO
492	! endif ! of if (planet_type=="earth")
493	!----------------
494
495	c Geopotentiel calcule par rapport a la surface locale:
496	c -----------------------------------------------------
497
498	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi)
499	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis)
500	DO l=1,llm
501	DO ig=1,ngridmx
502	zphi(ig,l)=zphi(ig,l)-zphis(ig)
503	ENDDO
504	ENDDO
505
506	if (flag_moyzon) then
507	tmpvar(:,1) = pphis(:,1)
508	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
509	zphisbar(1) = tmpvarbar(1)
510	tmpvar(:,:) = pphi(:,1,:)
511	call moyzon(llm,tmpvar,tmpvarbar)
512	zphibar(1,:) = tmpvarbar
513	ig0 = 2
514	do j = 2, jjm
515	tmpvar(:,1) = pphis(:,j)
516	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
517	zphisbar(ig0) = tmpvarbar(1)
518	tmpvar(:,:) = pphi(:,j,:)
519	call moyzon(llm,tmpvar,tmpvarbar)
520	zphibar(ig0,:) = tmpvarbar
521	ig0 = ig0+1
522	do i=2,iim
523	zphisbar(ig0) = zphisbar(ig0-1)
524	zphibar(ig0,:) = zphibar(ig0-1,:)
525	ig0 = ig0+1
526	enddo
527	enddo
528	tmpvar(:,1) = pphis(:,jjp1)
529	call moyzon(1,tmpvar(:,1),tmpvarbar(1))
530	zphisbar(ngridmx) = tmpvarbar(1)
531	tmpvar(:,:) = pphi(:,jjp1,:)
532	call moyzon(llm,tmpvar,tmpvarbar)
533	zphibar(ngridmx,:) = tmpvarbar
534	endif
535
536	c .... Calcul de la vitesse verticale ( en Pams ou Kg/s ) ....
537	c JG : ancien calcule de omega utilise dans physiq.F. Maintenant le flux
538	c de masse est calclue dans advtrac.F
539	c DO l=1,llm
540	c pvervel(1,l)=pw(1,1,l) * g /apoln
541	c ig0=2
542	c DO j=2,jjm
543	c DO i = 1, iim
544	c pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
545	c ig0 = ig0 + 1
546	c ENDDO
547	c ENDDO
548	c pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
549	c ENDDO
550
551	c
552	c 45. champ u:
553	c ------------
554
555	DO 50 l=1,llm
556
557	DO 25 j=2,jjm
558	ig0 = 1+(j-2)*iim
559	zufi(ig0+1,l)= 0.5 *
560	$ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) )
561	! pcvgu(ig0+1,l)= 0.5 *
562	! $ ( pducov(iim,j,l)/cu(iim,j) + pducov(1,j,l)/cu(1,j) )
563	DO 10 i=2,iim
564	zufi(ig0+i,l)= 0.5 *
565	$ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) )
566	! pcvgu(ig0+i,l)= 0.5 *
567	! $ ( pducov(i-1,j,l)/cu(i-1,j) + pducov(i,j,l)/cu(i,j) )
568	10 CONTINUE
569	25 CONTINUE
570
571	50 CONTINUE
572
573
574	c 46.champ v:
575	c -----------
576
577	DO l=1,llm
578	DO j=2,jjm
579	ig0=1+(j-2)*iim
580	DO i=1,iim
581	zvfi(ig0+i,l)= 0.5 *
582	$ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) )
583	c pcvgv(ig0+i,l)= 0.5 *
584	c $ ( pdvcov(i,j-1,l)/cv(i,j-1) + pdvcov(i,j,l)/cv(i,j) )
585	ENDDO
586	ENDDO
587	ENDDO
588
589
590	c 47. champs de vents aux pole nord
591	c ------------------------------
592	c U = 1 / pi * integrale [ v * cos(long) * d long ]
593	c V = 1 / pi * integrale [ v * sin(long) * d long ]
594
595	DO l=1,llm
596
597	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,1,l)/cv(1,1)
598	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,1,l)/cv(1,1)
599	DO i=2,iim
600	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1)
601	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1)
602	ENDDO
603
604	DO i=1,iim
605	zcos(i) = COS(rlonv(i))*z1(i)
606	zcosbis(i)= COS(rlonv(i))*z1bis(i)
607	zsin(i) = SIN(rlonv(i))*z1(i)
608	zsinbis(i)= SIN(rlonv(i))*z1bis(i)
609	ENDDO
610
611	zufi(1,l) = SSUM(iim,zcos,1)/pi
612	! pcvgu(1,l) = SSUM(iim,zcosbis,1)/pi
613	zvfi(1,l) = SSUM(iim,zsin,1)/pi
614	! pcvgv(1,l) = SSUM(iim,zsinbis,1)/pi
615
616	ENDDO
617
618
619	c 48. champs de vents aux pole sud:
620	c ---------------------------------
621	c U = 1 / pi * integrale [ v * cos(long) * d long ]
622	c V = 1 / pi * integrale [ v * sin(long) * d long ]
623
624	DO l=1,llm
625
626	z1(1) =(rlonu(1)-rlonu(iim)+2.pi)pvcov(1,jjm,l)/cv(1,jjm)
627	z1bis(1)=(rlonu(1)-rlonu(iim)+2.pi)pdvcov(1,jjm,l)/cv(1,jjm)
628	DO i=2,iim
629	z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm)
630	z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm)
631	ENDDO
632
633	DO i=1,iim
634	zcos(i) = COS(rlonv(i))*z1(i)
635	zcosbis(i) = COS(rlonv(i))*z1bis(i)
636	zsin(i) = SIN(rlonv(i))*z1(i)
637	zsinbis(i) = SIN(rlonv(i))*z1bis(i)
638	ENDDO
639
640	zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi
641	! pcvgu(ngridmx,l) = SSUM(iim,zcosbis,1)/pi
642	zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi
643	! pcvgv(ngridmx,l) = SSUM(iim,zsinbis,1)/pi
644
645	ENDDO
646	c
647	c On change de grille, dynamique vers physiq, pour le flux de masse verticale
648	CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
649
650	c-----------------------------------------------------------------------
651	c Appel de la physique:
652	c ---------------------
653
654	! Appel de la physique: pose probleme quand on tourne
655	! SANS physique, car physiq.F est dans le repertoire phy[]...
656	! Il faut une cle CPP_PHYS
657
658	! Le fait que les arguments de physiq soient differents selon les planetes
659	! ne pose pas de probleme a priori.
660
661	! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
662	zdt_split=dtphys/nsplit_phys
663	zdufic(:,:)=0.
664	zdvfic(:,:)=0.
665	zdtfic(:,:)=0.
666	zdqfic(:,:,:)=0.
667
668	#ifdef CPP_PHYS
669
670	do isplit=1,nsplit_phys
671
672	jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
673	debut_split=debut.and.isplit==1
674	lafin_split=lafin.and.isplit==nsplit_phys
675
676	if (planet_type.eq."earth") then
677	CALL physiq (ngridmx,
678	. llm,
679	. debut_split,
680	. lafin_split,
681	. jD_cur,
682	. jH_cur_split,
683	. zdt_split,
684	. zplev,
685	. zplay,
686	. zphi,
687	. zphis,
688	. presnivs,
689	. zufi,
690	. zvfi,
691	. ztfi,
692	. zqfi,
693	. flxwfi,
694	. zdufi,
695	. zdvfi,
696	. zdtfi,
697	. zdqfi,
698	. zdpsrf,
699	. pducov)
700
701	else if ( planet_type=="generic" ) then
702
703	CALL physiq (ngridmx, !! ngrid
704	. llm, !! nlayer
705	. nqtot, !! nq
706	. tname, !! tracer names from dynamical core (given in infotrac)
707	. debut_split, !! firstcall
708	. lafin_split, !! lastcall
709	. jD_cur, !! pday. see leapfrog
710	. jH_cur_split, !! ptime "fraction of day"
711	. zdt_split, !! ptimestep
712	. zplev, !! pplev
713	. zplay, !! pplay
714	. zphi, !! pphi
715	. zufi, !! pu
716	. zvfi, !! pv
717	. ztfi, !! pt
718	. zqfi, !! pq
719	. flxwfi, !! pw !! or 0. anyway this is for diagnostic. not used in physiq.
720	. zdufi, !! pdu
721	. zdvfi, !! pdv
722	. zdtfi, !! pdt
723	. zdqfi, !! pdq
724	. zdpsrf, !! pdpsrf
725	. tracerdyn) !! tracerdyn <-- utilite ???
726
727	else if ( planet_type=="mars" ) then
728
729	CALL physiq (ngridmx, ! ngrid
730	. llm, ! nlayer
731	. nqtot, ! nq
732	. debut_split, ! firstcall
733	. lafin_split, ! lastcall
734	. jD_cur, ! pday
735	. jH_cur_split, ! ptime
736	. zdt_split, ! ptimestep
737	. zplev, ! pplev
738	. zplay, ! pplay
739	. zphi, ! pphi
740	. zufi, ! pu
741	. zvfi, ! pv
742	. ztfi, ! pt
743	. zqfi, ! pq
744	. flxwfi, ! pw
745	. zdufi, ! pdu
746	. zdvfi, ! pdv
747	. zdtfi, ! pdt
748	. zdqfi, ! pdq
749	. zdpsrf, ! pdpsrf
750	. tracerdyn) ! tracerdyn (somewhat obsolete)
751
752	else if ((planet_type=="titan").or.(planet_type=="venus")) then
753
754	CALL physiq (ngridmx,
755	. llm,
756	. nqtot,
757	. debut_split,
758	. lafin_split,
759	. jD_cur,
760	. jH_cur_split,
761	. zdt_split,
762	. zplev,
763	. zplay,
764	. zpk,
765	. zphi,
766	. zphis,
767	. presnivs,
768	. zufi,
769	. zvfi,
770	. ztfi,
771	. zqfi,
772	. flxwfi,
773	. zdufi,
774	. zdvfi,
775	. zdtfi,
776	. zdqfi,
777	. zdpsrf)
778
779	else ! unknown "planet_type"
780
781	write(lunout,*) "calfis_p: error, unknown planet_type: ",
782	& trim(planet_type)
783	stop
784
785	endif ! planet_type
786
787	zufi(:,:)=zufi(:,:)+zdufi(:,:)*zdt_split
788	zvfi(:,:)=zvfi(:,:)+zdvfi(:,:)*zdt_split
789	ztfi(:,:)=ztfi(:,:)+zdtfi(:,:)*zdt_split
790	zqfi(:,:,:)=zqfi(:,:,:)+zdqfi(:,:,:)*zdt_split
791
792	zdufic(:,:)=zdufic(:,:)+zdufi(:,:)
793	zdvfic(:,:)=zdvfic(:,:)+zdvfi(:,:)
794	zdtfic(:,:)=zdtfic(:,:)+zdtfi(:,:)
795	zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
796
797	enddo ! of do isplit=1,nsplit_phys
798
799	! ATTENTION...
800	if (flag_moyzon.and.(nsplit_phys.ne.1)) then
801	print*,"WARNING ! flag_moyzon + nsplit_phys"
802	print*,"zqfibar n'est pas implemente au cours des iterations"
803	print*,"Donc a revoir..."
804	stop
805	endif
806
807	#endif
808	! #endif of #ifdef CPP_PHYS
809
810	zdufi(:,:)=zdufic(:,:)/nsplit_phys
811	zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
812	zdtfi(:,:)=zdtfic(:,:)/nsplit_phys
813	zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
814
815
816	500 CONTINUE
817
818	c-----------------------------------------------------------------------
819	c transformation des tendances physiques en tendances dynamiques:
820	c ---------------------------------------------------------------
821
822	c tendance sur la pression :
823	c -----------------------------------
824
825	CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi)
826	c
827	c 62. enthalpie potentielle
828	c ---------------------
829
830	! ADAPTATION GCM POUR CP(T)
831	call t2tpot(ngridmx*llm,ztfi,zteta,zpk)
832
833	DO i=1,iip1
834	pdhfi(i,1,:) = (zteta(1,:) - pteta(i,1,:))/dtphys
835	pdhfi(i,jjp1,:) = (zteta(ngridmx,:) - pteta(i,jjp1,:))/dtphys
836	ENDDO
837
838	DO j=2,jjm
839	ig0=1+(j-2)*iim
840	DO i=1,iim
841	pdhfi(i,j,:) = (zteta(ig0+i,:) - pteta(i,j,:))/dtphys
842	ENDDO
843	pdhfi(iip1,j,:) = pdhfi(1,j,:)
844	ENDDO
845
846
847	c 62. humidite specifique
848	c ---------------------
849	! Ehouarn: removed this useless bit: was overwritten at step 63 anyways
850	! DO iq=1,nqtot
851	! DO l=1,llm
852	! DO i=1,iip1
853	! pdqfi(i,1,l,iq) = zdqfi(1,l,iq)
854	! pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq)
855	! ENDDO
856	! DO j=2,jjm
857	! ig0=1+(j-2)*iim
858	! DO i=1,iim
859	! pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq)
860	! ENDDO
861	! pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq)
862	! ENDDO
863	! ENDDO
864	! ENDDO
865
866	c 63. traceurs (tous en intensifs)
867	c ------------
868	C initialisation des tendances
869	pdqfi(:,:,:,:)=0.
870	C
871	DO iq=1,nqtot
872	iiq=niadv(iq)
873	DO l=1,llm
874	DO i=1,iip1
875	pdqfi(i,1,l,iiq) = zdqfi(1,l,iq)
876	pdqfi(i,jjp1,l,iiq) = zdqfi(ngridmx,l,iq)
877	ENDDO
878	DO j=2,jjm
879	ig0=1+(j-2)*iim
880	DO i=1,iim
881	pdqfi(i,j,l,iiq) = zdqfi(ig0+i,l,iq)
882	ENDDO
883	pdqfi(iip1,j,l,iiq) = pdqfi(1,j,l,iq)
884	ENDDO
885	ENDDO
886	ENDDO
887
888	c 65. champ u:
889	c ------------
890
891	DO l=1,llm
892
893	DO i=1,iip1
894	pdufi(i,1,l) = 0.
895	pdufi(i,jjp1,l) = 0.
896	ENDDO
897
898	DO j=2,jjm
899	ig0=1+(j-2)*iim
900	DO i=1,iim-1
901	pdufi(i,j,l)=
902	$ 0.5(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))cu(i,j)
903	ENDDO
904	pdufi(iim,j,l)=
905	$ 0.5(zdufi(ig0+1,l)+zdufi(ig0+iim,l))cu(iim,j)
906	pdufi(iip1,j,l)=pdufi(1,j,l)
907	ENDDO
908
909	ENDDO
910
911
912	c 67. champ v:
913	c ------------
914
915	DO l=1,llm
916
917	DO j=2,jjm-1
918	ig0=1+(j-2)*iim
919	DO i=1,iim
920	pdvfi(i,j,l)=
921	$ 0.5(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))cv(i,j)
922	ENDDO
923	pdvfi(iip1,j,l) = pdvfi(1,j,l)
924	ENDDO
925	ENDDO
926
927
928	c 68. champ v pres des poles:
929	c ---------------------------
930	c v = U * cos(long) + V * SIN(long)
931
932	DO l=1,llm
933
934	DO i=1,iim
935	pdvfi(i,1,l)=
936	$ zdufi(1,l)COS(rlonv(i))+zdvfi(1,l)SIN(rlonv(i))
937	pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i))
938	$ +zdvfi(ngridmx,l)*SIN(rlonv(i))
939	pdvfi(i,1,l)=
940	$ 0.5(pdvfi(i,1,l)+zdvfi(i+1,l))cv(i,1)
941	pdvfi(i,jjm,l)=
942	$ 0.5(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))cv(i,jjm)
943	ENDDO
944
945	pdvfi(iip1,1,l) = pdvfi(1,1,l)
946	pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l)
947
948	ENDDO
949
950	c-----------------------------------------------------------------------
951
952	700 CONTINUE
953
954	firstcal = .FALSE.
955
956	#endif
957	! of #ifndef CPP_PARA
958	END

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

Download in other formats:

Original Format