/[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

trunk/libf/dyn3d/bilan_dyn.f90 revision 57 by guez, Mon Jan 30 12:54:02 2012 UTC trunk/dyn3d/bilan_dyn.f revision 91 by guez, Wed Mar 26 17:18:58 2014 UTC
# Line 18  contains Line 18  contains
18      USE dimens_m, ONLY: iim, jjm, llm      USE dimens_m, ONLY: iim, jjm, llm
19      USE histwrite_m, ONLY: histwrite      USE histwrite_m, ONLY: histwrite
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 massbar_m, only: massbar
22      USE paramet_m, ONLY: iip1, jjp1      USE paramet_m, ONLY: iip1, jjp1
23    
     ! Arguments:  
   
24      real, intent(in):: ps(iip1, jjp1)      real, intent(in):: ps(iip1, jjp1)
25      real, intent(in):: masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)      real, intent(in):: masse(iip1, jjp1, llm), pk(iip1, jjp1, llm)
26      real, intent(in):: flux_u(iip1, jjp1, llm)      real, intent(in):: flux_u(iip1, jjp1, llm)
27      real, intent(in):: flux_v(iip1, jjm, llm)      real, intent(in):: flux_v(iip1, jjm, llm)
28      real, intent(in):: teta(iip1, jjp1, llm)      real, intent(in):: teta(iip1, jjp1, llm)
29      real, intent(in):: phi(iip1, jjp1, llm)      real, intent(in):: phi(iip1, jjp1, llm)
30      real, intent(in):: ucov(iip1, jjp1, llm)      real, intent(in):: ucov(:, :, :) ! (iip1, jjp1, llm)
31      real, intent(in):: vcov(iip1, jjm, llm)      real, intent(in):: vcov(iip1, jjm, llm)
32      real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)      real, intent(in):: trac(:, :, :) ! (iim + 1, jjm + 1, llm)
33    
# Line 36  contains Line 35  contains
35    
36      integer:: icum  = 0      integer:: icum  = 0
37      integer:: itau = 0      integer:: itau = 0
38      real zqy, zfactv(jjm, llm)      real qy, factv(jjm, llm)
   
     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 w(iip1, jjp1, llm), ecin(iip1, jjp1, llm), convm(iip1, jjp1, llm)      REAL ecin(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 57  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 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:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5      integer, parameter:: iave = 1, itot = 2, immc = 3, itrs = 4, istn = 5
61    
62      real zvQ(jjm, llm, ntr, nQ), zvQtmp(jjm, llm)      real vq(jjm, llm, ntr, nQ), vqtmp(jjm, llm)
63      real zavQ(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)      real avq(jjm, 2: ntr, nQ), psiQ(jjm, llm + 1, nQ)
64      real zmasse(jjm, llm)      real zmasse(jjm, llm)
65      real zv(jjm, llm), psi(jjm, llm + 1)      real v(jjm, llm), psi(jjm, llm + 1)
66      integer i, j, l, iQ      integer i, j, l, iQ
67    
68      !-----------------------------------------------------------------      !-----------------------------------------------------------------
# Line 79  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(:, :, :, 1) = teta * pk / cpp      Q(:, :, :, 1) = teta * pk / cpp
84      Q(:, :, :, 2) = phi      Q(:, :, :, 2) = phi
# Line 112  contains Line 108  contains
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      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 122  contains
122           = flux_vQ_cum(:, j, :, iQ) &           = flux_vQ_cum(:, j, :, iQ) &
123           + flux_v(:, j, :) * 0.5 * (Q(:, j, :, iQ) + Q(:, j + 1, :, iQ))           + flux_v(:, j, :) * 0.5 * (Q(:, j, :, iQ) + Q(:, j + 1, :, iQ))
124    
     ! 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  
   
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
           Q_cum(:, :, :, iQ) = Q_cum(:, :, :, iQ)/masse_cum  
        enddo  
128         ps_cum = ps_cum / ncum         ps_cum = ps_cum / ncum
129         masse_cum = masse_cum / ncum         masse_cum = masse_cum / ncum
130         flux_u_cum = flux_u_cum / ncum         flux_u_cum = flux_u_cum / ncum
131         flux_v_cum = flux_v_cum / ncum         flux_v_cum = flux_v_cum / ncum
132         flux_uQ_cum = flux_uQ_cum / ncum         flux_uQ_cum = flux_uQ_cum / ncum
133         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  
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         do iQ = 1, nQ         do iQ = 1, nQ
198            do itr = 1, ntr            do itr = 1, ntr
199               call histwrite(fileid, znom(itr, iQ), itau, zvQ(:, :, itr, iQ))               call histwrite(fileid, znom(itr, iQ), itau, vq(:, :, itr, iQ))
200            enddo            enddo
201            call histwrite(fileid, 'psi'//nom(iQ), itau, psiQ(:, :llm, iQ))            call histwrite(fileid, 'psi' // nom(iQ), itau, psiQ(:, :llm, iQ))
202         enddo         enddo
203    
204         call histwrite(fileid, 'masse', itau, zmasse)         call histwrite(fileid, 'masse', itau, zmasse)
205         call histwrite(fileid, 'v', itau, zv)         call histwrite(fileid, 'v', itau, v)
206         psi = psi*1.e-9         psi = psi * 1e-9
207         call histwrite(fileid, 'psi', itau, psi(:, :llm))         call histwrite(fileid, 'psi', itau, psi(:, :llm))
208    
209         ! Intégrale verticale         ! Intégrale verticale
210    
211         forall (iQ = 1: nQ, itr = 2: ntr) zavQ(:, itr, iQ) &         forall (iQ = 1: nQ, itr = 2: ntr) avq(:, itr, iQ) &
212              = sum(zvQ(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)              = sum(vq(:, :, itr, iQ) * zmasse, dim=2) / cv_2d(1, :)
213    
214         do iQ = 1, nQ         do iQ = 1, nQ
215            do itr = 2, ntr            do itr = 2, ntr
216               call histwrite(fileid, 'a'//znom(itr, iQ), itau, zavQ(:, itr, iQ))               call histwrite(fileid, 'a' // znom(itr, iQ), itau, avq(:, itr, iQ))
217            enddo            enddo
218         enddo         enddo
219    
        ! On doit pouvoir tracer systematiquement la fonction de courant.  
220         icum = 0         icum = 0
221      endif writing_step      endif writing_step
222    
Legend:
Removed from v.57  
changed lines
  Added in v.91
  ViewVC Help
Powered by ViewVC 1.1.21