Context Navigation

bilan_dyn_p.F @ 1352

Last change on this file since 1352 was 1279, checked in by Laurent Fairhead, 15 years ago
Merged LMDZ4-dev branch changes r1241:1278 into the trunk Running trunk and LMDZ4-dev in LMDZOR configuration on local machine (sequential) and SX8 (4-proc) yields identical results (restart and restartphy are identical binarily) Log history from r1241 to r1278 is available by switching to source:LMDZ4/branches/LMDZ4-dev-20091210
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 19.9 KB

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

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

Download in other formats:

Original Format