source: trunk/LMDZ.TITAN/Tools/epflux.F90 @ 819

Last change on this file since 819 was 816, checked in by slebonnois, 12 years ago

SL: tools for postprocessing (Veznus and Titan); see DOC/documentation/vt-tools.pdf

File size: 14.2 KB
Line 
1subroutine  epflux(iip1,jjp1,llm,indefini,rlatu,rbar &
2                  ,teta,u3d,v3d,w3d,press &
3                  ,epy2d,epz2d,divep2d,vtem2d,wtem2d,acc_meridien2d &
4!                 ,vpupbar,wpupbar,vpvpbar,wpvpbar,vptetapbar,wptetapbar &
5                  )
6
7      IMPLICIT NONE
8!=======================================================================
9!
10!   Audrey Crespin  Sept 2007
11!
12!   source: Francois Forget    Avril 1996
13!
14!
15! MODIF SLebonnois nov 2009
16!
17! Cette subroutine Calcule Les flux d'Eliassen Palm a partir
18! des composantes INTERPOLE EN NIVEAU DE PRESSION
19! Toutes les variables sont sur la grille de pression press(1:llm)
20! On suit la methode de Peixoto et Oort
21!
22! On est dans le plan meridien
23!
24! Calcule la divergence & composantes du flux d'EP
25! Calcule la deviation a la circulation meridienne moyenne
26! Calcule les TEM
27! Calcule les termes de Reynolds
28!
29!=======================================================================
30!-----------------------------------------------------------------------
31!   declarations:
32!   -------------
33
34include "planet.h"
35
36!  --------
37!  ARGUMENTS
38!  ---------
39
40! Inputs:
41     
42      integer :: iip1,jjp1,llm
43      real :: indefini,rlatu(jjp1),rbar(jjp1,llm)
44      REAL :: teta(iip1,jjp1,llm)
45      REAL :: u3d(iip1,jjp1,llm),v3d(iip1,jjp1,llm)
46      REAL :: w3d(iip1,jjp1,llm)
47      REAL :: press(llm)
48
49! Outputs:
50
51      REAL :: epy2d(jjp1,llm),epz2d(jjp1,llm),divep2d(jjp1,llm)
52      REAL :: vtem2d(jjp1,llm),wtem2d(jjp1,llm),acc_meridien2d(jjp1,llm)
53      REAL :: vpupbar(jjp1,llm),wpupbar(jjp1,llm)
54      REAL :: wpvpbar(jjp1,llm),vpvpbar(jjp1,llm)
55      REAL :: vptetapbar(jjp1,llm),wptetapbar(jjp1,llm)
56
57! Variables non zonales (vpup signifie : v^prim * u^prim)
58! -------------------
59      REAL :: vpup(iip1,jjp1,llm), wpup(iip1,jjp1,llm)
60      REAL :: vpvp(iip1,jjp1,llm), wpvp(iip1,jjp1,llm)
61      REAL :: vptetap(iip1,jjp1,llm), wptetap(iip1,jjp1,llm)
62
63!  Moyennes zonales
64!  -------------
65      REAL :: ubar(jjp1,llm),vbar(jjp1,llm),tetabar(jjp1,llm)
66      REAL :: wbar(jjp1,llm)
67
68
69      REAL :: dtetabardp(jjp1,llm),dubardp(jjp1,llm)
70      REAL :: vz(jjp1,llm),wlat(jjp1,llm),dvzdp(jjp1,llm)
71      REAL :: dwlatdlat(jjp1,llm)
72
73      REAL :: depzdp(jjp1,llm)
74
75!  Autres
76!  ------
77      real :: pi
78      REAL :: tetas(llm)
79      logical :: peixoto
80
81      REAL :: f(jjp1), ducosdlat(jjp1,llm)
82      REAL :: depycosdlat
83      INTEGER :: i,j,l, n, iim, jjm
84
85!-----------------------------------------------------------------------
86
87      iim = iip1-1
88      jjm = jjp1-1
89      pi = 2.*asin(1.)
90
91!     Calcul de moyennes zonales
92!     ---------------------------
93      call moyzon(iim,jjp1,llm,indefini,u3d,ubar)
94      call moyzon(iim,jjp1,llm,indefini,v3d,vbar)
95      call moyzon(iim,jjp1,llm,indefini,teta,tetabar)
96      call moyzon(iim,jjp1,llm,indefini,w3d,wbar)
97
98
99! --*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--
100      peixoto = .false.
101
102      if (peixoto) then
103
104!     MODIF speciale Peixoto: on utilise la moyenne
105!     globale de teta: tetas(l) a la place de tetabar
106      do l=1,llm
107         tetas(l) = 0
108         n= 0
109         do j=2,jjm
110            if (tetabar(j,l).lt.indefini) then
111               tetas(l) = tetas(l) + tetabar(j,l)
112               n = n+1
113            end if
114         end do
115         if (n.eq.0) stop 'bug dans elias '
116        tetas(l) = tetas(l) / float(n)
117
118        do j=1,jjp1
119             tetabar(j,l) = tetas(l)
120        end do
121      end do
122!       write (*,*) 'tetas(l) ' , tetas
123!       write(*,*)
124
125      endif
126! --*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--
127
128! coriolis
129! --------
130      do j=1,jjp1
131            f(j) = 2*omega*sin(rlatu(j))
132      enddo
133
134
135!     Calcul des termes non zonaux
136!     ----------------------------
137      do l=1,llm
138         do j=1,jjp1
139           do i=1,iip1
140             if ((  v3d(i,j,l).eq.indefini).or. &
141                 (   vbar(j,l).eq.indefini).or. &
142                 ( teta(i,j,l).eq.indefini).or. &
143                 (tetabar(j,l).eq.indefini)) then
144               vptetap(i,j,l)= indefini
145             else
146               vptetap(i,j,l)=(v3d(i,j,l)-vbar(j,l)) &
147                  *(teta(i,j,l)-tetabar(j,l))
148             end if
149             if ((  w3d(i,j,l).eq.indefini).or. &
150                 (   wbar(j,l).eq.indefini).or. &
151                 ( teta(i,j,l).eq.indefini).or. &
152                 (tetabar(j,l).eq.indefini)) then
153               wptetap(i,j,l)= indefini
154             else
155               wptetap(i,j,l)=(w3d(i,j,l)-wbar(j,l)) &
156                  *(teta(i,j,l)-tetabar(j,l))
157             end if
158             if ((v3d(i,j,l).eq.indefini).or. &
159                 ( vbar(j,l).eq.indefini).or. &
160                 (u3d(i,j,l).eq.indefini).or. &
161                 ( ubar(j,l).eq.indefini)) then
162               vpup(i,j,l)= indefini
163             else
164               vpup(i,j,l)=(v3d(i,j,l)-vbar(j,l))*(u3d(i,j,l)-ubar(j,l))
165             end if
166             if ((w3d(i,j,l).eq.indefini).or. &
167                 ( wbar(j,l).eq.indefini).or. &
168                 (u3d(i,j,l).eq.indefini).or. &
169                 ( ubar(j,l).eq.indefini)) then
170               wpup(i,j,l)= indefini
171             else
172               wpup(i,j,l)=(w3d(i,j,l)-wbar(j,l))*(u3d(i,j,l)-ubar(j,l))
173             end if
174             if ((v3d(i,j,l).eq.indefini).or. &
175                 ( vbar(j,l).eq.indefini)) then
176               vpvp(i,j,l)= indefini
177             else
178               vpvp(i,j,l)=(v3d(i,j,l)-vbar(j,l))*(v3d(i,j,l)-vbar(j,l))
179             end if
180             if ((w3d(i,j,l).eq.indefini).or. &
181                 ( wbar(j,l).eq.indefini).or. &
182                 (v3d(i,j,l).eq.indefini).or. &
183                 ( vbar(j,l).eq.indefini)) then
184               wpvp(i,j,l)= indefini
185             else
186               wpvp(i,j,l)=(w3d(i,j,l)-wbar(j,l))*(v3d(i,j,l)-vbar(j,l))
187             end if
188           end do
189         end do
190      end do
191
192!     Moyennes zonales des termes non zonaux
193!     --------------------------------------
194
195!     Termes de Reynolds
196
197!     flux zonaux de quantite de mouvement (vpupbar et wpupbar)
198!     flux meridiens de quantite de mouvement (vpvpbar et wpvpbar)
199!     flux meridien et vertical de chaleur (vptetapbar et wptetapbar)
200
201      call moyzon(iim,jjp1,llm,indefini,vptetap,vptetapbar)
202      call moyzon(iim,jjp1,llm,indefini,wptetap,wptetapbar)
203      call moyzon(iim,jjp1,llm,indefini,vpup,vpupbar)
204      call moyzon(iim,jjp1,llm,indefini,wpup,wpupbar)
205      call moyzon(iim,jjp1,llm,indefini,vpvp,vpvpbar)
206      call moyzon(iim,jjp1,llm,indefini,wpvp,wpvpbar)
207
208!     write(*,*) 'vptetapbar(2,23) =' , vptetapbar(2,23)
209!     write(*,*) 'wptetapbar(2,23) =' , wptetapbar(2,23)
210!     write(*,*) 'vpupbar(2,23) =' , vpupbar(2,23)
211
212 
213!     Derivees d/dp des moyennes zonales
214!     --------------------------------------
215
216      call dx_dp(jjp1,llm,indefini,press,ubar,dubardp)
217      call dx_dp(jjp1,llm,indefini,press,tetabar,dtetabardp)
218!     write (*,*) 'dtetabardp(l) (K/Pa)',(dtetabardp(6,l),l=1,llm)
219!       write(*,*)
220
221!     Controle (on triche !) de dtetabardp
222      do l=1,llm
223         do j=1,jjp1
224            if ((dtetabardp(j,l).gt.-1.e-3).and. &
225               (dtetabardp(j,l).ne.indefini)) then
226               dtetabardp(j,l) = -1.e-3
227!              write(*,*) 'profil presque instable en j = ',j,' l= ',l
228            end if
229         end do
230      end do
231!     write(*,*)'dtetabardp(l) (K/Pa) corr ',(dtetabardp(6,l),l=1,llm)
232!       write(*,*)
233
234
235!     calculs intermediaires
236!     ----------------------
237      do l=1,llm
238
239       if ((ubar(2,l).lt.indefini).and.(ubar(1,l).lt.indefini))then
240      ducosdlat(1,l) = ( ubar(2,l)*cos(rlatu(2))-ubar(1,l)*cos(rlatu(1)) ) &
241                      /( rlatu(2) - rlatu(1))
242       else
243      ducosdlat(1,l) = indefini
244       end if
245        do j=2,jjm
246       if ((ubar(j+1,l).lt.indefini).and.(ubar(j-1,l).lt.indefini))then
247      ducosdlat(j,l) = ( ubar(j+1,l)*cos(rlatu(j+1))-ubar(j-1,l)*cos(rlatu(j-1))) &
248                      /( rlatu(j+1) - rlatu(j-1))
249       else
250      ducosdlat(j,l) = indefini
251       end if
252        enddo
253       if ((ubar(jjp1,l).lt.indefini).and.(ubar(jjm,l).lt.indefini))then
254      ducosdlat(jjp1,l) = ( ubar(jjp1,l)*cos(rlatu(jjp1))  &
255                               -ubar(jjm,l)*cos(rlatu(jjm))  ) &
256                         /( rlatu(jjp1) - rlatu(jjm))
257       else
258      ducosdlat(jjp1,l) = indefini
259       end if
260           
261        do j=1,jjp1
262       if ((vptetapbar(j,l).lt.indefini).and. &
263           (dtetabardp(j,l).lt.indefini)) then
264           vz(j,l) = vptetapbar(j,l)/dtetabardp(j,l)
265         wlat(j,l) = cos(rlatu(j))*vz(j,l)
266       else
267           vz(j,l) = indefini
268         wlat(j,l) = indefini
269       end if
270        enddo
271
272      enddo
273             
274!      Calcul des vecteurs flux Eliassen Palm et divergence
275!      -----------------------------------------------------
276! ref: Read 1986
277
278      do l=1,llm
279        epy2d(1,l) = indefini
280        epz2d(1,l) = indefini
281        do j=2,jjm
282 
283!     Composante y
284!     ------------
285
286            if (   (rbar(j,l).lt.indefini).and. &
287                (vpupbar(j,l).lt.indefini).and. &
288                (dubardp(j,l).lt.indefini).and. &
289                     (vz(j,l).lt.indefini)) then
290           
291     epy2d(j,l) = rbar(j,l) * cos(rlatu(j))* ( vpupbar(j,l) &
292                                   - dubardp(j,l) * vz(j,l) & ! terme non geo
293                                             )       
294
295!     if ((j.eq.1).and.(l.eq.23))write(*,*) 'epy2d(1,23) =' , epy2d(1,23)
296!     if ((j.eq.2).and.(l.eq.23))write(*,*) 'epy2d(2,23) =' , epy2d(2,23)
297!     if ((j.eq.3).and.(l.eq.23))write(*,*) 'epy2d(3,23) =' , epy2d(3,23)
298           else
299     epy2d(j,l) = indefini
300            end if
301
302!     Composante z
303!     ------------
304
305            if (   (rbar(j,l).lt.indefini).and. &
306                (wpupbar(j,l).lt.indefini).and. &
307              (ducosdlat(j,l).lt.indefini).and. &
308                     (vz(j,l).lt.indefini)) then
309
310     epz2d(j,l) = rbar(j,l) * cos(rlatu(j))*                    &
311          (  vz(j,l) * ducosdlat(j,l)/(rbar(j,l)*cos(rlatu(j))) & ! terme non geo
312           - vz(j,l) * f(j)                                     &
313           + wpupbar(j,l)                                       &
314          )
315           else
316     epz2d(j,l) = indefini
317           end if
318
319!       if ((j.eq.2).and.(l.eq.23))write(*,*) 'epz2d(2,23) =' , epz2d(2,23)
320        end do
321        epy2d(jjp1,l) = indefini
322        epz2d(jjp1,l) = indefini
323      end do
324
325!    Moyennes euleriennes transformees (vtem2d et wtem2d)
326!    -----------------------------------------------
327     
328      call dx_dp(jjp1,llm,indefini,press,vz,dvzdp)
329
330      do l=1,llm
331
332          if ( (wlat(2,l).lt.indefini) .and. &
333               (wlat(1,l).lt.indefini) ) then
334               dwlatdlat(1,l)=(wlat(2,l)-wlat(1,l)) &
335                     /(rlatu(2)-rlatu(1))
336          else
337               dwlatdlat(1,l)=indefini
338          end if
339        do j=2,jjm
340          if ( (wlat(j+1,l).lt.indefini) .and. &
341               (wlat(j-1,l).lt.indefini) ) then
342               dwlatdlat(j,l)=(wlat(j+1,l)-wlat(j-1,l)) &
343                     /(rlatu(j+1)-rlatu(j-1))
344          else
345               dwlatdlat(j,l)=indefini
346          end if
347        enddo
348          if ( (wlat(jjp1,l).lt.indefini) .and. &
349               (wlat(jjm, l).lt.indefini) ) then
350               dwlatdlat(jjp1,l)=(wlat(jjp1,l)-wlat(jjm,l)) &
351                     /(rlatu(jjp1)-rlatu(jjm))
352          else
353               dwlatdlat(jjp1,l)=indefini
354          end if
355
356        do j=1,jjp1
357          if ( (dvzdp(j,l).lt.indefini) &
358           .and.(vbar(j,l).lt.indefini)) then
359            vtem2d(j,l) = vbar(j,l)-dvzdp(j,l)
360          else
361            vtem2d(j,l) = indefini
362          endif
363        enddo
364
365        do j=1,jjp1
366          if (   (rbar(j,l).lt.indefini).and. &
367            (dwlatdlat(j,l).lt.indefini).and. &
368                 (wbar(j,l).lt.indefini)) then
369           wtem2d(j,l) = wbar(j,l)+dwlatdlat(j,l)/(rbar(j,l)*cos(rlatu(j)))
370          else
371           wtem2d(j,l) = indefini
372          endif
373        enddo
374
375      enddo
376
377!      write(*,*) 'vtem2d(2,23) =' , vtem2d(2,23)
378!      write(*,*) 'wtem2d(2,23) =' , wtem2d(2,23)
379
380!      print*,"OK"
381
382!    Deviation par rapport a la circulation meridienne moyenne
383!    --------------------------------------------------------
384!                    U*.Nabla ubar
385      do l=1,llm
386        do j=1,jjp1
387          if (   (rbar(j,l).lt.indefini) &
388          .and.(vtem2d(j,l).lt.indefini) &
389          .and.(wtem2d(j,l).lt.indefini) &
390       .and.(ducosdlat(j,l).lt.indefini) &
391         .and.(dubardp(j,l).lt.indefini)) then
392       acc_meridien2d(j,l) =  &
393           vtem2d(j,l)*(ducosdlat(j,l)/(rbar(j,l)*cos(rlatu(j))) &
394                       -f(j)) &
395         + wtem2d(j,l)*dubardp(j,l)     
396          else
397       acc_meridien2d(j,l) = indefini
398          endif
399        enddo
400      enddo       
401
402!      write(*,*) 'acc_meridien2d(2,23) =' , acc_meridien2d(2,23)
403
404!     Divergence du flux
405!     ------------------
406      call dx_dp(jjp1,llm,indefini,press,epz2d,depzdp)
407!     write(*,*) 'depzdp(2,23)',depzdp(j,l)
408
409      do l=1,llm
410        divep2d(1,l) =0
411        do j=2,jjm
412          if (    (rbar(j,l).lt.indefini) .and. &
413               (epy2d(j+1,l).lt.indefini) .and. &
414               (epy2d(j-1,l).lt.indefini) .and. &
415                (depzdp(j,l).lt.indefini)    ) then
416
417           depycosdlat=(epy2d(j+1,l)*cos(rlatu(j+1))   &
418                      - epy2d(j-1,l)*cos(rlatu(j-1)) ) &
419                       /(rlatu(j+1) - rlatu(j-1))
420
421!      DIVERGENCE DU FLUX :
422
423           divep2d(j,l) = depycosdlat/(rbar(j,l)*cos(rlatu(j)))+depzdp(j,l)
424
425!          if ((j.eq.2).and.(l.eq.23)) then
426!              write(*,*) 'depycosdlat(2,23)', depycosdlat
427!              write(*,*) 'depzdp(2,23)',depzdp(j,l)
428!              write(*,*) 'divergence(2,23)',divep2d(j,l)
429!          end if
430
431! Wave driving (Du/Dt et non Dm/dt) :
432           divep2d(j,l) = divep2d(j,l)/(rbar(j,l)*cos(rlatu(j)))
433          else
434            divep2d(j,l) = indefini
435          end if
436        end do
437        divep2d(jjp1,l) =0
438      end do
439 
440!     write(*,*) ' fin :divep2d(2,23) =' , divep2d(2,23)
441
442! Preparation pour sortie graphique (flux ``chapeau'' Peixoto et Oort p 391)
443! A VOIR...
444!      do l=1,llm
445!        do j=1,jjp1
446!           if (epy2d(j,l).ne.indefini) &
447!           epy2d(j,l) = 2*pi*a0    *(1./g0)*cos(rlatu(j))*epy2d(j,l)
448!           if (epz2d(j,l).ne.indefini) &
449!           epz2d(j,l) = 2*pi*a0**2 *(1./g0)*cos(rlatu(j))*epz2d(j,l)
450!        end do
451!      end do
452
453      return
454      end
Note: See TracBrowser for help on using the repository browser.