/[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 56 by guez, Tue Jan 10 19:02: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, dt_app)
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    
     USE histcom, ONLY: histbeg_totreg, histdef, histend, histvert  
16      USE calendar, ONLY: ymds2ju      USE calendar, ONLY: ymds2ju
17      USE histwrite_m, ONLY: histwrite      USE conf_gcm_m, ONLY: day_step, iperiod, periodav
     USE dimens_m, ONLY: iim, jjm, llm  
     USE paramet_m, ONLY: iip1, jjp1  
18      USE comconst, ONLY: cpp      USE comconst, ONLY: cpp
19      USE comvert, ONLY: presnivs      USE comvert, ONLY: presnivs
20      USE comgeom, ONLY: constang_2d, cu_2d, cv_2d, rlatv      USE comgeom, ONLY: constang_2d, cu_2d, cv_2d, rlatv
21      USE temps, ONLY: annee_ref, day_ref, itau_dyn      USE dimens_m, ONLY: iim, jjm, llm
22      USE inigrads_m, ONLY: inigrads      USE histcom, ONLY: histbeg_totreg, histdef, histend, histvert
23        USE histwrite_m, ONLY: histwrite
24      USE nr_util, ONLY: pi      USE nr_util, ONLY: pi
25        USE paramet_m, ONLY: iip1, jjp1
26        USE temps, ONLY: annee_ref, day_ref, itau_dyn
27    
28      ! Arguments:      ! Arguments:
29    
30      real, intent(in):: dt_app, dt_cum      real, intent(in):: ps(iip1, jjp1)
31      real ps(iip1, jjp1)      real, intent(in):: masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)
32      real masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)      real, intent(in):: flux_u(iip1, jjp1, llm)
33      real flux_u(iip1, jjp1, llm)      real, intent(in):: flux_v(iip1, jjm, llm)
     real flux_v(iip1, jjm, llm)  
34      real, intent(in):: teta(iip1, jjp1, llm)      real, intent(in):: teta(iip1, jjp1, llm)
35      real phi(iip1, jjp1, llm)      real, intent(in):: phi(iip1, jjp1, llm)
36      real ucov(iip1, jjp1, llm)      real, intent(in):: ucov(iip1, jjp1, llm)
37      real vcov(iip1, jjm, llm)      real, intent(in):: vcov(iip1, jjm, llm)
38      real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)      real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)
39        real, intent(in):: dt_app
40    
41      ! Local:      ! Local:
42    
43      integer, save:: icum, ncum      real dt_cum
44        integer:: icum  = 0
45        integer, save:: ncum
46      logical:: first = .true.      logical:: first = .true.
47      real zz, zqy, zfactv(jjm, llm)      real 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             'u   ', 'ovap', 'un  '/)
52        character(len=5), parameter:: unites(nQ) = (/'K    ', 'm2/s2', 'm2/s2', &
53             'ang  ', 'm/s  ', 'kg/kg', 'un   '/)
54    
     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.  
55      integer:: itau = 0      integer:: itau = 0
   
     data itemp, igeop, iecin, iang, iu, iovap, iun/1, 2, 3, 4, 5, 6, 7/  
   
56      real ww      real ww
57    
58      ! Variables dynamiques intermédiaires      ! Variables dynamiques intermédiaires
59      REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)      REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)
60      REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)      REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)
61      REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)      REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
     REAL vorpot(iip1, jjm, llm)  
62      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)  
63    
64      ! Champ contenant les scalaires advectés      ! Champ contenant les scalaires advectés
65      real Q(iip1, jjp1, llm, nQ)      real Q(iip1, jjp1, llm, nQ)
# Line 82  contains Line 72  contains
72      real, save:: Q_cum(iip1, jjp1, llm, nQ)      real, save:: Q_cum(iip1, jjp1, llm, nQ)
73      real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)      real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)
74      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)  
75      real dQ(iip1, jjp1, llm, nQ)      real dQ(iip1, jjp1, llm, nQ)
76    
77      ! champs de tansport en moyenne zonale      ! champs de tansport en moyenne zonale
78      integer itr      integer itr
79      integer, parameter:: ntr=5      integer, parameter:: ntr = 5
80    
81      character(len=10), save:: znom(ntr, nQ)      character(len=10), save:: znom(ntr, nQ)
82      character(len=20), save:: znoml(ntr, nQ)      character(len=26), save:: znoml(ntr, nQ)
83      character(len=10), save:: zunites(ntr, nQ)      character(len=12), save:: zunites(ntr, nQ)
84    
85      integer iave, itot, immc, itrs, istn      integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5
86      data iave, itot, immc, itrs, istn/1, 2, 3, 4, 5/      character(len=3), parameter:: ctrs(ntr) = (/'   ', 'TOT', 'MMC', 'TRS', &
87      character(len=3) ctrs(ntr)           'STN'/)
     data ctrs/' ', 'TOT', 'MMC', 'TRS', 'STN'/  
88    
89      real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)      real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)
90      real zavQ(jjm, ntr, nQ), psiQ(jjm, llm+1, nQ)      real zavQ(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)
91      real zmasse(jjm, llm), zamasse(jjm)      real zmasse(jjm, llm)
92        real zv(jjm, llm), psi(jjm, llm + 1)
     real zv(jjm, llm), psi(jjm, llm+1)  
   
93      integer i, j, l, iQ      integer i, j, l, iQ
94    
95      ! Initialisation du fichier contenant les moyennes zonales.      ! Initialisation du fichier contenant les moyennes zonales.
96    
97      integer, save:: fileid      integer, save:: fileid
98      integer thoriid, zvertiid      integer thoriid, zvertiid
     integer ndex3d(jjm*llm)  
   
     ! Variables locales  
99    
100      real zjulian      real zjulian
     character(len=3) str  
     character(len=10) ctrac  
     integer ii, jj  
101      integer zan, dayref      integer zan, dayref
102    
103      real rlong(jjm), rlatg(jjm)      real rlong(jjm), rlatg(jjm)
# Line 126  contains Line 106  contains
106    
107      !!print *, "Call sequence information: bilan_dyn"      !!print *, "Call sequence information: bilan_dyn"
108    
109      ! Initialisation      first_call: if (first) then
   
     time=time+dt_app  
     itau=itau+1  
   
     if (first) then  
        icum=0  
110         ! initialisation des fichiers         ! initialisation des fichiers
111         first=.false.         first = .false.
112         ! ncum est la frequence de stokage en pas de temps         ! ncum est la frequence de stokage en pas de temps
113         ncum=dt_cum/dt_app         ncum = day_step / iperiod * periodav
114         if (abs(ncum * dt_app - dt_cum) > 1e-5 * dt_app) then         dt_cum = ncum * dt_app
           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'  
115    
116         ! Initialisation du fichier contenant les moyennes zonales         ! Initialisation du fichier contenant les moyennes zonales
117    
# Line 172  contains Line 119  contains
119         dayref = day_ref         dayref = day_ref
120         CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)         CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
121    
122         rlong=0.         rlong = 0.
123         rlatg=rlatv*180./pi         rlatg = rlatv*180./pi
124    
125         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, &
126              zjulian, dt_cum, thoriid, fileid)              zjulian, dt_cum, thoriid, fileid)
# Line 184  contains Line 131  contains
131              zvertiid)              zvertiid)
132    
133         ! Appels à histdef pour la définition des variables à sauvegarder         ! Appels à histdef pour la définition des variables à sauvegarder
134         do iQ=1, nQ         do iQ = 1, nQ
135            do itr=1, ntr            do itr = 1, ntr
136               if(itr == 1) then               if (itr == 1) then
137                  znom(itr, iQ)=nom(iQ)                  znom(itr, iQ) = nom(iQ)
138                  znoml(itr, iQ)=nom(iQ)                  znoml(itr, iQ) = nom(iQ)
139                  zunites(itr, iQ)=unites(iQ)                  zunites(itr, iQ) = unites(iQ)
140               else               else
141                  znom(itr, iQ)=ctrs(itr)//'v'//nom(iQ)                  znom(itr, iQ) = ctrs(itr)//'v'//nom(iQ)
142                  znoml(itr, iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)                  znoml(itr, iQ) = 'transport : v * '//nom(iQ)//' '//ctrs(itr)
143                  zunites(itr, iQ)='m/s * '//unites(iQ)                  zunites(itr, iQ) = 'm/s * '//unites(iQ)
144               endif               endif
145            enddo            enddo
146         enddo         enddo
147    
148         ! Déclarations des champs avec dimension verticale         ! Déclarations des champs avec dimension verticale
149         do iQ=1, nQ         do iQ = 1, nQ
150            do itr=1, ntr            do itr = 1, ntr
151               call histdef(fileid, znom(itr, iQ), znoml(itr, iQ), &               call histdef(fileid, znom(itr, iQ), znoml(itr, iQ), &
152                    zunites(itr, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &                    zunites(itr, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &
153                    'ave(X)', dt_cum, dt_cum)                    'ave(X)', dt_cum, dt_cum)
154            enddo            enddo
155            ! Declarations pour les fonctions de courant            ! Déclarations pour les fonctions de courant
156            call histdef(fileid, 'psi'//nom(iQ), 'stream fn. '//znoml(itot, iQ), &            call histdef(fileid, 'psi'//nom(iQ), 'stream fn. '//znoml(itot, iQ), &
157                 zunites(itot, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &                 zunites(itot, iQ), 1, jjm, thoriid, llm, 1, llm, zvertiid, &
158                 'ave(X)', dt_cum, dt_cum)                 'ave(X)', dt_cum, dt_cum)
159         enddo         enddo
160    
161         ! Declarations pour les champs de transport d'air         ! Déclarations pour les champs de transport d'air
162         call histdef(fileid, 'masse', 'masse', &         call histdef(fileid, 'masse', 'masse', &
163              'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid, &              'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
164              'ave(X)', dt_cum, dt_cum)              'ave(X)', dt_cum, dt_cum)
165         call histdef(fileid, 'v', 'v', &         call histdef(fileid, 'v', 'v', &
166              'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid, &              'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid, &
167              'ave(X)', dt_cum, dt_cum)              'ave(X)', dt_cum, dt_cum)
168         ! Declarations pour les fonctions de courant         ! Déclarations pour les fonctions de courant
169         call histdef(fileid, 'psi', 'stream fn. MMC ', 'mega t/s', &         call histdef(fileid, 'psi', 'stream fn. MMC ', 'mega t/s', &
170              1, jjm, thoriid, llm, 1, llm, zvertiid, &              1, jjm, thoriid, llm, 1, llm, zvertiid, &
171              'ave(X)', dt_cum, dt_cum)              'ave(X)', dt_cum, dt_cum)
172    
173         ! Declaration des champs 1D de transport en latitude         ! Déclaration des champs 1D de transport en latitude
174         do iQ=1, nQ         do iQ = 1, nQ
175            do itr=2, ntr            do itr = 2, ntr
176               call histdef(fileid, 'a'//znom(itr, iQ), znoml(itr, iQ), &               call histdef(fileid, 'a'//znom(itr, iQ), znoml(itr, iQ), &
177                    zunites(itr, iQ), 1, jjm, thoriid, 1, 1, 1, -99, &                    zunites(itr, iQ), 1, jjm, thoriid, 1, 1, 1, -99, &
178                    'ave(X)', dt_cum, dt_cum)                    'ave(X)', dt_cum, dt_cum)
# Line 233  contains Line 180  contains
180         enddo         enddo
181    
182         CALL histend(fileid)         CALL histend(fileid)
183      endif      endif first_call
184    
185        itau = itau + 1
186    
187      ! Calcul des champs dynamiques      ! Calcul des champs dynamiques
188    
# Line 243  contains Line 192  contains
192      CALL enercin(vcov, ucov, vcont, ucont, ecin)      CALL enercin(vcov, ucov, vcont, ucont, ecin)
193    
194      ! moment cinétique      ! moment cinétique
195      do l=1, llm      do l = 1, llm
196         ang(:, :, l)=ucov(:, :, l)+constang_2d         ang(:, :, l) = ucov(:, :, l) + constang_2d
197         unat(:, :, l)=ucont(:, :, l)*cu_2d         unat(:, :, l) = ucont(:, :, l)*cu_2d
198      enddo      enddo
199    
200      Q(:, :, :, itemp)=teta*pk/cpp      Q(:, :, :, 1) = teta * pk / cpp
201      Q(:, :, :, igeop)=phi      Q(:, :, :, 2) = phi
202      Q(:, :, :, iecin)=ecin      Q(:, :, :, 3) = ecin
203      Q(:, :, :, iang)=ang      Q(:, :, :, 4) = ang
204      Q(:, :, :, iu)=unat      Q(:, :, :, 5) = unat
205      Q(:, :, :, iovap)=trac      Q(:, :, :, 6) = trac
206      Q(:, :, :, iun)=1.      Q(:, :, :, 7) = 1.
207    
208      ! Cumul      ! Cumul
209    
210      if(icum == 0) then      if (icum == 0) then
211         ps_cum=0.         ps_cum = 0.
212         masse_cum=0.         masse_cum = 0.
213         flux_u_cum=0.         flux_u_cum = 0.
214         flux_v_cum=0.         flux_v_cum = 0.
215         Q_cum=0.         Q_cum = 0.
216         flux_vQ_cum=0.         flux_vQ_cum = 0.
217         flux_uQ_cum=0.         flux_uQ_cum = 0.
218      endif      endif
219    
220      icum=icum+1      icum = icum + 1
221    
222      ! Accumulation des flux de masse horizontaux      ! Accumulation des flux de masse horizontaux
223      ps_cum=ps_cum+ps      ps_cum = ps_cum + ps
224      masse_cum=masse_cum+masse      masse_cum = masse_cum + masse
225      flux_u_cum=flux_u_cum+flux_u      flux_u_cum = flux_u_cum + flux_u
226      flux_v_cum=flux_v_cum+flux_v      flux_v_cum = flux_v_cum + flux_v
227      do iQ=1, nQ      do iQ = 1, nQ
228         Q_cum(:, :, :, iQ)=Q_cum(:, :, :, iQ)+Q(:, :, :, iQ)*masse         Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) + Q(:, :, :, iQ)*masse
229      enddo      enddo
230    
231      ! FLUX ET TENDANCES      ! FLUX ET TENDANCES
232    
233      ! Flux longitudinal      ! Flux longitudinal
234      do iQ=1, nQ      forall (iQ = 1: nQ, i = 1: iim) flux_uQ_cum(i, :, :, iQ) &
235         do l=1, llm           = flux_uQ_cum(i, :, :, iQ) &
236            do j=1, jjp1           + flux_u(i, :, :) * 0.5 * (Q(i, :, :, iQ) + Q(i + 1, :, :, iQ))
237               do i=1, iim      flux_uQ_cum(iip1, :, :, :) = flux_uQ_cum(1, :, :, :)
238                  flux_uQ_cum(i, j, l, iQ)=flux_uQ_cum(i, j, l, iQ) &  
239                       +flux_u(i, j, l)*0.5*(Q(i, j, l, iQ)+Q(i+1, j, l, iQ))      ! Flux méridien
240               enddo      forall (iQ = 1: nQ, j = 1: jjm) flux_vQ_cum(:, j, :, iQ) &
241               flux_uQ_cum(iip1, j, l, iQ)=flux_uQ_cum(1, j, l, iQ)           = flux_vQ_cum(:, j, :, iQ) &
242            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  
243    
244      ! tendances      ! tendances
245    
# Line 315  contains Line 250  contains
250      call convmas(flux_u_cum, flux_v_cum, convm)      call convmas(flux_u_cum, flux_v_cum, convm)
251      CALL vitvert(convm, w)      CALL vitvert(convm, w)
252    
253      do iQ=1, nQ      do iQ = 1, nQ
254         do l=1, llm-1         do l = 1, llm-1
255            do j=1, jjp1            do j = 1, jjp1
256               do i=1, iip1               do i = 1, iip1
257                  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))
258                  dQ(i, j, l , iQ)=dQ(i, j, l , iQ)-ww                  dQ(i, j, l, iQ) = dQ(i, j, l, iQ)-ww
259                  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
260               enddo               enddo
261            enddo            enddo
262         enddo         enddo
# Line 331  contains Line 266  contains
266    
267      writing_step: if (icum == ncum) then      writing_step: if (icum == ncum) then
268         ! Normalisation         ! Normalisation
269         do iQ=1, nQ         do iQ = 1, nQ
270            Q_cum(:, :, :, iQ)=Q_cum(:, :, :, iQ)/masse_cum            Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ)/masse_cum
271         enddo         enddo
272         zz=1./float(ncum)         ps_cum = ps_cum / ncum
273         ps_cum=ps_cum*zz         masse_cum = masse_cum / ncum
274         masse_cum=masse_cum*zz         flux_u_cum = flux_u_cum / ncum
275         flux_u_cum=flux_u_cum*zz         flux_v_cum = flux_v_cum / ncum
276         flux_v_cum=flux_v_cum*zz         flux_uQ_cum = flux_uQ_cum / ncum
277         flux_uQ_cum=flux_uQ_cum*zz         flux_vQ_cum = flux_vQ_cum / ncum
278         flux_vQ_cum=flux_vQ_cum*zz         dQ = dQ / ncum
        dQ=dQ*zz  
279    
280         ! A retravailler eventuellement         ! A retravailler eventuellement
281         ! division de dQ par la masse pour revenir aux bonnes grandeurs         ! division de dQ par la masse pour revenir aux bonnes grandeurs
282         do iQ=1, nQ         do iQ = 1, nQ
283            dQ(:, :, :, iQ)=dQ(:, :, :, iQ)/masse_cum            dQ(:, :, :, iQ) = dQ(:, :, :, iQ)/masse_cum
284         enddo         enddo
285    
286         ! Transport méridien         ! Transport méridien
287    
288         ! cumul zonal des masses des mailles         ! cumul zonal des masses des mailles
289    
290         zv=0.         zv = 0.
291         zmasse=0.         zmasse = 0.
292         call massbar(masse_cum, massebx, masseby)         call massbar(masse_cum, massebx, masseby)
293         do l=1, llm         do l = 1, llm
294            do j=1, jjm            do j = 1, jjm
295               do i=1, iim               do i = 1, iim
296                  zmasse(j, l)=zmasse(j, l)+masseby(i, j, l)                  zmasse(j, l) = zmasse(j, l) + masseby(i, j, l)
297                  zv(j, l)=zv(j, l)+flux_v_cum(i, j, l)                  zv(j, l) = zv(j, l) + flux_v_cum(i, j, l)
298               enddo               enddo
299               zfactv(j, l)=cv_2d(1, j)/zmasse(j, l)               zfactv(j, l) = cv_2d(1, j)/zmasse(j, l)
300            enddo            enddo
301         enddo         enddo
302    
303         ! Transport dans le plan latitude-altitude         ! Transport dans le plan latitude-altitude
304    
305         zvQ=0.         zvQ = 0.
306         psiQ=0.         psiQ = 0.
307         do iQ=1, nQ         do iQ = 1, nQ
308            zvQtmp=0.            zvQtmp = 0.
309            do l=1, llm            do l = 1, llm
310               do j=1, jjm               do j = 1, jjm
311                  ! Calcul des moyennes zonales du transort total et de zvQtmp                  ! Calcul des moyennes zonales du transort total et de zvQtmp
312                  do i=1, iim                  do i = 1, iim
313                     zvQ(j, l, itot, iQ)=zvQ(j, l, itot, iQ) &                     zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ) &
314                          +flux_vQ_cum(i, j, l, iQ)                          + flux_vQ_cum(i, j, l, iQ)
315                     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) &
316                          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))
317                     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 &
318                          /(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)))
319                     zvQ(j, l, iave, iQ)=zvQ(j, l, iave, iQ)+zqy                     zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ) + zqy
320                  enddo                  enddo
321                  ! Decomposition                  ! Decomposition
322                  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)
323                  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)
324                  zvQtmp(j, l)=zvQtmp(j, l)*zfactv(j, l)                  zvQtmp(j, l) = zvQtmp(j, l)*zfactv(j, l)
325                  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)
326                  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)
327                  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)
328               enddo               enddo
329            enddo            enddo
330            ! fonction de courant meridienne pour la quantite Q            ! fonction de courant meridienne pour la quantite Q
331            do l=llm, 1, -1            do l = llm, 1, -1
332               do j=1, jjm               do j = 1, jjm
333                  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)
334               enddo               enddo
335            enddo            enddo
336         enddo         enddo
337    
338         ! fonction de courant pour la circulation meridienne moyenne         ! fonction de courant pour la circulation meridienne moyenne
339         psi=0.         psi = 0.
340         do l=llm, 1, -1         do l = llm, 1, -1
341            do j=1, jjm            do j = 1, jjm
342               psi(j, l)=psi(j, l+1)+zv(j, l)               psi(j, l) = psi(j, l + 1) + zv(j, l)
343               zv(j, l)=zv(j, l)*zfactv(j, l)               zv(j, l) = zv(j, l)*zfactv(j, l)
344            enddo            enddo
345         enddo         enddo
346    
347         ! sorties proprement dites         ! sorties proprement dites
348         if (i_sortie == 1) then         do iQ = 1, nQ
349            do iQ=1, nQ            do itr = 1, ntr
350               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))  
351            enddo            enddo
352              call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, :llm, iQ))
353           enddo
354    
355            call histwrite(fileid, 'masse', itau, zmasse)         call histwrite(fileid, 'masse', itau, zmasse)
356            call histwrite(fileid, 'v', itau, zv)         call histwrite(fileid, 'v', itau, zv)
357            psi=psi*1.e-9         psi = psi*1.e-9
358            call histwrite(fileid, 'psi', itau, psi(:, 1:llm))         call histwrite(fileid, 'psi', itau, psi(:, :llm))
359         endif  
360           ! Intégrale verticale
361         ! Moyenne verticale  
362           forall (iQ = 1: nQ, itr = 2: ntr) zavQ(:, itr, iQ) &
363         zamasse=0.              = sum(zvQ(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)
364         do l=1, llm  
365            zamasse(:)=zamasse(:)+zmasse(:, l)         do iQ = 1, nQ
366         enddo            do itr = 2, ntr
        zavQ=0.  
        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(:)  
367               call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))               call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))
368            enddo            enddo
369         enddo         enddo
370    
371         ! on doit pouvoir tracer systematiquement la fonction de courant.         ! On doit pouvoir tracer systematiquement la fonction de courant.
372         icum=0         icum = 0
373      endif writing_step      endif writing_step
374    
375    end SUBROUTINE bilan_dyn    end SUBROUTINE bilan_dyn
Legend:
Removed from v.44  
changed lines
  Added in v.56
  ViewVC Help
Powered by ViewVC 1.1.21