Context Navigation

bilan_dyn.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: 27.6 KB

Line
1	!
2	! $Id: bilan_dyn.F 1403 2010-07-01 09:02:53Z fairhead $
3	!
4	SUBROUTINE bilan_dyn (dt_app,dt_cum,
5	s ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,
6	s ducovdyn,ducovdis,ducovspg,ducovphy)
7	c si besoin des traceurs:
8	c SUBROUTINE bilan_dyn (ntrac,dt_app,dt_cum,
9	c s ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac,
10	c s ducovdyn,ducovdis,ducovspg,ducovphy)
11
12	c AFAIRE
13	c Prevoir en champ nq+1 le diagnostique de l'energie
14	c en faisant Qzon=Cv T + L * ...
15	c vQ..A=Cp T + L * ...
16
17	#ifdef CPP_IOIPSL
18	USE IOIPSL
19	#endif
20
21	USE control_mod, ONLY: planet_type
22	USE cpdet_mod, only: tpot2t
23
24	IMPLICIT NONE
25
26	#include "dimensions.h"
27	#include "paramet.h"
28	#include "comconst.h"
29	#include "comvert.h"
30	#include "comgeom2.h"
31	#include "temps.h"
32	#include "iniprint.h"
33
34	c====================================================================
35	c
36	c Sous-programme consacre à des diagnostics dynamiques de base
37	c
38	c
39	c De facon generale, les moyennes des scalaires Q sont ponderees par
40	c la masse.
41	c
42	c Les flux de masse sont eux simplement moyennes.
43	c
44	c====================================================================
45
46	c Arguments :
47	c ===========
48
49	c integer ntrac
50	real dt_app,dt_cum
51	real ps(iip1,jjp1)
52	real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
53	real flux_u(iip1,jjp1,llm)
54	real flux_v(iip1,jjm,llm)
55	real teta(iip1,jjp1,llm)
56	real phi(iip1,jjp1,llm)
57	real ucov(iip1,jjp1,llm)
58	real vcov(iip1,jjm,llm)
59	c real trac(iip1,jjp1,llm,ntrac)
60	c Tendances en m/s2 :
61	real ducovdyn(iip1,jjp1,llm)
62	real ducovdis(iip1,jjp1,llm)
63	real ducovspg(iip1,jjp1,llm)
64	real ducovphy(iip1,jjp1,llm)
65
66	c Local :
67	c =======
68
69	integer icum,ncum
70	save icum,ncum
71
72	integer i,j,l,iQ,num
73	real zz,zqy,zfactv(jjm,llm),zfactw(jjm,llm)
74	character*2 strd2
75	real ww
76
77	logical first
78	save first
79	data first/.true./
80
81	integer i_sortie
82	save i_sortie
83	data i_sortie/1/
84
85	real time
86	integer itau
87	save time,itau
88	data time,itau/0.,0/
89
90	! facteur = -1. pour Venus
91	real fact_geovenus
92	save fact_geovenus
93
94	c variables dynamiques intermédiaires
95	c -----------------------------------
96	REAL vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
97	REAL ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
98	REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
99	REAL vorpot(iip1,jjm,llm)
100	REAL w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm)
101	real temp(iip1,jjp1,llm)
102	real dudyn(iip1,jjp1,llm)
103	real dudis(iip1,jjp1,llm)
104	real duspg(iip1,jjp1,llm)
105	real duphy(iip1,jjp1,llm)
106
107	c CHAMPS SCALAIRES Q ADVECTES
108	c ----------------------------
109	integer nQ
110	c avec tous les composes, ca fait trop.... Je les enleve
111	c parameter (nQ=6+nqmx)
112	parameter (nQ=6)
113
114	character*6,save :: nom(nQ)
115	character*6,save :: unites(nQ)
116
117	integer itemp,igeop,iecin,iang,iu,iun
118	save itemp,igeop,iecin,iang,iu,iun
119	data itemp,igeop,iecin,iang,iu,iun/1,2,3,4,5,6/
120
121	c champ contenant les scalaires advectés.
122	real Q(iip1,jjp1,llm,nQ)
123
124	c champs cumulés
125	real ps_cum(iip1,jjp1)
126	real masse_cum(iip1,jjp1,llm)
127	real flux_u_cum(iip1,jjp1,llm)
128	real flux_v_cum(iip1,jjm,llm)
129	real flux_w_cum(iip1,jjp1,llm)
130	real Q_cum(iip1,jjp1,llm,nQ)
131	real flux_uQ_cum(iip1,jjp1,llm,nQ)
132	real flux_vQ_cum(iip1,jjm,llm,nQ)
133	real flux_wQ_cum(iip1,jjp1,llm,nQ)
134	real dQ(iip1,jjp1,llm,nQ)
135
136	save ps_cum,masse_cum,flux_u_cum,flux_v_cum
137	save Q_cum,flux_uQ_cum,flux_vQ_cum
138	save flux_w_cum,flux_wQ_cum
139
140	c champs de transport en moyenne zonale
141	integer ntr,itr
142	parameter (ntr=5)
143
144	character*10,save :: znom(ntr,nQ)
145	character*20,save :: znoml(ntr,nQ)
146	character*10,save :: zunites(ntr,nQ)
147	character*10,save :: znom2(ntr,nQ)
148	character*20,save :: znom2l(ntr,nQ)
149	character*10,save :: zunites2(ntr,nQ)
150	character*10,save :: znom3(nQ)
151	character*20,save :: znom3l(nQ)
152	character*10,save :: zunites3(nQ)
153
154	integer iave,itot,immc,itrs,istn
155	data iave,itot,immc,itrs,istn/1,2,3,4,5/
156	character*3 ctrs(ntr)
157	data ctrs/' ','TOT','MMC','TRS','STN'/
158
159	real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
160	real zwQ(jjm,llm,ntr,nQ),zwQtmp(jjm,llm)
161	real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
162	real zawQ(jjm,llm,ntr,nQ)
163	real zdQ(jjm,llm,nQ)
164	real zmasse(jjm,llm),zavmasse(jjm),zawmasse(llm)
165	real psbar(jjm)
166
167	real zv(jjm,llm),zw(jjp1,llm),psi(jjm,llm+1)
168
169	c TENDANCES POUR MOMENT CINETIQUE
170	c -------------------------------
171
172	integer ntdc,itdc
173	parameter (ntdc=4)
174
175	integer itdcdyn,itdcdis,itdcspg,itdcphy
176	data itdcdyn,itdcdis,itdcspg,itdcphy/1,2,3,4/
177
178	character*6,save :: nomtdc(ntdc)
179
180	c champ contenant les tendances du moment cinetique.
181	real tdc(iip1,jjp1,llm,ntdc)
182	real ztdc(jjm,llm,ntdc) ! moyenne zonale
183
184	c champs cumulés
185	real tdc_cum(iip1,jjp1,llm,ntdc)
186	save tdc_cum
187
188	c integrations completes
189	real mtot,mctot,dmctot(ntdc)
190
191	c Initialisation du fichier contenant les moyennes zonales.
192	c ---------------------------------------------------------
193
194	character*10 infile
195
196	integer fileid
197	integer thoriid, zvertiid
198	save fileid
199
200	integer ndex3d(jjm*llm)
201	real ztmp3d(jjm,llm)
202
203	C Variables locales
204	C
205	integer tau0
206	real zjulian
207	integer zan, dayref
208	C
209	real rlong(jjm),rlatg(jjm)
210
211
212
213	c=====================================================================
214	c Initialisation
215	c=====================================================================
216
217	ndex3d=0
218
219	if (first) then
220
221	if (planet_type.eq."venus") then
222	fact_geovenus = -1.
223	else
224	fact_geovenus = 1.
225	endif
226
227	icum=0
228	c initialisation des fichiers
229	first=.false.
230	c ncum est la frequence de stokage en pas de temps
231	ncum=dt_cum/dt_app
232	if (abs(ncumdt_app-dt_cum).gt.1.e-2dt_app) then
233	if (abs((ncum+1)dt_app-dt_cum).lt.1.e-2dt_app) then
234	ncum = ncum+1
235	elseif (abs((ncum-1)dt_app-dt_cum).lt.1.e-2dt_app) then
236	ncum = ncum-1
237	else
238	WRITE(lunout,*)
239	. 'Pb : le pas de cumule doit etre multiple du pas'
240	WRITE(lunout,*)'dt_app=',dt_app
241	WRITE(lunout,*)'dt_cum=',dt_cum
242	WRITE(lunout,)'ncumdt_app=',ncum*dt_app
243	WRITE(lunout,*)'ncum=',ncum
244	stop
245	endif
246	endif
247
248	c VARIABLES ADVECTEES:
249
250	nom(itemp)='temp'
251	nom(igeop)='gz'
252	nom(iecin)='ecin'
253	nom(iang)='ang'
254	nom(iu)='u'
255	nom(iun)='un'
256
257	unites(itemp)='K'
258	unites(igeop)='m2/s2'
259	unites(iecin)='m2/s2'
260	unites(iang)='ang'
261	unites(iu)='m/s'
262	unites(iun)='un'
263
264	c avec tous les composes, ca fait trop.... Je les enleve
265	c do num=1,ntrac
266	c write(strd2,'(i2.2)') num
267	c nom(6+num)='trac'//strd2
268	c unites(6+num)='kg/kg'
269	c enddo
270
271	c TENDANCES MOMENT CIN:
272
273	nomtdc(itdcdyn) ='dmcdyn'
274	nomtdc(itdcdis) ='dmcdis'
275	nomtdc(itdcspg) ='dmcspg'
276	nomtdc(itdcphy) ='dmcphy'
277
278	c Initialisation du fichier contenant les moyennes zonales.
279	c ---------------------------------------------------------
280
281	infile='dynzon'
282
283	zan = annee_ref
284	dayref = day_ref
285	CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
286	c tau0 = itau_dyn
287	tau0 = 0
288	itau = tau0
289
290	rlong=0.
291	rlatg=rlatv180./pifact_geovenus
292
293	call histbeg(infile, 1, rlong, jjm, rlatg,
294	. 1, 1, 1, jjm,
295	. tau0, zjulian, dt_cum, thoriid, fileid)
296
297	C
298	C Appel a histvert pour la grille verticale
299	C
300	call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
301	. llm, presnivs, zvertiid)
302	C
303	C Appels a histdef pour la definition des variables a sauvegarder
304
305	do iQ=1,nQ
306	do itr=1,ntr
307	if(itr.eq.1) then
308	znom(itr,iQ) =nom(iQ)
309	znoml(itr,iQ) =nom(iQ)
310	zunites(itr,iQ) =unites(iQ)
311	else
312	znom(itr,iQ) =ctrs(itr)//'v'//nom(iQ)
313	znoml(itr,iQ) ='transport : v * '//nom(iQ)//' '//ctrs(itr)
314	zunites(itr,iQ) ='m/s * '//unites(iQ)
315	znom2(itr,iQ) =ctrs(itr)//'w'//nom(iQ)
316	znom2l(itr,iQ) ='transport: w * '//nom(iQ)//' '//ctrs(itr)
317	zunites2(itr,iQ)='Pa/s * '//unites(iQ)
318	endif
319	enddo
320	znom3(iQ)='d'//nom(iQ)
321	znom3l(iQ)='convergence: '//nom(iQ)
322	zunites3(iQ)=unites(iQ)//' /s'
323	c print*,'DEBUG:',znom3(iQ),znom3l(iQ),zunites3(iQ)
324	enddo
325
326	c Declarations des champs avec dimension verticale
327
328	if (1.eq.0) then ! on les sort, ou pas...
329
330	c do iQ=1,nQ
331	c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
332	do iQ=1,4,3
333	do itr=1,ntr
334	IF (prt_level > 5)
335	. WRITE(lunout,*)'var ',itr,iQ
336	. ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
337	call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
338	. zunites(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
339	. 32,'ave(X)',dt_cum,dt_cum)
340	enddo
341	c transport vertical:
342	do itr=2,ntr
343	IF (prt_level > 5)
344	. WRITE(lunout,*)'var ',itr,iQ
345	. ,znom2(itr,iQ),znom2l(itr,iQ),zunites2(itr,iQ)
346	call histdef(fileid,znom2(itr,iQ),znom2l(itr,iQ),
347	. zunites2(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
348	. 32,'ins(X)',dt_cum,dt_cum)
349	enddo
350
351	c Declarations pour convergences
352	IF (prt_level > 5)
353	. WRITE(lunout,*)'var ',iQ
354	. ,znom3(iQ),znom3l(iQ),zunites3(iQ)
355	call histdef(fileid,znom3(iQ),znom3l(iQ),
356	. zunites3(iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
357	. 32,'ins(X)',dt_cum,dt_cum)
358
359	c Declarations pour les fonctions de courant
360	c Non sorties ici...
361	c call histdef(fileid,'psi'//nom(iQ)
362	c . ,'stream fn. '//znoml(itot,iQ),
363	c . zunites(itot,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
364	c . 32,'ave(X)',dt_cum,dt_cum)
365
366	enddo ! iQ
367
368	endif ! 1=1 sortie ou non...
369
370	c Declarations pour les champs de transport d'air
371	call histdef(fileid, 'masse', 'masse',
372	. 'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid,
373	. 32, 'ave(X)', dt_cum, dt_cum)
374	call histdef(fileid, 'v', 'v',
375	. 'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
376	. 32, 'ave(X)', dt_cum, dt_cum)
377	call histdef(fileid, 'w', 'w',
378	. 'Pa/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
379	. 32, 'ins(X)', dt_cum, dt_cum)
380
381	c Declarations pour la fonction de courant
382	call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
383	. 1,jjm,thoriid,llm,1,llm,zvertiid,
384	. 32,'ave(X)',dt_cum,dt_cum)
385
386
387	c Declarations pour les tendances de moment cinetique
388	do itdc=1,ntdc
389	call histdef(fileid, nomtdc(itdc), nomtdc(itdc),
390	. 'ang/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
391	. 32, 'ins(X)', dt_cum, dt_cum)
392	enddo
393
394	c Declaration des champs 1D de transport en latitude
395	c do iQ=1,nQ
396	c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
397	do iQ=1,4,3
398	do itr=2,ntr
399	call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
400	. zunites(itr,iQ),1,jjm,thoriid,1,1,1,-99,
401	. 32,'ave(X)',dt_cum,dt_cum)
402	c JE VIRE LE VERTICAL POUR L'INSTANT
403	c call histdef(fileid,'a'//znom2(itr,iQ),znom2l(itr,iQ),
404	c . zunites2(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
405	c . 32,'ins(X)',dt_cum,dt_cum)
406	enddo
407	enddo
408
409	CALL histend(fileid)
410
411
412	endif ! first
413
414
415	c=====================================================================
416	c Calcul des champs dynamiques
417	c ----------------------------
418
419	c énergie cinétique
420	ucont(:,:,:)=0
421	CALL covcont(llm,ucov,vcov,ucont,vcont)
422	CALL enercin(vcov,ucov,vcont,ucont,ecin)
423
424	c moment cinétique et tendances
425	dudyn = 0.
426	dudis = 0.
427	duspg = 0.
428	duphy = 0.
429	do l=1,llm
430	unat(:,:,l)=ucont(:,:,l)*cu(:,:)
431	dudyn(:,2:jjm,l) = ducovdyn(:,2:jjm,l)/cu(:,2:jjm)
432	dudis(:,2:jjm,l) = ducovdis(:,2:jjm,l)/cu(:,2:jjm)
433	duspg(:,2:jjm,l) = ducovspg(:,2:jjm,l)/cu(:,2:jjm)
434	duphy(:,2:jjm,l) = ducovphy(:,2:jjm,l)/cu(:,2:jjm)
435	do j=1,jjp1
436	ang(:,j,l)= radcos(rlatu(j))
437	. ( unat(:,j,l) + radcos(rlatu(j))omeg )
438	tdc(:,j,l,1) = radcos(rlatu(j))dudyn(:,j,l)
439	tdc(:,j,l,2) = radcos(rlatu(j))dudis(:,j,l)
440	tdc(:,j,l,3) = radcos(rlatu(j))duspg(:,j,l)
441	tdc(:,j,l,4) = radcos(rlatu(j))duphy(:,j,l)
442	enddo
443	enddo
444	c Normalisation:
445	ang = ang / (2./3. radrad*omeg)
446	do itdc=1,ntdc
447	tdc(:,:,:,itdc)=tdc(:,:,:,itdc) / (2./3. radrad*omeg)
448	enddo
449
450	! ADAPTATION GCM POUR CP(T)
451	call tpot2t(ip1jmp1*llm,teta,temp,pk)
452	Q(:,:,:,itemp) = temp(:,:,:)
453	Q(:,:,:,igeop) =phi(:,:,:)
454	Q(:,:,:,iecin) =ecin(:,:,:)
455	Q(:,:,:,iang) =ang(:,:,:)
456	Q(:,:,:,iu) =unat(:,:,:)
457	Q(:,:,:,iun) =1.
458	c avec tous les composes, ca fait trop.... Je les enleve
459	c do num=1,ntrac
460	c Q(:,:,:,6+num)=trac(:,:,:,num)
461	c enddo
462
463	c calcul du flux de masse vertical (+ vers le bas)
464	call convmas(flux_u,flux_v,convm)
465	CALL vitvert(convm,w)
466
467	c=====================================================================
468	c Cumul
469	c=====================================================================
470	c
471	if(icum.EQ.0) then
472	ps_cum = 0.
473	masse_cum = 0.
474	flux_u_cum = 0.
475	flux_v_cum = 0.
476	flux_w_cum = 0.
477	Q_cum = 0.
478	flux_vQ_cum = 0.
479	flux_uQ_cum = 0.
480	flux_wQ_cum = 0.
481	tdc_cum = 0.
482	endif
483
484	IF (prt_level > 5)
485	. WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
486	icum=icum+1
487
488	c accumulation des flux de masse horizontaux
489	ps_cum = ps_cum + ps
490	masse_cum = masse_cum + masse
491	flux_u_cum = flux_u_cum + flux_u
492	flux_v_cum = flux_v_cum + flux_v
493	flux_w_cum = flux_w_cum + w
494	do iQ=1,nQ
495	Q_cum(:,:,:,iQ) = Q_cum(:,:,:,iQ) + Q(:,:,:,iQ)*masse(:,:,:)
496	enddo
497	do itdc=1,ntdc
498	tdc_cum(:,:,:,itdc) =
499	. tdc_cum(:,:,:,itdc) + tdc(:,:,:,itdc)*masse(:,:,:)
500	enddo
501
502	c=====================================================================
503	c FLUX ET TENDANCES
504	c=====================================================================
505
506	c Flux longitudinal
507	c -----------------
508	do iQ=1,nQ
509	do l=1,llm
510	do j=1,jjp1
511	do i=1,iim
512	flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
513	s +flux_u(i,j,l)0.5(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
514	enddo
515	flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
516	enddo
517	enddo
518	enddo
519
520	c flux méridien
521	c -------------
522	do iQ=1,nQ
523	do l=1,llm
524	do j=1,jjm
525	do i=1,iip1
526	flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
527	s +flux_v(i,j,l)0.5(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
528	enddo
529	enddo
530	enddo
531	enddo
532
533	c Flux vertical
534	c -------------
535	do iQ=1,nQ
536	do l=2,llm
537	do j=1,jjp1
538	do i=1,iip1
539	flux_wQ_cum(i,j,l,iQ)=flux_wQ_cum(i,j,l,iQ)
540	s +w(i,j,l)0.5(Q(i,j,l-1,iQ)+Q(i,j,l,iQ))
541	enddo
542	enddo
543	enddo
544	flux_wQ_cum(:,:,1,iQ)=0.0e0
545	enddo
546
547	c tendances
548	c ---------
549
550	c convergence horizontale
551	call convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
552
553	c calcul de la vitesse verticale
554	call convmas(flux_u_cum,flux_v_cum,convm)
555	CALL vitvert(convm,w)
556
557	c ajustement tendances (vitesse verticale)
558	do iQ=1,nQ
559	do l=1,llm-1
560	do j=1,jjp1
561	do i=1,iip1
562	ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
563	dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
564	dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
565	enddo
566	enddo
567	enddo
568	enddo
569	IF (prt_level > 5)
570	. WRITE(lunout,*)'Apres les calculs fait a chaque pas'
571
572	c=====================================================================
573	c PAS DE TEMPS D'ECRITURE
574	c=====================================================================
575	if (icum.eq.ncum) then
576	c=====================================================================
577
578	time=time+dt_cum
579	itau=itau+1
580
581	IF (prt_level > 5)
582	. WRITE(lunout,*)'Pas d ecriture'
583
584	c Normalisation
585	do iQ=1,nQ
586	Q_cum(:,:,:,iQ) = Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
587	dQ(:,:,:,iQ) = dQ(:,:,:,iQ) /masse_cum(:,:,:)
588	enddo
589	do itdc=1,ntdc
590	tdc_cum(:,:,:,itdc) = tdc_cum(:,:,:,itdc)/masse_cum(:,:,:)
591	enddo
592
593	zz=1./REAL(ncum)
594	ps_cum = ps_cum *zz
595	masse_cum = masse_cum *zz
596	flux_u_cum = flux_u_cum *zz
597	flux_v_cum = flux_v_cum *zz
598	flux_w_cum = flux_w_cum *zz
599	flux_uQ_cum = flux_uQ_cum *zz
600	flux_vQ_cum = flux_vQ_cum *zz
601	flux_wQ_cum = flux_wQ_cum *zz
602
603	c Integration complete
604	mtot = 0.
605	mctot = 0.
606	dmctot = 0.
607	do l=1,llm
608	do j=2,jjm
609	do i=1,iim
610	mtot = mtot + masse_cum(i,j,l)
611	mctot = mctot + Q_cum(i,j,l,iang)*masse_cum(i,j,l)
612	enddo
613	enddo
614	enddo
615	mctot = mctot/mtot
616	do itdc=1,ntdc
617	do l=1,llm
618	do j=2,jjm
619	do i=1,iim
620	dmctot(itdc) = dmctot(itdc)
621	. + tdc_cum(i,j,l,itdc)*masse_cum(i,j,l)/mtot
622	enddo
623	enddo
624	enddo
625	enddo
626
627	c=====================================================================
628	c Transport méridien
629	c=====================================================================
630
631	c cumul zonal des masses des mailles
632	c ----------------------------------
633	zv=0.
634	zw=0.
635	zmasse=0.
636	call massbar(masse_cum,massebx,masseby)
637
638	c moy zonale de la ps cumulee
639	do j=1,jjm
640	psbar(j)=0.
641	do i=1,iim
642	psbar(j)=psbar(j)+ps_cum(i,j)/iim
643	enddo
644	enddo
645
646	do l=1,llm
647	do j=1,jjm
648	do i=1,iim
649	zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
650	zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
651	zw(j,l)=zw(j,l)+flux_w_cum(i,j,l)
652	enddo
653	zfactv(j,l)=cv(1,j)/zmasse(j,l)
654	zfactw(j,l)=((ap(l)-ap(l+1))+psbar(j)*(bp(l)-bp(l+1)))
655	. /zmasse(j,l)
656	enddo
657	do i=1,iim
658	zw(jjp1,l)=zw(jjp1,l)+flux_w_cum(i,jjp1,l)
659	enddo
660	enddo
661
662	c print*,'3OK'
663	c --------------------------------------------------------------
664	c calcul de la moyenne zonale du transport :
665	c ------------------------------------------
666	c
667	c --
668	c TOT : la circulation totale [ vq ]
669	c
670	c - -
671	c MMC : mean meridional circulation [ v ] [ q ]
672	c
673	c ---- -- - -
674	c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
675	c
676	c - * - * - - - -
677	c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
678	c
679	c - -
680	c on utilise aussi l'intermediaire TMP : [ v q ]
681	c
682	c la variable zfactv transforme un transport meridien cumule
683	c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
684	c la variable zfactw transforme un transport vertical cumule
685	c en kg/s * unte-du-champ-transporte en Pa/s * unite-du-champ-transporte
686	c
687	c --------------------------------------------------------------
688
689
690	c ----------------------------------------
691	c Transport dans le plan latitude-altitude
692	c ----------------------------------------
693
694	zvQ=0.
695	zwQ=0.
696	zdQ=0.
697	psiQ=0.
698
699	do iQ=1,nQ
700
701	c transport meridien
702	zvQtmp=0.
703	do l=1,llm
704	do j=1,jjm
705	c print*,'j,l,iQ=',j,l,iQ
706	c Calcul des moyennes zonales du transport total et de zvQtmp
707	do i=1,iim
708	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
709	s +flux_vQ_cum(i,j,l,iQ)
710	zqy= 0.5(Q_cum(i,j,l,iQ)masse_cum(i,j,l)+
711	s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
712	zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
713	s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
714	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
715	enddo
716	c print*,'aOK'
717	c Decomposition
718	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
719	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
720	zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
721	zvQ(j,l,immc,iQ)=zv(j,l)zvQ(j,l,iave,iQ)zfactv(j,l)
722	zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
723	zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
724	enddo
725	enddo
726	c fonction de courant meridienne pour la quantite Q
727	do l=llm,1,-1
728	do j=1,jjm
729	psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)/zfactv(j,l)
730	enddo
731	enddo
732	!! enddo
733
734	c transport vertical
735	zwQtmp=0.
736	do l=1,llm
737	do j=1,jjm
738	c print*,'j,l,iQ=',j,l,iQ
739	c Calcul des moyennes zonales du transport vertical total et de zwQtmp
740	do i=1,iim
741	zwQ(j,l,itot,iQ)=zwQ(j,l,itot,iQ)
742	s +flux_wQ_cum(i,j,l,iQ)
743	zqy= 0.5(Q_cum(i,j,l,iQ)masse_cum(i,j,l)+
744	s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
745	zwQtmp(j,l)=zwQtmp(j,l)+flux_w_cum(i,j,l)*zqy
746	s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
747	zwQ(j,l,iave,iQ)=zwQ(j,l,iave,iQ)+zqy
748	enddo
749	c Decomposition
750	zwQ(j,l,iave,iQ)=zwQ(j,l,iave,iQ)/zmasse(j,l)
751	zwQ(j,l,itot,iQ)=zwQ(j,l,itot,iQ)*zfactw(j,l)
752	zwQtmp(j,l)=zwQtmp(j,l)*zfactw(j,l)
753	zwQ(j,l,immc,iQ)=zw(j,l)zwQ(j,l,iave,iQ)zfactw(j,l)
754	zwQ(j,l,itrs,iQ)=zwQ(j,l,itot,iQ)-zwQtmp(j,l)
755	zwQ(j,l,istn,iQ)=zwQtmp(j,l)-zwQ(j,l,immc,iQ)
756	enddo
757	enddo
758
759	c convergence
760	c Calcul moyenne zonale de la convergence totale
761	do l=1,llm
762	do j=1,jjm
763	c print*,'j,l,iQ=',j,l,iQ
764	do i=1,iim
765	zdQ(j,l,iQ)=zdQ(j,l,iQ) +
766	. ( dQ(i,j,l,iQ) * masse_cum(i,j,l)
767	. + dQ(i,j+1,l,iQ) * masse_cum(i,j+1,l))
768	. / ( masse_cum(i,j,l)+masse_cum(i,j+1,l))
769	enddo
770	enddo
771	enddo
772	enddo ! of do iQ=1,nQ
773
774	c fonction de courant pour la circulation meridienne moyenne
775	psi=0.
776	do l=llm,1,-1
777	do j=1,jjm
778	psi(j,l)= psi(j,l+1)+zv(j,l)
779	zv(j,l) = zv(j,l)*zfactv(j,l)
780	zw(j,l) = 0.5(zw(j,l)+zw(j+1,l))zfactw(j,l)
781	enddo
782	enddo
783
784	c Calcul moyenne zonale des tendances moment cin.
785	ztdc=0.
786	do itdc=1,ntdc
787	do l=1,llm
788	do j=1,jjm
789	do i=1,iim
790	ztdc(j,l,itdc)=ztdc(j,l,itdc) +
791	. ( tdc_cum(i,j,l,itdc) * masse_cum(i,j,l)
792	. + tdc_cum(i,j+1,l,itdc) * masse_cum(i,j+1,l))
793	. / ( masse_cum(i,j,l)+masse_cum(i,j+1,l))
794	enddo
795	enddo
796	enddo
797	enddo
798
799	c print*,'4OK'
800
801	c--------------------------------------
802	c--------------------------------------
803	c sorties proprement dites
804	c--------------------------------------
805	c--------------------------------------
806
807	if (i_sortie.eq.1) then
808
809	c sortie des integrations completes dans le listing
810	write(*,'(A12,5(1PE11.4,X))') "BILANMCDYN ",mctot,dmctot
811
812	c sorties dans fichier dynzon
813
814	if (1.eq.0) then ! on les sort, ou pas...
815
816	c avec tous les composes, ca fait trop.... Je les enleve
817	c do iQ=1,nQ
818	c do iQ=1,6
819	c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
820	do iQ=1,4,3
821
822	ztmp3d(:,:)= zvQ(:,:,1,iQ) ! valeur moyenne
823	call histwrite(fileid,znom(1,iQ),itau,ztmp3d
824	s ,jjm*llm,ndex3d)
825	do itr=2,ntr
826	ztmp3d(:,:)= zvQ(:,:,itr,iQ)*fact_geovenus ! transport horizontal
827	call histwrite(fileid,znom(itr,iQ),itau,ztmp3d
828	s ,jjm*llm,ndex3d)
829	enddo
830
831	do itr=2,ntr
832	ztmp3d(:,:)=zwQ(:,:,itr,iQ)
833	call histwrite(fileid,znom2(itr,iQ),itau,ztmp3d
834	s ,jjm*llm,ndex3d)
835	enddo
836
837	ztmp3d(:,:)= zdQ(:,:,iQ)
838	call histwrite(fileid,znom3(iQ),itau,ztmp3d
839	s ,jjm*llm,ndex3d)
840
841	c ztmp3d(:,:)= psiQ(:,1:llm,iQ)*fact_geovenus
842	c call histwrite(fileid,'psi'//nom(iQ),itau,ztmp3d
843	c s ,jjm*llm,ndex3d)
844	enddo
845
846	endif ! 1=1 sortie ou non...
847
848	ztmp3d=zmasse
849	call histwrite(fileid,'masse',itau,ztmp3d
850	s ,jjm*llm,ndex3d)
851
852	ztmp3d= zv*fact_geovenus
853	call histwrite(fileid,'v',itau,ztmp3d
854	s ,jjm*llm,ndex3d)
855	ztmp3d(:,:)=zw(1:jjm,:)
856	call histwrite(fileid,'w',itau,ztmp3d
857	s ,jjm*llm,ndex3d)
858	ztmp3d= psi(:,1:llm)1.e-9fact_geovenus
859	call histwrite(fileid,'psi',itau,ztmp3d,jjm*llm,ndex3d)
860
861	do itdc=1,ntdc
862	ztmp3d(:,:)= ztdc(:,:,itdc)
863	call histwrite(fileid,nomtdc(itdc),itau,ztmp3d
864	s ,jjm*llm,ndex3d)
865	enddo
866
867	endif ! i_sortie
868
869
870	c -----------------
871	c Moyenne verticale
872	c -----------------
873
874	zavmasse=0.
875	do l=1,llm
876	zavmasse(:)=zavmasse(:)+zmasse(:,l)
877	enddo
878	zavQ=0.
879
880	c avec tous les composes, ca fait trop.... Je les enleve
881	c do iQ=1,nQ
882	c do iQ=1,6
883	c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
884	do iQ=1,4,3
885	do itr=2,ntr
886	do l=1,llm
887	zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
888	enddo
889	zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zavmasse(:)
890	if (i_sortie.eq.1) then
891	ztmp3d=0.0
892	ztmp3d(:,1)= zavQ(:,itr,iQ)*fact_geovenus
893	call histwrite(fileid,'a'//znom(itr,iQ),itau,ztmp3d
894	. ,jjm*llm,ndex3d)
895	endif
896	enddo
897	enddo
898
899	c ------------------
900	c Moyenne meridienne
901	c ------------------
902
903	zawmasse=0.
904	do j=1,jjm
905	do l=1,llm
906	zawmasse(l)=zawmasse(l)+zmasse(j,l)
907	enddo
908	enddo
909	zawQ=0.
910
911	c avec tous les composes, ca fait trop.... Je les enleve
912	c do iQ=1,nQ
913	c do iQ=1,6
914	c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
915	do iQ=1,4,3
916	do itr=2,ntr
917	do l=1,llm
918	do j=1,jjm
919	zawQ(1,l,itr,iQ)=zawQ(1,l,itr,iQ)+zwQ(j,l,itr,iQ)*zmasse(j,l)
920	enddo
921	zawQ(1,l,itr,iQ)=zawQ(1,l,itr,iQ)/zawmasse(l)
922	enddo
923	if (i_sortie.eq.1) then
924	ztmp3d=0.0
925	do l=1,llm
926	ztmp3d(1,l)=zawQ(1,l,itr,iQ)
927	enddo
928	c JE VIRE LE VERTICAL POUR L'INSTANT
929	c call histwrite(fileid,'a'//znom2(itr,iQ),itau,ztmp3d
930	c . ,jjm*llm,ndex3d)
931	endif
932	enddo
933	enddo
934
935	call histsync(fileid)
936
937	c=====================================================================
938	c/////////////////////////////////////////////////////////////////////
939	icum=0 !///////////////////////////////////////
940	endif ! icum.eq.ncum !///////////////////////////////////////
941	c/////////////////////////////////////////////////////////////////////
942	c=====================================================================
943
944	return
945	end

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

Download in other formats:

Original Format