source: trunk/LMDZ.GENERIC/libf/dyn3d/amont_qsat.F @ 487

Last change on this file since 487 was 135, checked in by aslmd, 14 years ago

CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET

File size: 4.3 KB
RevLine 
[135]1      SUBROUTINE amont_qsat (nq,iq,q,teta,pk,w, pbaru, pbarv,dq)
2      IMPLICIT NONE
3
4c=======================================================================
5c
6c   Auteur:  P. Le Van / F.Forget
7c   -------
8c
9c   ********************************************************************
10c   Transport d'un traceur q par shema amont 3D  avec contrainte thermi-
11c           -que
12c   ********************************************************************
13c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
14c     dq               sont des arguments de sortie pour le s-pg ....
15c
16c=======================================================================
17
18
19#include "dimensions.h"
20#include "paramet.h"
21#include "logic.h"
22#include "comvert.h"
23#include "comconst.h"
24#include "comgeom.h"
25#include "serre.h"
26
27c   Arguments:
28c   ----------
29      INTEGER nq,iq 
30      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
31      REAL q(ip1jmp1,llm,nq), dq( ip1jmp1,llm,nq )
32      REAL teta(ip1jmp1,llm),pk(ip1jmp1,llm)
33      REAL w(ip1jmp1,llm)
34
35c   Local:
36c   ------
37
38      INTEGER   l,ij
39
40      REAL qbyv( ip1jm,llm ), qbxu( ip1jmp1,llm )
41      REAL ww
42      REAL dqh( ip1jmp1,llm)
43
44      REAL ptarg,pdelarg,foeew,zdelta
45      REAL rtt,retv,r3les,r3ies,r4les,r4ies,unskappa,play
46      REAL tempe(ip1jmp1),qsat(ip1jmp1,llm),r2es
47
48      EXTERNAL     convflu
49
50       FOEEW ( PTARG,PDELARG ) = EXP (
51     *          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT)
52     * / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
53C
54
55
56c     ----------------------------------------------------------------
57c     --        calcul transport horiz. pour  dq                    --
58c     -                                                              -
59c     -         shema ' amont '                                      -
60c     ----------------------------------------------------------------
61 
62        r2es  = 380.11733
63        r3les = 17.269
64        r3ies = 21.875
65        r4les = 35.86
66        r4ies = 7.66
67        retv = 0.6077667
68        rtt  = 273.16
69 
70        unskappa =  1./ kappa
71
72        DO l = 1, llm
73         DO ij = 1, ip1jmp1
74          tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
75         ENDDO
76         DO ij = 1, ip1jmp1
77          zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
78          play   =  preff * (pk(ij,l)/cpp) ** unskappa
79          qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
80          qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
81         ENDDO
82        ENDDO
83
84        DO 10 l = 1, llm
85
86        DO 6  ij     = iip2, ip1jm - 1
87         IF( pbaru(ij,l).GT.0  )  THEN
88          qbxu( ij,l ) = pbaru( ij,l ) * MIN( q(ij,l,iq),qsat(ij +1,l))
89         ELSE
90          qbxu( ij,l ) = pbaru( ij,l ) * MIN( q(ij+1,l,iq),qsat(ij,l) )
91         ENDIF
92       
93   6    CONTINUE
94 
95c     ..... correction  pour  qbxu(iip1,j,l)   .....
96c     ...   qbxu(iip1,j,l)= qbxu(1,j,l)  ...
97 
98CDIR$ IVDEP
99        DO  7   ij   = iip1 +iip1, ip1jm, iip1
100        qbxu( ij,l ) = qbxu( ij - iim, l )
101   7    CONTINUE
102 
103        DO  8  ij     = 1, ip1jm
104         IF( pbarv(ij,l).GT.0  )  THEN
105          qbyv(ij,l) = pbarv( ij,l ) * MIN( q(ij+iip1,l,iq),qsat(ij,l) )
106         ELSE
107          qbyv(ij,l) = pbarv( ij,l ) * MIN( q(ij,l,iq),qsat(ij+iip1,l) )
108         ENDIF
109   8    CONTINUE
110c
111  10    CONTINUE
112 
113c     stockage dans  dqh de la convergence horiz.du flux d'humidite  .
114c     -------------------------------------------------------------
115c
116                  CALL convflu(qbxu, qbyv, llm, dqh )
117
118
119c ---------------------------------------------------------------
120c   .... calcul des termes d'advection vertic.pour q. Shema amont
121c ---------------------------------------------------------------
122
123c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dqh pour calculer dq
124c
125
126       DO 20 l = 1,llmm1
127c
128         DO 11 ij = 1,ip1jmp1
129           IF( w(ij,l+1).LT.0. )   THEN
130             ww =  -w(ij,l+1) * MIN( q(ij, l ,iq),qsat(ij ,l+1) )
131           ELSE
132             ww =  -w(ij,l+1) * MIN( q(ij,l+1,iq),qsat(ij , l ) )
133           ENDIF
134
135cc     .... Modif . P. Le Van  (26/09/95) ....
136
137          dq (ij, l ,iq ) = dqh(ij, l )   -  ww
138          dqh(ij,l+1    ) = dqh(ij,l+1)   +  ww
139  11     CONTINUE
140c
141  20   CONTINUE
142
143c
144c       special dq (ij,llm) (qui n'a pas ete rempli ! )
145c
146        DO  ij = 1,ip1jmp1
147          dq( ij,llm,iq ) = dqh( ij,llm )
148        END DO
149           
150      RETURN
151      END
Note: See TracBrowser for help on using the repository browser.