Context Navigation

calfis.F @ 1549

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

All GCMs:
Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation (up to rev r2420 of LMDZ5)

all physics packages:
added module callphysiq_mod.F90 in dynphy_lonlat/phy* which contains the routine "call_physiq" which is called by calfis* and calls the physics. This way different "physiq" routine from different physics packages may be called: The calfis* routines now exposes all available fields that might be transmitted to physiq but which is actually send (ie: expected/needed by physiq) is decided in call_physiq.
turned "physiq.F[90]" into module "physiq_mod.F[90]" for better control of "physiq" arguments. for phyvenus/phytitan, extracted gr_fi_ecrit from physiq.F as gr_fi_ecrit.F90 (note that it can only work in serial).

misc:
updated wxios.F90 to keep up with LMDZ5 modifications.

dyn3d_common:
infotrac.F90 keep up with LMDZ5 modifications (cosmetics)

dyn3d:
gcm.F90: cosmetic cleanup.
leapfrog.F90: fix computation of date as function of itau.

dyn3dpar:
gcm.F: cosmetic cleanup.
leapfrog_p.F90: fix computation of date as function of itau.

NB: physics are given the date corresponding to the end of the
physics step.

dynphy_lonlat:
calfis.F : added computation of relative wind vorticity.
calfis_p.F: added computation of relative wind vorticity (input required by Earth physics)

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

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

Download in other formats:

Original Format