source: LMDZ4/trunk/libf/phylmd/thermcell_env.F90 @ 1080

Last change on this file since 1080 was 970, checked in by Laurent Fairhead, 16 years ago

Corrections petits bugs divers
LF

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