/[lmdze]/trunk/dyn3d/bilan_dyn.f
ViewVC logotype

Diff of /trunk/dyn3d/bilan_dyn.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 44 by guez, Wed Apr 13 12:29:18 2011 UTC revision 57 by guez, Mon Jan 30 12:54:02 2012 UTC
# Line 5  module bilan_dyn_m Line 5  module bilan_dyn_m
5  contains  contains
6    
7    SUBROUTINE bilan_dyn(ps, masse, pk, flux_u, flux_v, teta, phi, ucov, vcov, &    SUBROUTINE bilan_dyn(ps, masse, pk, flux_u, flux_v, teta, phi, ucov, vcov, &
8         trac, dt_app, dt_cum)         trac)
9    
10      ! From LMDZ4/libf/dyn3d/bilan_dyn.F, version 1.5 2005/03/16      ! From LMDZ4/libf/dyn3d/bilan_dyn.F, version 1.5 2005/03/16 10:12:17
     ! 10:12:17 fairhead  
11    
12      ! Sous-programme consacré à des diagnostics dynamiques de base      ! Sous-programme consacré à des diagnostics dynamiques de base.
13      ! De façon générale, les moyennes des scalaires Q sont pondérées par      ! De façon générale, les moyennes des scalaires Q sont pondérées
14      ! la masse. Les flux de masse sont eux simplement moyennés.      ! par la masse. Les flux de masse sont, eux, simplement moyennés.
15    
16      USE histcom, ONLY: histbeg_totreg, histdef, histend, histvert      USE comconst, ONLY: cpp
17      USE calendar, ONLY: ymds2ju      USE comgeom, ONLY: constang_2d, cu_2d, cv_2d
     USE histwrite_m, ONLY: histwrite  
18      USE dimens_m, ONLY: iim, jjm, llm      USE dimens_m, ONLY: iim, jjm, llm
19        USE histwrite_m, ONLY: histwrite
20        use init_dynzon_m, only: ncum, fileid, znom, ntr, nq, nom
21      USE paramet_m, ONLY: iip1, jjp1      USE paramet_m, ONLY: iip1, jjp1
     USE comconst, ONLY: cpp  
     USE comvert, ONLY: presnivs  
     USE comgeom, ONLY: constang_2d, cu_2d, cv_2d, rlatv  
     USE temps, ONLY: annee_ref, day_ref, itau_dyn  
     USE inigrads_m, ONLY: inigrads  
     USE nr_util, ONLY: pi  
22    
23      ! Arguments:      ! Arguments:
24    
25      real, intent(in):: dt_app, dt_cum      real, intent(in):: ps(iip1, jjp1)
26      real ps(iip1, jjp1)      real, intent(in):: masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)
27      real masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)      real, intent(in):: flux_u(iip1, jjp1, llm)
28      real flux_u(iip1, jjp1, llm)      real, intent(in):: flux_v(iip1, jjm, llm)
     real flux_v(iip1, jjm, llm)  
29      real, intent(in):: teta(iip1, jjp1, llm)      real, intent(in):: teta(iip1, jjp1, llm)
30      real phi(iip1, jjp1, llm)      real, intent(in):: phi(iip1, jjp1, llm)
31      real ucov(iip1, jjp1, llm)      real, intent(in):: ucov(iip1, jjp1, llm)
32      real vcov(iip1, jjm, llm)      real, intent(in):: vcov(iip1, jjm, llm)
33      real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)      real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)
34    
35      ! Local:      ! Local:
36    
37      integer, save:: icum, ncum      integer:: icum  = 0
     logical:: first = .true.  
     real zz, zqy, zfactv(jjm, llm)  
   
     integer, parameter:: nQ=7  
   
     character(len=6), save:: nom(nQ)  
     character(len=6), save:: unites(nQ)  
   
     character(len=10) file  
     integer, parameter:: ifile=4  
   
     integer itemp, igeop, iecin, iang, iu, iovap, iun  
     integer:: i_sortie = 1  
   
     real:: time = 0.  
38      integer:: itau = 0      integer:: itau = 0
39        real zqy, zfactv(jjm, llm)
     data itemp, igeop, iecin, iang, iu, iovap, iun/1, 2, 3, 4, 5, 6, 7/  
40    
41      real ww      real ww
42    
# Line 67  contains Line 44  contains
44      REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)      REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)
45      REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)      REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)
46      REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)      REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
     REAL vorpot(iip1, jjm, llm)  
47      REAL w(iip1, jjp1, llm), ecin(iip1, jjp1, llm), convm(iip1, jjp1, llm)      REAL w(iip1, jjp1, llm), ecin(iip1, jjp1, llm), convm(iip1, jjp1, llm)
     REAL bern(iip1, jjp1, llm)  
48    
49      ! Champ contenant les scalaires advectés      ! Champ contenant les scalaires advectés
50      real Q(iip1, jjp1, llm, nQ)      real Q(iip1, jjp1, llm, nQ)
# Line 82  contains Line 57  contains
57      real, save:: Q_cum(iip1, jjp1, llm, nQ)      real, save:: Q_cum(iip1, jjp1, llm, nQ)
58      real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)      real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)
59      real, save:: flux_vQ_cum(iip1, jjm, llm, nQ)      real, save:: flux_vQ_cum(iip1, jjm, llm, nQ)
     real flux_wQ_cum(iip1, jjp1, llm, nQ)  
60      real dQ(iip1, jjp1, llm, nQ)      real dQ(iip1, jjp1, llm, nQ)
61    
62      ! champs de tansport en moyenne zonale      ! champs de tansport en moyenne zonale
63      integer itr      integer itr
64      integer, parameter:: ntr=5      integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5
   
     character(len=10), save:: znom(ntr, nQ)  
     character(len=20), save:: znoml(ntr, nQ)  
     character(len=10), save:: zunites(ntr, nQ)  
   
     integer iave, itot, immc, itrs, istn  
     data iave, itot, immc, itrs, istn/1, 2, 3, 4, 5/  
     character(len=3) ctrs(ntr)  
     data ctrs/' ', 'TOT', 'MMC', 'TRS', 'STN'/  
65    
66      real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)      real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)
67      real zavQ(jjm, ntr, nQ), psiQ(jjm, llm+1, nQ)      real zavQ(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)
68      real zmasse(jjm, llm), zamasse(jjm)      real zmasse(jjm, llm)
69        real zv(jjm, llm), psi(jjm, llm + 1)
     real zv(jjm, llm), psi(jjm, llm+1)  
   
70      integer i, j, l, iQ      integer i, j, l, iQ
71    
     ! Initialisation du fichier contenant les moyennes zonales.  
   
     integer, save:: fileid  
     integer thoriid, zvertiid  
     integer ndex3d(jjm*llm)  
   
     ! Variables locales  
   
     real zjulian  
     character(len=3) str  
     character(len=10) ctrac  
     integer ii, jj  
     integer zan, dayref  
   
     real rlong(jjm), rlatg(jjm)  
   
72      !-----------------------------------------------------------------      !-----------------------------------------------------------------
73    
     !!print *, "Call sequence information: bilan_dyn"  
   
     ! Initialisation  
   
     time=time+dt_app  
     itau=itau+1  
   
     if (first) then  
        icum=0  
        ! initialisation des fichiers  
        first=.false.  
        ! ncum est la frequence de stokage en pas de temps  
        ncum=dt_cum/dt_app  
        if (abs(ncum * dt_app - dt_cum) > 1e-5 * dt_app) then  
           print *, 'Problème : le pas de cumul doit être multiple du pas'  
           print *, 'dt_app=', dt_app  
           print *, 'dt_cum=', dt_cum  
           stop 1  
        endif  
   
        if (i_sortie == 1) then  
           file='dynzon'  
           call inigrads(ifile , (/0./), 180./pi, 0., 0., rlatv, -90., 90., &  
                180./pi , presnivs, 1. , dt_cum, file, 'dyn_zon ')  
        endif  
   
        nom(itemp)='T'  
        nom(igeop)='gz'  
        nom(iecin)='K'  
        nom(iang)='ang'  
        nom(iu)='u'  
        nom(iovap)='ovap'  
        nom(iun)='un'  
   
        unites(itemp)='K'  
        unites(igeop)='m2/s2'  
        unites(iecin)='m2/s2'  
        unites(iang)='ang'  
        unites(iu)='m/s'  
        unites(iovap)='kg/kg'  
        unites(iun)='un'  
   
        ! Initialisation du fichier contenant les moyennes zonales  
   
        zan = annee_ref  
        dayref = day_ref  
        CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)  
   
        rlong=0.  
        rlatg=rlatv*180./pi  
   
        call histbeg_totreg('dynzon', rlong(:1), rlatg, 1, 1, 1, jjm, itau_dyn, &  
             zjulian, dt_cum, thoriid, fileid)  
   
        ! Appel à histvert pour la grille verticale  
   
        call histvert(fileid, 'presnivs', 'Niveaux sigma', 'mb', llm, presnivs, &  
             zvertiid)  
   
        ! Appels à histdef pour la définition des variables à sauvegarder  
        do iQ=1, nQ  
           do itr=1, ntr  
              if(itr == 1) then  
                 znom(itr, iQ)=nom(iQ)  
                 znoml(itr, iQ)=nom(iQ)  
                 zunites(itr, iQ)=unites(iQ)  
              else  
                 znom(itr, iQ)=ctrs(itr)//'v'//nom(iQ)  
                 znoml(itr, iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)  
                 zunites(itr, iQ)='m/s * '//unites(iQ)  
              endif  
           enddo  
        enddo  
   
        ! Déclarations des champs avec dimension verticale  
        do iQ=1, nQ  
           do itr=1, ntr  
              call histdef(fileid, znom(itr, iQ), znoml(itr, iQ), &  
                   zunites(itr, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &  
                   'ave(X)', dt_cum, dt_cum)  
           enddo  
           ! Declarations pour les fonctions de courant  
           call histdef(fileid, 'psi'//nom(iQ), 'stream fn. '//znoml(itot, iQ), &  
                zunites(itot, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &  
                'ave(X)', dt_cum, dt_cum)  
        enddo  
   
        ! Declarations pour les champs de transport d'air  
        call histdef(fileid, 'masse', 'masse', &  
             'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid, &  
             'ave(X)', dt_cum, dt_cum)  
        call histdef(fileid, 'v', 'v', &  
             'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid, &  
             'ave(X)', dt_cum, dt_cum)  
        ! Declarations pour les fonctions de courant  
        call histdef(fileid, 'psi', 'stream fn. MMC ', 'mega t/s', &  
             1, jjm, thoriid, llm, 1, llm, zvertiid, &  
             'ave(X)', dt_cum, dt_cum)  
   
        ! Declaration des champs 1D de transport en latitude  
        do iQ=1, nQ  
           do itr=2, ntr  
              call histdef(fileid, 'a'//znom(itr, iQ), znoml(itr, iQ), &  
                   zunites(itr, iQ), 1, jjm, thoriid, 1, 1, 1, -99, &  
                   'ave(X)', dt_cum, dt_cum)  
           enddo  
        enddo  
   
        CALL histend(fileid)  
     endif  
   
74      ! Calcul des champs dynamiques      ! Calcul des champs dynamiques
75    
76      ! Énergie cinétique      ! Énergie cinétique
# Line 243  contains Line 79  contains
79      CALL enercin(vcov, ucov, vcont, ucont, ecin)      CALL enercin(vcov, ucov, vcont, ucont, ecin)
80    
81      ! moment cinétique      ! moment cinétique
82      do l=1, llm      do l = 1, llm
83         ang(:, :, l)=ucov(:, :, l)+constang_2d         ang(:, :, l) = ucov(:, :, l) + constang_2d
84         unat(:, :, l)=ucont(:, :, l)*cu_2d         unat(:, :, l) = ucont(:, :, l)*cu_2d
85      enddo      enddo
86    
87      Q(:, :, :, itemp)=teta*pk/cpp      Q(:, :, :, 1) = teta * pk / cpp
88      Q(:, :, :, igeop)=phi      Q(:, :, :, 2) = phi
89      Q(:, :, :, iecin)=ecin      Q(:, :, :, 3) = ecin
90      Q(:, :, :, iang)=ang      Q(:, :, :, 4) = ang
91      Q(:, :, :, iu)=unat      Q(:, :, :, 5) = unat
92      Q(:, :, :, iovap)=trac      Q(:, :, :, 6) = trac
93      Q(:, :, :, iun)=1.      Q(:, :, :, 7) = 1.
94    
95      ! Cumul      ! Cumul
96    
97      if(icum == 0) then      if (icum == 0) then
98         ps_cum=0.         ps_cum = 0.
99         masse_cum=0.         masse_cum = 0.
100         flux_u_cum=0.         flux_u_cum = 0.
101         flux_v_cum=0.         flux_v_cum = 0.
102         Q_cum=0.         Q_cum = 0.
103         flux_vQ_cum=0.         flux_vQ_cum = 0.
104         flux_uQ_cum=0.         flux_uQ_cum = 0.
105      endif      endif
106    
107      icum=icum+1      itau = itau + 1
108        icum = icum + 1
109    
110      ! Accumulation des flux de masse horizontaux      ! Accumulation des flux de masse horizontaux
111      ps_cum=ps_cum+ps      ps_cum = ps_cum + ps
112      masse_cum=masse_cum+masse      masse_cum = masse_cum + masse
113      flux_u_cum=flux_u_cum+flux_u      flux_u_cum = flux_u_cum + flux_u
114      flux_v_cum=flux_v_cum+flux_v      flux_v_cum = flux_v_cum + flux_v
115      do iQ=1, nQ      do iQ = 1, nQ
116         Q_cum(:, :, :, iQ)=Q_cum(:, :, :, iQ)+Q(:, :, :, iQ)*masse         Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) + Q(:, :, :, iQ)*masse
117      enddo      enddo
118    
119      ! FLUX ET TENDANCES      ! FLUX ET TENDANCES
120    
121      ! Flux longitudinal      ! Flux longitudinal
122      do iQ=1, nQ      forall (iQ = 1: nQ, i = 1: iim) flux_uQ_cum(i, :, :, iQ) &
123         do l=1, llm           = flux_uQ_cum(i, :, :, iQ) &
124            do j=1, jjp1           + flux_u(i, :, :) * 0.5 * (Q(i, :, :, iQ) + Q(i + 1, :, :, iQ))
125               do i=1, iim      flux_uQ_cum(iip1, :, :, :) = flux_uQ_cum(1, :, :, :)
126                  flux_uQ_cum(i, j, l, iQ)=flux_uQ_cum(i, j, l, iQ) &  
127                       +flux_u(i, j, l)*0.5*(Q(i, j, l, iQ)+Q(i+1, j, l, iQ))      ! Flux méridien
128               enddo      forall (iQ = 1: nQ, j = 1: jjm) flux_vQ_cum(:, j, :, iQ) &
129               flux_uQ_cum(iip1, j, l, iQ)=flux_uQ_cum(1, j, l, iQ)           = flux_vQ_cum(:, j, :, iQ) &
130            enddo           + flux_v(:, j, :) * 0.5 * (Q(:, j, :, iQ) + Q(:, j + 1, :, iQ))
        enddo  
     enddo  
   
     ! flux méridien  
     do iQ=1, nQ  
        do l=1, llm  
           do j=1, jjm  
              do i=1, iip1  
                 flux_vQ_cum(i, j, l, iQ)=flux_vQ_cum(i, j, l, iQ) &  
                      +flux_v(i, j, l)*0.5*(Q(i, j, l, iQ)+Q(i, j+1, l, iQ))  
              enddo  
           enddo  
        enddo  
     enddo  
131    
132      ! tendances      ! tendances
133    
# Line 315  contains Line 138  contains
138      call convmas(flux_u_cum, flux_v_cum, convm)      call convmas(flux_u_cum, flux_v_cum, convm)
139      CALL vitvert(convm, w)      CALL vitvert(convm, w)
140    
141      do iQ=1, nQ      do iQ = 1, nQ
142         do l=1, llm-1         do l = 1, llm-1
143            do j=1, jjp1            do j = 1, jjp1
144               do i=1, iip1               do i = 1, iip1
145                  ww=-0.5*w(i, j, l+1)*(Q(i, j, l, iQ)+Q(i, j, l+1, iQ))                  ww = -0.5*w(i, j, l + 1)*(Q(i, j, l, iQ) + Q(i, j, l + 1, iQ))
146                  dQ(i, j, l , iQ)=dQ(i, j, l , iQ)-ww                  dQ(i, j, l, iQ) = dQ(i, j, l, iQ)-ww
147                  dQ(i, j, l+1, iQ)=dQ(i, j, l+1, iQ)+ww                  dQ(i, j, l + 1, iQ) = dQ(i, j, l + 1, iQ) + ww
148               enddo               enddo
149            enddo            enddo
150         enddo         enddo
# Line 331  contains Line 154  contains
154    
155      writing_step: if (icum == ncum) then      writing_step: if (icum == ncum) then
156         ! Normalisation         ! Normalisation
157         do iQ=1, nQ         do iQ = 1, nQ
158            Q_cum(:, :, :, iQ)=Q_cum(:, :, :, iQ)/masse_cum            Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ)/masse_cum
159         enddo         enddo
160         zz=1./float(ncum)         ps_cum = ps_cum / ncum
161         ps_cum=ps_cum*zz         masse_cum = masse_cum / ncum
162         masse_cum=masse_cum*zz         flux_u_cum = flux_u_cum / ncum
163         flux_u_cum=flux_u_cum*zz         flux_v_cum = flux_v_cum / ncum
164         flux_v_cum=flux_v_cum*zz         flux_uQ_cum = flux_uQ_cum / ncum
165         flux_uQ_cum=flux_uQ_cum*zz         flux_vQ_cum = flux_vQ_cum / ncum
166         flux_vQ_cum=flux_vQ_cum*zz         dQ = dQ / ncum
        dQ=dQ*zz  
167    
168         ! A retravailler eventuellement         ! A retravailler eventuellement
169         ! division de dQ par la masse pour revenir aux bonnes grandeurs         ! division de dQ par la masse pour revenir aux bonnes grandeurs
170         do iQ=1, nQ         do iQ = 1, nQ
171            dQ(:, :, :, iQ)=dQ(:, :, :, iQ)/masse_cum            dQ(:, :, :, iQ) = dQ(:, :, :, iQ)/masse_cum
172         enddo         enddo
173    
174         ! Transport méridien         ! Transport méridien
175    
176         ! cumul zonal des masses des mailles         ! cumul zonal des masses des mailles
177    
178         zv=0.         zv = 0.
179         zmasse=0.         zmasse = 0.
180         call massbar(masse_cum, massebx, masseby)         call massbar(masse_cum, massebx, masseby)
181         do l=1, llm         do l = 1, llm
182            do j=1, jjm            do j = 1, jjm
183               do i=1, iim               do i = 1, iim
184                  zmasse(j, l)=zmasse(j, l)+masseby(i, j, l)                  zmasse(j, l) = zmasse(j, l) + masseby(i, j, l)
185                  zv(j, l)=zv(j, l)+flux_v_cum(i, j, l)                  zv(j, l) = zv(j, l) + flux_v_cum(i, j, l)
186               enddo               enddo
187               zfactv(j, l)=cv_2d(1, j)/zmasse(j, l)               zfactv(j, l) = cv_2d(1, j)/zmasse(j, l)
188            enddo            enddo
189         enddo         enddo
190    
191         ! Transport dans le plan latitude-altitude         ! Transport dans le plan latitude-altitude
192    
193         zvQ=0.         zvQ = 0.
194         psiQ=0.         psiQ = 0.
195         do iQ=1, nQ         do iQ = 1, nQ
196            zvQtmp=0.            zvQtmp = 0.
197            do l=1, llm            do l = 1, llm
198               do j=1, jjm               do j = 1, jjm
199                  ! Calcul des moyennes zonales du transort total et de zvQtmp                  ! Calcul des moyennes zonales du transort total et de zvQtmp
200                  do i=1, iim                  do i = 1, iim
201                     zvQ(j, l, itot, iQ)=zvQ(j, l, itot, iQ) &                     zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ) &
202                          +flux_vQ_cum(i, j, l, iQ)                          + flux_vQ_cum(i, j, l, iQ)
203                     zqy= 0.5*(Q_cum(i, j, l, iQ)*masse_cum(i, j, l)+ &                     zqy =  0.5 * (Q_cum(i, j, l, iQ) * masse_cum(i, j, l) &
204                          Q_cum(i, j+1, l, iQ)*masse_cum(i, j+1, l))                          + Q_cum(i, j + 1, l, iQ) * masse_cum(i, j + 1, l))
205                     zvQtmp(j, l)=zvQtmp(j, l)+flux_v_cum(i, j, l)*zqy &                     zvQtmp(j, l) = zvQtmp(j, l) + flux_v_cum(i, j, l) * zqy &
206                          /(0.5*(masse_cum(i, j, l)+masse_cum(i, j+1, l)))                          / (0.5 * (masse_cum(i, j, l) + masse_cum(i, j + 1, l)))
207                     zvQ(j, l, iave, iQ)=zvQ(j, l, iave, iQ)+zqy                     zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ) + zqy
208                  enddo                  enddo
209                  ! Decomposition                  ! Decomposition
210                  zvQ(j, l, iave, iQ)=zvQ(j, l, iave, iQ)/zmasse(j, l)                  zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ)/zmasse(j, l)
211                  zvQ(j, l, itot, iQ)=zvQ(j, l, itot, iQ)*zfactv(j, l)                  zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ)*zfactv(j, l)
212                  zvQtmp(j, l)=zvQtmp(j, l)*zfactv(j, l)                  zvQtmp(j, l) = zvQtmp(j, l)*zfactv(j, l)
213                  zvQ(j, l, immc, iQ)=zv(j, l)*zvQ(j, l, iave, iQ)*zfactv(j, l)                  zvQ(j, l, immc, iQ) = zv(j, l)*zvQ(j, l, iave, iQ)*zfactv(j, l)
214                  zvQ(j, l, itrs, iQ)=zvQ(j, l, itot, iQ)-zvQtmp(j, l)                  zvQ(j, l, itrs, iQ) = zvQ(j, l, itot, iQ)-zvQtmp(j, l)
215                  zvQ(j, l, istn, iQ)=zvQtmp(j, l)-zvQ(j, l, immc, iQ)                  zvQ(j, l, istn, iQ) = zvQtmp(j, l)-zvQ(j, l, immc, iQ)
216               enddo               enddo
217            enddo            enddo
218            ! fonction de courant meridienne pour la quantite Q            ! fonction de courant meridienne pour la quantite Q
219            do l=llm, 1, -1            do l = llm, 1, -1
220               do j=1, jjm               do j = 1, jjm
221                  psiQ(j, l, iQ)=psiQ(j, l+1, iQ)+zvQ(j, l, itot, iQ)                  psiQ(j, l, iQ) = psiQ(j, l + 1, iQ) + zvQ(j, l, itot, iQ)
222               enddo               enddo
223            enddo            enddo
224         enddo         enddo
225    
226         ! fonction de courant pour la circulation meridienne moyenne         ! fonction de courant pour la circulation meridienne moyenne
227         psi=0.         psi = 0.
228         do l=llm, 1, -1         do l = llm, 1, -1
229            do j=1, jjm            do j = 1, jjm
230               psi(j, l)=psi(j, l+1)+zv(j, l)               psi(j, l) = psi(j, l + 1) + zv(j, l)
231               zv(j, l)=zv(j, l)*zfactv(j, l)               zv(j, l) = zv(j, l)*zfactv(j, l)
232            enddo            enddo
233         enddo         enddo
234    
235         ! sorties proprement dites         ! sorties proprement dites
236         if (i_sortie == 1) then         do iQ = 1, nQ
237            do iQ=1, nQ            do itr = 1, ntr
238               do itr=1, ntr               call histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ))
                 call histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ))  
              enddo  
              call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, 1:llm, iQ))  
239            enddo            enddo
240              call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, :llm, iQ))
           call histwrite(fileid, 'masse', itau, zmasse)  
           call histwrite(fileid, 'v', itau, zv)  
           psi=psi*1.e-9  
           call histwrite(fileid, 'psi', itau, psi(:, 1:llm))  
        endif  
   
        ! Moyenne verticale  
   
        zamasse=0.  
        do l=1, llm  
           zamasse(:)=zamasse(:)+zmasse(:, l)  
241         enddo         enddo
242         zavQ=0.  
243         do iQ=1, nQ         call histwrite(fileid, 'masse', itau, zmasse)
244            do itr=2, ntr         call histwrite(fileid, 'v', itau, zv)
245               do l=1, llm         psi = psi*1.e-9
246                  zavQ(:, itr, iQ) = zavQ(:, itr, iQ) &         call histwrite(fileid, 'psi', itau, psi(:, :llm))
247                       + zvQ(:, l, itr, iQ) * zmasse(:, l)  
248               enddo         ! Intégrale verticale
249               zavQ(:, itr, iQ)=zavQ(:, itr, iQ)/zamasse(:)  
250           forall (iQ = 1: nQ, itr = 2: ntr) zavQ(:, itr, iQ) &
251                = sum(zvQ(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)
252    
253           do iQ = 1, nQ
254              do itr = 2, ntr
255               call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))               call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))
256            enddo            enddo
257         enddo         enddo
258    
259         ! on doit pouvoir tracer systematiquement la fonction de courant.         ! On doit pouvoir tracer systematiquement la fonction de courant.
260         icum=0         icum = 0
261      endif writing_step      endif writing_step
262    
263    end SUBROUTINE bilan_dyn    end SUBROUTINE bilan_dyn
Legend:
Removed from v.44  
changed lines
  Added in v.57
  ViewVC Help
Powered by ViewVC 1.1.21