/[lmdze]/trunk/dyn3d/vlspltqs.f
ViewVC logotype

Diff of /trunk/dyn3d/vlspltqs.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 43 by guez, Fri Apr 8 12:43:31 2011 UTC revision 44 by guez, Wed Apr 13 12:29:18 2011 UTC
# Line 1  Line 1 
1  !  SUBROUTINE vlspltqs(q, pente_max, masse, w, pbaru, pbarv, pdt, p, pk, teta)
 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/vlspltqs.F,v 1.2 2005/02/24 12:16:57 fairhead Exp $  
 !  
        SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt, &  
                                         p,pk,teta                 )  
 !  
 !     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron  
 !  
 !    ********************************************************************  
 !          Shema  d'advection " pseudo amont " .  
 !      + test sur humidite specifique: Q advecte< Qsat aval  
 !                   (F. Codron, 10/99)  
 !    ********************************************************************  
 !     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....  
 !  
 !     pente_max facteur de limitation des pentes: 2 en general  
 !                                                0 pour un schema amont  
 !     pbaru,pbarv,w flux de masse en u ,v ,w  
 !     pdt pas de temps  
 !  
 !     teta temperature potentielle, p pression aux interfaces,  
 !     pk exner au milieu des couches necessaire pour calculer Qsat  
 !   --------------------------------------------------------------------  
       use dimens_m  
       use paramet_m  
       use comconst  
       use comvert  
       use logic  
       IMPLICIT NONE  
 !  
   
 !  
 !   Arguments:  
 !   ----------  
       REAL masse(ip1jmp1,llm),pente_max  
       REAL, intent(in):: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)  
       REAL q(ip1jmp1,llm)  
       REAL w(ip1jmp1,llm)  
       real, intent(in):: pdt  
       REAL, intent(in):: p(ip1jmp1,llmp1)  
       real teta(ip1jmp1,llm),pk(ip1jmp1,llm)  
 !  
 !      Local  
 !   ---------  
 !  
       INTEGER i,ij,l,j,ii  
 !  
       REAL qsat(ip1jmp1,llm)  
       REAL zm(ip1jmp1,llm)  
       REAL mu(ip1jmp1,llm)  
       REAL mv(ip1jm,llm)  
       REAL mw(ip1jmp1,llm+1)  
       REAL zq(ip1jmp1,llm)  
       REAL temps1,temps2,temps3  
       REAL zzpbar, zzw  
       LOGICAL testcpu  
       SAVE testcpu  
       SAVE temps1,temps2,temps3  
   
       REAL qmin,qmax  
       DATA qmin,qmax/0.,1.e33/  
       DATA testcpu/.false./  
       DATA temps1,temps2,temps3/0.,0.,0./  
   
 !--pour rapport de melange saturant--  
   
       REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play  
       REAL ptarg,pdelarg,foeew,zdelta  
       REAL tempe(ip1jmp1)  
   
 !    fonction psat(T)  
   
        FOEEW ( PTARG,PDELARG ) = EXP ( &  
                 (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) &  
        / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )  
   
         r2es  = 380.11733  
         r3les = 17.269  
         r3ies = 21.875  
         r4les = 35.86  
         r4ies = 7.66  
         retv = 0.6077667  
         rtt  = 273.16  
   
 !-- Calcul de Qsat en chaque point  
 !-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2  
 !   pour eviter une exponentielle.  
         DO l = 1, llm  
          DO ij = 1, ip1jmp1  
           tempe(ij) = teta(ij,l) * pk(ij,l) /cpp  
          ENDDO  
          DO ij = 1, ip1jmp1  
           zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )  
           play   = 0.5*(p(ij,l)+p(ij,l+1))  
           qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )  
           qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )  
          ENDDO  
         ENDDO  
   
         zzpbar = 0.5 * pdt  
         zzw    = pdt  
       DO l=1,llm  
         DO ij = iip2,ip1jm  
             mu(ij,l)=pbaru(ij,l) * zzpbar  
          ENDDO  
          DO ij=1,ip1jm  
             mv(ij,l)=pbarv(ij,l) * zzpbar  
          ENDDO  
          DO ij=1,ip1jmp1  
             mw(ij,l)=w(ij,l) * zzw  
          ENDDO  
       ENDDO  
   
       DO ij=1,ip1jmp1  
          mw(ij,llm+1)=0.  
       ENDDO  
   
       CALL SCOPY(ijp1llm,q,1,zq,1)  
       CALL SCOPY(ijp1llm,masse,1,zm,1)  
   
 !      call minmaxq(zq,qmin,qmax,'avant vlxqs     ')  
       call vlxqs(zq,pente_max,zm,mu,qsat)  
   
   
 !     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')  
   
       call vlyqs(zq,pente_max,zm,mv,qsat)  
   
   
 !      call minmaxq(zq,qmin,qmax,'avant vlz     ')  
   
       call vlz(zq,pente_max,zm,mw)  
   
   
 !     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')  
 !     call minmaxq(zm,qmin,qmax,'M avant vlyqs     ')  
   
       call vlyqs(zq,pente_max,zm,mv,qsat)  
   
   
 !     call minmaxq(zq,qmin,qmax,'avant vlxqs     ')  
 !     call minmaxq(zm,qmin,qmax,'M avant vlxqs     ')  
   
       call vlxqs(zq,pente_max,zm,mu,qsat)  
   
 !     call minmaxq(zq,qmin,qmax,'apres vlxqs     ')  
 !     call minmaxq(zm,qmin,qmax,'M apres vlxqs     ')  
   
   
       DO l=1,llm  
          DO ij=1,ip1jmp1  
            q(ij,l)=zq(ij,l)  
          ENDDO  
          DO ij=1,ip1jm+1,iip1  
             q(ij+iim,l)=q(ij,l)  
          ENDDO  
       ENDDO  
2    
3        RETURN    ! From LMDZ4/libf/dyn3d/vlspltqs.F, version 1.2 2005/02/24 12:16:57 fairhead
4        END  
5      ! Authors: P. Le Van, F. Hourdin, F. Forget, F. Codron
6    
7      ! Schéma d'advection "pseudo amont"
8      ! + test sur humidité spécifique : Q advecté < Qsat aval
9      ! (F. Codron, 10/99)
10    
11      ! q, pbaru, pbarv, w sont des arguments d'entree pour le sous-programme
12    
13      ! pente_max facteur de limitation des pentes: 2 en général
14      ! 0 pour un schéma amont
15      ! pbaru, pbarv, w flux de masse en u , v , w
16      ! pdt pas de temps
17    
18      ! teta température potentielle, p pression aux interfaces,
19      ! pk exner au milieu des couches nécessaire pour calculer Qsat
20    
21      USE dimens_m, ONLY : iim, llm
22      USE paramet_m, ONLY : iip1, iip2, ijp1llm, ip1jm, ip1jmp1, llmp1
23      USE comconst, ONLY : cpp
24    
25      IMPLICIT NONE
26    
27      ! Arguments:
28    
29      REAL masse(ip1jmp1, llm), pente_max
30      REAL, intent(in):: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
31      REAL q(ip1jmp1, llm)
32      REAL w(ip1jmp1, llm)
33      real, intent(in):: pdt
34      REAL, intent(in):: p(ip1jmp1, llmp1)
35      real, intent(in):: teta(ip1jmp1, llm)
36      real pk(ip1jmp1, llm)
37    
38      ! Local
39    
40      INTEGER i, ij, l, j, ii
41    
42      REAL qsat(ip1jmp1, llm)
43      REAL zm(ip1jmp1, llm)
44      REAL mu(ip1jmp1, llm)
45      REAL mv(ip1jm, llm)
46      REAL mw(ip1jmp1, llm+1)
47      REAL zq(ip1jmp1, llm)
48      REAL temps1, temps2, temps3
49      REAL zzpbar, zzw
50      LOGICAL testcpu
51      SAVE testcpu
52      SAVE temps1, temps2, temps3
53    
54      REAL qmin, qmax
55      DATA qmin, qmax/0., 1.e33/
56      DATA testcpu/.false./
57      DATA temps1, temps2, temps3/0., 0., 0./
58    
59      !--pour rapport de melange saturant--
60    
61      REAL rtt, retv, r2es, r3les, r3ies, r4les, r4ies, play
62      REAL ptarg, pdelarg, foeew, zdelta
63      REAL tempe(ip1jmp1)
64    
65      ! fonction psat(T)
66      FOEEW(PTARG, PDELARG) = EXP ((R3LES*(1.-PDELARG)+R3IES*PDELARG) &
67           * (PTARG-RTT) / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)))
68    
69      !------------------------------------------------------------------
70    
71      r2es = 380.11733
72      r3les = 17.269
73      r3ies = 21.875
74      r4les = 35.86
75      r4ies = 7.66
76      retv = 0.6077667
77      rtt = 273.16
78    
79      !-- Calcul de Qsat en chaque point
80      !-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
81      ! pour eviter une exponentielle.
82      DO l = 1, llm
83         DO ij = 1, ip1jmp1
84            tempe(ij) = teta(ij, l) * pk(ij, l) /cpp
85         ENDDO
86         DO ij = 1, ip1jmp1
87            zdelta = MAX(0., SIGN(1., rtt - tempe(ij)))
88            play = 0.5*(p(ij, l)+p(ij, l+1))
89            qsat(ij, l) = MIN(0.5, r2es* FOEEW(tempe(ij), zdelta) / play)
90            qsat(ij, l) = qsat(ij, l) / (1. - retv * qsat(ij, l))
91         ENDDO
92      ENDDO
93    
94      zzpbar = 0.5 * pdt
95      zzw = pdt
96      DO l=1, llm
97         DO ij = iip2, ip1jm
98            mu(ij, l)=pbaru(ij, l) * zzpbar
99         ENDDO
100         DO ij=1, ip1jm
101            mv(ij, l)=pbarv(ij, l) * zzpbar
102         ENDDO
103         DO ij=1, ip1jmp1
104            mw(ij, l)=w(ij, l) * zzw
105         ENDDO
106      ENDDO
107    
108      DO ij=1, ip1jmp1
109         mw(ij, llm+1)=0.
110      ENDDO
111    
112      CALL SCOPY(ijp1llm, q, 1, zq, 1)
113      CALL SCOPY(ijp1llm, masse, 1, zm, 1)
114    
115      ! call minmaxq(zq, qmin, qmax, 'avant vlxqs ')
116      call vlxqs(zq, pente_max, zm, mu, qsat)
117    
118      ! call minmaxq(zq, qmin, qmax, 'avant vlyqs ')
119    
120      call vlyqs(zq, pente_max, zm, mv, qsat)
121    
122      ! call minmaxq(zq, qmin, qmax, 'avant vlz ')
123    
124      call vlz(zq, pente_max, zm, mw)
125    
126      ! call minmaxq(zq, qmin, qmax, 'avant vlyqs ')
127      ! call minmaxq(zm, qmin, qmax, 'M avant vlyqs ')
128    
129      call vlyqs(zq, pente_max, zm, mv, qsat)
130    
131      ! call minmaxq(zq, qmin, qmax, 'avant vlxqs ')
132      ! call minmaxq(zm, qmin, qmax, 'M avant vlxqs ')
133    
134      call vlxqs(zq, pente_max, zm, mu, qsat)
135    
136      ! call minmaxq(zq, qmin, qmax, 'apres vlxqs ')
137      ! call minmaxq(zm, qmin, qmax, 'M apres vlxqs ')
138    
139      DO l=1, llm
140         DO ij=1, ip1jmp1
141            q(ij, l)=zq(ij, l)
142         ENDDO
143         DO ij=1, ip1jm+1, iip1
144            q(ij+iim, l)=q(ij, l)
145         ENDDO
146      ENDDO
147    
148    END SUBROUTINE vlspltqs
Legend:
Removed from v.43  
changed lines
  Added in v.44
  ViewVC Help
Powered by ViewVC 1.1.21