source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/tlift.F90 @ 2404

Last change on this file since 2404 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

  • 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
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.7 KB
Line 
1
2! $Header$
3
4SUBROUTINE tlift(p, t, rr, rs, gz, plcl, icb, nk, tvp, tpk, clw, nd, nl, &
5    dtvpdt1, dtvpdq1)
6
7  ! Argument NK ajoute (jyg) = Niveau de depart de la
8  ! convection
9
10  PARAMETER (na=60)
11  REAL gz(nd), tpk(nd), clw(nd)
12  REAL t(nd), rr(nd), rs(nd), tvp(nd), p(nd)
13  REAL dtvpdt1(nd), dtvpdq1(nd) ! Derivatives of parcel virtual
14  ! temperature wrt T1 and Q1
15
16  REAL clw_new(na), qi(na)
17  REAL dtpdt1(na), dtpdq1(na) ! Derivatives of parcel temperature
18  ! wrt T1 and Q1
19
20
21  LOGICAL ice_conv
22
23  ! ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***
24
25  ! sb        CPD=1005.7
26  ! sb      CPV=1870.0
27  ! sb        CL=4190.0
28  ! sb        CPVMCL=2320.0
29  ! sb        RV=461.5
30  ! sb        RD=287.04
31  ! sb        EPS=RD/RV
32  ! sb        ALV0=2.501E6
33  ! cccccccccccccccccccccc
34  ! constantes coherentes avec le modele du Centre Europeen
35  ! sb      RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644
36  ! sb      RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153
37  ! sb      CPD = 3.5 * RD
38  ! sb      CPV = 4.0 * RV
39  ! sb      CL = 4218.0
40  ! sb      CI=2090.0
41  ! sb      CPVMCL=CL-CPV
42  ! sb      CLMCI=CL-CI
43  ! sb      EPS=RD/RV
44  ! sb      ALV0=2.5008E+06
45  ! sb      ALF0=3.34E+05
46
47  ! ccccccccccc
48  ! on utilise les constantes thermo du Centre Europeen: (SB)
49
50  include "YOMCST.h"
51  gravity = rg !sb: Pr que gravite ne devienne pas humidite!
52
53  cpd = rcpd
54  cpv = rcpv
55  cl = rcw
56  ci = rcs
57  cpvmcl = cl - cpv
58  clmci = cl - ci
59  eps = rd/rv
60  alv0 = rlvtt
61  alf0 = rlmlt ! (ALF0 = RLSTT-RLVTT)
62
63  ! ccccccccccccccccccccc
64
65  ! ***  CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY   ***
66
67  icb1 = max(icb, 2)
68  icb1 = min(icb, nl)
69
70  ! jyg1
71  ! C      CPP=CPD*(1.-RR(1))+RR(1)*CPV
72  cpp = cpd*(1.-rr(nk)) + rr(nk)*cpv
73  ! jyg2
74  cpinv = 1./cpp
75  ! jyg1
76  ! ICB may be below condensation level
77  ! CC        DO 100 I=1,ICB1-1
78  ! CC         TPK(I)=T(1)-GZ(I)*CPINV
79  ! CC         TVP(I)=TPK(I)*(1.+RR(1)/EPS)
80  DO i = 1, icb1
81    clw(i) = 0.0
82  END DO
83
84  DO i = nk, icb1
85    tpk(i) = t(nk) - (gz(i)-gz(nk))*cpinv
86    ! jyg1
87    ! CC         TVP(I)=TPK(I)*(1.+RR(NK)/EPS)
88    tvp(i) = tpk(i)*(1.+rr(nk)/eps-rr(nk))
89    ! jyg2
90    dtvpdt1(i) = 1. + rr(nk)/eps - rr(nk)
91    dtvpdq1(i) = tpk(i)*(1./eps-1.)
92
93    ! jyg2
94
95  END DO
96
97
98  ! ***  FIND LIFTED PARCEL TEMPERATURE AND MIXING RATIO    ***
99
100  ! jyg1
101  ! C        AH0=(CPD*(1.-RR(1))+CL*RR(1))*T(1)
102  ! C     $     +RR(1)*(ALV0-CPVMCL*(T(1)-273.15))
103  ah0 = (cpd*(1.-rr(nk))+cl*rr(nk))*t(nk) + rr(nk)*(alv0-cpvmcl*(t(nk)-273.15 &
104    )) + gz(nk)
105  ! jyg2
106
107  ! jyg1
108  imin = icb1
109  ! If ICB is below LCL, start loop at ICB+1
110  IF (plcl<p(icb1)) imin = min(imin+1, nl)
111
112  ! CC        DO 300 I=ICB1,NL
113  DO i = imin, nl
114    ! jyg2
115    alv = alv0 - cpvmcl*(t(i)-273.15)
116    alf = alf0 + clmci*(t(i)-273.15)
117
118    rg = rs(i)
119    tg = t(i)
120    ! S=CPD+ALV*ALV*RG/(RV*T(I)*T(I))
121    ! jyg1
122    ! C        S=CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*RG/(RV*T(I)*T(I))
123    s = cpd*(1.-rr(nk)) + cl*rr(nk) + alv*alv*rg/(rv*t(i)*t(i))
124    ! jyg2
125    s = 1./s
126
127    DO j = 1, 2
128      ! jyg1
129      ! C         AHG=CPD*TG+(CL-CPD)*RR(1)*TG+ALV*RG+GZ(I)
130      ahg = cpd*tg + (cl-cpd)*rr(nk)*tg + alv*rg + gz(i)
131      ! jyg2
132      tg = tg + s*(ah0-ahg)
133      tc = tg - 273.15
134      denom = 243.5 + tc
135      denom = max(denom, 1.0)
136
137      ! FORMULE DE BOLTON POUR PSAT
138
139      es = 6.112*exp(17.67*tc/denom)
140      rg = eps*es/(p(i)-es*(1.-eps))
141
142
143    END DO
144
145    ! jyg1
146    ! C        TPK(I)=(AH0-GZ(I)-ALV*RG)/(CPD+(CL-CPD)*RR(1))
147    tpk(i) = (ah0-gz(i)-alv*rg)/(cpd+(cl-cpd)*rr(nk))
148    ! jyg2
149    ! TPK(I)=(AH0-GZ(I)-ALV*RG-(CL-CPD)*T(I)*RR(1))/CPD
150
151    ! jyg1
152    ! C        CLW(I)=RR(1)-RG
153    clw(i) = rr(nk) - rg
154    ! jyg2
155    clw(i) = max(0.0, clw(i))
156    ! jyg1
157    ! CC        TVP(I)=TPK(I)*(1.+RG/EPS)
158    tvp(i) = tpk(i)*(1.+rg/eps-rr(nk))
159    ! jyg2
160
161    ! jyg1       Derivatives
162
163    dtpdt1(i) = cpd*s
164    dtpdq1(i) = alv*s
165
166    dtvpdt1(i) = dtpdt1(i)*(1.+rg/eps-rr(nk)+alv*rg/(rd*tpk(i)))
167    dtvpdq1(i) = dtpdq1(i)*(1.+rg/eps-rr(nk)+alv*rg/(rd*tpk(i))) - tpk(i)
168
169    ! jyg2
170
171  END DO
172
173  ice_conv = .FALSE.
174
175  IF (ice_conv) THEN
176
177    ! JAM
178    ! RAJOUT DE LA PROCEDURE ICEFRAC
179
180    ! sb        CALL ICEFRAC(T,CLW,CLW_NEW,QI,ND,NL)
181
182    DO i = icb1, nl
183      IF (t(i)<263.15) THEN
184        tg = tpk(i)
185        tc = tpk(i) - 273.15
186        denom = 243.5 + tc
187        es = 6.112*exp(17.67*tc/denom)
188        alv = alv0 - cpvmcl*(t(i)-273.15)
189        alf = alf0 + clmci*(t(i)-273.15)
190
191        DO j = 1, 4
192          esi = exp(23.33086-(6111.72784/tpk(i))+0.15215*log(tpk(i)))
193          qsat_new = eps*esi/(p(i)-esi*(1.-eps))
194          ! CC        SNEW=
195          ! CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*QSAT_NEW/(RV*TPK(I)*TPK(I))
196          snew = cpd*(1.-rr(nk)) + cl*rr(nk) + alv*alv*qsat_new/(rv*tpk(i)* &
197            tpk(i))
198
199          snew = 1./snew
200          tpk(i) = tg + (alf*qi(i)+alv*rg*(1.-(esi/es)))*snew
201          ! @$$        PRINT*,'################################'
202          ! @$$        PRINT*,TPK(I)
203          ! @$$        PRINT*,(ALF*QI(I)+ALV*RG*(1.-(ESI/ES)))*SNEW
204        END DO
205        ! CC        CLW(I)=RR(1)-QSAT_NEW
206        clw(i) = rr(nk) - qsat_new
207        clw(i) = max(0.0, clw(i))
208        ! jyg1
209        ! CC        TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS)
210        tvp(i) = tpk(i)*(1.+qsat_new/eps-rr(nk))
211        ! jyg2
212      ELSE
213        CONTINUE
214      END IF
215
216    END DO
217
218  END IF
219
220
221  ! *****************************************************
222  ! * BK :  RAJOUT DE LA TEMPERATURE DES ASCENDANCES
223  ! *   NON DILUES AU  NIVEAU KLEV = ND
224  ! *   POSONS LE ENVIRON EGAL A CELUI DE KLEV-1
225  ! *******************************************************
226
227  tpk(nl+1) = tpk(nl)
228
229  ! ******************************************************
230
231  rg = gravity ! RG redevient la gravite de YOMCST (sb)
232
233
234  RETURN
235END SUBROUTINE tlift
236
237
238
239
240
241
242
243
Note: See TracBrowser for help on using the repository browser.