Context Navigation

bilan_dyn.F @ 1035

Last change on this file since 1035 was 524, checked in by lmdzadmin, 20 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 15.8 KB

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

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

Download in other formats:

Original Format