/[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 57 by guez, Mon Jan 30 12:54:02 2012 UTC revision 62 by guez, Thu Jul 26 14:37:37 2012 UTC
# Line 20  contains Line 20  contains
20      use init_dynzon_m, only: ncum, fileid, znom, ntr, nq, nom      use init_dynzon_m, only: ncum, fileid, znom, ntr, nq, nom
21      USE paramet_m, ONLY: iip1, jjp1      USE paramet_m, ONLY: iip1, jjp1
22    
     ! Arguments:  
   
23      real, intent(in):: ps(iip1, jjp1)      real, intent(in):: ps(iip1, jjp1)
24      real, intent(in):: masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)      real, intent(in):: masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)
25      real, intent(in):: flux_u(iip1, jjp1, llm)      real, intent(in):: flux_u(iip1, jjp1, llm)
26      real, intent(in):: flux_v(iip1, jjm, llm)      real, intent(in):: flux_v(iip1, jjm, llm)
27      real, intent(in):: teta(iip1, jjp1, llm)      real, intent(in):: teta(iip1, jjp1, llm)
28      real, intent(in):: phi(iip1, jjp1, llm)      real, intent(in):: phi(iip1, jjp1, llm)
29      real, intent(in):: ucov(iip1, jjp1, llm)      real, intent(in):: ucov(:, :, :) ! (iip1, jjp1, llm)
30      real, intent(in):: vcov(iip1, jjm, llm)      real, intent(in):: vcov(iip1, jjm, llm)
31      real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)      real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)
32    
# Line 36  contains Line 34  contains
34    
35      integer:: icum  = 0      integer:: icum  = 0
36      integer:: itau = 0      integer:: itau = 0
37      real zqy, zfactv(jjm, llm)      real qy, factv(jjm, llm)
   
     real ww  
38    
39      ! Variables dynamiques intermédiaires      ! Variables dynamiques intermédiaires
40      REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)      REAL vcont(iip1, jjm, llm), ucont(iip1, jjp1, llm)
41      REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)      REAL ang(iip1, jjp1, llm), unat(iip1, jjp1, llm)
42      REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)      REAL massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
43      REAL w(iip1, jjp1, llm), ecin(iip1, jjp1, llm), convm(iip1, jjp1, llm)      REAL ecin(iip1, jjp1, llm)
44    
45      ! Champ contenant les scalaires advectés      ! Champ contenant les scalaires advectés
46      real Q(iip1, jjp1, llm, nQ)      real Q(iip1, jjp1, llm, nQ)
# Line 57  contains Line 53  contains
53      real, save:: Q_cum(iip1, jjp1, llm, nQ)      real, save:: Q_cum(iip1, jjp1, llm, nQ)
54      real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)      real, save:: flux_uQ_cum(iip1, jjp1, llm, nQ)
55      real, save:: flux_vQ_cum(iip1, jjm, llm, nQ)      real, save:: flux_vQ_cum(iip1, jjm, llm, nQ)
     real dQ(iip1, jjp1, llm, nQ)  
56    
57      ! champs de tansport en moyenne zonale      ! champs de tansport en moyenne zonale
58      integer itr      integer itr
59      integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5      integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5
60    
61      real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)      real vq(jjm, llm, ntr, nQ), vqtmp(jjm, llm)
62      real zavQ(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)      real avq(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)
63      real zmasse(jjm, llm)      real zmasse(jjm, llm)
64      real zv(jjm, llm), psi(jjm, llm + 1)      real v(jjm, llm), psi(jjm, llm + 1)
65      integer i, j, l, iQ      integer i, j, l, iQ
66    
67      !-----------------------------------------------------------------      !-----------------------------------------------------------------
# Line 79  contains Line 74  contains
74      CALL enercin(vcov, ucov, vcont, ucont, ecin)      CALL enercin(vcov, ucov, vcont, ucont, ecin)
75    
76      ! moment cinétique      ! moment cinétique
77      do l = 1, llm      forall (l = 1: llm)
78         ang(:, :, l) = ucov(:, :, l) + constang_2d         ang(:, :, l) = ucov(:, :, l) + constang_2d
79         unat(:, :, l) = ucont(:, :, l)*cu_2d         unat(:, :, l) = ucont(:, :, l) * cu_2d
80      enddo      end forall
81    
82      Q(:, :, :, 1) = teta * pk / cpp      Q(:, :, :, 1) = teta * pk / cpp
83      Q(:, :, :, 2) = phi      Q(:, :, :, 2) = phi
# Line 112  contains Line 107  contains
107      masse_cum = masse_cum + masse      masse_cum = masse_cum + masse
108      flux_u_cum = flux_u_cum + flux_u      flux_u_cum = flux_u_cum + flux_u
109      flux_v_cum = flux_v_cum + flux_v      flux_v_cum = flux_v_cum + flux_v
110      do iQ = 1, nQ      forall (iQ = 1: nQ) Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) &
111         Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) + Q(:, :, :, iQ)*masse           + Q(:, :, :, iQ) * masse
     enddo  
   
     ! FLUX ET TENDANCES  
112    
113      ! Flux longitudinal      ! Flux longitudinal
114      forall (iQ = 1: nQ, i = 1: iim) flux_uQ_cum(i, :, :, iQ) &      forall (iQ = 1: nQ, i = 1: iim) flux_uQ_cum(i, :, :, iQ) &
# Line 129  contains Line 121  contains
121           = flux_vQ_cum(:, j, :, iQ) &           = flux_vQ_cum(:, j, :, iQ) &
122           + flux_v(:, j, :) * 0.5 * (Q(:, j, :, iQ) + Q(:, j + 1, :, iQ))           + flux_v(:, j, :) * 0.5 * (Q(:, j, :, iQ) + Q(:, j + 1, :, iQ))
123    
     ! 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      writing_step: if (icum == ncum) then      writing_step: if (icum == ncum) then
125         ! Normalisation         ! Normalisation
126         do iQ = 1, nQ         forall (iQ = 1: nQ) Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ) / masse_cum
           Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ)/masse_cum  
        enddo  
127         ps_cum = ps_cum / ncum         ps_cum = ps_cum / ncum
128         masse_cum = masse_cum / ncum         masse_cum = masse_cum / ncum
129         flux_u_cum = flux_u_cum / ncum         flux_u_cum = flux_u_cum / ncum
130         flux_v_cum = flux_v_cum / ncum         flux_v_cum = flux_v_cum / ncum
131         flux_uQ_cum = flux_uQ_cum / ncum         flux_uQ_cum = flux_uQ_cum / ncum
132         flux_vQ_cum = flux_vQ_cum / ncum         flux_vQ_cum = flux_vQ_cum / ncum
        dQ = dQ / ncum  
   
        ! 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  
133    
134         ! Transport méridien         ! Transport méridien
135    
136         ! cumul zonal des masses des mailles         ! Cumul zonal des masses des mailles
137    
138         zv = 0.         v = 0.
139         zmasse = 0.         zmasse = 0.
140         call massbar(masse_cum, massebx, masseby)         call massbar(masse_cum, massebx, masseby)
141         do l = 1, llm         do l = 1, llm
142            do j = 1, jjm            do j = 1, jjm
143               do i = 1, iim               do i = 1, iim
144                  zmasse(j, l) = zmasse(j, l) + masseby(i, j, l)                  zmasse(j, l) = zmasse(j, l) + masseby(i, j, l)
145                  zv(j, l) = zv(j, l) + flux_v_cum(i, j, l)                  v(j, l) = v(j, l) + flux_v_cum(i, j, l)
146               enddo               enddo
147               zfactv(j, l) = cv_2d(1, j)/zmasse(j, l)               factv(j, l) = cv_2d(1, j) / zmasse(j, l)
148            enddo            enddo
149         enddo         enddo
150    
151         ! Transport dans le plan latitude-altitude         ! Transport dans le plan latitude-altitude
152    
153         zvQ = 0.         vq = 0.
154         psiQ = 0.         psiQ = 0.
155         do iQ = 1, nQ         do iQ = 1, nQ
156            zvQtmp = 0.            vqtmp = 0.
157            do l = 1, llm            do l = 1, llm
158               do j = 1, jjm               do j = 1, jjm
159                  ! Calcul des moyennes zonales du transort total et de zvQtmp                  ! Calcul des moyennes zonales du transport total et de vqtmp
160                  do i = 1, iim                  do i = 1, iim
161                     zvQ(j, l, itot, iQ) = zvQ(j, l, itot, iQ) &                     vq(j, l, itot, iQ) = vq(j, l, itot, iQ) &
162                          + flux_vQ_cum(i, j, l, iQ)                          + flux_vQ_cum(i, j, l, iQ)
163                     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) &
164                          + 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))
165                     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 &
166                          / (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)))
167                     zvQ(j, l, iave, iQ) = zvQ(j, l, iave, iQ) + zqy                     vq(j, l, iave, iQ) = vq(j, l, iave, iQ) + qy
168                  enddo                  enddo
169                  ! Decomposition                  ! Decomposition
170                  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)
171                  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)
172                  zvQtmp(j, l) = zvQtmp(j, l)*zfactv(j, l)                  vqtmp(j, l) = vqtmp(j, l) * factv(j, l)
173                  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)
174                  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)
175                  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)
176               enddo               enddo
177            enddo            enddo
178            ! fonction de courant meridienne pour la quantite Q            ! Fonction de courant méridienne pour la quantité Q
179            do l = llm, 1, -1            do l = llm, 1, -1
180               do j = 1, jjm               do j = 1, jjm
181                  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)
182               enddo               enddo
183            enddo            enddo
184         enddo         enddo
185    
186         ! fonction de courant pour la circulation meridienne moyenne         ! Fonction de courant pour la circulation méridienne moyenne
187         psi = 0.         psi = 0.
188         do l = llm, 1, -1         do l = llm, 1, -1
189            do j = 1, jjm            do j = 1, jjm
190               psi(j, l) = psi(j, l + 1) + zv(j, l)               psi(j, l) = psi(j, l + 1) + v(j, l)
191               zv(j, l) = zv(j, l)*zfactv(j, l)               v(j, l) = v(j, l) * factv(j, l)
192            enddo            enddo
193         enddo         enddo
194    
195         ! sorties proprement dites         ! Sorties proprement dites
196         do iQ = 1, nQ         do iQ = 1, nQ
197            do itr = 1, ntr            do itr = 1, ntr
198               call histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ))               call histwrite(fileid, znom(itr, iQ), itau, vq(:, :, itr, iQ))
199            enddo            enddo
200            call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, :llm, iQ))            call histwrite(fileid, 'psi' // nom(iQ), itau, psiQ(:, :llm, iQ))
201         enddo         enddo
202    
203         call histwrite(fileid, 'masse', itau, zmasse)         call histwrite(fileid, 'masse', itau, zmasse)
204         call histwrite(fileid, 'v', itau, zv)         call histwrite(fileid, 'v', itau, v)
205         psi = psi*1.e-9         psi = psi * 1e-9
206         call histwrite(fileid, 'psi', itau, psi(:, :llm))         call histwrite(fileid, 'psi', itau, psi(:, :llm))
207    
208         ! Intégrale verticale         ! Intégrale verticale
209    
210         forall (iQ = 1: nQ, itr = 2: ntr) zavQ(:, itr, iQ) &         forall (iQ = 1: nQ, itr = 2: ntr) avq(:, itr, iQ) &
211              = sum(zvQ(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)              = sum(vq(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)
212    
213         do iQ = 1, nQ         do iQ = 1, nQ
214            do itr = 2, ntr            do itr = 2, ntr
215               call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))               call histwrite(fileid, 'a' // znom(itr, iQ), itau, avq(:, itr, iQ))
216            enddo            enddo
217         enddo         enddo
218    
        ! On doit pouvoir tracer systematiquement la fonction de courant.  
219         icum = 0         icum = 0
220      endif writing_step      endif writing_step
221    
Legend:
Removed from v.57  
changed lines
  Added in v.62
  ViewVC Help
Powered by ViewVC 1.1.21