Context Navigation

bilan_dyn.F @ 541

Last change on this file since 541 was 541, checked in by lmdzadmin, 20 years ago

Convergence avec la version d'Olivia Coindreau incluant:

le offline
les thermiques
mellor & yamada dans la couche limite

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

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

Download in other formats:

Original Format