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

Diff of /trunk/Sources/dyn3d/bilan_dyn.f

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

trunk/libf/dyn3d/bilan_dyn.f90 revision 40 by guez, Tue Feb 22 13:49:36 2011 UTC trunk/Sources/dyn3d/bilan_dyn.f revision 134 by guez, Wed Apr 29 15:47:56 2015 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 massbar_m, only: massbar
22      USE paramet_m, ONLY: iip1, jjp1      USE paramet_m, ONLY: iip1, jjp1
23      USE comconst, ONLY: cpp  
24      USE comvert, ONLY: presnivs      real, intent(in):: ps(iip1, jjp1)
25      USE comgeom, ONLY: constang_2d, cu_2d, cv_2d, rlatv      real, intent(in):: masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)
26      USE temps, ONLY: annee_ref, day_ref, itau_dyn      real, intent(in):: flux_u(iip1, jjp1, llm)
27      USE inigrads_m, ONLY: inigrads      real, intent(in):: flux_v(iip1, jjm, llm)
28      USE nr_util, ONLY: pi      real, intent(in):: teta(iip1, jjp1, llm)
29        real, intent(in):: phi(iip1, jjp1, llm)
30      ! Arguments:      real, intent(in):: ucov(:, :, :) ! (iip1, jjp1, llm)
31        real, intent(in):: vcov(iip1, jjm, llm)
     real, intent(in):: dt_app, dt_cum  
     real ps(iip1, jjp1)  
     real masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)  
     real flux_u(iip1, jjp1, llm)  
     real flux_v(iip1, jjm, llm)  
     real teta(iip1, jjp1, llm)  
     real phi(iip1, jjp1, llm)  
     real ucov(iip1, jjp1, llm)  
     real vcov(iip1, jjm, llm)  
32      real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)      real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)
33    
34      ! Local:      ! Local:
35    
36      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.  
37      integer:: itau = 0      integer:: itau = 0
38        real qy, factv(jjm, llm)
     data itemp, igeop, iecin, iang, iu, iovap, iun/1, 2, 3, 4, 5, 6, 7/  
   
     real ww  
39    
40      ! Variables dynamiques intermédiaires      ! Variables dynamiques intermédiaires
41      REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)      REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)
42      REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)      REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)
43      REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)      REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
44      REAL vorpot(iip1, jjm, llm)      REAL ecin(iip1, jjp1, llm)
     REAL w(iip1, jjp1, llm), ecin(iip1, jjp1, llm), convm(iip1, jjp1, llm)  
     REAL bern(iip1, jjp1, llm)  
45    
46      ! Champ contenant les scalaires advectés      ! Champ contenant les scalaires advectés
47      real Q(iip1, jjp1, llm, nQ)      real Q(iip1, jjp1, llm, nQ)
# Line 82  contains Line 54  contains
54      real, save:: Q_cum(iip1, jjp1, llm, nQ)      real, save:: Q_cum(iip1, jjp1, llm, nQ)
55      real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)      real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)
56      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)  
     real dQ(iip1, jjp1, llm, nQ)  
57    
58      ! champs de tansport en moyenne zonale      ! champs de tansport en moyenne zonale
59      integer itr      integer itr
60      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'/  
   
     real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)  
     real zavQ(jjm, ntr, nQ), psiQ(jjm, llm+1, nQ)  
     real zmasse(jjm, llm), zamasse(jjm)  
   
     real zv(jjm, llm), psi(jjm, llm+1)  
61    
62        real vq(jjm, llm, ntr, nQ), vqtmp(jjm, llm)
63        real avq(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)
64        real zmasse(jjm, llm)
65        real v(jjm, llm), psi(jjm, llm + 1)
66      integer i, j, l, iQ      integer i, j, l, iQ
67    
     ! 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)  
   
68      !-----------------------------------------------------------------      !-----------------------------------------------------------------
69    
     !!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  
   
70      ! Calcul des champs dynamiques      ! Calcul des champs dynamiques
71    
72      ! Énergie cinétique      ! Énergie cinétique
# Line 243  contains Line 75  contains
75      CALL enercin(vcov, ucov, vcont, ucont, ecin)      CALL enercin(vcov, ucov, vcont, ucont, ecin)
76    
77      ! moment cinétique      ! moment cinétique
78      do l=1, llm      forall (l = 1: llm)
79         ang(:, :, l)=ucov(:, :, l)+constang_2d         ang(:, :, l) = ucov(:, :, l) + constang_2d
80         unat(:, :, l)=ucont(:, :, l)*cu_2d         unat(:, :, l) = ucont(:, :, l) * cu_2d
81      enddo      end forall
82    
83      Q(:, :, :, itemp)=teta*pk/cpp      Q(:, :, :, 1) = teta * pk / cpp
84      Q(:, :, :, igeop)=phi      Q(:, :, :, 2) = phi
85      Q(:, :, :, iecin)=ecin      Q(:, :, :, 3) = ecin
86      Q(:, :, :, iang)=ang      Q(:, :, :, 4) = ang
87      Q(:, :, :, iu)=unat      Q(:, :, :, 5) = unat
88      Q(:, :, :, iovap)=trac      Q(:, :, :, 6) = trac
89      Q(:, :, :, iun)=1.      Q(:, :, :, 7) = 1.
90    
91      ! Cumul      ! Cumul
92    
93      if(icum == 0) then      if (icum == 0) then
94         ps_cum=0.         ps_cum = 0.
95         masse_cum=0.         masse_cum = 0.
96         flux_u_cum=0.         flux_u_cum = 0.
97         flux_v_cum=0.         flux_v_cum = 0.
98         Q_cum=0.         Q_cum = 0.
99         flux_vQ_cum=0.         flux_vQ_cum = 0.
100         flux_uQ_cum=0.         flux_uQ_cum = 0.
101      endif      endif
102    
103      icum=icum+1      itau = itau + 1
104        icum = icum + 1
105    
106      ! Accumulation des flux de masse horizontaux      ! Accumulation des flux de masse horizontaux
107      ps_cum=ps_cum+ps      ps_cum = ps_cum + ps
108      masse_cum=masse_cum+masse      masse_cum = masse_cum + masse
109      flux_u_cum=flux_u_cum+flux_u      flux_u_cum = flux_u_cum + flux_u
110      flux_v_cum=flux_v_cum+flux_v      flux_v_cum = flux_v_cum + flux_v
111      do iQ=1, nQ      forall (iQ = 1: nQ) Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) &
112         Q_cum(:, :, :, iQ)=Q_cum(:, :, :, iQ)+Q(:, :, :, iQ)*masse           + Q(:, :, :, iQ) * masse
     enddo  
   
     ! FLUX ET TENDANCES  
113    
114      ! Flux longitudinal      ! Flux longitudinal
115      do iQ=1, nQ      forall (iQ = 1: nQ, i = 1: iim) flux_uQ_cum(i, :, :, iQ) &
116         do l=1, llm           = flux_uQ_cum(i, :, :, iQ) &
117            do j=1, jjp1           + flux_u(i, :, :) * 0.5 * (Q(i, :, :, iQ) + Q(i + 1, :, :, iQ))
118               do i=1, iim      flux_uQ_cum(iip1, :, :, :) = flux_uQ_cum(1, :, :, :)
119                  flux_uQ_cum(i, j, l, iQ)=flux_uQ_cum(i, j, l, iQ) &  
120                       +flux_u(i, j, l)*0.5*(Q(i, j, l, iQ)+Q(i+1, j, l, iQ))      ! Flux méridien
121               enddo      forall (iQ = 1: nQ, j = 1: jjm) flux_vQ_cum(:, j, :, iQ) &
122               flux_uQ_cum(iip1, j, l, iQ)=flux_uQ_cum(1, j, l, iQ)           = flux_vQ_cum(:, j, :, iQ) &
123            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  
   
     ! tendances  
   
     ! convergence horizontale  
     call convflu(flux_uQ_cum, flux_vQ_cum, llm*nQ, dQ)  
   
     ! calcul de la vitesse verticale  
     call convmas(flux_u_cum, flux_v_cum, convm)  
     CALL vitvert(convm, w)  
   
     do iQ=1, nQ  
        do l=1, llm-1  
           do j=1, jjp1  
              do i=1, iip1  
                 ww=-0.5*w(i, j, l+1)*(Q(i, j, l, iQ)+Q(i, j, l+1, iQ))  
                 dQ(i, j, l , iQ)=dQ(i, j, l , iQ)-ww  
                 dQ(i, j, l+1, iQ)=dQ(i, j, l+1, iQ)+ww  
              enddo  
           enddo  
        enddo  
     enddo  
   
     ! PAS DE TEMPS D'ECRITURE  
124    
125      writing_step: if (icum == ncum) then      writing_step: if (icum == ncum) then
126         ! Normalisation         ! Normalisation
127         do iQ=1, nQ         forall (iQ = 1: nQ) Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) / masse_cum
128            Q_cum(:, :, :, iQ)=Q_cum(:, :, :, iQ)/masse_cum         ps_cum = ps_cum / ncum
129         enddo         masse_cum = masse_cum / ncum
130         zz=1./float(ncum)         flux_u_cum = flux_u_cum / ncum
131         ps_cum=ps_cum*zz         flux_v_cum = flux_v_cum / ncum
132         masse_cum=masse_cum*zz         flux_uQ_cum = flux_uQ_cum / ncum
133         flux_u_cum=flux_u_cum*zz         flux_vQ_cum = flux_vQ_cum / ncum
        flux_v_cum=flux_v_cum*zz  
        flux_uQ_cum=flux_uQ_cum*zz  
        flux_vQ_cum=flux_vQ_cum*zz  
        dQ=dQ*zz  
   
        ! A retravailler eventuellement  
        ! division de dQ par la masse pour revenir aux bonnes grandeurs  
        do iQ=1, nQ  
           dQ(:, :, :, iQ)=dQ(:, :, :, iQ)/masse_cum  
        enddo  
134    
135         ! Transport méridien         ! Transport méridien
136    
137         ! cumul zonal des masses des mailles         ! Cumul zonal des masses des mailles
138    
139         zv=0.         v = 0.
140         zmasse=0.         zmasse = 0.
141         call massbar(masse_cum, massebx, masseby)         call massbar(masse_cum, massebx, masseby)
142         do l=1, llm         do l = 1, llm
143            do j=1, jjm            do j = 1, jjm
144               do i=1, iim               do i = 1, iim
145                  zmasse(j, l)=zmasse(j, l)+masseby(i, j, l)                  zmasse(j, l) = zmasse(j, l) + masseby(i, j, l)
146                  zv(j, l)=zv(j, l)+flux_v_cum(i, j, l)                  v(j, l) = v(j, l) + flux_v_cum(i, j, l)
147               enddo               enddo
148               zfactv(j, l)=cv_2d(1, j)/zmasse(j, l)               factv(j, l) = cv_2d(1, j) / zmasse(j, l)
149            enddo            enddo
150         enddo         enddo
151    
152         ! Transport dans le plan latitude-altitude         ! Transport dans le plan latitude-altitude
153    
154         zvQ=0.         vq = 0.
155         psiQ=0.         psiQ = 0.
156         do iQ=1, nQ         do iQ = 1, nQ
157            zvQtmp=0.            vqtmp = 0.
158            do l=1, llm            do l = 1, llm
159               do j=1, jjm               do j = 1, jjm
160                  ! Calcul des moyennes zonales du transort total et de zvQtmp                  ! Calcul des moyennes zonales du transport total et de vqtmp
161                  do i=1, iim                  do i = 1, iim
162                     zvQ(j, l, itot, iQ)=zvQ(j, l, itot, iQ) &                     vq(j, l, itot, iQ) = vq(j, l, itot, iQ) &
163                          +flux_vQ_cum(i, j, l, iQ)                          + flux_vQ_cum(i, j, l, iQ)
164                     zqy= 0.5*(Q_cum(i, j, l, iQ)*masse_cum(i, j, l)+ &                     qy =  0.5 * (Q_cum(i, j, l, iQ) * masse_cum(i, j, l) &
165                          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))
166                     zvQtmp(j, l)=zvQtmp(j, l)+flux_v_cum(i, j, l)*zqy &                     vqtmp(j, l) = vqtmp(j, l) + flux_v_cum(i, j, l) * qy &
167                          /(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)))
168                     zvQ(j, l, iave, iQ)=zvQ(j, l, iave, iQ)+zqy                     vq(j, l, iave, iQ) = vq(j, l, iave, iQ) + qy
169                  enddo                  enddo
170                  ! Decomposition                  ! Decomposition
171                  zvQ(j, l, iave, iQ)=zvQ(j, l, iave, iQ)/zmasse(j, l)                  vq(j, l, iave, iQ) = vq(j, l, iave, iQ) / zmasse(j, l)
172                  zvQ(j, l, itot, iQ)=zvQ(j, l, itot, iQ)*zfactv(j, l)                  vq(j, l, itot, iQ) = vq(j, l, itot, iQ) * factv(j, l)
173                  zvQtmp(j, l)=zvQtmp(j, l)*zfactv(j, l)                  vqtmp(j, l) = vqtmp(j, l) * factv(j, l)
174                  zvQ(j, l, immc, iQ)=zv(j, l)*zvQ(j, l, iave, iQ)*zfactv(j, l)                  vq(j, l, immc, iQ) = v(j, l) * vq(j, l, iave, iQ) * factv(j, l)
175                  zvQ(j, l, itrs, iQ)=zvQ(j, l, itot, iQ)-zvQtmp(j, l)                  vq(j, l, itrs, iQ) = vq(j, l, itot, iQ) - vqtmp(j, l)
176                  zvQ(j, l, istn, iQ)=zvQtmp(j, l)-zvQ(j, l, immc, iQ)                  vq(j, l, istn, iQ) = vqtmp(j, l) - vq(j, l, immc, iQ)
177               enddo               enddo
178            enddo            enddo
179            ! fonction de courant meridienne pour la quantite Q            ! Fonction de courant méridienne pour la quantité Q
180            do l=llm, 1, -1            do l = llm, 1, -1
181               do j=1, jjm               do j = 1, jjm
182                  psiQ(j, l, iQ)=psiQ(j, l+1, iQ)+zvQ(j, l, itot, iQ)                  psiQ(j, l, iQ) = psiQ(j, l + 1, iQ) + vq(j, l, itot, iQ)
183               enddo               enddo
184            enddo            enddo
185         enddo         enddo
186    
187         ! fonction de courant pour la circulation meridienne moyenne         ! Fonction de courant pour la circulation méridienne moyenne
188         psi=0.         psi = 0.
189         do l=llm, 1, -1         do l = llm, 1, -1
190            do j=1, jjm            do j = 1, jjm
191               psi(j, l)=psi(j, l+1)+zv(j, l)               psi(j, l) = psi(j, l + 1) + v(j, l)
192               zv(j, l)=zv(j, l)*zfactv(j, l)               v(j, l) = v(j, l) * factv(j, l)
193            enddo            enddo
194         enddo         enddo
195    
196         ! sorties proprement dites         ! Sorties proprement dites
197         if (i_sortie == 1) then         do iQ = 1, nQ
198            do iQ=1, nQ            do itr = 1, ntr
199               do itr=1, ntr               call histwrite(fileid, znom(itr, iQ), itau, vq(:, :, itr, iQ))
                 call histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ))  
              enddo  
              call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, 1:llm, iQ))  
200            enddo            enddo
201              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)  
202         enddo         enddo
203         zavQ=0.  
204         do iQ=1, nQ         call histwrite(fileid, 'masse', itau, zmasse)
205            do itr=2, ntr         call histwrite(fileid, 'v', itau, v)
206               do l=1, llm         psi = psi * 1e-9
207                  zavQ(:, itr, iQ) = zavQ(:, itr, iQ) &         call histwrite(fileid, 'psi', itau, psi(:, :llm))
208                       + zvQ(:, l, itr, iQ) * zmasse(:, l)  
209               enddo         ! Intégrale verticale
210               zavQ(:, itr, iQ)=zavQ(:, itr, iQ)/zamasse(:)  
211               call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))         forall (iQ = 1: nQ, itr = 2: ntr) avq(:, itr, iQ) &
212                = sum(vq(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)
213    
214           do iQ = 1, nQ
215              do itr = 2, ntr
216                 call histwrite(fileid, 'a' // znom(itr, iQ), itau, avq(:, itr, iQ))
217            enddo            enddo
218         enddo         enddo
219    
220         ! on doit pouvoir tracer systematiquement la fonction de courant.         icum = 0
        icum=0  
221      endif writing_step      endif writing_step
222    
223    end SUBROUTINE bilan_dyn    end SUBROUTINE bilan_dyn
Legend:
Removed from v.40  
changed lines
  Added in v.134
  ViewVC Help
Powered by ViewVC 1.1.21