source: LMDZ6/branches/blowing_snow/libf/phylmd/tlift.F90 @ 5018

Last change on this file since 5018 was 2197, checked in by Ehouarn Millour, 10 years ago

Added 'implicit none' statements and proper variable definitions where they were missing.
EM

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