Context Navigation

bilan_dyn.F @ 431

Last change on this file since 431 was 431, checked in by lmdzadmin, 21 years ago

Convergence avec la version de Ionela Musat

bilan_dyn.F : correction bug dimension (PF)
conf_gcm.F : remize a zero de la date initiale (LF)
control.h : remize a zero de la date initiale (LF)
diagedyn.F : correction bug dimension (IM)
etat0_netcdf.F : lecture de co2_ppm & solaire ds physiq.def
gcm.F : remize a zero de la date initiale (LF)

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 15.7 KB

Line
1	SUBROUTINE bilan_dyn (ntrac,dt_app,dt_cum,
2	s ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
3
4	c AFAIRE
5	c Prevoir en champ nq+1 le diagnostique de l'energie
6	c en faisant Qzon=Cv T + L * ...
7	c vQ..A=Cp T + L * ...
8
9	USE IOIPSL
10
11	IMPLICIT NONE
12
13	#include "dimensions.h"
14	#include "paramet.h"
15	#include "comconst.h"
16	#include "comvert.h"
17	#include "comgeom2.h"
18
19	c====================================================================
20	c
21	c Sous-programme consacre à des diagnostics dynamiques de base
22	c
23	c
24	c De facon generale, les moyennes des scalaires Q sont ponderees par
25	c la masse.
26	c
27	c Les flux de masse sont eux simplement moyennes.
28	c
29	c====================================================================
30
31	c Arguments :
32	c ===========
33
34	integer ntrac
35	real dt_app,dt_cum
36	real ps(iip1,jjp1)
37	real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
38	real flux_u(iip1,jjp1,llm)
39	real flux_v(iip1,jjm,llm)
40	real teta(iip1,jjp1,llm)
41	real phi(iip1,jjp1,llm)
42	real ucov(iip1,jjp1,llm)
43	real vcov(iip1,jjm,llm)
44	real trac(iip1,jjp1,llm,ntrac)
45
46	c Local :
47	c =======
48
49	integer icum,ncum
50	logical first
51	real zz,zqy,zfactv(jjm,llm)
52
53	integer nQ
54	parameter (nQ=7)
55
56
57	character*6 nom(nQ)
58	character*6 unites(nQ)
59	character*10 file
60	integer ifile
61	parameter (ifile=4)
62
63	integer itemp,igeop,iecin,iang,iu,iovap,iun
64	integer i_sortie
65
66	save first,icum,ncum
67	save itemp,igeop,iecin,iang,iu,iovap,iun
68	save i_sortie
69
70	real time
71	integer itau
72	save time,itau
73	data time,itau/0.,0/
74
75	data first/.true./
76	data itemp,igeop,iecin,iang,iu,iovap,iun/1,2,3,4,5,6,7/
77	data i_sortie/1/
78
79	real ww
80
81	c variables dynamiques intermédiaires
82	REAL vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
83	REAL ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
84	REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
85	REAL vorpot(iip1,jjm,llm)
86	REAL w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm)
87	REAL bern(iip1,jjp1,llm)
88
89	c champ contenant les scalaires advectés.
90	real Q(iip1,jjp1,llm,nQ)
91
92	c champs cumulés
93	real ps_cum(iip1,jjp1)
94	real masse_cum(iip1,jjp1,llm)
95	real flux_u_cum(iip1,jjp1,llm)
96	real flux_v_cum(iip1,jjm,llm)
97	real Q_cum(iip1,jjp1,llm,nQ)
98	real flux_uQ_cum(iip1,jjp1,llm,nQ)
99	real flux_vQ_cum(iip1,jjm,llm,nQ)
100	real flux_wQ_cum(iip1,jjp1,llm,nQ)
101	real dQ(iip1,jjp1,llm,nQ)
102
103	save ps_cum,masse_cum,flux_u_cum,flux_v_cum
104	save Q_cum,flux_uQ_cum,flux_vQ_cum
105
106	c champs de tansport en moyenne zonale
107	integer ntr,itr
108	parameter (ntr=5)
109
110	character*10 znom(ntr,nQ)
111	character*20 znoml(ntr,nQ)
112	character*10 zunites(ntr,nQ)
113	integer iave,itot,immc,itrs,istn
114	data iave,itot,immc,itrs,istn/1,2,3,4,5/
115	character*3 ctrs(ntr)
116	data ctrs/' ','TOT','MMC','TRS','STN'/
117
118	real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
119	real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
120	real zmasse(jjm,llm),zamasse(jjm)
121
122	real zv(jjm,llm),psi(jjm,llm+1)
123
124	integer i,j,l,iQ
125
126
127	c Initialisation du fichier contenant les moyennes zonales.
128	c ---------------------------------------------------------
129
130	character*10 infile
131
132	integer*4 day0, anne0
133	integer fileid
134	integer thoriid, zvertiid
135	save fileid
136
137	integer ndex3d(jjm*llm)
138
139	C Variables locales
140	C
141	integer tau0
142	real zjulian
143	character*3 str
144	character*10 ctrac
145	integer ii,jj
146	integer zan
147	C
148	real rlong(jjm),rlatg(jjm)
149
150
151
152	c=====================================================================
153	c Initialisation
154	c=====================================================================
155
156	time=time+dt_app
157	itau=itau+1
158
159	if (first) then
160
161
162	icum=0
163	c initialisation des fichiers
164	first=.false.
165	c ncum est la frequence de stokage en pas de temps
166	ncum=dt_cum/dt_app
167	if (abs(ncumdt_app-dt_cum).gt.1.e-5dt_app) then
168	print*,'Pb : le pas de cumule doit etre multiple du pas'
169	print*,'dt_app=',dt_app
170	print*,'dt_cum=',dt_cum
171	stop
172	endif
173
174	if (i_sortie.eq.1) then
175	file='dynzon'
176	call inigrads(ifile,1
177	s ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi
178	s ,llm,presnivs,1.
179	s ,dt_cum,file,'dyn_zon ')
180	endif
181
182	nom(itemp)='T'
183	nom(igeop)='gz'
184	nom(iecin)='K'
185	nom(iang)='ang'
186	nom(iu)='u'
187	nom(iovap)='ovap'
188	nom(iun)='un'
189
190	unites(itemp)='K'
191	unites(igeop)='m2/s2'
192	unites(iecin)='m2/s2'
193	unites(iang)='ang'
194	unites(iu)='m/s'
195	unites(iovap)='kg/kg'
196	unites(iun)='un'
197
198
199	c Initialisation du fichier contenant les moyennes zonales.
200	c ---------------------------------------------------------
201
202	infile='dynzon'
203
204	day0=0
205	zan = 0.
206
207	CALL ymds2ju(zan, 1, 1, 0.0, zjulian)
208	zjulian = zjulian + day0
209	tau0 = 0
210
211	rlong=0.
212	rlatg=rlatv*180./pi
213
214	call histbeg(infile, 1, rlong, jjm, rlatg,
215	. 1, 1, 1, jjm,
216	. tau0, zjulian, dt_cum, thoriid, fileid)
217
218	C
219	C Appel a histvert pour la grille verticale
220	C
221	call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
222	. llm, presnivs, zvertiid)
223	C
224	C Appels a histdef pour la definition des variables a sauvegarder
225	do iQ=1,nQ
226	do itr=1,ntr
227	if(itr.eq.1) then
228	znom(itr,iQ)=nom(iQ)
229	znoml(itr,iQ)=nom(iQ)
230	zunites(itr,iQ)=unites(iQ)
231	else
232	znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
233	znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
234	zunites(itr,iQ)='m/s * '//unites(iQ)
235	endif
236	enddo
237	enddo
238
239	c Declarations des champs avec dimension verticale
240	print*,'1HISTDEF'
241	do iQ=1,nQ
242	do itr=1,ntr
243	print*,'var ',itr,iQ
244	. ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
245	call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
246	. zunites(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
247	. 32,'ave(X)',dt_cum,dt_cum)
248	enddo
249	c Declarations pour les fonctions de courant
250	print*,'2HISTDEF'
251	call histdef(fileid,'psi'//nom(iQ)
252	. ,'stream fn. '//znoml(itot,iQ),
253	. zunites(itot,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
254	. 32,'ave(X)',dt_cum,dt_cum)
255	enddo
256
257
258	c Declarations pour les champs de transport d'air
259	print*,'3HISTDEF'
260	call histdef(fileid, 'masse', 'masse',
261	. 'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid,
262	. 32, 'ave(X)', dt_cum, dt_cum)
263	call histdef(fileid, 'v', 'v',
264	. 'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
265	. 32, 'ave(X)', dt_cum, dt_cum)
266	c Declarations pour les fonctions de courant
267	print*,'4HISTDEF'
268	call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
269	. 1,jjm,thoriid,llm,1,llm,zvertiid,
270	. 32,'ave(X)',dt_cum,dt_cum)
271
272
273	c Declaration des champs 1D de transport en latitude
274	print*,'5HISTDEF'
275	do iQ=1,nQ
276	do itr=2,ntr
277	call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
278	. zunites(itr,iQ),1,jjm,thoriid,1,1,1,-99,
279	. 32,'ave(X)',dt_cum,dt_cum)
280	enddo
281	enddo
282
283
284	print*,'8HISTDEF'
285	CALL histend(fileid)
286
287
288	endif
289
290
291	c=====================================================================
292	c Calcul des champs dynamiques
293	c ----------------------------
294
295	c énergie cinétique
296	CALL covcont(llm,ucov,vcov,ucont,vcont)
297	CALL enercin(vcov,ucov,vcont,ucont,ecin)
298
299	c moment cinétique
300	do l=1,llm
301	ang(:,:,l)=ucov(:,:,l)+constang(:,:)
302	unat(:,:,l)=ucont(:,:,l)*cu(:,:)
303	enddo
304
305	Q(:,:,:,itemp)=teta(:,:,:)*pk(:,:,:)/cpp
306	Q(:,:,:,igeop)=phi(:,:,:)
307	Q(:,:,:,iecin)=ecin(:,:,:)
308	Q(:,:,:,iang)=ang(:,:,:)
309	Q(:,:,:,iu)=unat(:,:,:)
310	Q(:,:,:,iovap)=q(:,:,:,1)
311	Q(:,:,:,iun)=1.
312
313
314	c=====================================================================
315	c Cumul
316	c=====================================================================
317	c
318	if(icum.EQ.0) then
319	ps_cum=0.
320	masse_cum=0.
321	flux_u_cum=0.
322	flux_v_cum=0.
323	Q_cum=0.
324	flux_vQ_cum=0.
325	flux_uQ_cum=0.
326	endif
327
328	print*,'dans bilan_dyn ',icum,'->',icum+1
329	icum=icum+1
330
331	c accumulation des flux de masse horizontaux
332	ps_cum=ps_cum+ps
333	masse_cum=masse_cum+masse
334	flux_u_cum=flux_u_cum+flux_u
335	flux_v_cum=flux_v_cum+flux_v
336	do iQ=1,nQ
337	Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)+Q(:,:,:,iQ)*masse(:,:,:)
338	enddo
339
340	c=====================================================================
341	c FLUX ET TENDANCES
342	c=====================================================================
343
344	c Flux longitudinal
345	c -----------------
346	do iQ=1,nQ
347	do l=1,llm
348	do j=1,jjp1
349	do i=1,iim
350	flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
351	s +flux_u(i,j,l)0.5(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
352	enddo
353	flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
354	enddo
355	enddo
356	enddo
357
358	c flux méridien
359	c -------------
360	do iQ=1,nQ
361	do l=1,llm
362	do j=1,jjm
363	do i=1,iip1
364	flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
365	s +flux_v(i,j,l)0.5(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
366	enddo
367	enddo
368	enddo
369	enddo
370
371
372	c tendances
373	c ---------
374
375	c convergence horizontale
376	call convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
377
378	c calcul de la vitesse verticale
379	call convmas(flux_u_cum,flux_v_cum,convm)
380	CALL vitvert(convm,w)
381
382	do iQ=1,nQ
383	do l=1,llm-1
384	do j=1,jjp1
385	do i=1,iip1
386	ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
387	dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
388	dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
389	enddo
390	enddo
391	enddo
392	enddo
393
394	print*,'Apres les calculs fait a chaque pas'
395	c=====================================================================
396	c PAS DE TEMPS D'ECRITURE
397	c=====================================================================
398	if (icum.eq.ncum) then
399	c=====================================================================
400
401	print*,'Pas d ecriture'
402
403	c Normalisation
404	do iQ=1,nQ
405	Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
406	enddo
407	zz=1./float(ncum)
408	ps_cum=ps_cum*zz
409	masse_cum=masse_cum*zz
410	flux_u_cum=flux_u_cum*zz
411	flux_v_cum=flux_v_cum*zz
412	flux_uQ_cum=flux_uQ_cum*zz
413	flux_vQ_cum=flux_vQ_cum*zz
414	dQ=dQ*zz
415
416	c print*,'1OK'
417
418	c A retravailler eventuellement
419	c division de dQ par la masse pour revenir aux bonnes grandeurs
420	do iQ=1,nQ
421	dQ(:,:,:,iQ)=dQ(:,:,:,iQ)/masse_cum(:,:,:)
422	enddo
423
424	c print*,'2OK'
425	c=====================================================================
426	c Transport méridien
427	c=====================================================================
428
429	c cumul zonal des masses des mailles
430	c ----------------------------------
431	zv=0.
432	zmasse=0.
433	call massbar(masse_cum,massebx,masseby)
434	do l=1,llm
435	do j=1,jjm
436	do i=1,iim
437	zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
438	zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
439	enddo
440	zfactv(j,l)=cv(1,j)/zmasse(j,l)
441	enddo
442	enddo
443
444	c print*,'3OK'
445	c --------------------------------------------------------------
446	c calcul de la moyenne zonale du transport :
447	c ------------------------------------------
448	c
449	c --
450	c TOT : la circulation totale [ vq ]
451	c
452	c - -
453	c MMC : mean meridional circulation [ v ] [ q ]
454	c
455	c ---- -- - -
456	c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
457	c
458	c - * - * - - - -
459	c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
460	c
461	c - -
462	c on utilise aussi l'intermediaire TMP : [ v q ]
463	c
464	c la variable zfactv transforme un transport meridien cumule
465	c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
466	c
467	c --------------------------------------------------------------
468
469
470	c ----------------------------------------
471	c Transport dans le plan latitude-altitude
472	c ----------------------------------------
473
474	zvQ=0.
475	psiQ=0.
476	do iQ=1,nQ
477	zvQtmp=0.
478	do l=1,llm
479	do j=1,jjm
480	c print*,'j,l,iQ=',j,l,iQ
481	c Calcul des moyennes zonales du transort total et de zvQtmp
482	do i=1,iim
483	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
484	s +flux_vQ_cum(i,j,l,iQ)
485	zqy= 0.5(Q_cum(i,j,l,iQ)masse_cum(i,j,l)+
486	s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
487	zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
488	s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
489	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
490	enddo
491	c print*,'aOK'
492	c Decomposition
493	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
494	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
495	zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
496	zvQ(j,l,immc,iQ)=zv(j,l)zvQ(j,l,iave,iQ)zfactv(j,l)
497	zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
498	zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
499	enddo
500	enddo
501	c fonction de courant meridienne pour la quantite Q
502	do l=llm,1,-1
503	do j=1,jjm
504	psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
505	enddo
506	enddo
507	enddo
508
509	c fonction de courant pour la circulation meridienne moyenne
510	psi=0.
511	do l=llm,1,-1
512	do j=1,jjm
513	psi(j,l)=psi(j,l+1)+zv(j,l)
514	zv(j,l)=zv(j,l)*zfactv(j,l)
515	enddo
516	enddo
517
518	c print*,'4OK'
519	c sorties proprement dites
520	if (i_sortie.eq.1) then
521	do iQ=1,nQ
522	do itr=1,ntr
523	call histwrite(fileid,znom(itr,iQ),itau,zvQ(:,:,itr,iQ)
524	s ,jjm*llm,ndex3d)
525	enddo
526	call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ)
527	s ,jjm*llm,ndex3d)
528	enddo
529
530	call histwrite(fileid,'masse',itau,zmasse
531	s ,jjm*llm,ndex3d)
532	call histwrite(fileid,'v',itau,zv
533	s ,jjm*llm,ndex3d)
534	psi=psi*1.e-9
535	call histwrite(fileid,'psi',itau,psi(:,1:llm),jjm*llm,ndex3d)
536
537	endif
538
539
540	c -----------------
541	c Moyenne verticale
542	c -----------------
543
544	zamasse=0.
545	do l=1,llm
546	zamasse(:)=zamasse(:)+zmasse(:,l)
547	enddo
548	zavQ=0.
549	do iQ=1,nQ
550	do itr=2,ntr
551	do l=1,llm
552	zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
553	enddo
554	zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zamasse(:)
555	call histwrite(fileid,'a'//znom(itr,iQ),itau,zavQ(:,itr,iQ)
556	s ,jjm*llm,ndex3d)
557	enddo
558	enddo
559
560	c on doit pouvoir tracer systematiquement la fonction de courant.
561
562	c=====================================================================
563	c/////////////////////////////////////////////////////////////////////
564	icum=0 !///////////////////////////////////////
565	endif ! icum.eq.ncum !///////////////////////////////////////
566	c/////////////////////////////////////////////////////////////////////
567	c=====================================================================
568
569	return
570	end

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

Download in other formats:

Original Format