source: trunk/LMDZ.COMMON/libf/dyn3d/bilan_dyn.F @ 1422

Last change on this file since 1422 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

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