source: trunk/libf/dyn3d/bilan_dyn.F @ 6

Last change on this file since 6 was 6, checked in by slebonnois, 14 years ago

cf commit_v6.log :

  • manipulation traceurs
  • homogeneisation .def
  • bilan_dyn
  • etats initiaux start.nc
  • appels specifiques pour physique
File size: 27.5 KB
Line 
1!
2! $Id: bilan_dyn.F 1403 2010-07-01 09:02:53Z fairhead $
3!
4      SUBROUTINE bilan_dyn (dt_app,dt_cum,
5     s  ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,
6     s  ducovdyn,ducovdis,ducovspg,ducovphy)
7c si besoin des traceurs:
8c      SUBROUTINE bilan_dyn (ntrac,dt_app,dt_cum,
9c     s  ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac,
10c     s  ducovdyn,ducovdis,ducovspg,ducovphy)
11
12c   AFAIRE
13c   Prevoir en champ nq+1 le diagnostique de l'energie
14c   en faisant Qzon=Cv T + L * ...
15c             vQ..A=Cp T + L * ...
16
17#ifdef CPP_IOIPSL
18      USE IOIPSL
19#endif
20
21      IMPLICIT NONE
22
23#include "dimensions.h"
24#include "paramet.h"
25#include "comconst.h"
26#include "comvert.h"
27#include "comgeom2.h"
28#include "temps.h"
29#include "iniprint.h"
30
31c====================================================================
32c
33c   Sous-programme consacre à des diagnostics dynamiques de base
34c
35c
36c   De facon generale, les moyennes des scalaires Q sont ponderees par
37c   la masse.
38c
39c   Les flux de masse sont eux simplement moyennes.
40c
41c====================================================================
42
43c   Arguments :
44c   ===========
45
46c      integer ntrac
47      real dt_app,dt_cum
48      real ps(iip1,jjp1)
49      real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
50      real flux_u(iip1,jjp1,llm)
51      real flux_v(iip1,jjm,llm)
52      real teta(iip1,jjp1,llm)
53      real phi(iip1,jjp1,llm)
54      real ucov(iip1,jjp1,llm)
55      real vcov(iip1,jjm,llm)
56c      real trac(iip1,jjp1,llm,ntrac)
57c Tendances en m/s2 :
58      real ducovdyn(iip1,jjp1,llm)
59      real ducovdis(iip1,jjp1,llm)
60      real ducovspg(iip1,jjp1,llm)
61      real ducovphy(iip1,jjp1,llm)
62
63c   Local :
64c   =======
65
66      integer icum,ncum
67      save icum,ncum
68
69      integer i,j,l,iQ,num
70      real zz,zqy,zfactv(jjm,llm),zfactw(jjm,llm)
71      character*2 strd2
72      real ww
73
74      logical first
75      save first
76      data first/.true./
77
78      integer i_sortie
79      save i_sortie
80      data i_sortie/1/
81
82      real time
83      integer itau
84      save time,itau
85      data time,itau/0.,0/
86
87! facteur = -1. pour Venus
88      real    fact_geovenus
89      save    fact_geovenus
90
91c   variables dynamiques intermédiaires
92c -----------------------------------
93      REAL vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
94      REAL ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
95      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
96      REAL vorpot(iip1,jjm,llm)
97      REAL w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm)
98      real temp(iip1,jjp1,llm)
99      real dudyn(iip1,jjp1,llm)
100      real dudis(iip1,jjp1,llm)
101      real duspg(iip1,jjp1,llm)
102      real duphy(iip1,jjp1,llm)
103
104c CHAMPS SCALAIRES Q ADVECTES
105c ----------------------------
106      integer nQ
107c avec tous les composes, ca fait trop.... Je les enleve
108c     parameter (nQ=6+nqmx)
109      parameter (nQ=6)
110
111      character*6,save :: nom(nQ)
112      character*6,save :: unites(nQ)
113
114      integer itemp,igeop,iecin,iang,iu,iun
115      save itemp,igeop,iecin,iang,iu,iun
116      data itemp,igeop,iecin,iang,iu,iun/1,2,3,4,5,6/
117
118c   champ contenant les scalaires advectés.
119      real Q(iip1,jjp1,llm,nQ)
120   
121c   champs cumulés
122      real ps_cum(iip1,jjp1)
123      real masse_cum(iip1,jjp1,llm)
124      real flux_u_cum(iip1,jjp1,llm)
125      real flux_v_cum(iip1,jjm,llm)
126      real flux_w_cum(iip1,jjp1,llm)
127      real Q_cum(iip1,jjp1,llm,nQ)
128      real flux_uQ_cum(iip1,jjp1,llm,nQ)
129      real flux_vQ_cum(iip1,jjm,llm,nQ)
130      real flux_wQ_cum(iip1,jjp1,llm,nQ)
131      real dQ(iip1,jjp1,llm,nQ)
132
133      save ps_cum,masse_cum,flux_u_cum,flux_v_cum
134      save Q_cum,flux_uQ_cum,flux_vQ_cum
135      save flux_w_cum,flux_wQ_cum
136
137c   champs de transport en moyenne zonale
138      integer ntr,itr
139      parameter (ntr=5)
140
141      character*10,save :: znom(ntr,nQ)
142      character*20,save :: znoml(ntr,nQ)
143      character*10,save :: zunites(ntr,nQ)
144      character*10,save :: znom2(ntr,nQ)
145      character*20,save :: znom2l(ntr,nQ)
146      character*10,save :: zunites2(ntr,nQ)
147      character*10,save :: znom3(ntr,nQ)
148      character*20,save :: znom3l(ntr,nQ)
149      character*10,save :: zunites3(ntr,nQ)
150
151      integer iave,itot,immc,itrs,istn
152      data iave,itot,immc,itrs,istn/1,2,3,4,5/
153      character*3 ctrs(ntr)
154      data ctrs/'  ','TOT','MMC','TRS','STN'/
155
156      real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
157      real zwQ(jjm,llm,ntr,nQ),zwQtmp(jjm,llm)
158      real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
159      real zawQ(jjm,llm,ntr,nQ)
160      real zdQ(jjm,llm,nQ)
161      real zmasse(jjm,llm),zavmasse(jjm),zawmasse(llm)
162      real psbar(jjm)
163
164      real zv(jjm,llm),zw(jjp1,llm),psi(jjm,llm+1)
165
166c TENDANCES POUR MOMENT CINETIQUE
167c -------------------------------
168
169      integer ntdc,itdc
170      parameter (ntdc=4)
171
172      integer itdcdyn,itdcdis,itdcspg,itdcphy
173      data    itdcdyn,itdcdis,itdcspg,itdcphy/1,2,3,4/
174
175      character*6,save :: nomtdc(ntdc)
176
177c   champ contenant les tendances du moment cinetique.
178      real    tdc(iip1,jjp1,llm,ntdc)
179      real    ztdc(jjm,llm,ntdc)   ! moyenne zonale
180   
181c   champs cumulés
182      real tdc_cum(iip1,jjp1,llm,ntdc)
183      save tdc_cum
184
185c   integrations completes
186      real mtot,mctot,dmctot(ntdc)
187
188c   Initialisation du fichier contenant les moyennes zonales.
189c   ---------------------------------------------------------
190
191      character*10 infile
192
193      integer fileid
194      integer thoriid, zvertiid
195      save fileid
196
197      integer ndex3d(jjm*llm)
198      real    ztmp3d(jjm,llm)
199
200C   Variables locales
201C
202      integer tau0
203      real zjulian
204      integer zan, dayref
205C
206      real rlong(jjm),rlatg(jjm)
207
208
209
210c=====================================================================
211c   Initialisation
212c=====================================================================
213
214      ndex3d=0
215
216      if (first) then
217
218        if (planet_type.eq."venus") then
219            fact_geovenus = -1.
220        else
221            fact_geovenus = 1.
222        endif
223
224        icum=0
225c       initialisation des fichiers
226        first=.false.
227c   ncum est la frequence de stokage en pas de temps
228        ncum=dt_cum/dt_app
229        if (abs(ncum*dt_app-dt_cum).gt.1.e-2*dt_app) then
230         if (abs((ncum+1)*dt_app-dt_cum).lt.1.e-2*dt_app) then
231           ncum = ncum+1
232         elseif (abs((ncum-1)*dt_app-dt_cum).lt.1.e-2*dt_app) then
233           ncum = ncum-1
234         else
235           WRITE(lunout,*)
236     .            'Pb : le pas de cumule doit etre multiple du pas'
237           WRITE(lunout,*)'dt_app=',dt_app
238           WRITE(lunout,*)'dt_cum=',dt_cum
239           WRITE(lunout,*)'ncum*dt_app=',ncum*dt_app
240           WRITE(lunout,*)'ncum=',ncum
241           stop
242         endif
243        endif
244
245c VARIABLES ADVECTEES:
246
247        nom(itemp)='temp'
248        nom(igeop)='gz'
249        nom(iecin)='ecin'
250        nom(iang)='ang'
251        nom(iu)='u'
252        nom(iun)='un'
253
254        unites(itemp)='K'
255        unites(igeop)='m2/s2'
256        unites(iecin)='m2/s2'
257        unites(iang)='ang'
258        unites(iu)='m/s'
259        unites(iun)='un'
260
261c avec tous les composes, ca fait trop.... Je les enleve
262c       do num=1,ntrac
263c        write(strd2,'(i2.2)') num
264c        nom(6+num)='trac'//strd2
265c        unites(6+num)='kg/kg'
266c       enddo
267
268c TENDANCES MOMENT CIN:
269       
270        nomtdc(itdcdyn) ='dmcdyn'
271        nomtdc(itdcdis) ='dmcdis'
272        nomtdc(itdcspg) ='dmcspg'
273        nomtdc(itdcphy) ='dmcphy'
274
275c   Initialisation du fichier contenant les moyennes zonales.
276c   ---------------------------------------------------------
277
278      infile='dynzon'
279
280      zan = annee_ref
281      dayref = day_ref
282      CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
283c     tau0 = itau_dyn
284      tau0 = 0
285      itau = tau0
286     
287      rlong=0.
288      rlatg=rlatv*180./pi*fact_geovenus
289       
290      call histbeg(infile, 1, rlong, jjm, rlatg,
291     .             1, 1, 1, jjm,
292     .             tau0, zjulian, dt_cum, thoriid, fileid)
293
294C
295C  Appel a histvert pour la grille verticale
296C
297      call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
298     .              llm, presnivs, zvertiid)
299C
300C  Appels a histdef pour la definition des variables a sauvegarder
301
302      do iQ=1,nQ
303         do itr=1,ntr
304            if(itr.eq.1) then
305               znom(itr,iQ)    =nom(iQ)
306               znoml(itr,iQ)   =nom(iQ)
307               zunites(itr,iQ) =unites(iQ)
308            else
309           znom(itr,iQ)    =ctrs(itr)//'v'//nom(iQ)
310           znoml(itr,iQ)   ='transport : v * '//nom(iQ)//' '//ctrs(itr)
311           zunites(itr,iQ) ='m/s * '//unites(iQ)
312           znom2(itr,iQ)   =ctrs(itr)//'w'//nom(iQ)
313           znom2l(itr,iQ)  ='transport: w * '//nom(iQ)//' '//ctrs(itr)
314           zunites2(itr,iQ)='Pa/s * '//unites(iQ)
315            endif
316         enddo
317               znom3(iQ)='d'//nom(iQ)
318               znom3l(iQ)='convergence: '//nom(iQ)
319               zunites3(iQ)=unites(iQ)//' /s'
320c          print*,'DEBUG:',znom3(iQ),znom3l(iQ),zunites3(iQ)
321      enddo
322
323c   Declarations des champs avec dimension verticale
324
325      if (1.eq.0) then  ! on les sort, ou pas...
326
327c     do iQ=1,nQ
328c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
329      do iQ=1,4,3
330         do itr=1,ntr
331      IF (prt_level > 5)
332     . WRITE(lunout,*)'var ',itr,iQ
333     .      ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
334            call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
335     .        zunites(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
336     .        32,'ave(X)',dt_cum,dt_cum)
337         enddo
338c transport vertical:
339         do itr=2,ntr
340      IF (prt_level > 5)
341     . WRITE(lunout,*)'var ',itr,iQ
342     .      ,znom2(itr,iQ),znom2l(itr,iQ),zunites2(itr,iQ)
343            call histdef(fileid,znom2(itr,iQ),znom2l(itr,iQ),
344     .        zunites2(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
345     .        32,'ins(X)',dt_cum,dt_cum)
346         enddo
347
348c Declarations pour convergences
349      IF (prt_level > 5)
350     . WRITE(lunout,*)'var ',iQ
351     .      ,znom3(iQ),znom3l(iQ),zunites3(iQ)
352            call histdef(fileid,znom3(iQ),znom3l(iQ),
353     .        zunites3(iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
354     .        32,'ins(X)',dt_cum,dt_cum)
355
356c   Declarations pour les fonctions de courant
357c   Non sorties ici...
358c          call histdef(fileid,'psi'//nom(iQ)
359c     .      ,'stream fn. '//znoml(itot,iQ),
360c     .      zunites(itot,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
361c     .      32,'ave(X)',dt_cum,dt_cum)
362
363      enddo ! iQ
364
365      endif ! 1=1 sortie ou non...
366
367c   Declarations pour les champs de transport d'air
368      call histdef(fileid, 'masse', 'masse',
369     .             'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid,
370     .             32, 'ave(X)', dt_cum, dt_cum)
371      call histdef(fileid, 'v', 'v',
372     .             'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
373     .             32, 'ave(X)', dt_cum, dt_cum)
374      call histdef(fileid, 'w', 'w',
375     .             'Pa/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
376     .             32, 'ins(X)', dt_cum, dt_cum)
377
378c   Declarations pour la fonction de courant
379          call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
380     .      1,jjm,thoriid,llm,1,llm,zvertiid,
381     .      32,'ave(X)',dt_cum,dt_cum)
382
383
384c   Declarations pour les tendances de moment cinetique
385      do itdc=1,ntdc
386      call histdef(fileid, nomtdc(itdc), nomtdc(itdc),
387     .             'ang/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
388     .             32, 'ins(X)', dt_cum, dt_cum)
389      enddo
390
391c   Declaration des champs 1D de transport en latitude
392c     do iQ=1,nQ
393c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
394      do iQ=1,4,3
395         do itr=2,ntr
396            call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
397     .        zunites(itr,iQ),1,jjm,thoriid,1,1,1,-99,
398     .        32,'ave(X)',dt_cum,dt_cum)
399c JE VIRE LE VERTICAL POUR L'INSTANT
400c           call histdef(fileid,'a'//znom2(itr,iQ),znom2l(itr,iQ),
401c    .        zunites2(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
402c    .        32,'ins(X)',dt_cum,dt_cum)
403         enddo
404      enddo
405
406               CALL histend(fileid)
407
408
409      endif  ! first
410
411
412c=====================================================================
413c   Calcul des champs dynamiques
414c   ----------------------------
415
416c   énergie cinétique
417      ucont(:,:,:)=0
418      CALL covcont(llm,ucov,vcov,ucont,vcont)
419      CALL enercin(vcov,ucov,vcont,ucont,ecin)
420
421c   moment cinétique et tendances
422      dudyn = 0.
423      dudis = 0.
424      duspg = 0.
425      duphy = 0.
426      do l=1,llm
427         unat(:,:,l)=ucont(:,:,l)*cu(:,:)
428         dudyn(:,2:jjm,l) = ducovdyn(:,2:jjm,l)/cu(:,2:jjm)
429         dudis(:,2:jjm,l) = ducovdis(:,2:jjm,l)/cu(:,2:jjm)
430         duspg(:,2:jjm,l) = ducovspg(:,2:jjm,l)/cu(:,2:jjm)
431         duphy(:,2:jjm,l) = ducovphy(:,2:jjm,l)/cu(:,2:jjm)
432         do j=1,jjp1
433          ang(:,j,l)= rad*cos(rlatu(j))*
434     .     ( unat(:,j,l) + rad*cos(rlatu(j))*omeg )
435          tdc(:,j,l,1) = rad*cos(rlatu(j))*dudyn(:,j,l)
436          tdc(:,j,l,2) = rad*cos(rlatu(j))*dudis(:,j,l)
437          tdc(:,j,l,3) = rad*cos(rlatu(j))*duspg(:,j,l)
438          tdc(:,j,l,4) = rad*cos(rlatu(j))*duphy(:,j,l)
439         enddo
440      enddo
441c Normalisation:
442      ang = ang / (2./3. *rad*rad*omeg)
443      do itdc=1,ntdc
444        tdc(:,:,:,itdc)=tdc(:,:,:,itdc) / (2./3. *rad*rad*omeg)
445      enddo
446
447! ADAPTATION GCM POUR CP(T)
448      call tpot2t(ip1jmp1*llm,teta,temp,pk)
449      Q(:,:,:,itemp) = temp(:,:,:)
450      Q(:,:,:,igeop) =phi(:,:,:)
451      Q(:,:,:,iecin) =ecin(:,:,:)
452      Q(:,:,:,iang)  =ang(:,:,:)
453      Q(:,:,:,iu)    =unat(:,:,:)
454      Q(:,:,:,iun)   =1.
455c avec tous les composes, ca fait trop.... Je les enleve
456c     do num=1,ntrac
457c      Q(:,:,:,6+num)=trac(:,:,:,num)
458c     enddo
459
460c   calcul du flux de masse vertical (+ vers le bas)
461      call convmas(flux_u,flux_v,convm)
462      CALL vitvert(convm,w)
463
464c=====================================================================
465c   Cumul
466c=====================================================================
467c
468      if(icum.EQ.0) then
469         ps_cum      = 0.
470         masse_cum   = 0.
471         flux_u_cum  = 0.
472         flux_v_cum  = 0.
473         flux_w_cum  = 0.
474         Q_cum       = 0.
475         flux_vQ_cum = 0.
476         flux_uQ_cum = 0.
477         flux_wQ_cum = 0.
478         tdc_cum     = 0.
479      endif
480
481      IF (prt_level > 5)
482     . WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
483      icum=icum+1
484
485c   accumulation des flux de masse horizontaux
486      ps_cum          = ps_cum     + ps
487      masse_cum       = masse_cum  + masse
488      flux_u_cum      = flux_u_cum + flux_u
489      flux_v_cum      = flux_v_cum + flux_v
490      flux_w_cum      = flux_w_cum + w
491      do iQ=1,nQ
492      Q_cum(:,:,:,iQ) = Q_cum(:,:,:,iQ) + Q(:,:,:,iQ)*masse(:,:,:)
493      enddo
494      do itdc=1,ntdc
495      tdc_cum(:,:,:,itdc) =
496     .       tdc_cum(:,:,:,itdc) + tdc(:,:,:,itdc)*masse(:,:,:)
497      enddo
498
499c=====================================================================
500c  FLUX ET TENDANCES
501c=====================================================================
502
503c   Flux longitudinal
504c   -----------------
505      do iQ=1,nQ
506         do l=1,llm
507            do j=1,jjp1
508               do i=1,iim
509                  flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
510     s            +flux_u(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
511               enddo
512               flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
513            enddo
514         enddo
515      enddo
516
517c    flux méridien
518c    -------------
519      do iQ=1,nQ
520         do l=1,llm
521            do j=1,jjm
522               do i=1,iip1
523                  flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
524     s            +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
525               enddo
526            enddo
527         enddo
528      enddo
529
530c   Flux vertical
531c   -------------
532      do iQ=1,nQ
533         do l=2,llm
534            do j=1,jjp1
535               do i=1,iip1
536                  flux_wQ_cum(i,j,l,iQ)=flux_wQ_cum(i,j,l,iQ)
537     s            +w(i,j,l)*0.5*(Q(i,j,l-1,iQ)+Q(i,j,l,iQ))
538               enddo
539            enddo
540         enddo
541         flux_wQ_cum(:,:,1,iQ)=0.0e0
542      enddo
543
544c    tendances
545c    ---------
546
547c   convergence horizontale
548      call  convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
549
550c   calcul de la vitesse verticale
551      call convmas(flux_u_cum,flux_v_cum,convm)
552      CALL vitvert(convm,w)
553
554c  ajustement tendances (vitesse verticale)
555      do iQ=1,nQ
556         do l=1,llm-1
557            do j=1,jjp1
558               do i=1,iip1
559                  ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
560                  dQ(i,j,l  ,iQ)=dQ(i,j,l  ,iQ)-ww
561                  dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
562               enddo
563            enddo
564         enddo
565      enddo
566      IF (prt_level > 5)
567     . WRITE(lunout,*)'Apres les calculs fait a chaque pas'
568
569c=====================================================================
570c   PAS DE TEMPS D'ECRITURE
571c=====================================================================
572      if (icum.eq.ncum) then
573c=====================================================================
574
575      time=time+dt_cum
576      itau=itau+1
577
578      IF (prt_level > 5)
579     . WRITE(lunout,*)'Pas d ecriture'
580
581c   Normalisation
582      do iQ=1,nQ
583         Q_cum(:,:,:,iQ) = Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
584         dQ(:,:,:,iQ)    = dQ(:,:,:,iQ)   /masse_cum(:,:,:)
585      enddo
586      do itdc=1,ntdc
587         tdc_cum(:,:,:,itdc) = tdc_cum(:,:,:,itdc)/masse_cum(:,:,:)
588      enddo
589
590      zz=1./REAL(ncum)
591      ps_cum      = ps_cum      *zz
592      masse_cum   = masse_cum   *zz
593      flux_u_cum  = flux_u_cum  *zz
594      flux_v_cum  = flux_v_cum  *zz
595      flux_w_cum  = flux_w_cum  *zz
596      flux_uQ_cum = flux_uQ_cum *zz
597      flux_vQ_cum = flux_vQ_cum *zz
598      flux_wQ_cum = flux_wQ_cum *zz
599
600c Integration complete
601      mtot  = 0.
602      mctot  = 0.
603      dmctot = 0.
604      do l=1,llm
605       do j=2,jjm
606        do i=1,iim
607          mtot  = mtot  + masse_cum(i,j,l)
608          mctot = mctot + Q_cum(i,j,l,iang)*masse_cum(i,j,l)
609        enddo
610       enddo
611      enddo
612      mctot = mctot/mtot
613      do itdc=1,ntdc
614      do l=1,llm
615       do j=2,jjm
616        do i=1,iim
617          dmctot(itdc) = dmctot(itdc)
618     .               + tdc_cum(i,j,l,itdc)*masse_cum(i,j,l)/mtot
619        enddo
620       enddo
621      enddo
622      enddo
623
624c=====================================================================
625c   Transport méridien
626c=====================================================================
627
628c   cumul zonal des masses des mailles
629c   ----------------------------------
630      zv=0.
631      zw=0.
632      zmasse=0.
633      call massbar(masse_cum,massebx,masseby)
634
635c moy zonale de la ps cumulee
636         do j=1,jjm
637            psbar(j)=0.
638            do i=1,iim
639               psbar(j)=psbar(j)+ps_cum(i,j)/iim
640            enddo
641         enddo
642
643      do l=1,llm
644         do j=1,jjm
645            do i=1,iim
646               zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
647               zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
648               zw(j,l)=zw(j,l)+flux_w_cum(i,j,l)
649            enddo
650            zfactv(j,l)=cv(1,j)/zmasse(j,l)
651            zfactw(j,l)=((ap(l)-ap(l+1))+psbar(j)*(bp(l)-bp(l+1)))
652     .                    /zmasse(j,l)
653         enddo
654            do i=1,iim
655               zw(jjp1,l)=zw(jjp1,l)+flux_w_cum(i,jjp1,l)
656            enddo
657      enddo
658
659c     print*,'3OK'
660c   --------------------------------------------------------------
661c   calcul de la moyenne zonale du transport :
662c   ------------------------------------------
663c
664c                                     --
665c TOT : la circulation totale       [ vq ]
666c
667c                                      -     -
668c MMC : mean meridional circulation [ v ] [ q ]
669c
670c                                     ----      --       - -
671c TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
672c
673c                                     - * - *       - -       -     -
674c STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
675c
676c                                              - -
677c    on utilise aussi l'intermediaire TMP :  [ v q ]
678c
679c    la variable zfactv transforme un transport meridien cumule
680c    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
681c    la variable zfactw transforme un transport vertical cumule
682c    en kg/s * unte-du-champ-transporte en Pa/s * unite-du-champ-transporte
683c
684c   --------------------------------------------------------------
685
686
687c   ----------------------------------------
688c   Transport dans le plan latitude-altitude
689c   ----------------------------------------
690
691      zvQ=0.
692      zwQ=0.
693      zdQ=0.
694      psiQ=0.
695
696      do iQ=1,nQ
697
698c   transport meridien
699         zvQtmp=0.
700         do l=1,llm
701            do j=1,jjm
702c              print*,'j,l,iQ=',j,l,iQ
703c   Calcul des moyennes zonales du transport total et de zvQtmp
704               do i=1,iim
705                  zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
706     s                            +flux_vQ_cum(i,j,l,iQ)
707                  zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+
708     s                           Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
709                  zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
710     s             /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
711                  zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
712               enddo
713c              print*,'aOK'
714c   Decomposition
715               zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
716               zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
717               zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
718               zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
719               zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
720               zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
721            enddo
722         enddo
723c   fonction de courant meridienne pour la quantite Q
724         do l=llm,1,-1
725            do j=1,jjm
726             psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)/zfactv(j,l)
727            enddo
728         enddo
729      enddo
730
731c   transport vertical
732         zwQtmp=0.
733         do l=1,llm
734            do j=1,jjm
735c              print*,'j,l,iQ=',j,l,iQ
736c   Calcul des moyennes zonales du transport vertical total et de zwQtmp
737               do i=1,iim
738                  zwQ(j,l,itot,iQ)=zwQ(j,l,itot,iQ)
739     s                            +flux_wQ_cum(i,j,l,iQ)
740                  zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+
741     s                           Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
742                  zwQtmp(j,l)=zwQtmp(j,l)+flux_w_cum(i,j,l)*zqy
743     s             /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
744                  zwQ(j,l,iave,iQ)=zwQ(j,l,iave,iQ)+zqy
745               enddo
746c   Decomposition
747               zwQ(j,l,iave,iQ)=zwQ(j,l,iave,iQ)/zmasse(j,l)
748               zwQ(j,l,itot,iQ)=zwQ(j,l,itot,iQ)*zfactw(j,l)
749               zwQtmp(j,l)=zwQtmp(j,l)*zfactw(j,l)
750               zwQ(j,l,immc,iQ)=zw(j,l)*zwQ(j,l,iave,iQ)*zfactw(j,l)
751               zwQ(j,l,itrs,iQ)=zwQ(j,l,itot,iQ)-zwQtmp(j,l)
752               zwQ(j,l,istn,iQ)=zwQtmp(j,l)-zwQ(j,l,immc,iQ)
753            enddo
754         enddo
755
756c   convergence
757c   Calcul moyenne zonale de la convergence totale
758         do l=1,llm
759            do j=1,jjm
760c              print*,'j,l,iQ=',j,l,iQ
761               do i=1,iim
762                  zdQ(j,l,iQ)=zdQ(j,l,iQ) +
763     .                   ( dQ(i,j,l,iQ)   * masse_cum(i,j,l)
764     .                   + dQ(i,j+1,l,iQ) * masse_cum(i,j+1,l))
765     .                 / ( masse_cum(i,j,l)+masse_cum(i,j+1,l))
766               enddo
767            enddo
768         enddo
769      enddo
770
771c   fonction de courant pour la circulation meridienne moyenne
772      psi=0.
773      do l=llm,1,-1
774         do j=1,jjm
775            psi(j,l)= psi(j,l+1)+zv(j,l)
776            zv(j,l) = zv(j,l)*zfactv(j,l)
777            zw(j,l) = 0.5*(zw(j,l)+zw(j+1,l))*zfactw(j,l)
778         enddo
779      enddo
780
781c   Calcul moyenne zonale des tendances moment cin.
782      ztdc=0.
783      do itdc=1,ntdc
784         do l=1,llm
785            do j=1,jjm
786               do i=1,iim
787                  ztdc(j,l,itdc)=ztdc(j,l,itdc) +
788     .            ( tdc_cum(i,j,l,itdc)   * masse_cum(i,j,l)
789     .            + tdc_cum(i,j+1,l,itdc) * masse_cum(i,j+1,l))
790     .          / ( masse_cum(i,j,l)+masse_cum(i,j+1,l))
791               enddo
792            enddo
793         enddo
794      enddo
795
796c     print*,'4OK'
797
798c--------------------------------------
799c--------------------------------------
800c   sorties proprement dites
801c--------------------------------------
802c--------------------------------------
803
804      if (i_sortie.eq.1) then
805
806c sortie des integrations completes dans le listing
807      write(*,'(A12,5(1PE11.4,X))') "BILANMCDYN  ",mctot,dmctot
808
809c sorties dans fichier dynzon
810
811      if (1.eq.0) then  ! on les sort, ou pas...
812
813c avec tous les composes, ca fait trop.... Je les enleve
814c      do iQ=1,nQ
815c      do iQ=1,6
816c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
817      do iQ=1,4,3
818
819         ztmp3d(:,:)= zvQ(:,:,1,iQ) ! valeur moyenne
820            call histwrite(fileid,znom(1,iQ),itau,ztmp3d
821     s      ,jjm*llm,ndex3d)
822       do itr=2,ntr
823         ztmp3d(:,:)= zvQ(:,:,itr,iQ)*fact_geovenus ! transport horizontal
824            call histwrite(fileid,znom(itr,iQ),itau,ztmp3d
825     s      ,jjm*llm,ndex3d)
826       enddo
827
828       do itr=2,ntr
829         ztmp3d(:,:)=zwQ(:,:,itr,iQ)
830            call histwrite(fileid,znom2(itr,iQ),itau,ztmp3d
831     s      ,jjm*llm,ndex3d)
832       enddo
833       
834         ztmp3d(:,:)= zdQ(:,:,iQ)
835            call histwrite(fileid,znom3(iQ),itau,ztmp3d
836     s      ,jjm*llm,ndex3d)
837
838c        ztmp3d(:,:)= psiQ(:,1:llm,iQ)*fact_geovenus
839c        call histwrite(fileid,'psi'//nom(iQ),itau,ztmp3d
840c    s      ,jjm*llm,ndex3d)
841      enddo
842
843      endif ! 1=1 sortie ou non...
844
845      ztmp3d=zmasse
846      call histwrite(fileid,'masse',itau,ztmp3d
847     s   ,jjm*llm,ndex3d)
848     
849      ztmp3d= zv*fact_geovenus
850      call histwrite(fileid,'v',itau,ztmp3d
851     s   ,jjm*llm,ndex3d)
852      ztmp3d(:,:)=zw(1:jjm,:)
853      call histwrite(fileid,'w',itau,ztmp3d
854     s   ,jjm*llm,ndex3d)
855      ztmp3d= psi(:,1:llm)*1.e-9*fact_geovenus
856      call histwrite(fileid,'psi',itau,ztmp3d,jjm*llm,ndex3d)
857
858      do itdc=1,ntdc
859         ztmp3d(:,:)= ztdc(:,:,itdc)
860         call histwrite(fileid,nomtdc(itdc),itau,ztmp3d
861     s    ,jjm*llm,ndex3d)
862      enddo
863
864      endif ! i_sortie
865
866
867c   -----------------
868c   Moyenne verticale
869c   -----------------
870
871      zavmasse=0.
872      do l=1,llm
873         zavmasse(:)=zavmasse(:)+zmasse(:,l)
874      enddo
875      zavQ=0.
876
877c avec tous les composes, ca fait trop.... Je les enleve
878c      do iQ=1,nQ
879c      do iQ=1,6
880c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
881      do iQ=1,4,3
882         do itr=2,ntr
883            do l=1,llm
884               zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
885            enddo
886            zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zavmasse(:)
887      if (i_sortie.eq.1) then
888         ztmp3d=0.0
889         ztmp3d(:,1)= zavQ(:,itr,iQ)*fact_geovenus
890         call histwrite(fileid,'a'//znom(itr,iQ),itau,ztmp3d
891     .      ,jjm*llm,ndex3d)     
892      endif
893         enddo
894      enddo
895
896c   ------------------
897c   Moyenne meridienne
898c   ------------------
899
900      zawmasse=0.
901      do j=1,jjm
902           do l=1,llm
903         zawmasse(l)=zawmasse(l)+zmasse(j,l)
904           enddo
905      enddo
906      zawQ=0.
907
908c avec tous les composes, ca fait trop.... Je les enleve
909c      do iQ=1,nQ
910c      do iQ=1,6
911c !!!! JE NE SORS ICI QUE temp et ang POUR CAUSE DE PLACE !
912      do iQ=1,4,3
913         do itr=2,ntr
914           do l=1,llm
915            do j=1,jjp1
916          zawQ(1,l,itr,iQ)=zawQ(1,l,itr,iQ)+zwQ(j,l,itr,iQ)*zmasse(j,l)
917            enddo
918            zawQ(1,l,itr,iQ)=zawQ(1,l,itr,iQ)/zawmasse(l)
919           enddo
920      if (i_sortie.eq.1) then
921         ztmp3d=0.0
922           do l=1,llm
923         ztmp3d(1,l)=zawQ(1,l,itr,iQ)
924           enddo
925c JE VIRE LE VERTICAL POUR L'INSTANT
926c        call histwrite(fileid,'a'//znom2(itr,iQ),itau,ztmp3d
927c    .      ,jjm*llm,ndex3d)     
928      endif
929         enddo
930      enddo
931
932      call histsync(fileid)
933
934c=====================================================================
935c/////////////////////////////////////////////////////////////////////
936      icum=0                  !///////////////////////////////////////
937      endif ! icum.eq.ncum    !///////////////////////////////////////
938c/////////////////////////////////////////////////////////////////////
939c=====================================================================
940
941      return
942      end
Note: See TracBrowser for help on using the repository browser.