/[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 54 by guez, Tue Dec 6 15:07:04 2011 UTC
# Line 12  contains Line 12  contains
12    
13      ! Sous-programme consacré à des diagnostics dynamiques de base      ! Sous-programme consacré à des diagnostics dynamiques de base
14      ! 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 par
15      ! la masse. Les flux de masse sont eux simplement moyennés.      ! la masse. Les flux de masse sont, eux, simplement moyennés.
16    
17      USE histcom, ONLY: histbeg_totreg, histdef, histend, histvert      USE histcom, ONLY: histbeg_totreg, histdef, histend, histvert
18      USE calendar, ONLY: ymds2ju      USE calendar, ONLY: ymds2ju
# Line 41  contains Line 41  contains
41    
42      ! Local:      ! Local:
43    
44      integer, save:: icum, ncum      integer:: icum  = 0
45        integer, save:: ncum
46      logical:: first = .true.      logical:: first = .true.
47      real zz, zqy, zfactv(jjm, llm)      real zz, zqy, zfactv(jjm, llm)
48    
49      integer, parameter:: nQ=7      integer, parameter:: nQ = 7
50        character(len=4), parameter:: nom(nQ) = (/'T   ', 'gz  ', 'K   ', 'ang ', &
51      character(len=6), save:: nom(nQ)           'u   ', 'ovap', 'un  '/)
52      character(len=6), save:: unites(nQ)      character(len=5), parameter:: unites(nQ) = (/'K    ', 'm2/s2', 'm2/s2', &
53             'ang  ', 'm/s  ', 'kg/kg', 'un   '/)
     character(len=10) file  
     integer, parameter:: ifile=4  
   
     integer itemp, igeop, iecin, iang, iu, iovap, iun  
     integer:: i_sortie = 1  
54    
55      real:: time = 0.      real:: time = 0.
56      integer:: itau = 0      integer:: itau = 0
   
     data itemp, igeop, iecin, iang, iu, iovap, iun/1, 2, 3, 4, 5, 6, 7/  
   
57      real ww      real ww
58    
59      ! Variables dynamiques intermédiaires      ! Variables dynamiques intermédiaires
60      REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)      REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)
61      REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)      REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)
62      REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)      REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
     REAL vorpot(iip1, jjm, llm)  
63      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)  
64    
65      ! Champ contenant les scalaires advectés      ! Champ contenant les scalaires advectés
66      real Q(iip1, jjp1, llm, nQ)      real Q(iip1, jjp1, llm, nQ)
# Line 82  contains Line 73  contains
73      real, save:: Q_cum(iip1, jjp1, llm, nQ)      real, save:: Q_cum(iip1, jjp1, llm, nQ)
74      real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)      real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)
75      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)  
76      real dQ(iip1, jjp1, llm, nQ)      real dQ(iip1, jjp1, llm, nQ)
77    
78      ! champs de tansport en moyenne zonale      ! champs de tansport en moyenne zonale
79      integer itr      integer itr
80      integer, parameter:: ntr=5      integer, parameter:: ntr = 5
81    
82      character(len=10), save:: znom(ntr, nQ)      character(len=10), save:: znom(ntr, nQ)
83      character(len=20), save:: znoml(ntr, nQ)      character(len=26), save:: znoml(ntr, nQ)
84      character(len=10), save:: zunites(ntr, nQ)      character(len=12), save:: zunites(ntr, nQ)
85    
86      integer iave, itot, immc, itrs, istn      integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5
87      data iave, itot, immc, itrs, istn/1, 2, 3, 4, 5/      character(len=3), parameter:: ctrs(ntr) = (/'   ', 'TOT', 'MMC', 'TRS', &
88      character(len=3) ctrs(ntr)           'STN'/)
     data ctrs/' ', 'TOT', 'MMC', 'TRS', 'STN'/  
89    
90      real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)      real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)
91      real zavQ(jjm, ntr, nQ), psiQ(jjm, llm+1, nQ)      real zavQ(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)
92      real zmasse(jjm, llm), zamasse(jjm)      real zmasse(jjm, llm)
93    
94      real zv(jjm, llm), psi(jjm, llm+1)      real zv(jjm, llm), psi(jjm, llm + 1)
95    
96      integer i, j, l, iQ      integer i, j, l, iQ
97    
# Line 110  contains Line 99  contains
99    
100      integer, save:: fileid      integer, save:: fileid
101      integer thoriid, zvertiid      integer thoriid, zvertiid
     integer ndex3d(jjm*llm)  
   
     ! Variables locales  
102    
103      real zjulian      real zjulian
     character(len=3) str  
     character(len=10) ctrac  
     integer ii, jj  
104      integer zan, dayref      integer zan, dayref
105    
106      real rlong(jjm), rlatg(jjm)      real rlong(jjm), rlatg(jjm)
# Line 128  contains Line 111  contains
111    
112      ! Initialisation      ! Initialisation
113    
114      time=time+dt_app      time = time + dt_app
115      itau=itau+1      itau = itau + 1
116    
117      if (first) then      first_call: if (first) then
        icum=0  
118         ! initialisation des fichiers         ! initialisation des fichiers
119         first=.false.         first = .false.
120         ! ncum est la frequence de stokage en pas de temps         ! ncum est la frequence de stokage en pas de temps
121         ncum=dt_cum/dt_app         ncum = dt_cum / dt_app
122         if (abs(ncum * dt_app - dt_cum) > 1e-5 * dt_app) then         if (abs(ncum * dt_app - dt_cum) > 1e-5 * dt_app) then
123            print *, 'Problème : le pas de cumul doit être multiple du pas'            print *, 'Problème : le pas de cumul doit être multiple du pas'
124            print *, 'dt_app=', dt_app            print *, 'dt_app = ', dt_app
125            print *, 'dt_cum=', dt_cum            print *, 'dt_cum = ', dt_cum
126            stop 1            stop 1
127         endif         endif
128    
129         if (i_sortie == 1) then         call inigrads(i_f=4, x=(/0./), fx=180./pi, xmin=0., xmax=0., y=rlatv, &
130            file='dynzon'              ymin=-90., ymax=90., fy=180./pi, z=presnivs, fz=1., dt=dt_cum, &
131            call inigrads(ifile , (/0./), 180./pi, 0., 0., rlatv, -90., 90., &              file='dynzon', titlel='dyn_zon ')
                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'  
132    
133         ! Initialisation du fichier contenant les moyennes zonales         ! Initialisation du fichier contenant les moyennes zonales
134    
# Line 172  contains Line 136  contains
136         dayref = day_ref         dayref = day_ref
137         CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)         CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
138    
139         rlong=0.         rlong = 0.
140         rlatg=rlatv*180./pi         rlatg = rlatv*180./pi
141    
142         call histbeg_totreg('dynzon', rlong(:1), rlatg, 1, 1, 1, jjm, itau_dyn, &         call histbeg_totreg('dynzon', rlong(:1), rlatg, 1, 1, 1, jjm, itau_dyn, &
143              zjulian, dt_cum, thoriid, fileid)              zjulian, dt_cum, thoriid, fileid)
# Line 184  contains Line 148  contains
148              zvertiid)              zvertiid)
149    
150         ! Appels à histdef pour la définition des variables à sauvegarder         ! Appels à histdef pour la définition des variables à sauvegarder
151         do iQ=1, nQ         do iQ = 1, nQ
152            do itr=1, ntr            do itr = 1, ntr
153               if(itr == 1) then               if (itr == 1) then
154                  znom(itr, iQ)=nom(iQ)                  znom(itr, iQ) = nom(iQ)
155                  znoml(itr, iQ)=nom(iQ)                  znoml(itr, iQ) = nom(iQ)
156                  zunites(itr, iQ)=unites(iQ)                  zunites(itr, iQ) = unites(iQ)
157               else               else
158                  znom(itr, iQ)=ctrs(itr)//'v'//nom(iQ)                  znom(itr, iQ) = ctrs(itr)//'v'//nom(iQ)
159                  znoml(itr, iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)                  znoml(itr, iQ) = 'transport : v * '//nom(iQ)//' '//ctrs(itr)
160                  zunites(itr, iQ)='m/s * '//unites(iQ)                  zunites(itr, iQ) = 'm/s * '//unites(iQ)
161               endif               endif
162            enddo            enddo
163         enddo         enddo
164    
165         ! Déclarations des champs avec dimension verticale         ! Déclarations des champs avec dimension verticale
166         do iQ=1, nQ         do iQ = 1, nQ
167            do itr=1, ntr            do itr = 1, ntr
168               call histdef(fileid, znom(itr, iQ), znoml(itr, iQ), &               call histdef(fileid, znom(itr, iQ), znoml(itr, iQ), &
169                    zunites(itr, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &                    zunites(itr, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &
170                    'ave(X)', dt_cum, dt_cum)                    'ave(X)', dt_cum, dt_cum)
# Line 224  contains Line 188  contains
188              'ave(X)', dt_cum, dt_cum)              'ave(X)', dt_cum, dt_cum)
189    
190         ! Declaration des champs 1D de transport en latitude         ! Declaration des champs 1D de transport en latitude
191         do iQ=1, nQ         do iQ = 1, nQ
192            do itr=2, ntr            do itr = 2, ntr
193               call histdef(fileid, 'a'//znom(itr, iQ), znoml(itr, iQ), &               call histdef(fileid, 'a'//znom(itr, iQ), znoml(itr, iQ), &
194                    zunites(itr, iQ), 1, jjm, thoriid, 1, 1, 1, -99, &                    zunites(itr, iQ), 1, jjm, thoriid, 1, 1, 1, -99, &
195                    'ave(X)', dt_cum, dt_cum)                    'ave(X)', dt_cum, dt_cum)
# Line 233  contains Line 197  contains
197         enddo         enddo
198    
199         CALL histend(fileid)         CALL histend(fileid)
200      endif      endif first_call
201    
202      ! Calcul des champs dynamiques      ! Calcul des champs dynamiques
203    
# Line 243  contains Line 207  contains
207      CALL enercin(vcov, ucov, vcont, ucont, ecin)      CALL enercin(vcov, ucov, vcont, ucont, ecin)
208    
209      ! moment cinétique      ! moment cinétique
210      do l=1, llm      do l = 1, llm
211         ang(:, :, l)=ucov(:, :, l)+constang_2d         ang(:, :, l) = ucov(:, :, l) + constang_2d
212         unat(:, :, l)=ucont(:, :, l)*cu_2d         unat(:, :, l) = ucont(:, :, l)*cu_2d
213      enddo      enddo
214    
215      Q(:, :, :, itemp)=teta*pk/cpp      Q(:, :, :, 1) = teta * pk / cpp
216      Q(:, :, :, igeop)=phi      Q(:, :, :, 2) = phi
217      Q(:, :, :, iecin)=ecin      Q(:, :, :, 3) = ecin
218      Q(:, :, :, iang)=ang      Q(:, :, :, 4) = ang
219      Q(:, :, :, iu)=unat      Q(:, :, :, 5) = unat
220      Q(:, :, :, iovap)=trac      Q(:, :, :, 6) = trac
221      Q(:, :, :, iun)=1.      Q(:, :, :, 7) = 1.
222    
223      ! Cumul      ! Cumul
224    
225      if(icum == 0) then      if (icum == 0) then
226         ps_cum=0.         ps_cum = 0.
227         masse_cum=0.         masse_cum = 0.
228         flux_u_cum=0.         flux_u_cum = 0.
229         flux_v_cum=0.         flux_v_cum = 0.
230         Q_cum=0.         Q_cum = 0.
231         flux_vQ_cum=0.         flux_vQ_cum = 0.
232         flux_uQ_cum=0.         flux_uQ_cum = 0.
233      endif      endif
234    
235      icum=icum+1      icum = icum + 1
236    
237      ! Accumulation des flux de masse horizontaux      ! Accumulation des flux de masse horizontaux
238      ps_cum=ps_cum+ps      ps_cum = ps_cum + ps
239      masse_cum=masse_cum+masse      masse_cum = masse_cum + masse
240      flux_u_cum=flux_u_cum+flux_u      flux_u_cum = flux_u_cum + flux_u
241      flux_v_cum=flux_v_cum+flux_v      flux_v_cum = flux_v_cum + flux_v
242      do iQ=1, nQ      do iQ = 1, nQ
243         Q_cum(:, :, :, iQ)=Q_cum(:, :, :, iQ)+Q(:, :, :, iQ)*masse         Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) + Q(:, :, :, iQ)*masse
244      enddo      enddo
245    
246      ! FLUX ET TENDANCES      ! FLUX ET TENDANCES
247    
248      ! Flux longitudinal      ! Flux longitudinal
249      do iQ=1, nQ      forall (iQ = 1: nQ, i = 1: iim) flux_uQ_cum(i, :, :, iQ) &
250         do l=1, llm           = flux_uQ_cum(i, :, :, iQ) &
251            do j=1, jjp1           + flux_u(i, :, :) * 0.5 * (Q(i, :, :, iQ) + Q(i + 1, :, :, iQ))
252               do i=1, iim      flux_uQ_cum(iip1, :, :, :) = flux_uQ_cum(1, :, :, :)
253                  flux_uQ_cum(i, j, l, iQ)=flux_uQ_cum(i, j, l, iQ) &  
254                       +flux_u(i, j, l)*0.5*(Q(i, j, l, iQ)+Q(i+1, j, l, iQ))      ! Flux méridien
255               enddo      forall (iQ = 1: nQ, j = 1: jjm) flux_vQ_cum(:, j, :, iQ) &
256               flux_uQ_cum(iip1, j, l, iQ)=flux_uQ_cum(1, j, l, iQ)           = flux_vQ_cum(:, j, :, iQ) &
257            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  
258    
259      ! tendances      ! tendances
260    
# Line 315  contains Line 265  contains
265      call convmas(flux_u_cum, flux_v_cum, convm)      call convmas(flux_u_cum, flux_v_cum, convm)
266      CALL vitvert(convm, w)      CALL vitvert(convm, w)
267    
268      do iQ=1, nQ      do iQ = 1, nQ
269         do l=1, llm-1         do l = 1, llm-1
270            do j=1, jjp1            do j = 1, jjp1
271               do i=1, iip1               do i = 1, iip1
272                  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))
273                  dQ(i, j, l , iQ)=dQ(i, j, l , iQ)-ww                  dQ(i, j, l, iQ) = dQ(i, j, l, iQ)-ww
274                  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
275               enddo               enddo
276            enddo            enddo
277         enddo         enddo
# Line 331  contains Line 281  contains
281    
282      writing_step: if (icum == ncum) then      writing_step: if (icum == ncum) then
283         ! Normalisation         ! Normalisation
284         do iQ=1, nQ         do iQ = 1, nQ
285            Q_cum(:, :, :, iQ)=Q_cum(:, :, :, iQ)/masse_cum            Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ)/masse_cum
286         enddo         enddo
287         zz=1./float(ncum)         zz = 1. / real(ncum)
288         ps_cum=ps_cum*zz         ps_cum = ps_cum*zz
289         masse_cum=masse_cum*zz         masse_cum = masse_cum*zz
290         flux_u_cum=flux_u_cum*zz         flux_u_cum = flux_u_cum*zz
291         flux_v_cum=flux_v_cum*zz         flux_v_cum = flux_v_cum*zz
292         flux_uQ_cum=flux_uQ_cum*zz         flux_uQ_cum = flux_uQ_cum*zz
293         flux_vQ_cum=flux_vQ_cum*zz         flux_vQ_cum = flux_vQ_cum*zz
294         dQ=dQ*zz         dQ = dQ*zz
295    
296         ! A retravailler eventuellement         ! A retravailler eventuellement
297         ! division de dQ par la masse pour revenir aux bonnes grandeurs         ! division de dQ par la masse pour revenir aux bonnes grandeurs
298         do iQ=1, nQ         do iQ = 1, nQ
299            dQ(:, :, :, iQ)=dQ(:, :, :, iQ)/masse_cum            dQ(:, :, :, iQ) = dQ(:, :, :, iQ)/masse_cum
300         enddo         enddo
301    
302         ! Transport méridien         ! Transport méridien
303    
304         ! cumul zonal des masses des mailles         ! cumul zonal des masses des mailles
305    
306         zv=0.         zv = 0.
307         zmasse=0.         zmasse = 0.
308         call massbar(masse_cum, massebx, masseby)         call massbar(masse_cum, massebx, masseby)
309         do l=1, llm         do l = 1, llm
310            do j=1, jjm            do j = 1, jjm
311               do i=1, iim               do i = 1, iim
312                  zmasse(j, l)=zmasse(j, l)+masseby(i, j, l)                  zmasse(j, l) = zmasse(j, l) + masseby(i, j, l)
313                  zv(j, l)=zv(j, l)+flux_v_cum(i, j, l)                  zv(j, l) = zv(j, l) + flux_v_cum(i, j, l)
314               enddo               enddo
315               zfactv(j, l)=cv_2d(1, j)/zmasse(j, l)               zfactv(j, l) = cv_2d(1, j)/zmasse(j, l)
316            enddo            enddo
317         enddo         enddo
318    
319         ! Transport dans le plan latitude-altitude         ! Transport dans le plan latitude-altitude
320    
321         zvQ=0.         zvQ = 0.
322         psiQ=0.         psiQ = 0.
323         do iQ=1, nQ         do iQ = 1, nQ
324            zvQtmp=0.            zvQtmp = 0.
325            do l=1, llm            do l = 1, llm
326               do j=1, jjm               do j = 1, jjm
327                  ! Calcul des moyennes zonales du transort total et de zvQtmp                  ! Calcul des moyennes zonales du transort total et de zvQtmp
328                  do i=1, iim                  do i = 1, iim
329                     zvQ(j, l, itot, iQ)=zvQ(j, l, itot, iQ) &                     zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ) &
330                          +flux_vQ_cum(i, j, l, iQ)                          + flux_vQ_cum(i, j, l, iQ)
331                     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) &
332                          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))
333                     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 &
334                          /(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)))
335                     zvQ(j, l, iave, iQ)=zvQ(j, l, iave, iQ)+zqy                     zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ) + zqy
336                  enddo                  enddo
337                  ! Decomposition                  ! Decomposition
338                  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)
339                  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)
340                  zvQtmp(j, l)=zvQtmp(j, l)*zfactv(j, l)                  zvQtmp(j, l) = zvQtmp(j, l)*zfactv(j, l)
341                  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)
342                  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)
343                  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)
344               enddo               enddo
345            enddo            enddo
346            ! fonction de courant meridienne pour la quantite Q            ! fonction de courant meridienne pour la quantite Q
347            do l=llm, 1, -1            do l = llm, 1, -1
348               do j=1, jjm               do j = 1, jjm
349                  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)
350               enddo               enddo
351            enddo            enddo
352         enddo         enddo
353    
354         ! fonction de courant pour la circulation meridienne moyenne         ! fonction de courant pour la circulation meridienne moyenne
355         psi=0.         psi = 0.
356         do l=llm, 1, -1         do l = llm, 1, -1
357            do j=1, jjm            do j = 1, jjm
358               psi(j, l)=psi(j, l+1)+zv(j, l)               psi(j, l) = psi(j, l + 1) + zv(j, l)
359               zv(j, l)=zv(j, l)*zfactv(j, l)               zv(j, l) = zv(j, l)*zfactv(j, l)
360            enddo            enddo
361         enddo         enddo
362    
363         ! sorties proprement dites         ! sorties proprement dites
364         if (i_sortie == 1) then         do iQ = 1, nQ
365            do iQ=1, nQ            do itr = 1, ntr
366               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))  
367            enddo            enddo
368              call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, :llm, iQ))
369           enddo
370    
371            call histwrite(fileid, 'masse', itau, zmasse)         call histwrite(fileid, 'masse', itau, zmasse)
372            call histwrite(fileid, 'v', itau, zv)         call histwrite(fileid, 'v', itau, zv)
373            psi=psi*1.e-9         psi = psi*1.e-9
374            call histwrite(fileid, 'psi', itau, psi(:, 1:llm))         call histwrite(fileid, 'psi', itau, psi(:, :llm))
        endif  
375    
376         ! Moyenne verticale         ! Moyenne verticale
377    
378         zamasse=0.         forall (iQ = 1: nQ, itr = 2: ntr) zavQ(:, itr, iQ) &
379         do l=1, llm              = sum(zvQ(:, :, itr, iQ) * zmasse, dim=2) / sum(zmasse, dim=2)
380            zamasse(:)=zamasse(:)+zmasse(:, l)  
381         enddo         do iQ = 1, nQ
382         zavQ=0.            do itr = 2, ntr
        do iQ=1, nQ  
           do itr=2, ntr  
              do l=1, llm  
                 zavQ(:, itr, iQ) = zavQ(:, itr, iQ) &  
                      + zvQ(:, l, itr, iQ) * zmasse(:, l)  
              enddo  
              zavQ(:, itr, iQ)=zavQ(:, itr, iQ)/zamasse(:)  
383               call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))               call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))
384            enddo            enddo
385         enddo         enddo
386    
387         ! on doit pouvoir tracer systematiquement la fonction de courant.         ! On doit pouvoir tracer systematiquement la fonction de courant.
388         icum=0         icum = 0
389      endif writing_step      endif writing_step
390    
391    end SUBROUTINE bilan_dyn    end SUBROUTINE bilan_dyn
Legend:
Removed from v.44  
changed lines
  Added in v.54
  ViewVC Help
Powered by ViewVC 1.1.21