Context Navigation

bilan_dyn.F @ 3356

Last change on this file since 3356 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

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

File size: 16.1 KB

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

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

Download in other formats:

Original Format