source: LMDZ5/branches/testing/libf/phylmd/thermcell_qsat.F90 @ 3946

Last change on this file since 3946 was 1910, checked in by Laurent Fairhead, 11 years ago

Merged trunk changes r1860:1909 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 2.7 KB
RevLine 
[1338]1subroutine thermcell_qsat(klon,active,pplev,ztemp,zqta,zqsat)
[1336]2implicit none
3
4#include "YOMCST.h"
5#include "YOETHF.h"
6#include "FCTTRE.h"
7
[1338]8
[1336]9!====================================================================
10! DECLARATIONS
11!====================================================================
12
13! Arguments
14INTEGER klon
[1338]15REAL zpspsk(klon),pplev(klon)
16REAL ztemp(klon),zqta(klon),zqsat(klon)
17LOGICAL active(klon)
[1336]18
19! Variables locales
20INTEGER ig,iter
21REAL Tbef(klon),DT(klon)
[1338]22REAL tdelta,qsatbef,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT
[1336]23logical Zsat
24REAL RLvCp
25REAL, SAVE :: DDT0=.01
26LOGICAL afaire(klon),tout_converge
27
28!====================================================================
29! INITIALISATIONS
30!====================================================================
31
32RLvCp = RLVTT/RCPD
33tout_converge=.false.
[1338]34afaire(:)=.false.
[1336]35DT(:)=0.
36
37
38!====================================================================
[1338]39! Routine a vectoriser en copiant active dans converge et en mettant
[1336]40! la boucle sur les iterations a l'exterieur est en mettant
41! converge= false des que la convergence est atteinte.
42!====================================================================
43
44do ig=1,klon
[1338]45   if (active(ig)) then
46               Tbef(ig)=ztemp(ig)
[1336]47               zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
[1338]48               qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
49               qsatbef=MIN(0.5,qsatbef)
50               zcor=1./(1.-retv*qsatbef)
51               qsatbef=qsatbef*zcor
52               qlbef=max(0.,zqta(ig)-qsatbef)
[1336]53               DT(ig) = 0.5*RLvCp*qlbef
[1338]54               zqsat(ig)=qsatbef
[1336]55     endif
56enddo
57
[1338]58! Traitement du cas ou il y a condensation mais faible
59! On ne condense pas mais on dit que le qsat est le qta
60do ig=1,klon
61   if (active(ig)) then
62      if (0.<abs(DT(ig)).and.abs(DT(ig))<=DDT0) then
63         zqsat(ig)=zqta(ig)
64      endif
65   endif
66enddo
[1336]67
68do iter=1,10
[1338]69    afaire(:)=abs(DT(:)).gt.DDT0
[1336]70    do ig=1,klon
71               if (afaire(ig)) then
72                 Tbef(ig)=Tbef(ig)+DT(ig)
73                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
[1338]74                 qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
75                 qsatbef=MIN(0.5,qsatbef)
76                 zcor=1./(1.-retv*qsatbef)
77                 qsatbef=qsatbef*zcor
78                 qlbef=zqta(ig)-qsatbef
[1336]79                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
80                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
[1338]81                 zcor=1./(1.-retv*qsatbef)
82                 dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef,zcor)
83                 num=-Tbef(ig)+ztemp(ig)+RLvCp*qlbef
[1336]84                 denom=1.+RLvCp*dqsat_dT
[1338]85                 zqsat(ig) = qsatbef
[1336]86                 DT(ig)=num/denom
87               endif
88    enddo
89enddo
90
91return
92end
Note: See TracBrowser for help on using the repository browser.