source: LMDZ4/trunk/libf/phytherm/thermcell_env.F90 @ 815

Last change on this file since 815 was 814, checked in by Laurent Fairhead, 17 years ago

Rajout de la physique utilisant les thermiques FH
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.0 KB
Line 
1      SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
2     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
3
4!--------------------------------------------------------------
5!thermcell_env: calcule les caracteristiques de l environnement
6!necessaires au calcul des proprietes dans le thermique
7!--------------------------------------------------------------
8
9      IMPLICIT NONE
10
11#include "YOMCST.h"
12#include "YOETHF.h"
13#include "FCTTRE.h"     
14
15      INTEGER ngrid,nlay
16      REAL po(ngrid,nlay)
17      REAL pt(ngrid,nlay)
18      REAL pu(ngrid,nlay)
19      REAL pv(ngrid,nlay)
20      REAL pplay(ngrid,nlay)
21      REAL pplev(ngrid,nlay+1)
22      integer lev_out                           ! niveau pour les print
23
24      REAL zo(ngrid,nlay)
25      REAL zl(ngrid,nlay)
26      REAL zh(ngrid,nlay)
27      REAL ztv(ngrid,nlay)
28      REAL zthl(ngrid,nlay)
29      REAL zpspsk(ngrid,nlay)
30      REAL zu(ngrid,nlay)
31      REAL zv(ngrid,nlay)
32      REAL zqsat(ngrid,nlay)
33
34      INTEGER ig,l,ll
35
36      real zcor,zdelta,zcvm5,qlbef
37      real Tbef,qsatbef
38      real dqsat_dT,DT,num,denom
39      REAL RLvCp,DDT0
40      PARAMETER (DDT0=.01)
41      LOGICAL Zsat
42
43      Zsat=.false.
44!
45! Pr Tprec=Tl calcul de qsat
46! Si qsat>qT T=Tl, q=qT
47! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
48! On cherche DDT < DDT0
49!
50! calcul des caracteristiques de l environnement
51       DO  ll=1,nlay
52         DO ig=1,ngrid
53            zo(ig,ll)=po(ig,ll)
54            zl(ig,ll)=0.
55            zh(ig,ll)=pt(ig,ll)
56            zqsat(ig,ll)=0.
57         EndDO
58       EndDO
59!
60!
61!recherche de saturation dans l environnement
62       DO ll=1,nlay
63! les points insatures sont definitifs
64         DO ig=1,ngrid
65            Tbef=pt(ig,ll)
66            zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
67            qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll)
68            qsatbef=MIN(0.5,qsatbef)
69            zcor=1./(1.-retv*qsatbef)
70            qsatbef=qsatbef*zcor
71            Zsat = (max(0.,po(ig,ll)-qsatbef) .gt. 1.e-10)
72            if (Zsat) then
73            qlbef=max(0.,po(ig,ll)-qsatbef)
74! si sature: ql est surestime, d'ou la sous-relax
75            DT = 0.5*RLvCp*qlbef
76! on pourra enchainer 2 ou 3 calculs sans Do while
77            do while (abs(DT).gt.DDT0)
78! il faut verifier si c,a conserve quand on repasse en insature ...
79              Tbef=Tbef+DT
80              zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
81              qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll)
82              qsatbef=MIN(0.5,qsatbef)
83              zcor=1./(1.-retv*qsatbef)
84              qsatbef=qsatbef*zcor
85! on veut le signe de qlbef
86              qlbef=po(ig,ll)-qsatbef
87!          dqsat_dT
88              zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
89              zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
90              zcor=1./(1.-retv*qsatbef)
91              dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
92              num=-Tbef+pt(ig,ll)+RLvCp*qlbef
93              denom=1.+RLvCp*dqsat_dT
94              if (denom.lt.1.e-10) then
95                  print*,'pb denom'
96              endif
97              DT=num/denom
98            enddo
99! on ecrit de maniere conservative (sat ou non)
100            zl(ig,ll) = max(0.,qlbef)
101!          T = Tl +Lv/Cp ql
102            zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)
103            zo(ig,ll) = po(ig,ll)-zl(ig,ll)
104           endif
105!on ecrit zqsat
106            zqsat(ig,ll)=qsatbef     
107         EndDO
108       EndDO
109!
110!
111!-----------------------------------------------------------------------
112!   incrementation eventuelle de tendances precedentes:
113!   ---------------------------------------------------
114
115      if (lev_out.ge.1) print*,'0 OK convect8'
116
117      DO 1010 l=1,nlay
118         DO 1015 ig=1,ngrid
119             zpspsk(ig,l)=(pplay(ig,l)/100000.)**RKAPPA
120             zu(ig,l)=pu(ig,l)
121             zv(ig,l)=pv(ig,l)
122!attention zh est maintenant le profil de T et plus le profil de theta !
123!
124!   T-> Theta
125            ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
126!Theta_v
127            ztv(ig,l)=ztv(ig,l)*(1.+RETV*(zo(ig,l))  &
128     &           -zl(ig,l))
129!Thetal
130            zthl(ig,l)=pt(ig,l)/zpspsk(ig,l)
131!           
1321015     CONTINUE
1331010  CONTINUE
134 
135      RETURN
136      END
Note: See TracBrowser for help on using the repository browser.