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

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

Serie de modifs SL pour homogeneisation des phytitan et phyvenus
Ca touche aussi aux liens phy/dyn (surtout a propos de clesphy0),
a verifier avec les autres, donc...

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