source: LMDZ4/tags/LMDZ4_V3_3/libf/dyn3dpar/bilan_dyn_p.F @ 919

Last change on this file since 919 was 919, checked in by (none), 16 years ago

This commit was manufactured by cvs2svn to create tag 'LMDZ4_V3_3'.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.9 KB
Line 
1!
2! $Header$
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
7c   AFAIRE
8c   Prevoir en champ nq+1 le diagnostique de l'energie
9c   en faisant Qzon=Cv T + L * ...
10c             vQ..A=Cp T + L * ...
11
12      USE IOIPSL
13      USE parallel
14      USE mod_hallo
15      use misc_mod
16      use write_field
17      USE Write_Field_p
18      IMPLICIT NONE
19
20#include "dimensions.h"
21#include "paramet.h"
22#include "comconst.h"
23#include "comvert.h"
24#include "comgeom2.h"
25#include "temps.h"
26#include "iniprint.h"
27
28c====================================================================
29c
30c   Sous-programme consacre à des diagnostics dynamiques de base
31c
32c
33c   De facon generale, les moyennes des scalaires Q sont ponderees par
34c   la masse.
35c
36c   Les flux de masse sont eux simplement moyennes.
37c
38c====================================================================
39
40c   Arguments :
41c   ===========
42
43      integer ntrac
44      real dt_app,dt_cum
45      real ps(iip1,jjp1)
46      real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
47      real flux_u(iip1,jjp1,llm)
48      real flux_v(iip1,jjm,llm)
49      real teta(iip1,jjp1,llm)
50      real phi(iip1,jjp1,llm)
51      real ucov(iip1,jjp1,llm)
52      real vcov(iip1,jjm,llm)
53      real trac(iip1,jjp1,llm,ntrac)
54
55c   Local :
56c   =======
57
58      integer icum,ncum
59      logical first
60      real zz,zqy,zfactv(jjm,llm)
61
62      integer nQ
63      parameter (nQ=7)
64
65
66cym      character*6 nom(nQ)
67cym      character*6 unites(nQ)
68      character*6,save :: nom(nQ)
69      character*6,save :: unites(nQ)
70
71      character*10 file
72      integer ifile
73      parameter (ifile=4)
74
75      integer itemp,igeop,iecin,iang,iu,iovap,iun
76      integer i_sortie
77
78      save first,icum,ncum
79      save itemp,igeop,iecin,iang,iu,iovap,iun
80      save i_sortie
81
82      real time
83      integer itau
84      save time,itau
85      data time,itau/0.,0/
86
87      data first/.true./
88      data itemp,igeop,iecin,iang,iu,iovap,iun/1,2,3,4,5,6,7/
89      data i_sortie/1/
90
91      real ww
92
93c   variables dynamiques intermédiaires
94      REAL vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
95      REAL ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
96      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
97      REAL vorpot(iip1,jjm,llm)
98      REAL w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm)
99      REAL bern(iip1,jjp1,llm)
100
101c   champ contenant les scalaires advectés.
102      real Q(iip1,jjp1,llm,nQ)
103   
104c   champs cumulés
105      real ps_cum(iip1,jjp1)
106      real masse_cum(iip1,jjp1,llm)
107      real flux_u_cum(iip1,jjp1,llm)
108      real flux_v_cum(iip1,jjm,llm)
109      real Q_cum(iip1,jjp1,llm,nQ)
110      real flux_uQ_cum(iip1,jjp1,llm,nQ)
111      real flux_vQ_cum(iip1,jjm,llm,nQ)
112      real flux_wQ_cum(iip1,jjp1,llm,nQ)
113      real dQ(iip1,jjp1,llm,nQ)
114
115      save ps_cum,masse_cum,flux_u_cum,flux_v_cum
116      save Q_cum,flux_uQ_cum,flux_vQ_cum
117
118c   champs de tansport en moyenne zonale
119      integer ntr,itr
120      parameter (ntr=5)
121
122cym      character*10 znom(ntr,nQ)
123cym      character*20 znoml(ntr,nQ)
124cym      character*10 zunites(ntr,nQ)
125      character*10,save :: znom(ntr,nQ)
126      character*20,save :: znoml(ntr,nQ)
127      character*10,save :: zunites(ntr,nQ)
128
129      integer iave,itot,immc,itrs,istn
130      data iave,itot,immc,itrs,istn/1,2,3,4,5/
131      character*3 ctrs(ntr)
132      data ctrs/'  ','TOT','MMC','TRS','STN'/
133
134      real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
135      real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
136      real zmasse(jjm,llm),zamasse(jjm)
137
138      real zv(jjm,llm),psi(jjm,llm+1)
139
140      integer i,j,l,iQ
141
142
143c   Initialisation du fichier contenant les moyennes zonales.
144c   ---------------------------------------------------------
145
146      character*10 infile
147
148      integer fileid
149      integer thoriid, zvertiid
150      save fileid
151
152      integer ndex3d(jjm*llm)
153
154C   Variables locales
155C
156      integer tau0
157      real zjulian
158      character*3 str
159      character*10 ctrac
160      integer ii,jj
161      integer zan, dayref
162C
163      real rlong(jjm),rlatg(jjm)
164      integer :: jjb,jje,jjn,ijb,ije
165      type(Request) :: Req
166
167! definition du domaine d'ecriture pour le rebuild
168
169      INTEGER,DIMENSION(1) :: ddid
170      INTEGER,DIMENSION(1) :: dsg
171      INTEGER,DIMENSION(1) :: dsl
172      INTEGER,DIMENSION(1) :: dpf
173      INTEGER,DIMENSION(1) :: dpl
174      INTEGER,DIMENSION(1) :: dhs
175      INTEGER,DIMENSION(1) :: dhe
176     
177      INTEGER :: bilan_dyn_domain_id
178
179
180c=====================================================================
181c   Initialisation
182c=====================================================================
183      ndex3d=0
184      if (adjust) return
185     
186      time=time+dt_app
187      itau=itau+1
188
189      if (first) then
190
191
192        icum=0
193c       initialisation des fichiers
194        first=.false.
195c   ncum est la frequence de stokage en pas de temps
196        ncum=dt_cum/dt_app
197        if (abs(ncum*dt_app-dt_cum).gt.1.e-5*dt_app) then
198           WRITE(lunout,*)
199     .            'Pb : le pas de cumule doit etre multiple du pas'
200           WRITE(lunout,*)'dt_app=',dt_app
201           WRITE(lunout,*)'dt_cum=',dt_cum
202           stop
203        endif
204
205        if (i_sortie.eq.1) then
206         file='dynzon'
207         if (mpi_rank==0) then
208         call inigrads(ifile,1
209     s  ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi
210     s  ,llm,presnivs,1.
211     s  ,dt_cum,file,'dyn_zon ')
212         endif
213        endif
214
215        nom(itemp)='T'
216        nom(igeop)='gz'
217        nom(iecin)='K'
218        nom(iang)='ang'
219        nom(iu)='u'
220        nom(iovap)='ovap'
221        nom(iun)='un'
222
223        unites(itemp)='K'
224        unites(igeop)='m2/s2'
225        unites(iecin)='m2/s2'
226        unites(iang)='ang'
227        unites(iu)='m/s'
228        unites(iovap)='kg/kg'
229        unites(iun)='un'
230
231
232c   Initialisation du fichier contenant les moyennes zonales.
233c   ---------------------------------------------------------
234
235      infile='dynzon'
236
237      zan = annee_ref
238      dayref = day_ref
239      CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
240      tau0 = itau_dyn
241     
242      rlong=0.
243      rlatg=rlatv*180./pi
244
245      jjb=jj_begin
246      jje=jj_end
247      jjn=jj_nb
248      IF (pole_sud) THEN
249        jjn=jj_nb-1
250        jje=jj_end-1
251      ENDIF
252
253      ddid=(/ 2 /)
254      dsg=(/ jjm /)
255      dsl=(/ jjn /)
256      dpf=(/ jjb /)
257      dpl=(/ jje /)
258      dhs=(/ 0 /)
259      dhe=(/ 0 /)
260
261      call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe,
262     .                 'box',bilan_dyn_domain_id)
263       
264      call histbeg(trim(infile),
265     .             1, rlong(jjb:jje), jjn, rlatg(jjb:jje),
266     .             1, 1, 1, jjn,
267     .             tau0, zjulian, dt_cum, thoriid, fileid,
268     .             bilan_dyn_domain_id)
269
270C
271C  Appel a histvert pour la grille verticale
272C
273      call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
274     .              llm, presnivs, zvertiid)
275C
276C  Appels a histdef pour la definition des variables a sauvegarder
277      do iQ=1,nQ
278         do itr=1,ntr
279            if(itr.eq.1) then
280               znom(itr,iQ)=nom(iQ)
281               znoml(itr,iQ)=nom(iQ)
282               zunites(itr,iQ)=unites(iQ)
283            else
284               znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
285               znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
286               zunites(itr,iQ)='m/s * '//unites(iQ)
287            endif
288         enddo
289      enddo
290
291c   Declarations des champs avec dimension verticale
292c      print*,'1HISTDEF'
293      do iQ=1,nQ
294         do itr=1,ntr
295      IF (prt_level > 5)
296     . WRITE(lunout,*)'var ',itr,iQ
297     .      ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
298            call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
299     .        zunites(itr,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
300     .        32,'ave(X)',dt_cum,dt_cum)
301         enddo
302c   Declarations pour les fonctions de courant
303c      print*,'2HISTDEF'
304          call histdef(fileid,'psi'//nom(iQ)
305     .      ,'stream fn. '//znoml(itot,iQ),
306     .      zunites(itot,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
307     .      32,'ave(X)',dt_cum,dt_cum)
308      enddo
309
310
311c   Declarations pour les champs de transport d'air
312c      print*,'3HISTDEF'
313      call histdef(fileid, 'masse', 'masse',
314     .             'kg', 1, jjn, thoriid, llm, 1, llm, zvertiid,
315     .             32, 'ave(X)', dt_cum, dt_cum)
316      call histdef(fileid, 'v', 'v',
317     .             'm/s', 1, jjn, thoriid, llm, 1, llm, zvertiid,
318     .             32, 'ave(X)', dt_cum, dt_cum)
319c   Declarations pour les fonctions de courant
320c      print*,'4HISTDEF'
321          call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
322     .      1,jjn,thoriid,llm,1,llm,zvertiid,
323     .      32,'ave(X)',dt_cum,dt_cum)
324
325
326c   Declaration des champs 1D de transport en latitude
327c      print*,'5HISTDEF'
328      do iQ=1,nQ
329         do itr=2,ntr
330            call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
331     .        zunites(itr,iQ),1,jjn,thoriid,1,1,1,-99,
332     .        32,'ave(X)',dt_cum,dt_cum)
333         enddo
334      enddo
335
336
337c      print*,'8HISTDEF'
338               CALL histend(fileid)
339
340
341      endif
342
343
344c=====================================================================
345c   Calcul des champs dynamiques
346c   ----------------------------
347
348      jjb=jj_begin
349      jje=jj_end
350   
351c   énergie cinétique
352      ucont(:,jjb:jje,:)=0
353
354      call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Req)
355      call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Req)
356      call SendRequest(Req)
357      call WaitRequest(Req)
358
359      CALL covcont_p(llm,ucov,vcov,ucont,vcont)
360      CALL enercin_p(vcov,ucov,vcont,ucont,ecin)
361
362c   moment cinétique
363      do l=1,llm
364         ang(:,jjb:jje,l)=ucov(:,jjb:jje,l)+constang(:,jjb:jje)
365         unat(:,jjb:jje,l)=ucont(:,jjb:jje,l)*cu(:,jjb:jje)
366      enddo
367
368      Q(:,jjb:jje,:,itemp)=teta(:,jjb:jje,:)*pk(:,jjb:jje,:)/cpp
369      Q(:,jjb:jje,:,igeop)=phi(:,jjb:jje,:)
370      Q(:,jjb:jje,:,iecin)=ecin(:,jjb:jje,:)
371      Q(:,jjb:jje,:,iang)=ang(:,jjb:jje,:)
372      Q(:,jjb:jje,:,iu)=unat(:,jjb:jje,:)
373      Q(:,jjb:jje,:,iovap)=trac(:,jjb:jje,:,1)
374      Q(:,jjb:jje,:,iun)=1.
375
376
377c=====================================================================
378c   Cumul
379c=====================================================================
380c
381      if(icum.EQ.0) then
382         jjb=jj_begin
383         jje=jj_end
384
385         ps_cum(:,jjb:jje)=0.
386         masse_cum(:,jjb:jje,:)=0.
387         flux_u_cum(:,jjb:jje,:)=0.
388         Q_cum(:,jjb:jje,:,:)=0.
389         flux_uQ_cum(:,jjb:jje,:,:)=0.
390         flux_v_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
400c   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
420c=====================================================================
421c  FLUX ET TENDANCES
422c=====================================================================
423
424c   Flux longitudinal
425c   -----------------
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
438c    flux méridien
439c    -------------
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
462c    tendances
463c    ---------
464
465c   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
473c   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'
498c=====================================================================
499c   PAS DE TEMPS D'ECRITURE
500c=====================================================================
501      if (icum.eq.ncum) then
502c=====================================================================
503
504      IF (prt_level > 5)
505     . WRITE(lunout,*)'Pas d ecriture'
506
507c   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./float(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
531c   A retravailler eventuellement
532c   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 
537c=====================================================================
538c   Transport méridien
539c=====================================================================
540
541c   cumul zonal des masses des mailles
542c   ----------------------------------
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
574c     print*,'3OK'
575c   --------------------------------------------------------------
576c   calcul de la moyenne zonale du transport :
577c   ------------------------------------------
578c
579c                                     --
580c TOT : la circulation totale       [ vq ]
581c
582c                                      -     -
583c MMC : mean meridional circulation [ v ] [ q ]
584c
585c                                     ----      --       - -
586c TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
587c
588c                                     - * - *       - -       -     -
589c STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
590c
591c                                              - -
592c    on utilise aussi l'intermediaire TMP :  [ v q ]
593c
594c    la variable zfactv transforme un transport meridien cumule
595c    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
596c
597c   --------------------------------------------------------------
598
599
600c   ----------------------------------------
601c   Transport dans le plan latitude-altitude
602c   ----------------------------------------
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
614c              print*,'j,l,iQ=',j,l,iQ
615c   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
625c              print*,'aOK'
626c   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
635c   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
643c   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
652c     print*,'4OK'
653c   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
683c   -----------------
684c   Moyenne verticale
685c   -----------------
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
706c     on doit pouvoir tracer systematiquement la fonction de courant.
707
708c=====================================================================
709c/////////////////////////////////////////////////////////////////////
710      icum=0                  !///////////////////////////////////////
711      endif ! icum.eq.ncum    !///////////////////////////////////////
712c/////////////////////////////////////////////////////////////////////
713c=====================================================================
714
715      return
716      end
Note: See TracBrowser for help on using the repository browser.