source: trunk/LMDZ.TITAN.old/Tools/epflux.F90 @ 3094

Last change on this file since 3094 was 1454, checked in by slebonnois, 10 years ago

SL: corrections in Titan and Venus tools

File size: 14.4 KB
Line 
1subroutine  epflux(iip1,jjp1,llm,indefini,latdeg,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,latdeg(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 :: rlatu(jjp1) ! lat in radians (beware of rounding effects at poles)
78      real :: pi
79      REAL :: tetas(llm)
80      logical :: peixoto
81
82      REAL :: f(jjp1), ducosdlat(jjp1,llm)
83      REAL :: depycosdlat
84      INTEGER :: i,j,l, n, iim, jjm
85
86!-----------------------------------------------------------------------
87
88      iim = iip1-1
89      jjm = jjp1-1
90      pi = 2.*asin(1.)
91
92! To avoid rounding effects at poles
93      rlatu(:)=latdeg(:)*pi/180.00001
94
95!     Calcul de moyennes zonales
96!     ---------------------------
97      call moyzon(iim,jjp1,llm,indefini,u3d,ubar)
98      call moyzon(iim,jjp1,llm,indefini,v3d,vbar)
99      call moyzon(iim,jjp1,llm,indefini,teta,tetabar)
100      call moyzon(iim,jjp1,llm,indefini,w3d,wbar)
101
102
103! --*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--
104      peixoto = .false.
105
106      if (peixoto) then
107
108!     MODIF speciale Peixoto: on utilise la moyenne
109!     globale de teta: tetas(l) a la place de tetabar
110      do l=1,llm
111         tetas(l) = 0
112         n= 0
113         do j=2,jjm
114            if (tetabar(j,l).ne.indefini) then
115               tetas(l) = tetas(l) + tetabar(j,l)
116               n = n+1
117            end if
118         end do
119         if (n.eq.0) stop 'bug dans elias '
120        tetas(l) = tetas(l) / float(n)
121
122        do j=1,jjp1
123             tetabar(j,l) = tetas(l)
124        end do
125      end do
126!       write (*,*) 'tetas(l) ' , tetas
127!       write(*,*)
128
129      endif
130! --*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--
131
132! coriolis
133! --------
134      do j=1,jjp1
135            f(j) = 2*omega*sin(rlatu(j))
136      enddo
137
138
139!     Calcul des termes non zonaux
140!     ----------------------------
141      do l=1,llm
142         do j=1,jjp1
143           do i=1,iip1
144             if ((  v3d(i,j,l).eq.indefini).or. &
145                 (   vbar(j,l).eq.indefini).or. &
146                 ( teta(i,j,l).eq.indefini).or. &
147                 (tetabar(j,l).eq.indefini)) then
148               vptetap(i,j,l)= indefini
149             else
150               vptetap(i,j,l)=(v3d(i,j,l)-vbar(j,l)) &
151                  *(teta(i,j,l)-tetabar(j,l))
152             end if
153             if ((  w3d(i,j,l).eq.indefini).or. &
154                 (   wbar(j,l).eq.indefini).or. &
155                 ( teta(i,j,l).eq.indefini).or. &
156                 (tetabar(j,l).eq.indefini)) then
157               wptetap(i,j,l)= indefini
158             else
159               wptetap(i,j,l)=(w3d(i,j,l)-wbar(j,l)) &
160                  *(teta(i,j,l)-tetabar(j,l))
161             end if
162             if ((v3d(i,j,l).eq.indefini).or. &
163                 ( vbar(j,l).eq.indefini).or. &
164                 (u3d(i,j,l).eq.indefini).or. &
165                 ( ubar(j,l).eq.indefini)) then
166               vpup(i,j,l)= indefini
167             else
168               vpup(i,j,l)=(v3d(i,j,l)-vbar(j,l))*(u3d(i,j,l)-ubar(j,l))
169             end if
170             if ((w3d(i,j,l).eq.indefini).or. &
171                 ( wbar(j,l).eq.indefini).or. &
172                 (u3d(i,j,l).eq.indefini).or. &
173                 ( ubar(j,l).eq.indefini)) then
174               wpup(i,j,l)= indefini
175             else
176               wpup(i,j,l)=(w3d(i,j,l)-wbar(j,l))*(u3d(i,j,l)-ubar(j,l))
177             end if
178             if ((v3d(i,j,l).eq.indefini).or. &
179                 ( vbar(j,l).eq.indefini)) then
180               vpvp(i,j,l)= indefini
181             else
182               vpvp(i,j,l)=(v3d(i,j,l)-vbar(j,l))*(v3d(i,j,l)-vbar(j,l))
183             end if
184             if ((w3d(i,j,l).eq.indefini).or. &
185                 ( wbar(j,l).eq.indefini).or. &
186                 (v3d(i,j,l).eq.indefini).or. &
187                 ( vbar(j,l).eq.indefini)) then
188               wpvp(i,j,l)= indefini
189             else
190               wpvp(i,j,l)=(w3d(i,j,l)-wbar(j,l))*(v3d(i,j,l)-vbar(j,l))
191             end if
192           end do
193         end do
194      end do
195
196!     Moyennes zonales des termes non zonaux
197!     --------------------------------------
198
199!     Termes de Reynolds
200
201!     flux zonaux de quantite de mouvement (vpupbar et wpupbar)
202!     flux meridiens de quantite de mouvement (vpvpbar et wpvpbar)
203!     flux meridien et vertical de chaleur (vptetapbar et wptetapbar)
204
205      call moyzon(iim,jjp1,llm,indefini,vptetap,vptetapbar)
206      call moyzon(iim,jjp1,llm,indefini,wptetap,wptetapbar)
207      call moyzon(iim,jjp1,llm,indefini,vpup,vpupbar)
208      call moyzon(iim,jjp1,llm,indefini,wpup,wpupbar)
209      call moyzon(iim,jjp1,llm,indefini,vpvp,vpvpbar)
210      call moyzon(iim,jjp1,llm,indefini,wpvp,wpvpbar)
211
212!     write(*,*) 'vptetapbar(2,23) =' , vptetapbar(2,23)
213!     write(*,*) 'wptetapbar(2,23) =' , wptetapbar(2,23)
214!     write(*,*) 'vpupbar(2,23) =' , vpupbar(2,23)
215
216 
217!     Derivees d/dp des moyennes zonales
218!     --------------------------------------
219
220      call dx_dp(jjp1,llm,indefini,press,ubar,dubardp)
221      call dx_dp(jjp1,llm,indefini,press,tetabar,dtetabardp)
222!     write (*,*) 'dtetabardp(l) (K/Pa)',(dtetabardp(6,l),l=1,llm)
223!       write(*,*)
224
225!     Controle (on triche !) de dtetabardp
226      do l=1,llm
227         do j=1,jjp1
228            if ((dtetabardp(j,l).gt.-1.e-3).and. &
229               (dtetabardp(j,l).ne.indefini)) then
230               dtetabardp(j,l) = -1.e-3
231!              write(*,*) 'profil presque instable en j = ',j,' l= ',l
232            end if
233         end do
234      end do
235!     write(*,*)'dtetabardp(l) (K/Pa) corr ',(dtetabardp(6,l),l=1,llm)
236!       write(*,*)
237
238
239!     calculs intermediaires
240!     ----------------------
241      do l=1,llm
242
243       if ((ubar(2,l).ne.indefini).and.(ubar(1,l).ne.indefini))then
244      ducosdlat(1,l) = ( ubar(2,l)*cos(rlatu(2))-ubar(1,l)*cos(rlatu(1)) ) &
245                      /( rlatu(2) - rlatu(1))
246       else
247      ducosdlat(1,l) = indefini
248       end if
249        do j=2,jjm
250       if ((ubar(j+1,l).ne.indefini).and.(ubar(j-1,l).ne.indefini))then
251      ducosdlat(j,l) = ( ubar(j+1,l)*cos(rlatu(j+1))-ubar(j-1,l)*cos(rlatu(j-1))) &
252                      /( rlatu(j+1) - rlatu(j-1))
253       else
254      ducosdlat(j,l) = indefini
255       end if
256        enddo
257       if ((ubar(jjp1,l).ne.indefini).and.(ubar(jjm,l).ne.indefini))then
258      ducosdlat(jjp1,l) = ( ubar(jjp1,l)*cos(rlatu(jjp1))  &
259                               -ubar(jjm,l)*cos(rlatu(jjm))  ) &
260                         /( rlatu(jjp1) - rlatu(jjm))
261       else
262      ducosdlat(jjp1,l) = indefini
263       end if
264           
265        do j=1,jjp1
266       if ((vptetapbar(j,l).ne.indefini).and. &
267           (dtetabardp(j,l).ne.indefini)) then
268           vz(j,l) = vptetapbar(j,l)/dtetabardp(j,l)
269         wlat(j,l) = cos(rlatu(j))*vz(j,l)
270       else
271           vz(j,l) = indefini
272         wlat(j,l) = indefini
273       end if
274        enddo
275
276      enddo
277             
278!      Calcul des vecteurs flux Eliassen Palm et divergence
279!      -----------------------------------------------------
280! ref: Read 1986
281
282      do l=1,llm
283        epy2d(1,l) = indefini
284        epz2d(1,l) = indefini
285        do j=2,jjm
286 
287!     Composante y
288!     ------------
289
290            if (   (rbar(j,l).ne.indefini).and. &
291                (vpupbar(j,l).ne.indefini).and. &
292                (dubardp(j,l).ne.indefini).and. &
293                     (vz(j,l).ne.indefini)) then
294           
295     epy2d(j,l) = rbar(j,l) * cos(rlatu(j))* ( vpupbar(j,l) &
296                                   - dubardp(j,l) * vz(j,l) & ! terme non geo
297                                             )       
298
299!     if ((j.eq.1).and.(l.eq.23))write(*,*) 'epy2d(1,23) =' , epy2d(1,23)
300!     if ((j.eq.2).and.(l.eq.23))write(*,*) 'epy2d(2,23) =' , epy2d(2,23)
301!     if ((j.eq.3).and.(l.eq.23))write(*,*) 'epy2d(3,23) =' , epy2d(3,23)
302           else
303     epy2d(j,l) = indefini
304            end if
305
306!     Composante z
307!     ------------
308
309            if (   (rbar(j,l).ne.indefini).and. &
310                (wpupbar(j,l).ne.indefini).and. &
311              (ducosdlat(j,l).ne.indefini).and. &
312                     (vz(j,l).ne.indefini)) then
313
314     epz2d(j,l) = rbar(j,l) * cos(rlatu(j))*                    &
315          (  vz(j,l) * ducosdlat(j,l)/(rbar(j,l)*cos(rlatu(j))) & ! terme non geo
316           - vz(j,l) * f(j)                                     &
317           + wpupbar(j,l)                                       &
318          )
319           else
320     epz2d(j,l) = indefini
321           end if
322
323!       if ((j.eq.2).and.(l.eq.23))write(*,*) 'epz2d(2,23) =' , epz2d(2,23)
324        end do
325        epy2d(jjp1,l) = indefini
326        epz2d(jjp1,l) = indefini
327      end do
328
329!    Moyennes euleriennes transformees (vtem2d et wtem2d)
330!    -----------------------------------------------
331     
332      call dx_dp(jjp1,llm,indefini,press,vz,dvzdp)
333
334      do l=1,llm
335
336          if ( (wlat(2,l).ne.indefini) .and. &
337               (wlat(1,l).ne.indefini) ) then
338               dwlatdlat(1,l)=(wlat(2,l)-wlat(1,l)) &
339                     /(rlatu(2)-rlatu(1))
340          else
341               dwlatdlat(1,l)=indefini
342          end if
343        do j=2,jjm
344          if ( (wlat(j+1,l).ne.indefini) .and. &
345               (wlat(j-1,l).ne.indefini) ) then
346               dwlatdlat(j,l)=(wlat(j+1,l)-wlat(j-1,l)) &
347                     /(rlatu(j+1)-rlatu(j-1))
348          else
349               dwlatdlat(j,l)=indefini
350          end if
351        enddo
352          if ( (wlat(jjp1,l).ne.indefini) .and. &
353               (wlat(jjm, l).ne.indefini) ) then
354               dwlatdlat(jjp1,l)=(wlat(jjp1,l)-wlat(jjm,l)) &
355                     /(rlatu(jjp1)-rlatu(jjm))
356          else
357               dwlatdlat(jjp1,l)=indefini
358          end if
359
360        do j=1,jjp1
361          if ( (dvzdp(j,l).ne.indefini) &
362           .and.(vbar(j,l).ne.indefini)) then
363            vtem2d(j,l) = vbar(j,l)-dvzdp(j,l)
364          else
365            vtem2d(j,l) = indefini
366          endif
367        enddo
368
369        do j=1,jjp1
370          if (   (rbar(j,l).ne.indefini).and. &
371            (dwlatdlat(j,l).ne.indefini).and. &
372                 (wbar(j,l).ne.indefini)) then
373           wtem2d(j,l) = wbar(j,l)+dwlatdlat(j,l)/(rbar(j,l)*cos(rlatu(j)))
374          else
375           wtem2d(j,l) = indefini
376          endif
377        enddo
378
379      enddo
380
381!      write(*,*) 'vtem2d(2,23) =' , vtem2d(2,23)
382!      write(*,*) 'wtem2d(2,23) =' , wtem2d(2,23)
383
384!      print*,"OK"
385
386!    Deviation par rapport a la circulation meridienne moyenne
387!    --------------------------------------------------------
388!                    U*.Nabla ubar
389      do l=1,llm
390        do j=1,jjp1
391          if (   (rbar(j,l).ne.indefini) &
392          .and.(vtem2d(j,l).ne.indefini) &
393          .and.(wtem2d(j,l).ne.indefini) &
394       .and.(ducosdlat(j,l).ne.indefini) &
395         .and.(dubardp(j,l).ne.indefini)) then
396       acc_meridien2d(j,l) =  &
397           vtem2d(j,l)*(ducosdlat(j,l)/(rbar(j,l)*cos(rlatu(j))) &
398                       -f(j)) &
399         + wtem2d(j,l)*dubardp(j,l)     
400          else
401       acc_meridien2d(j,l) = indefini
402          endif
403        enddo
404      enddo       
405
406!      write(*,*) 'acc_meridien2d(2,23) =' , acc_meridien2d(2,23)
407
408!     Divergence du flux
409!     ------------------
410      call dx_dp(jjp1,llm,indefini,press,epz2d,depzdp)
411!     write(*,*) 'depzdp(2,23)',depzdp(j,l)
412
413      do l=1,llm
414        divep2d(1,l) =0
415        do j=2,jjm
416          if (    (rbar(j,l).ne.indefini) .and. &
417               (epy2d(j+1,l).ne.indefini) .and. &
418               (epy2d(j-1,l).ne.indefini) .and. &
419                (depzdp(j,l).ne.indefini)    ) then
420
421           depycosdlat=(epy2d(j+1,l)*cos(rlatu(j+1))   &
422                      - epy2d(j-1,l)*cos(rlatu(j-1)) ) &
423                       /(rlatu(j+1) - rlatu(j-1))
424
425!      DIVERGENCE DU FLUX :
426
427           divep2d(j,l) = depycosdlat/(rbar(j,l)*cos(rlatu(j)))+depzdp(j,l)
428
429!          if ((j.eq.2).and.(l.eq.23)) then
430!              write(*,*) 'depycosdlat(2,23)', depycosdlat
431!              write(*,*) 'depzdp(2,23)',depzdp(j,l)
432!              write(*,*) 'divergence(2,23)',divep2d(j,l)
433!          end if
434
435! Wave driving (Du/Dt et non Dm/dt) :
436           divep2d(j,l) = divep2d(j,l)/(rbar(j,l)*cos(rlatu(j)))
437          else
438            divep2d(j,l) = indefini
439          end if
440        end do
441        divep2d(jjp1,l) =0
442      end do
443 
444!     write(*,*) ' fin :divep2d(2,23) =' , divep2d(2,23)
445
446! Preparation pour sortie graphique (flux ``chapeau'' Peixoto et Oort p 391)
447! A VOIR...
448!      do l=1,llm
449!        do j=1,jjp1
450!           if (epy2d(j,l).ne.indefini) &
451!           epy2d(j,l) = 2*pi*a0    *(1./g0)*cos(rlatu(j))*epy2d(j,l)
452!           if (epz2d(j,l).ne.indefini) &
453!           epz2d(j,l) = 2*pi*a0**2 *(1./g0)*cos(rlatu(j))*epz2d(j,l)
454!        end do
455!      end do
456
457      return
458      end
Note: See TracBrowser for help on using the repository browser.