[92fab65] | 1 | ! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt |
---|
| 2 | ! SPDX-License-Identifier: GPL-3.0-or-later |
---|
[332fbbd] | 3 | |
---|
[e200b7a] | 4 | !************************************************************************** |
---|
| 5 | !**** SUBROUTINE CONVECT ***** |
---|
| 6 | !**** VERSION 4.3c ***** |
---|
| 7 | !**** 20 May, 2002 ***** |
---|
| 8 | !**** Kerry Emanuel ***** |
---|
| 9 | !************************************************************************** |
---|
| 10 | ! |
---|
| 11 | SUBROUTINE CONVECT & |
---|
| 12 | (ND, NL, DELT, IFLAG, & |
---|
| 13 | PRECIP, WD, TPRIME, QPRIME, CBMF ) |
---|
| 14 | ! |
---|
| 15 | !-cv ************************************************************************* |
---|
| 16 | !-cv C. Forster, November 2003 - May 2004: |
---|
| 17 | !-cv |
---|
| 18 | !-cv The subroutine has been downloaded from Kerry Emanuel's homepage, |
---|
| 19 | !-cv where further infos on the convection scheme can be found |
---|
| 20 | !-cv http://www-paoc.mit.edu/~emanuel/home.html |
---|
| 21 | !-cv |
---|
| 22 | !-cv The following changes have been made to integrate this subroutine |
---|
| 23 | !-cv into FLEXPART |
---|
| 24 | !-cv |
---|
| 25 | !-cv Putting most of the variables in a new common block |
---|
| 26 | !-cv renaming eps to eps0 because there is some eps already in includepar |
---|
| 27 | !-cv |
---|
| 28 | !-cv removing the arrays U,V,TRA and related arrays |
---|
| 29 | !-cv |
---|
| 30 | !-cv renaming the original arrays T,Q,QS,P,PH to |
---|
| 31 | !-cv TCONV,QCONV,QSCONV,PCONV_HPA,PHCONV_HPA |
---|
| 32 | !-cv |
---|
| 33 | !-cv Initialization of variables has been put into parameter statements |
---|
| 34 | !-cv instead of assignment of values at each call, in order to save |
---|
| 35 | !-cv computation time. |
---|
| 36 | !*************************************************************************** |
---|
| 37 | ! |
---|
| 38 | !----------------------------------------------------------------------------- |
---|
| 39 | ! *** On input: *** |
---|
| 40 | ! |
---|
| 41 | !T: Array of absolute temperature (K) of dimension ND, with first |
---|
| 42 | ! index corresponding to lowest model level. Note that this array |
---|
| 43 | ! will be altered by the subroutine if dry convective adjustment |
---|
| 44 | ! occurs and if IPBL is not equal to 0. |
---|
| 45 | ! |
---|
| 46 | !Q: Array of specific humidity (gm/gm) of dimension ND, with first |
---|
| 47 | ! index corresponding to lowest model level. Must be defined |
---|
| 48 | ! at same grid levels as T. Note that this array will be altered |
---|
| 49 | ! if dry convective adjustment occurs and if IPBL is not equal to 0. |
---|
| 50 | ! |
---|
| 51 | !QS: Array of saturation specific humidity of dimension ND, with first |
---|
| 52 | ! index corresponding to lowest model level. Must be defined |
---|
| 53 | ! at same grid levels as T. Note that this array will be altered |
---|
| 54 | ! if dry convective adjustment occurs and if IPBL is not equal to 0. |
---|
| 55 | ! |
---|
| 56 | !U: Array of zonal wind velocity (m/s) of dimension ND, witth first |
---|
| 57 | ! index corresponding with the lowest model level. Defined at |
---|
| 58 | ! same levels as T. Note that this array will be altered if |
---|
| 59 | ! dry convective adjustment occurs and if IPBL is not equal to 0. |
---|
| 60 | ! |
---|
| 61 | !V: Same as U but for meridional velocity. |
---|
| 62 | ! |
---|
| 63 | !TRA: Array of passive tracer mixing ratio, of dimensions (ND,NTRA), |
---|
| 64 | ! where NTRA is the number of different tracers. If no |
---|
| 65 | ! convective tracer transport is needed, define a dummy |
---|
| 66 | ! input array of dimension (ND,1). Tracers are defined at |
---|
| 67 | ! same vertical levels as T. Note that this array will be altered |
---|
| 68 | ! if dry convective adjustment occurs and if IPBL is not equal to 0. |
---|
| 69 | ! |
---|
| 70 | !P: Array of pressure (mb) of dimension ND, with first |
---|
| 71 | ! index corresponding to lowest model level. Must be defined |
---|
| 72 | ! at same grid levels as T. |
---|
| 73 | ! |
---|
| 74 | !PH: Array of pressure (mb) of dimension ND+1, with first index |
---|
| 75 | ! corresponding to lowest level. These pressures are defined at |
---|
| 76 | ! levels intermediate between those of P, T, Q and QS. The first |
---|
| 77 | ! value of PH should be greater than (i.e. at a lower level than) |
---|
| 78 | ! the first value of the array P. |
---|
| 79 | ! |
---|
| 80 | !ND: The dimension of the arrays T,Q,QS,P,PH,FT and FQ |
---|
| 81 | ! |
---|
| 82 | !NL: The maximum number of levels to which convection can |
---|
| 83 | ! penetrate, plus 1. |
---|
| 84 | ! NL MUST be less than or equal to ND-1. |
---|
| 85 | ! |
---|
| 86 | !NTRA:The number of different tracers. If no tracer transport |
---|
| 87 | ! is needed, set this equal to 1. (On most compilers, setting |
---|
| 88 | ! NTRA to 0 will bypass tracer calculation, saving some CPU.) |
---|
| 89 | ! |
---|
| 90 | !DELT: The model time step (sec) between calls to CONVECT |
---|
| 91 | ! |
---|
| 92 | !---------------------------------------------------------------------------- |
---|
| 93 | ! *** On Output: *** |
---|
| 94 | ! |
---|
| 95 | !IFLAG: An output integer whose value denotes the following: |
---|
| 96 | ! |
---|
| 97 | ! VALUE INTERPRETATION |
---|
| 98 | ! ----- -------------- |
---|
| 99 | ! 0 No moist convection; atmosphere is not |
---|
| 100 | ! unstable, or surface temperature is less |
---|
| 101 | ! than 250 K or surface specific humidity |
---|
| 102 | ! is non-positive. |
---|
| 103 | ! |
---|
| 104 | ! 1 Moist convection occurs. |
---|
| 105 | ! |
---|
| 106 | ! 2 No moist convection: lifted condensation |
---|
| 107 | ! level is above the 200 mb level. |
---|
| 108 | ! |
---|
| 109 | ! 3 No moist convection: cloud base is higher |
---|
| 110 | ! then the level NL-1. |
---|
| 111 | ! |
---|
| 112 | ! 4 Moist convection occurs, but a CFL condition |
---|
| 113 | ! on the subsidence warming is violated. This |
---|
| 114 | ! does not cause the scheme to terminate. |
---|
| 115 | ! |
---|
| 116 | !FT: Array of temperature tendency (K/s) of dimension ND, defined at same |
---|
| 117 | ! grid levels as T, Q, QS and P. |
---|
| 118 | ! |
---|
| 119 | !FQ: Array of specific humidity tendencies ((gm/gm)/s) of dimension ND, |
---|
| 120 | ! defined at same grid levels as T, Q, QS and P. |
---|
| 121 | ! |
---|
| 122 | !FU: Array of forcing of zonal velocity (m/s^2) of dimension ND, |
---|
| 123 | ! defined at same grid levels as T. |
---|
| 124 | ! |
---|
| 125 | !FV: Same as FU, but for forcing of meridional velocity. |
---|
| 126 | ! |
---|
| 127 | !FTRA: Array of forcing of tracer content, in tracer mixing ratio per |
---|
| 128 | ! second, defined at same levels as T. Dimensioned (ND,NTRA). |
---|
| 129 | ! |
---|
| 130 | !PRECIP: Scalar convective precipitation rate (mm/day). |
---|
| 131 | ! |
---|
| 132 | !WD: A convective downdraft velocity scale. For use in surface |
---|
| 133 | ! flux parameterizations. See convect.ps file for details. |
---|
| 134 | ! |
---|
| 135 | !TPRIME: A convective downdraft temperature perturbation scale (K). |
---|
| 136 | ! For use in surface flux parameterizations. See convect.ps |
---|
| 137 | ! file for details. |
---|
| 138 | ! |
---|
| 139 | !QPRIME: A convective downdraft specific humidity |
---|
| 140 | ! perturbation scale (gm/gm). |
---|
| 141 | ! For use in surface flux parameterizations. See convect.ps |
---|
| 142 | ! file for details. |
---|
| 143 | ! |
---|
| 144 | !CBMF: The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST |
---|
| 145 | ! BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT |
---|
| 146 | ! ITS NEXT CALL. That is, the value of CBMF must be "remembered" |
---|
| 147 | ! by the calling program between calls to CONVECT. |
---|
| 148 | ! |
---|
| 149 | !----------------------------------------------------------------------------- |
---|
| 150 | ! |
---|
| 151 | ! *** THE PARAMETER NA SHOULD IN GENERAL BE GREATER THAN *** |
---|
| 152 | ! *** OR EQUAL TO ND + 1 *** |
---|
| 153 | ! |
---|
| 154 | ! |
---|
| 155 | use par_mod |
---|
| 156 | use conv_mod |
---|
| 157 | |
---|
| 158 | implicit none |
---|
| 159 | ! |
---|
| 160 | !-cv====>Begin Module CONVECT File convect.f Undeclared variables |
---|
| 161 | ! |
---|
| 162 | !Argument variables |
---|
| 163 | ! |
---|
| 164 | integer :: iflag, nd, nl |
---|
| 165 | ! |
---|
| 166 | real :: cbmf, delt, precip, qprime, tprime, wd |
---|
| 167 | ! |
---|
| 168 | !Local variables |
---|
| 169 | ! |
---|
| 170 | integer :: i, icb, ihmin, inb, inb1, j, jtt, k |
---|
| 171 | integer :: nk |
---|
| 172 | ! |
---|
| 173 | real :: ad, afac, ahmax, ahmin, alt, altem |
---|
| 174 | real :: am, amp1, anum, asij, awat, b6, bf2, bsum, by |
---|
| 175 | real :: byp, c6, cape, capem, cbmfold, chi, coeff |
---|
| 176 | real :: cpinv, cwat, damps, dbo, dbosum |
---|
| 177 | real :: defrac, dei, delm, delp, delt0, delti, denom, dhdp |
---|
| 178 | real :: dpinv, dtma, dtmin, dtpbl, elacrit, ents |
---|
| 179 | real :: epmax, fac, fqold, frac, ftold |
---|
| 180 | real :: plcl, qp1, qsm, qstm, qti, rat |
---|
| 181 | real :: rdcp, revap, rh, scrit, sigt, sjmax |
---|
| 182 | real :: sjmin, smid, smin, stemp, tca |
---|
| 183 | real :: tvaplcl, tvpplcl, tvx, tvy, wdtrain |
---|
| 184 | |
---|
| 185 | !integer jc,jn |
---|
| 186 | !real alvnew,a2,ahm,alv,rm,sum,qnew,dphinv,tc,thbar,tnew,x |
---|
| 187 | |
---|
| 188 | real :: FUP(NA),FDOWN(NA) |
---|
| 189 | ! |
---|
| 190 | !-cv====>End Module CONVECT File convect.f |
---|
| 191 | |
---|
| 192 | INTEGER :: NENT(NA) |
---|
| 193 | REAL :: M(NA),MP(NA),MENT(NA,NA),QENT(NA,NA),ELIJ(NA,NA) |
---|
| 194 | REAL :: SIJ(NA,NA),TVP(NA),TV(NA),WATER(NA) |
---|
| 195 | REAL :: QP(NA),EP(NA),TH(NA),WT(NA),EVAP(NA),CLW(NA) |
---|
| 196 | REAL :: SIGP(NA),TP(NA),CPN(NA) |
---|
| 197 | REAL :: LV(NA),LVCP(NA),H(NA),HP(NA),GZ(NA),HM(NA) |
---|
| 198 | !REAL TOLD(NA) |
---|
| 199 | ! |
---|
| 200 | ! ----------------------------------------------------------------------- |
---|
| 201 | ! |
---|
| 202 | ! *** Specify Switches *** |
---|
| 203 | ! |
---|
| 204 | ! *** IPBL: Set to zero to bypass dry adiabatic adjustment *** |
---|
| 205 | ! *** Any other value results in dry adiabatic adjustment *** |
---|
| 206 | ! *** (Zero value recommended for use in models with *** |
---|
| 207 | ! *** boundary layer schemes) *** |
---|
| 208 | ! |
---|
| 209 | ! *** MINORIG: Lowest level from which convection may originate *** |
---|
| 210 | ! *** (Should be first model level at which T is defined *** |
---|
| 211 | ! *** for models using bulk PBL schemes; otherwise, it should *** |
---|
| 212 | ! *** be the first model level at which T is defined above *** |
---|
| 213 | ! *** the surface layer) *** |
---|
| 214 | ! |
---|
| 215 | INTEGER,PARAMETER :: IPBL=0 |
---|
| 216 | INTEGER,PARAMETER :: MINORIG=1 |
---|
| 217 | ! |
---|
| 218 | !------------------------------------------------------------------------------ |
---|
| 219 | ! |
---|
| 220 | ! *** SPECIFY PARAMETERS *** |
---|
| 221 | ! |
---|
| 222 | ! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) *** |
---|
| 223 | ! *** TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- *** |
---|
| 224 | ! *** CONVERSION THRESHOLD IS ASSUMED TO BE ZERO *** |
---|
| 225 | ! *** (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY *** |
---|
| 226 | ! *** BETWEEN 0 C AND TLCRIT) *** |
---|
| 227 | ! *** ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT *** |
---|
| 228 | ! *** FORMULATION *** |
---|
| 229 | ! *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** |
---|
| 230 | ! *** SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** |
---|
| 231 | ! *** OF CLOUD *** |
---|
| 232 | ! *** OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN *** |
---|
| 233 | ! *** OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW *** |
---|
| 234 | ! *** COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** |
---|
| 235 | ! *** OF RAIN *** |
---|
| 236 | ! *** COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** |
---|
| 237 | ! *** OF SNOW *** |
---|
| 238 | ! *** CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM *** |
---|
| 239 | ! *** TRANSPORT *** |
---|
| 240 | ! *** DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION *** |
---|
| 241 | ! *** A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC *** |
---|
| 242 | ! *** ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF *** |
---|
| 243 | ! *** APPROACH TO QUASI-EQUILIBRIUM *** |
---|
| 244 | ! *** (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) *** |
---|
| 245 | ! *** (DAMP MUST BE LESS THAN 1) *** |
---|
| 246 | ! |
---|
| 247 | REAL,PARAMETER :: ELCRIT=.0011 |
---|
| 248 | REAL,PARAMETER :: TLCRIT=-55.0 |
---|
| 249 | REAL,PARAMETER :: ENTP=1.5 |
---|
| 250 | REAL,PARAMETER :: SIGD=0.05 |
---|
| 251 | REAL,PARAMETER :: SIGS=0.12 |
---|
| 252 | REAL,PARAMETER :: OMTRAIN=50.0 |
---|
| 253 | REAL,PARAMETER :: OMTSNOW=5.5 |
---|
| 254 | REAL,PARAMETER :: COEFFR=1.0 |
---|
| 255 | REAL,PARAMETER :: COEFFS=0.8 |
---|
| 256 | REAL,PARAMETER :: CU=0.7 |
---|
| 257 | REAL,PARAMETER :: BETA=10.0 |
---|
| 258 | REAL,PARAMETER :: DTMAX=0.9 |
---|
| 259 | REAL,PARAMETER :: ALPHA=0.025 !original 0.2 |
---|
| 260 | REAL,PARAMETER :: DAMP=0.1 |
---|
| 261 | ! |
---|
| 262 | ! *** ASSIGN VALUES OF THERMODYNAMIC CONSTANTS, *** |
---|
| 263 | ! *** GRAVITY, AND LIQUID WATER DENSITY. *** |
---|
| 264 | ! *** THESE SHOULD BE CONSISTENT WITH *** |
---|
| 265 | ! *** THOSE USED IN CALLING PROGRAM *** |
---|
| 266 | ! *** NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT *** |
---|
| 267 | ! |
---|
| 268 | REAL,PARAMETER :: CPD=1005.7 |
---|
| 269 | REAL,PARAMETER :: CPV=1870.0 |
---|
| 270 | REAL,PARAMETER :: CL=2500.0 |
---|
| 271 | REAL,PARAMETER :: RV=461.5 |
---|
| 272 | REAL,PARAMETER :: RD=287.04 |
---|
| 273 | REAL,PARAMETER :: LV0=2.501E6 |
---|
| 274 | REAL,PARAMETER :: G=9.81 |
---|
| 275 | REAL,PARAMETER :: ROWL=1000.0 |
---|
| 276 | ! |
---|
| 277 | REAL,PARAMETER :: CPVMCL=CL-CPV |
---|
| 278 | REAL,PARAMETER :: EPS0=RD/RV |
---|
| 279 | REAL,PARAMETER :: EPSI=1./EPS0 |
---|
| 280 | REAL,PARAMETER :: GINV=1.0/G |
---|
| 281 | REAL,PARAMETER :: EPSILON=1.e-20 |
---|
| 282 | |
---|
| 283 | ! EPSILON IS A SMALL NUMBER USED TO EXCLUDE MASS FLUXES OF ZERO |
---|
| 284 | ! |
---|
| 285 | DELTI=1.0/DELT |
---|
| 286 | ! |
---|
| 287 | ! *** INITIALIZE OUTPUT ARRAYS AND PARAMETERS *** |
---|
| 288 | ! |
---|
| 289 | |
---|
| 290 | DO I=1,NL+1 |
---|
| 291 | FT(I)=0.0 |
---|
| 292 | FQ(I)=0.0 |
---|
| 293 | FDOWN(I)=0.0 |
---|
| 294 | SUB(I)=0.0 |
---|
| 295 | FUP(I)=0.0 |
---|
| 296 | M(I)=0.0 |
---|
| 297 | MP(I)=0.0 |
---|
| 298 | DO J=1,NL+1 |
---|
| 299 | FMASS(I,J)=0.0 |
---|
| 300 | MENT(I,J)=0.0 |
---|
| 301 | END DO |
---|
| 302 | END DO |
---|
| 303 | DO I=1,NL+1 |
---|
| 304 | RDCP=(RD*(1.-QCONV(I))+QCONV(I)*RV)/ & |
---|
| 305 | (CPD*(1.-QCONV(I))+QCONV(I)*CPV) |
---|
| 306 | TH(I)=TCONV(I)*(1000.0/PCONV_HPA(I))**RDCP |
---|
| 307 | END DO |
---|
| 308 | PRECIP=0.0 |
---|
| 309 | WD=0.0 |
---|
| 310 | TPRIME=0.0 |
---|
| 311 | QPRIME=0.0 |
---|
| 312 | IFLAG=0 |
---|
| 313 | ! |
---|
| 314 | ! IF(IPBL.NE.0)THEN |
---|
| 315 | ! |
---|
| 316 | !*** PERFORM DRY ADIABATIC ADJUSTMENT *** |
---|
| 317 | ! |
---|
| 318 | ! JC=0 |
---|
| 319 | ! DO 30 I=NL-1,1,-1 |
---|
| 320 | ! JN=0 |
---|
| 321 | ! SUM=TH(I)*(1.+QCONV(I)*EPSI-QCONV(I)) |
---|
| 322 | ! DO 10 J=I+1,NL |
---|
| 323 | ! SUM=SUM+TH(J)*(1.+QCONV(J)*EPSI-QCONV(J)) |
---|
| 324 | ! THBAR=SUM/REAL(J+1-I) |
---|
| 325 | ! IF((TH(J)*(1.+QCONV(J)*EPSI-QCONV(J))).LT.THBAR)JN=J |
---|
| 326 | ! 10 CONTINUE |
---|
| 327 | ! IF(I.EQ.1)JN=MAX(JN,2) |
---|
| 328 | ! IF(JN.EQ.0)GOTO 30 |
---|
| 329 | ! 12 CONTINUE |
---|
| 330 | ! AHM=0.0 |
---|
| 331 | ! RM=0.0 |
---|
| 332 | ! DO 15 J=I,JN |
---|
| 333 | ! AHM=AHM+(CPD*(1.-QCONV(J))+QCONV(J)*CPV)*TCONV(J)* |
---|
| 334 | ! + (PHCONV_HPA(J)-PHCONV_HPA(J+1)) |
---|
| 335 | ! RM=RM+QCONV(J)*(PHCONV_HPA(J)-PHCONV_HPA(J+1)) |
---|
| 336 | ! 15 CONTINUE |
---|
| 337 | ! DPHINV=1./(PHCONV_HPA(I)-PHCONV_HPA(JN+1)) |
---|
| 338 | ! RM=RM*DPHINV |
---|
| 339 | ! A2=0.0 |
---|
| 340 | ! DO 20 J=I,JN |
---|
| 341 | ! QCONV(J)=RM |
---|
| 342 | ! RDCP=(RD*(1.-QCONV(J))+QCONV(J)*RV)/ |
---|
| 343 | ! 1 (CPD*(1.-QCONV(J))+QCONV(J)*CPV) |
---|
| 344 | ! X=(0.001*PCONV_HPA(J))**RDCP |
---|
| 345 | ! TOLD(J)=TCONV(J) |
---|
| 346 | ! TCONV(J)=X |
---|
| 347 | ! A2=A2+(CPD*(1.-QCONV(J))+QCONV(J)*CPV)*X* |
---|
| 348 | ! 1 (PHCONV_HPA(J)-PHCONV_HPA(J+1)) |
---|
| 349 | ! 20 CONTINUE |
---|
| 350 | ! DO 25 J=I,JN |
---|
| 351 | ! TH(J)=AHM/A2 |
---|
| 352 | ! TCONV(J)=TCONV(J)*TH(J) |
---|
| 353 | ! TC=TOLD(J)-273.15 |
---|
| 354 | ! ALV=LV0-CPVMCL*TC |
---|
| 355 | ! QSCONV(J)=QSCONV(J)+QSCONV(J)*(1.+QSCONV(J)*(EPSI-1.))*ALV* |
---|
| 356 | ! 1 (TCONV(J)- TOLD(J))/(RV*TOLD(J)*TOLD(J)) |
---|
| 357 | ! if (qslev(j) .lt. 0.) then |
---|
| 358 | ! write(*,*) 'qslev.lt.0 ',j,qslev |
---|
| 359 | ! endif |
---|
| 360 | ! 25 CONTINUE |
---|
| 361 | ! IF((TH(JN+1)*(1.+QCONV(JN+1)*EPSI-QCONV(JN+1))).LT. |
---|
| 362 | ! 1 (TH(JN)*(1.+QCONV(JN)*EPSI-QCONV(JN))))THEN |
---|
| 363 | ! JN=JN+1 |
---|
| 364 | ! GOTO 12 |
---|
| 365 | ! END IF |
---|
| 366 | ! IF(I.EQ.1)JC=JN |
---|
| 367 | ! 30 CONTINUE |
---|
| 368 | ! |
---|
| 369 | ! *** Remove any supersaturation that results from adjustment *** |
---|
| 370 | ! |
---|
| 371 | !IF(JC.GT.1)THEN |
---|
| 372 | ! DO 38 J=1,JC |
---|
| 373 | ! IF(QSCONV(J).LT.QCONV(J))THEN |
---|
| 374 | ! ALV=LV0-CPVMCL*(TCONV(J)-273.15) |
---|
| 375 | ! TNEW=TCONV(J)+ALV*(QCONV(J)-QSCONV(J))/(CPD*(1.-QCONV(J))+ |
---|
| 376 | ! 1 CL*QCONV(J)+QSCONV(J)*(CPV-CL+ALV*ALV/(RV*TCONV(J)*TCONV(J)))) |
---|
| 377 | ! ALVNEW=LV0-CPVMCL*(TNEW-273.15) |
---|
| 378 | ! QNEW=(ALV*QCONV(J)-(TNEW-TCONV(J))*(CPD*(1.-QCONV(J)) |
---|
| 379 | ! 1 +CL*QCONV(J)))/ALVNEW |
---|
| 380 | ! PRECIP=PRECIP+24.*3600.*1.0E5*(PHCONV_HPA(J)-PHCONV_HPA(J+1))* |
---|
| 381 | ! 1 (QCONV(J)-QNEW)/(G*DELT*ROWL) |
---|
| 382 | ! TCONV(J)=TNEW |
---|
| 383 | ! QCONV(J)=QNEW |
---|
| 384 | ! QSCONV(J)=QNEW |
---|
| 385 | ! END IF |
---|
| 386 | ! 38 CONTINUE |
---|
| 387 | !END IF |
---|
| 388 | ! |
---|
| 389 | !END IF |
---|
| 390 | ! |
---|
| 391 | ! *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY |
---|
| 392 | ! |
---|
| 393 | GZ(1)=0.0 |
---|
| 394 | CPN(1)=CPD*(1.-QCONV(1))+QCONV(1)*CPV |
---|
| 395 | H(1)=TCONV(1)*CPN(1) |
---|
| 396 | LV(1)=LV0-CPVMCL*(TCONV(1)-273.15) |
---|
| 397 | HM(1)=LV(1)*QCONV(1) |
---|
| 398 | TV(1)=TCONV(1)*(1.+QCONV(1)*EPSI-QCONV(1)) |
---|
| 399 | AHMIN=1.0E12 |
---|
| 400 | IHMIN=NL |
---|
| 401 | DO I=2,NL+1 |
---|
| 402 | TVX=TCONV(I)*(1.+QCONV(I)*EPSI-QCONV(I)) |
---|
| 403 | TVY=TCONV(I-1)*(1.+QCONV(I-1)*EPSI-QCONV(I-1)) |
---|
| 404 | GZ(I)=GZ(I-1)+0.5*RD*(TVX+TVY)*(PCONV_HPA(I-1)-PCONV_HPA(I))/ & |
---|
| 405 | PHCONV_HPA(I) |
---|
| 406 | CPN(I)=CPD*(1.-QCONV(I))+CPV*QCONV(I) |
---|
| 407 | H(I)=TCONV(I)*CPN(I)+GZ(I) |
---|
| 408 | LV(I)=LV0-CPVMCL*(TCONV(I)-273.15) |
---|
| 409 | HM(I)=(CPD*(1.-QCONV(I))+CL*QCONV(I))*(TCONV(I)-TCONV(1))+ & |
---|
| 410 | LV(I)*QCONV(I)+GZ(I) |
---|
| 411 | TV(I)=TCONV(I)*(1.+QCONV(I)*EPSI-QCONV(I)) |
---|
| 412 | ! |
---|
| 413 | ! *** Find level of minimum moist static energy *** |
---|
| 414 | ! |
---|
| 415 | IF(I.GE.MINORIG.AND.HM(I).LT.AHMIN.AND.HM(I).LT.HM(I-1))THEN |
---|
| 416 | AHMIN=HM(I) |
---|
| 417 | IHMIN=I |
---|
| 418 | END IF |
---|
| 419 | END DO |
---|
| 420 | IHMIN=MIN(IHMIN, NL-1) |
---|
| 421 | ! |
---|
| 422 | ! *** Find that model level below the level of minimum moist *** |
---|
| 423 | ! *** static energy that has the maximum value of moist static energy *** |
---|
| 424 | ! |
---|
| 425 | AHMAX=0.0 |
---|
| 426 | ! *** bug fixed: need to assign an initial value to NK |
---|
| 427 | ! HSO, 05.08.2009 |
---|
| 428 | NK=MINORIG |
---|
| 429 | DO I=MINORIG,IHMIN |
---|
| 430 | IF(HM(I).GT.AHMAX)THEN |
---|
| 431 | NK=I |
---|
| 432 | AHMAX=HM(I) |
---|
| 433 | END IF |
---|
| 434 | END DO |
---|
| 435 | ! |
---|
| 436 | ! *** CHECK WHETHER PARCEL LEVEL TEMPERATURE AND SPECIFIC HUMIDITY *** |
---|
| 437 | ! *** ARE REASONABLE *** |
---|
| 438 | ! *** Skip convection if HM increases monotonically upward *** |
---|
| 439 | ! |
---|
| 440 | IF(TCONV(NK).LT.250.0.OR.QCONV(NK).LE.0.0.OR.IHMIN.EQ.(NL-1)) & |
---|
| 441 | THEN |
---|
| 442 | IFLAG=0 |
---|
| 443 | CBMF=0.0 |
---|
| 444 | RETURN |
---|
| 445 | END IF |
---|
| 446 | ! |
---|
| 447 | ! *** CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT PARCEL ORIGIN LEVEL *** |
---|
| 448 | ! *** (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980) *** |
---|
| 449 | ! |
---|
| 450 | RH=QCONV(NK)/QSCONV(NK) |
---|
| 451 | CHI=TCONV(NK)/(1669.0-122.0*RH-TCONV(NK)) |
---|
| 452 | PLCL=PCONV_HPA(NK)*(RH**CHI) |
---|
| 453 | IF(PLCL.LT.200.0.OR.PLCL.GE.2000.0)THEN |
---|
| 454 | IFLAG=2 |
---|
| 455 | CBMF=0.0 |
---|
| 456 | RETURN |
---|
| 457 | END IF |
---|
| 458 | ! |
---|
| 459 | ! *** CALCULATE FIRST LEVEL ABOVE LCL (=ICB) *** |
---|
| 460 | ! |
---|
| 461 | ICB=NL-1 |
---|
| 462 | DO I=NK+1,NL |
---|
| 463 | IF(PCONV_HPA(I).LT.PLCL)THEN |
---|
| 464 | ICB=MIN(ICB,I) |
---|
| 465 | END IF |
---|
| 466 | END DO |
---|
| 467 | IF(ICB.GE.(NL-1))THEN |
---|
| 468 | IFLAG=3 |
---|
| 469 | CBMF=0.0 |
---|
| 470 | RETURN |
---|
| 471 | END IF |
---|
| 472 | ! |
---|
| 473 | ! *** FIND TEMPERATURE UP THROUGH ICB AND TEST FOR INSTABILITY *** |
---|
| 474 | ! |
---|
| 475 | ! *** SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL *** |
---|
| 476 | ! *** TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC *** |
---|
| 477 | ! *** LIQUID WATER CONTENT *** |
---|
| 478 | ! |
---|
| 479 | CALL TLIFT(GZ,ICB,NK,TVP,TP,CLW,ND,NL,1) |
---|
| 480 | DO I=NK,ICB |
---|
| 481 | TVP(I)=TVP(I)-TP(I)*QCONV(NK) |
---|
| 482 | END DO |
---|
| 483 | ! |
---|
| 484 | ! *** If there was no convection at last time step and parcel *** |
---|
| 485 | ! *** is stable at ICB then skip rest of calculation *** |
---|
| 486 | ! |
---|
| 487 | IF(CBMF.EQ.0.0.AND.TVP(ICB).LE.(TV(ICB)-DTMAX))THEN |
---|
| 488 | IFLAG=0 |
---|
| 489 | RETURN |
---|
| 490 | END IF |
---|
| 491 | ! |
---|
| 492 | ! *** IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY *** |
---|
| 493 | ! |
---|
| 494 | IF(IFLAG.NE.4)IFLAG=1 |
---|
| 495 | ! |
---|
| 496 | ! *** FIND THE REST OF THE LIFTED PARCEL TEMPERATURES *** |
---|
| 497 | ! |
---|
| 498 | CALL TLIFT(GZ,ICB,NK,TVP,TP,CLW,ND,NL,2) |
---|
| 499 | ! |
---|
| 500 | ! *** SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF *** |
---|
| 501 | ! *** PRECIPITATION FALLING OUTSIDE OF CLOUD *** |
---|
| 502 | ! *** THESE MAY BE FUNCTIONS OF TP(I), PCONV_HPA(I) AND CLW(I) *** |
---|
| 503 | ! |
---|
| 504 | DO I=1,NK |
---|
| 505 | EP(I)=0.0 |
---|
| 506 | SIGP(I)=SIGS |
---|
| 507 | END DO |
---|
| 508 | DO I=NK+1,NL |
---|
| 509 | TCA=TP(I)-273.15 |
---|
| 510 | IF(TCA.GE.0.0)THEN |
---|
| 511 | ELACRIT=ELCRIT |
---|
| 512 | ELSE |
---|
| 513 | ELACRIT=ELCRIT*(1.0-TCA/TLCRIT) |
---|
| 514 | END IF |
---|
| 515 | ELACRIT=MAX(ELACRIT,0.0) |
---|
| 516 | EPMAX=0.999 |
---|
| 517 | EP(I)=EPMAX*(1.0-ELACRIT/MAX(CLW(I),1.0E-8)) |
---|
| 518 | EP(I)=MAX(EP(I),0.0) |
---|
| 519 | EP(I)=MIN(EP(I),EPMAX) |
---|
| 520 | SIGP(I)=SIGS |
---|
| 521 | END DO |
---|
| 522 | ! |
---|
| 523 | ! *** CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL *** |
---|
| 524 | ! *** VIRTUAL TEMPERATURE *** |
---|
| 525 | ! |
---|
| 526 | DO I=ICB+1,NL |
---|
| 527 | TVP(I)=TVP(I)-TP(I)*QCONV(NK) |
---|
| 528 | END DO |
---|
| 529 | TVP(NL+1)=TVP(NL)-(GZ(NL+1)-GZ(NL))/CPD |
---|
| 530 | ! |
---|
| 531 | ! *** NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS *** |
---|
| 532 | ! |
---|
| 533 | DO I=1,NL+1 |
---|
| 534 | HP(I)=H(I) |
---|
| 535 | NENT(I)=0 |
---|
| 536 | WATER(I)=0.0 |
---|
| 537 | EVAP(I)=0.0 |
---|
| 538 | WT(I)=OMTSNOW |
---|
| 539 | LVCP(I)=LV(I)/CPN(I) |
---|
| 540 | DO J=1,NL+1 |
---|
| 541 | QENT(I,J)=QCONV(J) |
---|
| 542 | ELIJ(I,J)=0.0 |
---|
| 543 | SIJ(I,J)=0.0 |
---|
| 544 | END DO |
---|
| 545 | END DO |
---|
| 546 | QP(1)=QCONV(1) |
---|
| 547 | DO I=2,NL+1 |
---|
| 548 | QP(I)=QCONV(I-1) |
---|
| 549 | END DO |
---|
| 550 | ! |
---|
| 551 | ! *** FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S *** |
---|
| 552 | ! *** HIGHEST LEVEL OF NEUTRAL BUOYANCY *** |
---|
| 553 | ! *** AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB) *** |
---|
| 554 | ! |
---|
| 555 | CAPE=0.0 |
---|
| 556 | CAPEM=0.0 |
---|
| 557 | INB=ICB+1 |
---|
| 558 | INB1=INB |
---|
| 559 | BYP=0.0 |
---|
| 560 | DO I=ICB+1,NL-1 |
---|
| 561 | BY=(TVP(I)-TV(I))*(PHCONV_HPA(I)-PHCONV_HPA(I+1))/PCONV_HPA(I) |
---|
| 562 | CAPE=CAPE+BY |
---|
| 563 | IF(BY.GE.0.0)INB1=I+1 |
---|
| 564 | IF(CAPE.GT.0.0)THEN |
---|
| 565 | INB=I+1 |
---|
| 566 | BYP=(TVP(I+1)-TV(I+1))*(PHCONV_HPA(I+1)-PHCONV_HPA(I+2))/ & |
---|
| 567 | PCONV_HPA(I+1) |
---|
| 568 | CAPEM=CAPE |
---|
| 569 | END IF |
---|
| 570 | END DO |
---|
| 571 | INB=MAX(INB,INB1) |
---|
| 572 | CAPE=CAPEM+BYP |
---|
| 573 | DEFRAC=CAPEM-CAPE |
---|
| 574 | DEFRAC=MAX(DEFRAC,0.001) |
---|
| 575 | FRAC=-CAPE/DEFRAC |
---|
| 576 | FRAC=MIN(FRAC,1.0) |
---|
| 577 | FRAC=MAX(FRAC,0.0) |
---|
| 578 | ! |
---|
| 579 | ! *** CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL *** |
---|
| 580 | ! |
---|
| 581 | DO I=ICB,INB |
---|
| 582 | HP(I)=H(NK)+(LV(I)+(CPD-CPV)*TCONV(I))*EP(I)*CLW(I) |
---|
| 583 | END DO |
---|
| 584 | ! |
---|
| 585 | ! *** CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I), *** |
---|
| 586 | ! *** AT EACH MODEL LEVEL *** |
---|
| 587 | ! |
---|
| 588 | DBOSUM=0.0 |
---|
| 589 | ! |
---|
| 590 | ! *** INTERPOLATE DIFFERENCE BETWEEN LIFTED PARCEL AND *** |
---|
| 591 | ! *** ENVIRONMENTAL TEMPERATURES TO LIFTED CONDENSATION LEVEL *** |
---|
| 592 | ! |
---|
| 593 | TVPPLCL=TVP(ICB-1)-RD*TVP(ICB-1)*(PCONV_HPA(ICB-1)-PLCL)/ & |
---|
| 594 | (CPN(ICB-1)*PCONV_HPA(ICB-1)) |
---|
| 595 | TVAPLCL=TV(ICB)+(TVP(ICB)-TVP(ICB+1))*(PLCL-PCONV_HPA(ICB))/ & |
---|
| 596 | (PCONV_HPA(ICB)-PCONV_HPA(ICB+1)) |
---|
| 597 | DTPBL=0.0 |
---|
| 598 | DO I=NK,ICB-1 |
---|
| 599 | DTPBL=DTPBL+(TVP(I)-TV(I))*(PHCONV_HPA(I)-PHCONV_HPA(I+1)) |
---|
| 600 | END DO |
---|
| 601 | DTPBL=DTPBL/(PHCONV_HPA(NK)-PHCONV_HPA(ICB)) |
---|
| 602 | DTMIN=TVPPLCL-TVAPLCL+DTMAX+DTPBL |
---|
| 603 | DTMA=DTMIN |
---|
| 604 | ! |
---|
| 605 | ! *** ADJUST CLOUD BASE MASS FLUX *** |
---|
| 606 | ! |
---|
| 607 | CBMFOLD=CBMF |
---|
| 608 | ! *** C. Forster: adjustment of CBMF is not allowed to depend on FLEXPART timestep |
---|
| 609 | DELT0=DELT/3. |
---|
| 610 | DAMPS=DAMP*DELT/DELT0 |
---|
| 611 | CBMF=(1.-DAMPS)*CBMF+0.1*ALPHA*DTMA |
---|
| 612 | CBMF=MAX(CBMF,0.0) |
---|
| 613 | ! |
---|
| 614 | ! *** If cloud base mass flux is zero, skip rest of calculation *** |
---|
| 615 | ! |
---|
| 616 | IF(CBMF.EQ.0.0.AND.CBMFOLD.EQ.0.0)THEN |
---|
| 617 | RETURN |
---|
| 618 | END IF |
---|
| 619 | |
---|
| 620 | ! |
---|
| 621 | ! *** CALCULATE RATES OF MIXING, M(I) *** |
---|
| 622 | ! |
---|
| 623 | M(ICB)=0.0 |
---|
| 624 | DO I=ICB+1,INB |
---|
| 625 | K=MIN(I,INB1) |
---|
| 626 | DBO=ABS(TV(K)-TVP(K))+ & |
---|
| 627 | ENTP*0.02*(PHCONV_HPA(K)-PHCONV_HPA(K+1)) |
---|
| 628 | DBOSUM=DBOSUM+DBO |
---|
| 629 | M(I)=CBMF*DBO |
---|
| 630 | END DO |
---|
| 631 | DO I=ICB+1,INB |
---|
| 632 | M(I)=M(I)/DBOSUM |
---|
| 633 | END DO |
---|
| 634 | ! |
---|
| 635 | ! *** CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING *** |
---|
| 636 | ! *** RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING *** |
---|
| 637 | ! *** FRACTION (SIJ) *** |
---|
| 638 | ! |
---|
| 639 | DO I=ICB+1,INB |
---|
| 640 | QTI=QCONV(NK)-EP(I)*CLW(I) |
---|
| 641 | DO J=ICB,INB |
---|
| 642 | BF2=1.+LV(J)*LV(J)*QSCONV(J)/(RV*TCONV(J)*TCONV(J)*CPD) |
---|
| 643 | ANUM=H(J)-HP(I)+(CPV-CPD)*TCONV(J)*(QTI-QCONV(J)) |
---|
| 644 | DENOM=H(I)-HP(I)+(CPD-CPV)*(QCONV(I)-QTI)*TCONV(J) |
---|
| 645 | DEI=DENOM |
---|
| 646 | IF(ABS(DEI).LT.0.01)DEI=0.01 |
---|
| 647 | SIJ(I,J)=ANUM/DEI |
---|
| 648 | SIJ(I,I)=1.0 |
---|
| 649 | ALTEM=SIJ(I,J)*QCONV(I)+(1.-SIJ(I,J))*QTI-QSCONV(J) |
---|
| 650 | ALTEM=ALTEM/BF2 |
---|
| 651 | CWAT=CLW(J)*(1.-EP(J)) |
---|
| 652 | STEMP=SIJ(I,J) |
---|
| 653 | IF((STEMP.LT.0.0.OR.STEMP.GT.1.0.OR. & |
---|
| 654 | ALTEM.GT.CWAT).AND.J.GT.I)THEN |
---|
| 655 | ANUM=ANUM-LV(J)*(QTI-QSCONV(J)-CWAT*BF2) |
---|
| 656 | DENOM=DENOM+LV(J)*(QCONV(I)-QTI) |
---|
| 657 | IF(ABS(DENOM).LT.0.01)DENOM=0.01 |
---|
| 658 | SIJ(I,J)=ANUM/DENOM |
---|
| 659 | ALTEM=SIJ(I,J)*QCONV(I)+(1.-SIJ(I,J))*QTI-QSCONV(J) |
---|
| 660 | ALTEM=ALTEM-(BF2-1.)*CWAT |
---|
| 661 | END IF |
---|
| 662 | IF(SIJ(I,J).GT.0.0.AND.SIJ(I,J).LT.0.9)THEN |
---|
| 663 | QENT(I,J)=SIJ(I,J)*QCONV(I)+(1.-SIJ(I,J))*QTI |
---|
| 664 | ELIJ(I,J)=ALTEM |
---|
| 665 | ELIJ(I,J)=MAX(0.0,ELIJ(I,J)) |
---|
| 666 | MENT(I,J)=M(I)/(1.-SIJ(I,J)) |
---|
| 667 | NENT(I)=NENT(I)+1 |
---|
| 668 | END IF |
---|
| 669 | SIJ(I,J)=MAX(0.0,SIJ(I,J)) |
---|
| 670 | SIJ(I,J)=MIN(1.0,SIJ(I,J)) |
---|
| 671 | END DO |
---|
| 672 | ! |
---|
| 673 | ! *** IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS *** |
---|
| 674 | ! *** AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES *** |
---|
| 675 | ! |
---|
| 676 | IF(NENT(I).EQ.0)THEN |
---|
| 677 | MENT(I,I)=M(I) |
---|
| 678 | QENT(I,I)=QCONV(NK)-EP(I)*CLW(I) |
---|
| 679 | ELIJ(I,I)=CLW(I) |
---|
| 680 | SIJ(I,I)=1.0 |
---|
| 681 | END IF |
---|
| 682 | END DO |
---|
| 683 | SIJ(INB,INB)=1.0 |
---|
| 684 | ! |
---|
| 685 | ! *** NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL *** |
---|
| 686 | ! *** PROBABILITIES OF MIXING *** |
---|
| 687 | ! |
---|
| 688 | DO I=ICB+1,INB |
---|
| 689 | IF(NENT(I).NE.0)THEN |
---|
| 690 | QP1=QCONV(NK)-EP(I)*CLW(I) |
---|
| 691 | ANUM=H(I)-HP(I)-LV(I)*(QP1-QSCONV(I)) |
---|
| 692 | DENOM=H(I)-HP(I)+LV(I)*(QCONV(I)-QP1) |
---|
| 693 | IF(ABS(DENOM).LT.0.01)DENOM=0.01 |
---|
| 694 | SCRIT=ANUM/DENOM |
---|
| 695 | ALT=QP1-QSCONV(I)+SCRIT*(QCONV(I)-QP1) |
---|
| 696 | IF(ALT.LT.0.0)SCRIT=1.0 |
---|
| 697 | SCRIT=MAX(SCRIT,0.0) |
---|
| 698 | ASIJ=0.0 |
---|
| 699 | SMIN=1.0 |
---|
| 700 | DO J=ICB,INB |
---|
| 701 | IF(SIJ(I,J).GT.0.0.AND.SIJ(I,J).LT.0.9)THEN |
---|
| 702 | IF(J.GT.I)THEN |
---|
| 703 | SMID=MIN(SIJ(I,J),SCRIT) |
---|
| 704 | SJMAX=SMID |
---|
| 705 | SJMIN=SMID |
---|
| 706 | IF(SMID.LT.SMIN.AND.SIJ(I,J+1).LT.SMID)THEN |
---|
| 707 | SMIN=SMID |
---|
| 708 | SJMAX=MIN(SIJ(I,J+1),SIJ(I,J),SCRIT) |
---|
| 709 | SJMIN=MAX(SIJ(I,J-1),SIJ(I,J)) |
---|
| 710 | SJMIN=MIN(SJMIN,SCRIT) |
---|
| 711 | END IF |
---|
| 712 | ELSE |
---|
| 713 | SJMAX=MAX(SIJ(I,J+1),SCRIT) |
---|
| 714 | SMID=MAX(SIJ(I,J),SCRIT) |
---|
| 715 | SJMIN=0.0 |
---|
| 716 | IF(J.GT.1)SJMIN=SIJ(I,J-1) |
---|
| 717 | SJMIN=MAX(SJMIN,SCRIT) |
---|
| 718 | END IF |
---|
| 719 | DELP=ABS(SJMAX-SMID) |
---|
| 720 | DELM=ABS(SJMIN-SMID) |
---|
| 721 | ASIJ=ASIJ+(DELP+DELM)*(PHCONV_HPA(J)-PHCONV_HPA(J+1)) |
---|
| 722 | MENT(I,J)=MENT(I,J)*(DELP+DELM)* & |
---|
| 723 | (PHCONV_HPA(J)-PHCONV_HPA(J+1)) |
---|
| 724 | END IF |
---|
| 725 | END DO |
---|
| 726 | ASIJ=MAX(1.0E-21,ASIJ) |
---|
| 727 | ASIJ=1.0/ASIJ |
---|
| 728 | DO J=ICB,INB |
---|
| 729 | MENT(I,J)=MENT(I,J)*ASIJ |
---|
| 730 | END DO |
---|
| 731 | BSUM=0.0 |
---|
| 732 | DO J=ICB,INB |
---|
| 733 | BSUM=BSUM+MENT(I,J) |
---|
| 734 | END DO |
---|
| 735 | IF(BSUM.LT.1.0E-18)THEN |
---|
| 736 | NENT(I)=0 |
---|
| 737 | MENT(I,I)=M(I) |
---|
| 738 | QENT(I,I)=QCONV(NK)-EP(I)*CLW(I) |
---|
| 739 | ELIJ(I,I)=CLW(I) |
---|
| 740 | SIJ(I,I)=1.0 |
---|
| 741 | END IF |
---|
| 742 | END IF |
---|
| 743 | END DO |
---|
| 744 | ! |
---|
| 745 | ! *** CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING *** |
---|
| 746 | ! *** DOWNDRAFT CALCULATION *** |
---|
| 747 | ! |
---|
| 748 | IF(EP(INB).LT.0.0001)GOTO 405 |
---|
| 749 | ! |
---|
| 750 | ! *** INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER *** |
---|
| 751 | ! *** AND CONDENSED WATER FLUX *** |
---|
| 752 | ! |
---|
| 753 | JTT=2 |
---|
| 754 | ! |
---|
| 755 | ! *** BEGIN DOWNDRAFT LOOP *** |
---|
| 756 | ! |
---|
| 757 | DO I=INB,1,-1 |
---|
| 758 | ! |
---|
| 759 | ! *** CALCULATE DETRAINED PRECIPITATION *** |
---|
| 760 | ! |
---|
| 761 | WDTRAIN=G*EP(I)*M(I)*CLW(I) |
---|
| 762 | IF(I.GT.1)THEN |
---|
| 763 | DO J=1,I-1 |
---|
| 764 | AWAT=ELIJ(J,I)-(1.-EP(I))*CLW(I) |
---|
| 765 | AWAT=MAX(0.0,AWAT) |
---|
| 766 | WDTRAIN=WDTRAIN+G*AWAT*MENT(J,I) |
---|
| 767 | END DO |
---|
| 768 | END IF |
---|
| 769 | ! |
---|
| 770 | ! *** FIND RAIN WATER AND EVAPORATION USING PROVISIONAL *** |
---|
| 771 | ! *** ESTIMATES OF QP(I)AND QP(I-1) *** |
---|
| 772 | ! |
---|
| 773 | ! |
---|
| 774 | ! *** Value of terminal velocity and coefficient of evaporation for snow *** |
---|
| 775 | ! |
---|
| 776 | COEFF=COEFFS |
---|
| 777 | WT(I)=OMTSNOW |
---|
| 778 | ! |
---|
| 779 | ! *** Value of terminal velocity and coefficient of evaporation for rain *** |
---|
| 780 | ! |
---|
| 781 | IF(TCONV(I).GT.273.0)THEN |
---|
| 782 | COEFF=COEFFR |
---|
| 783 | WT(I)=OMTRAIN |
---|
| 784 | END IF |
---|
| 785 | QSM=0.5*(QCONV(I)+QP(I+1)) |
---|
| 786 | AFAC=COEFF*PHCONV_HPA(I)*(QSCONV(I)-QSM)/ & |
---|
| 787 | (1.0E4+2.0E3*PHCONV_HPA(I)*QSCONV(I)) |
---|
| 788 | AFAC=MAX(AFAC,0.0) |
---|
| 789 | SIGT=SIGP(I) |
---|
| 790 | SIGT=MAX(0.0,SIGT) |
---|
| 791 | SIGT=MIN(1.0,SIGT) |
---|
| 792 | B6=100.*(PHCONV_HPA(I)-PHCONV_HPA(I+1))*SIGT*AFAC/WT(I) |
---|
| 793 | C6=(WATER(I+1)*WT(I+1)+WDTRAIN/SIGD)/WT(I) |
---|
| 794 | REVAP=0.5*(-B6+SQRT(B6*B6+4.*C6)) |
---|
| 795 | EVAP(I)=SIGT*AFAC*REVAP |
---|
| 796 | WATER(I)=REVAP*REVAP |
---|
| 797 | ! |
---|
| 798 | ! *** CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER *** |
---|
| 799 | ! *** HYDROSTATIC APPROXIMATION *** |
---|
| 800 | ! |
---|
| 801 | IF(I.EQ.1)GOTO 360 |
---|
| 802 | DHDP=(H(I)-H(I-1))/(PCONV_HPA(I-1)-PCONV_HPA(I)) |
---|
| 803 | DHDP=MAX(DHDP,10.0) |
---|
| 804 | MP(I)=100.*GINV*LV(I)*SIGD*EVAP(I)/DHDP |
---|
| 805 | MP(I)=MAX(MP(I),0.0) |
---|
| 806 | ! |
---|
| 807 | ! *** ADD SMALL AMOUNT OF INERTIA TO DOWNDRAFT *** |
---|
| 808 | ! |
---|
| 809 | FAC=20.0/(PHCONV_HPA(I-1)-PHCONV_HPA(I)) |
---|
| 810 | MP(I)=(FAC*MP(I+1)+MP(I))/(1.+FAC) |
---|
| 811 | ! |
---|
| 812 | ! *** FORCE MP TO DECREASE LINEARLY TO ZERO *** |
---|
| 813 | ! *** BETWEEN ABOUT 950 MB AND THE SURFACE *** |
---|
| 814 | ! |
---|
| 815 | IF(PCONV_HPA(I).GT.(0.949*PCONV_HPA(1)))THEN |
---|
| 816 | JTT=MAX(JTT,I) |
---|
| 817 | MP(I)=MP(JTT)*(PCONV_HPA(1)-PCONV_HPA(I))/(PCONV_HPA(1)- & |
---|
| 818 | PCONV_HPA(JTT)) |
---|
| 819 | END IF |
---|
| 820 | 360 CONTINUE |
---|
| 821 | ! |
---|
| 822 | ! *** FIND MIXING RATIO OF PRECIPITATING DOWNDRAFT *** |
---|
| 823 | ! |
---|
| 824 | IF(I.EQ.INB)GOTO 400 |
---|
| 825 | IF(I.EQ.1)THEN |
---|
| 826 | QSTM=QSCONV(1) |
---|
| 827 | ELSE |
---|
| 828 | QSTM=QSCONV(I-1) |
---|
| 829 | END IF |
---|
| 830 | IF(MP(I).GT.MP(I+1))THEN |
---|
| 831 | RAT=MP(I+1)/MP(I) |
---|
| 832 | QP(I)=QP(I+1)*RAT+QCONV(I)*(1.0-RAT)+100.*GINV* & |
---|
| 833 | SIGD*(PHCONV_HPA(I)-PHCONV_HPA(I+1))*(EVAP(I)/MP(I)) |
---|
| 834 | ELSE |
---|
| 835 | IF(MP(I+1).GT.0.0)THEN |
---|
| 836 | QP(I)=(GZ(I+1)-GZ(I)+QP(I+1)*(LV(I+1)+TCONV(I+1)*(CL-CPD))+ & |
---|
| 837 | CPD*(TCONV(I+1)-TCONV(I)))/(LV(I)+TCONV(I)*(CL-CPD)) |
---|
| 838 | END IF |
---|
| 839 | END IF |
---|
| 840 | QP(I)=MIN(QP(I),QSTM) |
---|
| 841 | QP(I)=MAX(QP(I),0.0) |
---|
| 842 | 400 CONTINUE |
---|
| 843 | END DO |
---|
| 844 | ! |
---|
| 845 | ! *** CALCULATE SURFACE PRECIPITATION IN MM/DAY *** |
---|
| 846 | ! |
---|
| 847 | PRECIP=PRECIP+WT(1)*SIGD*WATER(1)*3600.*24000./(ROWL*G) |
---|
| 848 | ! |
---|
| 849 | 405 CONTINUE |
---|
| 850 | ! |
---|
| 851 | ! *** CALCULATE DOWNDRAFT VELOCITY SCALE AND SURFACE TEMPERATURE AND *** |
---|
| 852 | ! *** WATER VAPOR FLUCTUATIONS *** |
---|
| 853 | ! |
---|
| 854 | WD=BETA*ABS(MP(ICB))*0.01*RD*TCONV(ICB)/(SIGD*PCONV_HPA(ICB)) |
---|
| 855 | QPRIME=0.5*(QP(1)-QCONV(1)) |
---|
| 856 | TPRIME=LV0*QPRIME/CPD |
---|
| 857 | ! |
---|
| 858 | ! *** CALCULATE TENDENCIES OF LOWEST LEVEL POTENTIAL TEMPERATURE *** |
---|
| 859 | ! *** AND MIXING RATIO *** |
---|
| 860 | ! |
---|
| 861 | DPINV=0.01/(PHCONV_HPA(1)-PHCONV_HPA(2)) |
---|
| 862 | AM=0.0 |
---|
| 863 | IF(NK.EQ.1)THEN |
---|
| 864 | DO K=2,INB |
---|
| 865 | AM=AM+M(K) |
---|
| 866 | END DO |
---|
| 867 | END IF |
---|
| 868 | ! save saturated upward mass flux for first level |
---|
| 869 | FUP(1)=AM |
---|
| 870 | IF((2.*G*DPINV*AM).GE.DELTI)IFLAG=4 |
---|
| 871 | FT(1)=FT(1)+G*DPINV*AM*(TCONV(2)-TCONV(1)+(GZ(2)-GZ(1))/CPN(1)) |
---|
| 872 | FT(1)=FT(1)-LVCP(1)*SIGD*EVAP(1) |
---|
| 873 | FT(1)=FT(1)+SIGD*WT(2)*(CL-CPD)*WATER(2)*(TCONV(2)- & |
---|
| 874 | TCONV(1))*DPINV/CPN(1) |
---|
| 875 | FQ(1)=FQ(1)+G*MP(2)*(QP(2)-QCONV(1))* & |
---|
| 876 | DPINV+SIGD*EVAP(1) |
---|
| 877 | FQ(1)=FQ(1)+G*AM*(QCONV(2)-QCONV(1))*DPINV |
---|
| 878 | DO J=2,INB |
---|
| 879 | FQ(1)=FQ(1)+G*DPINV*MENT(J,1)*(QENT(J,1)-QCONV(1)) |
---|
| 880 | END DO |
---|
| 881 | ! |
---|
| 882 | ! *** CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO *** |
---|
| 883 | ! *** AT LEVELS ABOVE THE LOWEST LEVEL *** |
---|
| 884 | ! |
---|
| 885 | ! *** FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES *** |
---|
| 886 | ! *** THROUGH EACH LEVEL *** |
---|
| 887 | ! |
---|
| 888 | DO I=2,INB |
---|
| 889 | DPINV=0.01/(PHCONV_HPA(I)-PHCONV_HPA(I+1)) |
---|
| 890 | CPINV=1.0/CPN(I) |
---|
| 891 | AMP1=0.0 |
---|
| 892 | AD=0.0 |
---|
| 893 | IF(I.GE.NK)THEN |
---|
| 894 | DO K=I+1,INB+1 |
---|
| 895 | AMP1=AMP1+M(K) |
---|
| 896 | END DO |
---|
| 897 | END IF |
---|
| 898 | DO K=1,I |
---|
| 899 | DO J=I+1,INB+1 |
---|
| 900 | AMP1=AMP1+MENT(K,J) |
---|
| 901 | END DO |
---|
| 902 | END DO |
---|
| 903 | ! save saturated upward mass flux |
---|
| 904 | FUP(I)=AMP1 |
---|
| 905 | IF((2.*G*DPINV*AMP1).GE.DELTI)IFLAG=4 |
---|
| 906 | DO K=1,I-1 |
---|
| 907 | DO J=I,INB |
---|
| 908 | AD=AD+MENT(J,K) |
---|
| 909 | END DO |
---|
| 910 | END DO |
---|
| 911 | ! save saturated downward mass flux |
---|
| 912 | FDOWN(I)=AD |
---|
| 913 | FT(I)=FT(I)+G*DPINV*(AMP1*(TCONV(I+1)-TCONV(I)+(GZ(I+1)-GZ(I))* & |
---|
| 914 | CPINV)-AD*(TCONV(I)-TCONV(I-1)+(GZ(I)-GZ(I-1))*CPINV)) & |
---|
| 915 | -SIGD*LVCP(I)*EVAP(I) |
---|
| 916 | FT(I)=FT(I)+G*DPINV*MENT(I,I)*(HP(I)-H(I)+ & |
---|
| 917 | TCONV(I)*(CPV-CPD)*(QCONV(I)-QENT(I,I)))*CPINV |
---|
| 918 | FT(I)=FT(I)+SIGD*WT(I+1)*(CL-CPD)*WATER(I+1)* & |
---|
| 919 | (TCONV(I+1)-TCONV(I))*DPINV*CPINV |
---|
| 920 | FQ(I)=FQ(I)+G*DPINV*(AMP1*(QCONV(I+1)-QCONV(I))- & |
---|
| 921 | AD*(QCONV(I)-QCONV(I-1))) |
---|
| 922 | DO K=1,I-1 |
---|
| 923 | AWAT=ELIJ(K,I)-(1.-EP(I))*CLW(I) |
---|
| 924 | AWAT=MAX(AWAT,0.0) |
---|
| 925 | FQ(I)=FQ(I)+G*DPINV*MENT(K,I)*(QENT(K,I)-AWAT-QCONV(I)) |
---|
| 926 | END DO |
---|
| 927 | DO K=I,INB |
---|
| 928 | FQ(I)=FQ(I)+G*DPINV*MENT(K,I)*(QENT(K,I)-QCONV(I)) |
---|
| 929 | END DO |
---|
| 930 | FQ(I)=FQ(I)+SIGD*EVAP(I)+G*(MP(I+1)* & |
---|
| 931 | (QP(I+1)-QCONV(I))-MP(I)*(QP(I)-QCONV(I-1)))*DPINV |
---|
| 932 | END DO |
---|
| 933 | ! |
---|
| 934 | ! *** Adjust tendencies at top of convection layer to reflect *** |
---|
| 935 | ! *** actual position of the level zero CAPE *** |
---|
| 936 | ! |
---|
| 937 | FQOLD=FQ(INB) |
---|
| 938 | FQ(INB)=FQ(INB)*(1.-FRAC) |
---|
| 939 | FQ(INB-1)=FQ(INB-1)+FRAC*FQOLD*((PHCONV_HPA(INB)- & |
---|
| 940 | PHCONV_HPA(INB+1))/ & |
---|
| 941 | (PHCONV_HPA(INB-1)-PHCONV_HPA(INB)))*LV(INB)/LV(INB-1) |
---|
| 942 | FTOLD=FT(INB) |
---|
| 943 | FT(INB)=FT(INB)*(1.-FRAC) |
---|
| 944 | FT(INB-1)=FT(INB-1)+FRAC*FTOLD*((PHCONV_HPA(INB)- & |
---|
| 945 | PHCONV_HPA(INB+1))/ & |
---|
| 946 | (PHCONV_HPA(INB-1)-PHCONV_HPA(INB)))*CPN(INB)/CPN(INB-1) |
---|
| 947 | ! |
---|
| 948 | ! *** Very slightly adjust tendencies to force exact *** |
---|
| 949 | ! *** enthalpy, momentum and tracer conservation *** |
---|
| 950 | ! |
---|
| 951 | ENTS=0.0 |
---|
| 952 | DO I=1,INB |
---|
| 953 | ENTS=ENTS+(CPN(I)*FT(I)+LV(I)*FQ(I))* & |
---|
[8a65cb0] | 954 | (PHCONV_HPA(I)-PHCONV_HPA(I+1)) |
---|
[e200b7a] | 955 | END DO |
---|
| 956 | ENTS=ENTS/(PHCONV_HPA(1)-PHCONV_HPA(INB+1)) |
---|
| 957 | DO I=1,INB |
---|
| 958 | FT(I)=FT(I)-ENTS/CPN(I) |
---|
| 959 | END DO |
---|
| 960 | |
---|
| 961 | ! ************************************************ |
---|
| 962 | ! **** DETERMINE MASS DISPLACEMENT MATRIX |
---|
| 963 | ! ***** AND COMPENSATING SUBSIDENCE |
---|
| 964 | ! ************************************************ |
---|
| 965 | |
---|
| 966 | ! mass displacement matrix due to saturated up-and downdrafts |
---|
| 967 | ! inside the cloud and determine compensating subsidence |
---|
| 968 | ! FUP(I) (saturated updrafts), FDOWN(I) (saturated downdrafts) are assumed to be |
---|
| 969 | ! balanced by compensating subsidence (SUB(I)) |
---|
| 970 | ! FDOWN(I) and SUB(I) defined positive downwards |
---|
| 971 | |
---|
| 972 | ! NCONVTOP IS THE TOP LEVEL AT WHICH CONVECTIVE MASS FLUXES ARE DIAGNOSED |
---|
| 973 | ! EPSILON IS A SMALL NUMBER |
---|
| 974 | |
---|
| 975 | SUB(1)=0. |
---|
| 976 | NCONVTOP=1 |
---|
| 977 | do i=1,INB+1 |
---|
| 978 | do j=1,INB+1 |
---|
| 979 | if (j.eq.NK) then |
---|
| 980 | FMASS(j,i)=FMASS(j,i)+M(i) |
---|
| 981 | endif |
---|
| 982 | FMASS(j,i)=FMASS(j,i)+MENT(j,i) |
---|
| 983 | IF (FMASS(J,I).GT.EPSILON) NCONVTOP=MAX(NCONVTOP,I,J) |
---|
| 984 | end do |
---|
| 985 | if (i.gt.1) then |
---|
| 986 | SUB(i)=FUP(i-1)-FDOWN(i) |
---|
| 987 | endif |
---|
| 988 | end do |
---|
| 989 | NCONVTOP=NCONVTOP+1 |
---|
| 990 | |
---|
| 991 | RETURN |
---|
| 992 | ! |
---|
| 993 | END SUBROUTINE CONVECT |
---|
| 994 | ! |
---|
| 995 | ! --------------------------------------------------------------------------- |
---|
| 996 | ! |
---|
| 997 | SUBROUTINE TLIFT(GZ,ICB,NK,TVP,TPK,CLW,ND,NL,KK) |
---|
| 998 | ! |
---|
| 999 | !-cv |
---|
| 1000 | use par_mod |
---|
| 1001 | use conv_mod |
---|
| 1002 | |
---|
| 1003 | implicit none |
---|
| 1004 | !-cv |
---|
| 1005 | !====>Begin Module TLIFT File convect.f Undeclared variables |
---|
| 1006 | ! |
---|
| 1007 | !Argument variables |
---|
| 1008 | ! |
---|
| 1009 | integer :: icb, kk, nd, nk, nl |
---|
| 1010 | ! |
---|
| 1011 | !Local variables |
---|
| 1012 | ! |
---|
| 1013 | integer :: i, j, nsb, nst |
---|
| 1014 | ! |
---|
| 1015 | real :: ah0, ahg, alv, cpinv, cpp, denom |
---|
| 1016 | real :: es, qg, rg, s, tc, tg |
---|
| 1017 | ! |
---|
| 1018 | !====>End Module TLIFT File convect.f |
---|
| 1019 | |
---|
| 1020 | REAL :: GZ(ND),TPK(ND),CLW(ND) |
---|
| 1021 | REAL :: TVP(ND) |
---|
| 1022 | ! |
---|
| 1023 | ! *** ASSIGN VALUES OF THERMODYNAMIC CONSTANTS *** |
---|
| 1024 | ! |
---|
| 1025 | REAL,PARAMETER :: CPD=1005.7 |
---|
| 1026 | REAL,PARAMETER :: CPV=1870.0 |
---|
| 1027 | REAL,PARAMETER :: CL=2500.0 |
---|
| 1028 | REAL,PARAMETER :: RV=461.5 |
---|
| 1029 | REAL,PARAMETER :: RD=287.04 |
---|
| 1030 | REAL,PARAMETER :: LV0=2.501E6 |
---|
| 1031 | ! |
---|
| 1032 | REAL,PARAMETER :: CPVMCL=CL-CPV |
---|
| 1033 | REAL,PARAMETER :: EPS0=RD/RV |
---|
| 1034 | REAL,PARAMETER :: EPSI=1./EPS0 |
---|
| 1035 | ! |
---|
| 1036 | ! *** CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY *** |
---|
| 1037 | ! |
---|
| 1038 | AH0=(CPD*(1.-QCONV(NK))+CL*QCONV(NK))*TCONV(NK)+QCONV(NK)* & |
---|
| 1039 | (LV0-CPVMCL*( & |
---|
| 1040 | TCONV(NK)-273.15))+GZ(NK) |
---|
| 1041 | CPP=CPD*(1.-QCONV(NK))+QCONV(NK)*CPV |
---|
| 1042 | CPINV=1./CPP |
---|
| 1043 | ! |
---|
| 1044 | IF(KK.EQ.1)THEN |
---|
| 1045 | ! |
---|
| 1046 | ! *** CALCULATE LIFTED PARCEL QUANTITIES BELOW CLOUD BASE *** |
---|
| 1047 | ! |
---|
| 1048 | DO I=1,ICB-1 |
---|
| 1049 | CLW(I)=0.0 |
---|
| 1050 | END DO |
---|
| 1051 | DO I=NK,ICB-1 |
---|
| 1052 | TPK(I)=TCONV(NK)-(GZ(I)-GZ(NK))*CPINV |
---|
| 1053 | TVP(I)=TPK(I)*(1.+QCONV(NK)*EPSI) |
---|
| 1054 | END DO |
---|
| 1055 | END IF |
---|
| 1056 | ! |
---|
| 1057 | ! *** FIND LIFTED PARCEL QUANTITIES ABOVE CLOUD BASE *** |
---|
| 1058 | ! |
---|
| 1059 | NST=ICB |
---|
| 1060 | NSB=ICB |
---|
| 1061 | IF(KK.EQ.2)THEN |
---|
| 1062 | NST=NL |
---|
| 1063 | NSB=ICB+1 |
---|
| 1064 | END IF |
---|
| 1065 | DO I=NSB,NST |
---|
| 1066 | TG=TCONV(I) |
---|
| 1067 | QG=QSCONV(I) |
---|
| 1068 | ALV=LV0-CPVMCL*(TCONV(I)-273.15) |
---|
| 1069 | DO J=1,2 |
---|
| 1070 | S=CPD+ALV*ALV*QG/(RV*TCONV(I)*TCONV(I)) |
---|
| 1071 | S=1./S |
---|
| 1072 | AHG=CPD*TG+(CL-CPD)*QCONV(NK)*TCONV(I)+ALV*QG+GZ(I) |
---|
| 1073 | TG=TG+S*(AH0-AHG) |
---|
| 1074 | TG=MAX(TG,35.0) |
---|
| 1075 | TC=TG-273.15 |
---|
| 1076 | DENOM=243.5+TC |
---|
| 1077 | IF(TC.GE.0.0)THEN |
---|
| 1078 | ES=6.112*EXP(17.67*TC/DENOM) |
---|
| 1079 | ELSE |
---|
| 1080 | ES=EXP(23.33086-6111.72784/TG+0.15215*LOG(TG)) |
---|
| 1081 | END IF |
---|
| 1082 | QG=EPS0*ES/(PCONV_HPA(I)-ES*(1.-EPS0)) |
---|
| 1083 | END DO |
---|
| 1084 | ALV=LV0-CPVMCL*(TCONV(I)-273.15) |
---|
| 1085 | TPK(I)=(AH0-(CL-CPD)*QCONV(NK)*TCONV(I)-GZ(I)-ALV*QG)/CPD |
---|
| 1086 | CLW(I)=QCONV(NK)-QG |
---|
| 1087 | CLW(I)=MAX(0.0,CLW(I)) |
---|
| 1088 | RG=QG/(1.-QCONV(NK)) |
---|
| 1089 | TVP(I)=TPK(I)*(1.+RG*EPSI) |
---|
| 1090 | END DO |
---|
| 1091 | RETURN |
---|
| 1092 | END SUBROUTINE TLIFT |
---|