My Project
 All Classes Files Functions Variables Macros
bilan_dyn_p.F
Go to the documentation of this file.
1 !
2 ! $Id: bilan_dyn_p.F 1550 2011-07-05 09:44:55Z lguez $
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(ncum*dt_app-dt_cum).gt.1.e-5*dt_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  if (pole_sud) jje=jj_end-1
392  flux_v_cum(:,jjb:jje,:)=0.
393  flux_vq_cum(:,jjb:jje,:,:)=0.
394  endif
395 
396  IF (prt_level > 5)
397  . WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
398  icum=icum+1
399 
400 c accumulation des flux de masse horizontaux
401  jjb=jj_begin
402  jje=jj_end
403 
404  ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)+ps(:,jjb:jje)
405  masse_cum(:,jjb:jje,:)=masse_cum(:,jjb:jje,:)+masse(:,jjb:jje,:)
406  flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:)
407  . +flux_u(:,jjb:jje,:)
408  if (pole_sud) jje=jj_end-1
409  flux_v_cum(:,jjb:jje,:)=flux_v_cum(:,jjb:jje,:)
410  . +flux_v(:,jjb:jje,:)
411 
412  jjb=jj_begin
413  jje=jj_end
414 
415  do iq=1,nq
416  q_cum(:,jjb:jje,:,iq)=q_cum(:,jjb:jje,:,iq)
417  . +q(:,jjb:jje,:,iq)*masse(:,jjb:jje,:)
418  enddo
419 
420 c=====================================================================
421 c FLUX ET TENDANCES
422 c=====================================================================
423 
424 c Flux longitudinal
425 c -----------------
426  do iq=1,nq
427  do l=1,llm
428  do j=jjb,jje
429  do i=1,iim
430  flux_uq_cum(i,j,l,iq)=flux_uq_cum(i,j,l,iq)
431  s +flux_u(i,j,l)*0.5*(q(i,j,l,iq)+q(i+1,j,l,iq))
432  enddo
433  flux_uq_cum(iip1,j,l,iq)=flux_uq_cum(1,j,l,iq)
434  enddo
435  enddo
436  enddo
437 
438 c flux méridien
439 c -------------
440  do iq=1,nq
441  call register_hallo(q(1,1,1,iq),ip1jmp1,llm,0,1,1,0,req)
442  enddo
443  call sendrequest(req)
444  call waitrequest(req)
445 
446  jjb=jj_begin
447  jje=jj_end
448  if (pole_sud) jje=jj_end-1
449 
450  do iq=1,nq
451  do l=1,llm
452  do j=jjb,jje
453  do i=1,iip1
454  flux_vq_cum(i,j,l,iq)=flux_vq_cum(i,j,l,iq)
455  s +flux_v(i,j,l)*0.5*(q(i,j,l,iq)+q(i,j+1,l,iq))
456  enddo
457  enddo
458  enddo
459  enddo
460 
461 
462 c tendances
463 c ---------
464 
465 c convergence horizontale
466  call register_hallo(flux_uq_cum,ip1jmp1,llm,2,2,2,2,req)
467  call register_hallo(flux_vq_cum,ip1jm,llm,2,2,2,2,req)
468  call sendrequest(req)
469  call waitrequest(req)
470 
471  call convflu_p(flux_uq_cum,flux_vq_cum,llm*nq,dq)
472 
473 c calcul de la vitesse verticale
474  call register_hallo(flux_u_cum,ip1jmp1,llm,2,2,2,2,req)
475  call register_hallo(flux_v_cum,ip1jm,llm,2,2,2,2,req)
476  call sendrequest(req)
477  call waitrequest(req)
478 
479  call convmas_p(flux_u_cum,flux_v_cum,convm)
480  CALL vitvert_p(convm,w)
481 
482  jjb=jj_begin
483  jje=jj_end
484 
485  do iq=1,nq
486  do l=1,llm-1
487  do j=jjb,jje
488  do i=1,iip1
489  ww=-0.5*w(i,j,l+1)*(q(i,j,l,iq)+q(i,j,l+1,iq))
490  dq(i,j,l ,iq)=dq(i,j,l ,iq)-ww
491  dq(i,j,l+1,iq)=dq(i,j,l+1,iq)+ww
492  enddo
493  enddo
494  enddo
495  enddo
496  IF (prt_level > 5)
497  . WRITE(lunout,*)'Apres les calculs fait a chaque pas'
498 c=====================================================================
499 c PAS DE TEMPS D'ECRITURE
500 c=====================================================================
501  if (icum.eq.ncum) then
502 c=====================================================================
503 
504  IF (prt_level > 5)
505  . WRITE(lunout,*)'Pas d ecriture'
506 
507 c Normalisation
508  do iq=1,nq
509  q_cum(:,jjb:jje,:,iq)=q_cum(:,jjb:jje,:,iq)
510  . /masse_cum(:,jjb:jje,:)
511  enddo
512  zz=1./REAL(ncum)
513 
514  jjb=jj_begin
515  jje=jj_end
516 
517  ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
518  masse_cum(:,jjb:jje,:)=masse_cum(:,jjb:jje,:)*zz
519  flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:)*zz
520  flux_uq_cum(:,jjb:jje,:,:)=flux_uq_cum(:,jjb:jje,:,:)*zz
521  dq(:,jjb:jje,:,:)=dq(:,jjb:jje,:,:)*zz
522 
523  IF (pole_sud) jje=jj_end-1
524  flux_v_cum(:,jjb:jje,:)=flux_v_cum(:,jjb:jje,:)*zz
525  flux_vq_cum(:,jjb:jje,:,:)=flux_vq_cum(:,jjb:jje,:,:)*zz
526 
527  jjb=jj_begin
528  jje=jj_end
529 
530 
531 c A retravailler eventuellement
532 c division de dQ par la masse pour revenir aux bonnes grandeurs
533  do iq=1,nq
534  dq(:,jjb:jje,:,iq)=dq(:,jjb:jje,:,iq)/masse_cum(:,jjb:jje,:)
535  enddo
536 
537 c=====================================================================
538 c Transport méridien
539 c=====================================================================
540 
541 c cumul zonal des masses des mailles
542 c ----------------------------------
543  jjb=jj_begin
544  jje=jj_end
545  if (pole_sud) jje=jj_end-1
546 
547  zv(jjb:jje,:)=0.
548  zmasse(jjb:jje,:)=0.
549 
550  call register_hallo(masse_cum,ip1jmp1,llm,1,1,1,1,req)
551  do iq=1,nq
552  call register_hallo(q_cum(1,1,1,iq),ip1jmp1,llm,0,1,1,0,req)
553  enddo
554 
555  call sendrequest(req)
556  call waitrequest(req)
557 
558  call massbar_p(masse_cum,massebx,masseby)
559 
560  jjb=jj_begin
561  jje=jj_end
562  if (pole_sud) jje=jj_end-1
563 
564  do l=1,llm
565  do j=jjb,jje
566  do i=1,iim
567  zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
568  zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
569  enddo
570  zfactv(j,l)=cv(1,j)/zmasse(j,l)
571  enddo
572  enddo
573 
574 c print*,'3OK'
575 c --------------------------------------------------------------
576 c calcul de la moyenne zonale du transport :
577 c ------------------------------------------
578 c
579 c --
580 c TOT : la circulation totale [ vq ]
581 c
582 c - -
583 c MMC : mean meridional circulation [ v ] [ q ]
584 c
585 c ---- -- - -
586 c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
587 c
588 c - * - * - - - -
589 c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
590 c
591 c - -
592 c on utilise aussi l'intermediaire TMP : [ v q ]
593 c
594 c la variable zfactv transforme un transport meridien cumule
595 c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
596 c
597 c --------------------------------------------------------------
598 
599 
600 c ----------------------------------------
601 c Transport dans le plan latitude-altitude
602 c ----------------------------------------
603 
604  jjb=jj_begin
605  jje=jj_end
606  if (pole_sud) jje=jj_end-1
607 
608  zvq=0.
609  psiq=0.
610  do iq=1,nq
611  zvqtmp=0.
612  do l=1,llm
613  do j=jjb,jje
614 c print*,'j,l,iQ=',j,l,iQ
615 c Calcul des moyennes zonales du transort total et de zvQtmp
616  do i=1,iim
617  zvq(j,l,itot,iq)=zvq(j,l,itot,iq)
618  s +flux_vq_cum(i,j,l,iq)
619  zqy= 0.5*(q_cum(i,j,l,iq)*masse_cum(i,j,l)+
620  s q_cum(i,j+1,l,iq)*masse_cum(i,j+1,l))
621  zvqtmp(j,l)=zvqtmp(j,l)+flux_v_cum(i,j,l)*zqy
622  s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
623  zvq(j,l,iave,iq)=zvq(j,l,iave,iq)+zqy
624  enddo
625 c print*,'aOK'
626 c Decomposition
627  zvq(j,l,iave,iq)=zvq(j,l,iave,iq)/zmasse(j,l)
628  zvq(j,l,itot,iq)=zvq(j,l,itot,iq)*zfactv(j,l)
629  zvqtmp(j,l)=zvqtmp(j,l)*zfactv(j,l)
630  zvq(j,l,immc,iq)=zv(j,l)*zvq(j,l,iave,iq)*zfactv(j,l)
631  zvq(j,l,itrs,iq)=zvq(j,l,itot,iq)-zvqtmp(j,l)
632  zvq(j,l,istn,iq)=zvqtmp(j,l)-zvq(j,l,immc,iq)
633  enddo
634  enddo
635 c fonction de courant meridienne pour la quantite Q
636  do l=llm,1,-1
637  do j=jjb,jje
638  psiq(j,l,iq)=psiq(j,l+1,iq)+zvq(j,l,itot,iq)
639  enddo
640  enddo
641  enddo
642 
643 c fonction de courant pour la circulation meridienne moyenne
644  psi(jjb:jje,:)=0.
645  do l=llm,1,-1
646  do j=jjb,jje
647  psi(j,l)=psi(j,l+1)+zv(j,l)
648  zv(j,l)=zv(j,l)*zfactv(j,l)
649  enddo
650  enddo
651 
652 c print*,'4OK'
653 c sorties proprement dites
654  if (i_sortie.eq.1) then
655  jjb=jj_begin
656  jje=jj_end
657  jjn=jj_nb
658  if (pole_sud) jje=jj_end-1
659  if (pole_sud) jjn=jj_nb-1
660 
661  do iq=1,nq
662  do itr=1,ntr
663  call histwrite(fileid,znom(itr,iq),itau,
664  s zvq(jjb:jje,:,itr,iq)
665  s ,jjn*llm,ndex3d)
666  enddo
667  call histwrite(fileid,'psi'//nom(iq),
668  s itau,psiq(jjb:jje,1:llm,iq)
669  s ,jjn*llm,ndex3d)
670  enddo
671 
672  call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
673  s ,jjn*llm,ndex3d)
674  call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
675  s ,jjn*llm,ndex3d)
676  psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
677  call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
678  s jjn*llm,ndex3d)
679 
680  endif
681 
682 
683 c -----------------
684 c Moyenne verticale
685 c -----------------
686 
687  zamasse(jjb:jje)=0.
688  do l=1,llm
689  zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
690  enddo
691 
692  zavq(jjb:jje,:,:)=0.
693  do iq=1,nq
694  do itr=2,ntr
695  do l=1,llm
696  zavq(jjb:jje,itr,iq)=zavq(jjb:jje,itr,iq)
697  s +zvq(jjb:jje,l,itr,iq)
698  s *zmasse(jjb:jje,l)
699  enddo
700  zavq(jjb:jje,itr,iq)=zavq(jjb:jje,itr,iq)/zamasse(jjb:jje)
701  call histwrite(fileid,'a'//znom(itr,iq),itau,
702  s zavq(jjb:jje,itr,iq),jjn*llm,ndex3d)
703  enddo
704  enddo
705 
706 c on doit pouvoir tracer systematiquement la fonction de courant.
707 
708 c=====================================================================
709 c/////////////////////////////////////////////////////////////////////
710  icum=0 !///////////////////////////////////////
711  endif ! icum.eq.ncum !///////////////////////////////////////
712 c/////////////////////////////////////////////////////////////////////
713 c=====================================================================
714 
715  return
716  end