source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/thermcell_qsat.F90 @ 4203

Last change on this file since 4203 was 3035, checked in by fhourdin, 7 years ago

Ajout de THREADPRIVATE (non indispensables)
Et un peu de nettoyage pour yamada4
  • 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
Line 
1subroutine thermcell_qsat(klon,active,pplev,ztemp,zqta,zqsat)
2implicit none
3
4#include "YOMCST.h"
5#include "YOETHF.h"
6#include "FCTTRE.h"
7
8
9!====================================================================
10! DECLARATIONS
11!====================================================================
12
13! Arguments
14INTEGER klon
15REAL zpspsk(klon),pplev(klon)
16REAL ztemp(klon),zqta(klon),zqsat(klon)
17LOGICAL active(klon)
18
19! Variables locales
20INTEGER ig,iter
21REAL Tbef(klon),DT(klon)
22REAL tdelta,qsatbef,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT
23logical Zsat
24REAL RLvCp
25
26REAL, SAVE :: DDT0=.01
27!$OMP THREADPRIVATE(DDT0)
28
29LOGICAL afaire(klon),tout_converge
30
31!====================================================================
32! INITIALISATIONS
33!====================================================================
34
35RLvCp = RLVTT/RCPD
36tout_converge=.false.
37afaire(:)=.false.
38DT(:)=0.
39
40
41!====================================================================
42! Routine a vectoriser en copiant active dans converge et en mettant
43! la boucle sur les iterations a l'exterieur est en mettant
44! converge= false des que la convergence est atteinte.
45!====================================================================
46
47do ig=1,klon
48   if (active(ig)) then
49               Tbef(ig)=ztemp(ig)
50               zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
51               qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
52               qsatbef=MIN(0.5,qsatbef)
53               zcor=1./(1.-retv*qsatbef)
54               qsatbef=qsatbef*zcor
55               qlbef=max(0.,zqta(ig)-qsatbef)
56               DT(ig) = 0.5*RLvCp*qlbef
57               zqsat(ig)=qsatbef
58     endif
59enddo
60
61! Traitement du cas ou il y a condensation mais faible
62! On ne condense pas mais on dit que le qsat est le qta
63do ig=1,klon
64   if (active(ig)) then
65      if (0.<abs(DT(ig)).and.abs(DT(ig))<=DDT0) then
66         zqsat(ig)=zqta(ig)
67      endif
68   endif
69enddo
70
71do iter=1,10
72    afaire(:)=abs(DT(:)).gt.DDT0
73    do ig=1,klon
74               if (afaire(ig)) then
75                 Tbef(ig)=Tbef(ig)+DT(ig)
76                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
77                 qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
78                 qsatbef=MIN(0.5,qsatbef)
79                 zcor=1./(1.-retv*qsatbef)
80                 qsatbef=qsatbef*zcor
81                 qlbef=zqta(ig)-qsatbef
82                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
83                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
84                 zcor=1./(1.-retv*qsatbef)
85                 dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef,zcor)
86                 num=-Tbef(ig)+ztemp(ig)+RLvCp*qlbef
87                 denom=1.+RLvCp*dqsat_dT
88                 zqsat(ig) = qsatbef
89                 DT(ig)=num/denom
90               endif
91    enddo
92enddo
93
94return
95end
Note: See TracBrowser for help on using the repository browser.