| 1 | !! Diagnostics: CAPE / CIN / LCL / LFC |
|---|
| 2 | !! Depend on 2D/3D flag send in |
|---|
| 3 | !! Code originally from RIP4 |
|---|
| 4 | |
|---|
| 5 | MODULE module_calc_cape |
|---|
| 6 | |
|---|
| 7 | real, dimension(150) :: buoy, zrel, benaccum |
|---|
| 8 | real, dimension(150) :: psadithte, psadiprs |
|---|
| 9 | real, dimension(150,150) :: psaditmk |
|---|
| 10 | |
|---|
| 11 | CONTAINS |
|---|
| 12 | SUBROUTINE calc_cape(SCR4, cname, cdesc, cunits, i3dflag) |
|---|
| 13 | |
|---|
| 14 | ! If i3dflag=1, this routine calculates CAPE and CIN (in m**2/s**2, |
|---|
| 15 | ! or J/kg) for every grid point in the entire 3D domain (treating |
|---|
| 16 | ! each grid point as a parcel). If i3dflag=0, then it |
|---|
| 17 | ! calculates CAPE and CIN only for the parcel with max theta-e in |
|---|
| 18 | ! the column, (i.e. something akin to Colman's MCAPE). By "parcel", |
|---|
| 19 | ! we mean a 500-m deep parcel, with actual temperature and moisture |
|---|
| 20 | ! averaged over that depth. |
|---|
| 21 | ! |
|---|
| 22 | ! In the case of i3dflag=0, |
|---|
| 23 | ! MCAPE, MCIN, LCL and LFC (2D fields are calculated) |
|---|
| 24 | |
|---|
| 25 | |
|---|
| 26 | USE constants_module |
|---|
| 27 | USE module_model_basics |
|---|
| 28 | |
|---|
| 29 | IMPLICIT NONE |
|---|
| 30 | |
|---|
| 31 | !Arguments |
|---|
| 32 | real, allocatable, dimension(:,:,:,:) :: SCR4 |
|---|
| 33 | character (len=128) :: cname, cdesc, cunits |
|---|
| 34 | integer :: i3dflag |
|---|
| 35 | |
|---|
| 36 | ! Local variables |
|---|
| 37 | integer :: i, j, k, kk, jt, ip, iustnlist |
|---|
| 38 | integer :: kpar, kpar1, kpar2, kmax, klev, kel |
|---|
| 39 | integer :: ilcl, klcl, klfc |
|---|
| 40 | integer :: nprs, nthte |
|---|
| 41 | real, dimension(west_east_dim,south_north_dim) :: ter |
|---|
| 42 | real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: prs, tmk, ght, qvp |
|---|
| 43 | real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: prsf, cape, cin |
|---|
| 44 | real :: ethpari, zlcl, tvenv |
|---|
| 45 | real :: p1, p2, pp1, pp2, pup, pdn |
|---|
| 46 | real :: totprs, totqvp, totthe |
|---|
| 47 | real :: eslift, ghtlift, qvplift, tmklift, tvlift |
|---|
| 48 | real :: ghtpari, prspari, qvppari, tmkpari |
|---|
| 49 | real :: tmkenv, qvpenv, tlcl |
|---|
| 50 | real :: fac1, fac2, facden, th, deltap |
|---|
| 51 | real :: benamin, davg, pavg, pressure, temp |
|---|
| 52 | real :: e, eth, ethmax, q, dz, cpm |
|---|
| 53 | character (len=20) :: fname |
|---|
| 54 | logical :: is_used |
|---|
| 55 | |
|---|
| 56 | |
|---|
| 57 | SCR4 = 0.0 |
|---|
| 58 | |
|---|
| 59 | |
|---|
| 60 | ! Open psadilookup.dat |
|---|
| 61 | DO iustnlist = 10,100 |
|---|
| 62 | INQUIRE(unit=iustnlist, opened=is_used) |
|---|
| 63 | if (.not. is_used) exit |
|---|
| 64 | END DO |
|---|
| 65 | fname = 'src/psadilookup.dat' |
|---|
| 66 | OPEN (unit=iustnlist,file=fname,form='formatted',status='old') |
|---|
| 67 | |
|---|
| 68 | DO i = 1,14 |
|---|
| 69 | READ (iustnlist,*) |
|---|
| 70 | ENDDO |
|---|
| 71 | READ (iustnlist,*) nthte,nprs |
|---|
| 72 | IF ( nthte.ne.150 .OR. nprs.ne.150 ) then |
|---|
| 73 | PRINT*, 'Number of pressure or theta_e levels in lookup table' |
|---|
| 74 | PRINT*, ' file not = 150. Check lookup table file.' |
|---|
| 75 | STOP |
|---|
| 76 | ENDIF |
|---|
| 77 | READ (iustnlist,173) (psadithte(jt),jt=1,nthte) |
|---|
| 78 | READ (iustnlist,173) (psadiprs(ip),ip=1,nprs) |
|---|
| 79 | READ (iustnlist,173) ((psaditmk(ip,jt),ip=1,nprs),jt=1,nthte) |
|---|
| 80 | 173 FORMAT (5e15.7) |
|---|
| 81 | CLOSE (iustnlist) |
|---|
| 82 | |
|---|
| 83 | |
|---|
| 84 | !! Get fields we want from the ones we have |
|---|
| 85 | ter = HGT |
|---|
| 86 | qvp = QV |
|---|
| 87 | prs = PRES * 0.01 ! pressure in hPa |
|---|
| 88 | tmk = TK ! temperature in K |
|---|
| 89 | ght = (PH+PHB) / G ! height in m |
|---|
| 90 | |
|---|
| 91 | |
|---|
| 92 | !! First calculate a pressure array on FULL SIGMA levels |
|---|
| 93 | !! Similar to the pfcalc.f routine from RIP4 |
|---|
| 94 | !! Top full sigma level is ommitted |
|---|
| 95 | prsf(:,:,1) = PSFC(:,:) !! Lowest full sigma set to surface pressure |
|---|
| 96 | DO k = 2, bottom_top_dim |
|---|
| 97 | prsf(:,:,k)=.5*(prs(:,:,k-1)+prs(:,:,k)) |
|---|
| 98 | END DO |
|---|
| 99 | |
|---|
| 100 | |
|---|
| 101 | cape = 0.0 |
|---|
| 102 | cin = 0.0 |
|---|
| 103 | |
|---|
| 104 | DO j = 1,south_north_dim !! BIG i/j loop |
|---|
| 105 | DO i = 1,west_east_dim |
|---|
| 106 | |
|---|
| 107 | IF ( i3dflag == 1 ) THEN !! 3D case |
|---|
| 108 | |
|---|
| 109 | kpar1=bottom_top_dim-1 |
|---|
| 110 | kpar2=1 |
|---|
| 111 | |
|---|
| 112 | ELSE !! 2D case |
|---|
| 113 | |
|---|
| 114 | ! Find parcel with max theta-e in lowest 3 km AGL. |
|---|
| 115 | ethmax = -1. |
|---|
| 116 | DO k = 1, bottom_top_dim |
|---|
| 117 | IF ( ght(i,j,k)-ter(i,j) .lt. 3000. ) THEN |
|---|
| 118 | q = max(qvp(i,j,k),1.e-15) |
|---|
| 119 | temp = tmk(i,j,k) |
|---|
| 120 | pressure = prs(i,j,k) |
|---|
| 121 | e = (q*pressure)/(EPS+q) |
|---|
| 122 | tlcl = TLCLC1/(log(temp**TLCLC2/e)-TLCLC3)+TLCLC4 |
|---|
| 123 | eth = temp*(1000./pressure)**( GAMMA_RIP*(1.+GAMMAMD*q) )* & |
|---|
| 124 | exp( (THTECON1/tlcl-THTECON2)*q*(1.+THTECON3*q) ) |
|---|
| 125 | IF ( eth .gt. ethmax ) THEN |
|---|
| 126 | klev=k |
|---|
| 127 | ethmax=eth |
|---|
| 128 | END IF |
|---|
| 129 | END IF |
|---|
| 130 | END DO |
|---|
| 131 | kpar1=klev |
|---|
| 132 | kpar2=klev |
|---|
| 133 | |
|---|
| 134 | ! Establish average properties of that parcel |
|---|
| 135 | ! (over depth of approximately davg meters) |
|---|
| 136 | davg = 500. |
|---|
| 137 | pavg = davg*prs(i,j,kpar1)*G / & |
|---|
| 138 | ( Rd*virtual(tmk(i,j,kpar1),qvp(i,j,kpar1)) ) |
|---|
| 139 | p2 = min(prs(i,j,kpar1)+.5*pavg,prsf(i,j,1)) |
|---|
| 140 | p1 = p2-pavg |
|---|
| 141 | totthe = 0. |
|---|
| 142 | totqvp = 0. |
|---|
| 143 | totprs = 0. |
|---|
| 144 | DO k = 1,bottom_top_dim-1 |
|---|
| 145 | IF ( prsf(i,j,k) .le. p1 ) GOTO 35 |
|---|
| 146 | IF ( prsf(i,j,k+1) .ge. p2 ) GOTO 34 |
|---|
| 147 | pressure = prs(i,j,k) |
|---|
| 148 | pup = prsf(i,j,k) |
|---|
| 149 | pdn = prsf(i,j,k+1) |
|---|
| 150 | q = max(qvp(i,j,k),1.e-15) |
|---|
| 151 | th = tmk(i,j,k)*(1000./pressure)**(GAMMA_RIP*(1.+GAMMAMD*q)) |
|---|
| 152 | pp1 = max(p1,pdn) |
|---|
| 153 | pp2 = min(p2,pup) |
|---|
| 154 | IF ( pp2 .gt. pp1 ) THEN |
|---|
| 155 | deltap = pp2-pp1 |
|---|
| 156 | totqvp = totqvp+q*deltap |
|---|
| 157 | totthe = totthe+th*deltap |
|---|
| 158 | totprs = totprs+deltap |
|---|
| 159 | END IF |
|---|
| 160 | 34 CONTINUE |
|---|
| 161 | END DO |
|---|
| 162 | 35 CONTINUE |
|---|
| 163 | qvppari = totqvp/totprs |
|---|
| 164 | tmkpari = (totthe/totprs)*(prs(i,j,kpar1)/1000.)** & |
|---|
| 165 | (GAMMA_RIP*(1.+GAMMAMD*qvp(i,j,kpar1))) |
|---|
| 166 | ENDIF |
|---|
| 167 | |
|---|
| 168 | !!! END of 2D / 3D specific bits |
|---|
| 169 | |
|---|
| 170 | |
|---|
| 171 | DO kpar = kpar1,kpar2,-1 !! This loop is done for both 2D / 3D |
|---|
| 172 | |
|---|
| 173 | ! Calculate temperature and moisture properties of parcel |
|---|
| 174 | |
|---|
| 175 | IF ( i3dflag == 1 ) THEN ! (Note, qvppari and tmkpari already calculated above for 2D case.) |
|---|
| 176 | qvppari = qvp(i,j,kpar) |
|---|
| 177 | tmkpari = tmk(i,j,kpar) |
|---|
| 178 | END IF |
|---|
| 179 | prspari = prs(i,j,kpar) |
|---|
| 180 | ghtpari = ght(i,j,kpar) |
|---|
| 181 | cpm = cp*(1.+CPMD*qvppari) |
|---|
| 182 | |
|---|
| 183 | e = max(1.e-20,qvppari*prspari/(EPS+qvppari)) |
|---|
| 184 | tlcl = TLCLC1/(log(tmkpari**TLCLC2/e)-TLCLC3)+TLCLC4 |
|---|
| 185 | ethpari = tmkpari*(1000./prspari)**(GAMMA_RIP*(1.+GAMMAMD*qvppari))* & |
|---|
| 186 | exp((THTECON1/tlcl-THTECON2)*qvppari* & |
|---|
| 187 | (1.+THTECON3*qvppari)) |
|---|
| 188 | zlcl = ghtpari+(tmkpari-tlcl)/(G/cpm) |
|---|
| 189 | |
|---|
| 190 | ! Calculate buoyancy and relative height of lifted parcel at |
|---|
| 191 | ! all levels, and store in bottom up arrays. Add a level at the LCL, |
|---|
| 192 | ! and at all points where buoyancy is zero. |
|---|
| 193 | |
|---|
| 194 | kk = 0 ! for arrays that go bottom to top |
|---|
| 195 | ilcl = 0 |
|---|
| 196 | IF ( ghtpari .ge. zlcl ) THEN |
|---|
| 197 | !! initial parcel already saturated or supersaturated. |
|---|
| 198 | ilcl = 2 |
|---|
| 199 | klcl = 1 |
|---|
| 200 | END IF |
|---|
| 201 | |
|---|
| 202 | DO k = kpar,bottom_top_dim |
|---|
| 203 | 33 kk=kk+1 ! for arrays that go bottom to top |
|---|
| 204 | |
|---|
| 205 | IF ( ght(i,j,k) .lt. zlcl ) THEN ! model level is below LCL |
|---|
| 206 | qvplift = qvppari |
|---|
| 207 | tmklift = tmkpari-G/cpm*(ght(i,j,k)-ghtpari) |
|---|
| 208 | tvenv = virtual(tmk(i,j,k),qvp(i,j,k)) |
|---|
| 209 | tvlift = virtual(tmklift,qvplift) |
|---|
| 210 | ghtlift = ght(i,j,k) |
|---|
| 211 | ELSE IF ( ght(i,j,k) .ge. zlcl .AND. ilcl .eq. 0 ) THEN |
|---|
| 212 | !! This model level and previous model level straddle the LCL, |
|---|
| 213 | !! so first create a new level in the bottom-up array, at the LCL. |
|---|
| 214 | tmklift = tlcl |
|---|
| 215 | qvplift = qvppari |
|---|
| 216 | facden = ght(i,j,k)-ght(i,j,k-1) |
|---|
| 217 | fac1 = (zlcl-ght(i,j,k-1))/facden |
|---|
| 218 | fac2 = (ght(i,j,k)-zlcl)/facden |
|---|
| 219 | tmkenv = tmk(i,j,k-1)*fac2+tmk(i,j,k)*fac1 |
|---|
| 220 | qvpenv = qvp(i,j,k-1)*fac2+qvp(i,j,k)*fac1 |
|---|
| 221 | tvenv = virtual(tmkenv,qvpenv) |
|---|
| 222 | tvlift = virtual(tmklift,qvplift) |
|---|
| 223 | ghtlift = zlcl |
|---|
| 224 | ilcl = 1 |
|---|
| 225 | ELSE |
|---|
| 226 | tmklift = tonpsadiabat(ethpari,prs(i,j,k)) |
|---|
| 227 | eslift = EZERO*exp(eslcon1*(tmklift-CELKEL)/ & |
|---|
| 228 | (tmklift-eslcon2)) |
|---|
| 229 | qvplift = EPS*eslift/(prs(i,j,k)-eslift) |
|---|
| 230 | tvenv = virtual(tmk(i,j,k),qvp(i,j,k)) |
|---|
| 231 | tvlift = virtual(tmklift,qvplift) |
|---|
| 232 | ghtlift = ght(i,j,k) |
|---|
| 233 | END IF |
|---|
| 234 | |
|---|
| 235 | buoy(kk) = G*(tvlift-tvenv)/tvenv ! buoyancy |
|---|
| 236 | zrel(kk) = ghtlift-ghtpari |
|---|
| 237 | |
|---|
| 238 | IF ( (buoy(kk)*buoy(kk-1).lt.0.0) .AND. (kk.gt.1) ) THEN |
|---|
| 239 | !! Parcel ascent curve crosses sounding curve, so create a new level |
|---|
| 240 | !! in the bottom-up array at the crossing. |
|---|
| 241 | kk = kk+1 |
|---|
| 242 | buoy(kk) = buoy(kk-1) |
|---|
| 243 | zrel(kk) = zrel(kk-1) |
|---|
| 244 | buoy(kk-1) = 0. |
|---|
| 245 | zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/(buoy(kk-2)-buoy(kk))* & |
|---|
| 246 | (zrel(kk)-zrel(kk-2)) |
|---|
| 247 | END IF |
|---|
| 248 | |
|---|
| 249 | IF (ilcl == 1) THEN |
|---|
| 250 | klcl = kk |
|---|
| 251 | ilcl = 2 |
|---|
| 252 | GOTO 33 |
|---|
| 253 | END IF |
|---|
| 254 | |
|---|
| 255 | END DO !! END DO k = kpar,bottom_top_dim |
|---|
| 256 | |
|---|
| 257 | kmax = kk |
|---|
| 258 | IF (kmax .gt. 150) THEN |
|---|
| 259 | PRINT*, 'in calc_cape: kmax got too big. kmax=',kmax |
|---|
| 260 | STOP |
|---|
| 261 | ENDIF |
|---|
| 262 | |
|---|
| 263 | |
|---|
| 264 | ! Get the accumulated buoyant energy from the parcel's starting |
|---|
| 265 | ! point, at all levels up to the top level. |
|---|
| 266 | |
|---|
| 267 | benaccum(1) = 0.0 |
|---|
| 268 | benamin = 9e9 |
|---|
| 269 | DO k = 2,kmax |
|---|
| 270 | dz = zrel(k)-zrel(k-1) |
|---|
| 271 | benaccum(k) = benaccum(k-1)+.5*dz*(buoy(k-1)+buoy(k)) |
|---|
| 272 | IF ( benaccum(k) .lt. benamin ) THEN |
|---|
| 273 | benamin = benaccum(k) |
|---|
| 274 | END IF |
|---|
| 275 | END DO |
|---|
| 276 | |
|---|
| 277 | |
|---|
| 278 | ! Determine equilibrium level (EL), which we define as the highest |
|---|
| 279 | ! level of non-negative buoyancy above the LCL. Note, this may be |
|---|
| 280 | ! the top level if the parcel is still buoyant there. |
|---|
| 281 | |
|---|
| 282 | DO k = kmax,klcl,-1 |
|---|
| 283 | IF ( buoy(k) .ge. 0. ) THEN |
|---|
| 284 | kel = k ! k of equilibrium level |
|---|
| 285 | GOTO 50 |
|---|
| 286 | END IF |
|---|
| 287 | END DO |
|---|
| 288 | |
|---|
| 289 | |
|---|
| 290 | ! If we got through that loop, then there is no non-negative |
|---|
| 291 | ! buoyancy above the LCL in the sounding. In these situations, |
|---|
| 292 | ! both CAPE and CIN will be set to -0.1 J/kg. Also, where CAPE is |
|---|
| 293 | ! non-zero, CAPE and CIN will be set to a minimum of +0.1 J/kg, so |
|---|
| 294 | ! that the zero contour in either the CIN or CAPE fields will |
|---|
| 295 | ! circumscribe regions of non-zero CAPE. |
|---|
| 296 | |
|---|
| 297 | cape(i,j,kpar) = -0.1 |
|---|
| 298 | cin(i,j,kpar) = -0.1 |
|---|
| 299 | klfc = kmax |
|---|
| 300 | |
|---|
| 301 | GOTO 102 |
|---|
| 302 | |
|---|
| 303 | 50 CONTINUE |
|---|
| 304 | |
|---|
| 305 | |
|---|
| 306 | |
|---|
| 307 | ! If there is an equilibrium level, then CAPE is positive. We'll |
|---|
| 308 | ! define the level of free convection (LFC) as the point below the |
|---|
| 309 | ! EL, but at or above the LCL, where accumulated buoyant energy is a |
|---|
| 310 | ! minimum. The net positive area (accumulated buoyant energy) from |
|---|
| 311 | ! the LFC up to the EL will be defined as the CAPE, and the net |
|---|
| 312 | ! negative area (negative of accumulated buoyant energy) from the |
|---|
| 313 | ! parcel starting point to the LFC will be defined as the convective |
|---|
| 314 | ! inhibition (CIN). |
|---|
| 315 | |
|---|
| 316 | ! First get the LFC according to the above definition. |
|---|
| 317 | |
|---|
| 318 | benamin = 9e9 |
|---|
| 319 | klfc = kmax |
|---|
| 320 | DO k = klcl,kel |
|---|
| 321 | IF ( benaccum(k) .lt. benamin ) THEN |
|---|
| 322 | benamin = benaccum(k) |
|---|
| 323 | klfc = k |
|---|
| 324 | END IF |
|---|
| 325 | END DO |
|---|
| 326 | |
|---|
| 327 | ! Now we can assign values to cape and cin |
|---|
| 328 | |
|---|
| 329 | cape(i,j,kpar) = max(benaccum(kel)-benamin,0.1) |
|---|
| 330 | cin(i,j,kpar) = max(-benamin,0.1) |
|---|
| 331 | |
|---|
| 332 | ! CIN is uninteresting when CAPE is small (< 100 J/kg), so set |
|---|
| 333 | ! CIN to -.1 in that case. |
|---|
| 334 | |
|---|
| 335 | IF ( cape(i,j,kpar) .lt. 100. ) cin(i,j,kpar) = -0.1 |
|---|
| 336 | |
|---|
| 337 | 102 CONTINUE |
|---|
| 338 | |
|---|
| 339 | ENDDO !! END of BIG 2D/3D loop |
|---|
| 340 | |
|---|
| 341 | |
|---|
| 342 | IF ( i3dflag == 0 ) THEN |
|---|
| 343 | SCR4(i,j,1,1) = cape(i,j,kpar1) |
|---|
| 344 | SCR4(i,j,1,2) = cin(i,j,kpar1) |
|---|
| 345 | SCR4(i,j,1,3) = zrel(klcl)+ghtpari-ter(i,j) ! meters AGL (LCL) |
|---|
| 346 | SCR4(i,j,1,4) = zrel(klfc)+ghtpari-ter(i,j) ! meters AGL (LFC) |
|---|
| 347 | ENDIF |
|---|
| 348 | |
|---|
| 349 | |
|---|
| 350 | END DO |
|---|
| 351 | END DO !! END BIG i/j loop |
|---|
| 352 | |
|---|
| 353 | |
|---|
| 354 | !! These will be set by module_diagnostics as we have more than 1 field |
|---|
| 355 | |
|---|
| 356 | IF ( i3dflag == 1 ) THEN |
|---|
| 357 | SCR4(:,:,:,1) = cape(:,:,:) |
|---|
| 358 | SCR4(:,:,:,2) = cin(:,:,:) |
|---|
| 359 | ENDIF |
|---|
| 360 | |
|---|
| 361 | cname = " " |
|---|
| 362 | cdesc = " " |
|---|
| 363 | cunits = " " |
|---|
| 364 | |
|---|
| 365 | |
|---|
| 366 | END SUBROUTINE calc_cape |
|---|
| 367 | |
|---|
| 368 | |
|---|
| 369 | !*********************************************************************c |
|---|
| 370 | !*********************************************************************c |
|---|
| 371 | FUNCTION tonpsadiabat (thte,prs) |
|---|
| 372 | |
|---|
| 373 | ! This function gives the temperature (in K) on a moist adiabat |
|---|
| 374 | ! (specified by thte in K) given pressure in hPa. It uses a |
|---|
| 375 | ! lookup table, with data that was generated by the Bolton (1980) |
|---|
| 376 | ! formula for theta_e. |
|---|
| 377 | |
|---|
| 378 | |
|---|
| 379 | ! First check if pressure is less than min pressure in lookup table. |
|---|
| 380 | ! If it is, assume parcel is so dry that the given theta-e value can |
|---|
| 381 | ! be interpretted as theta, and get temperature from the simple dry |
|---|
| 382 | ! theta formula. |
|---|
| 383 | |
|---|
| 384 | IF ( prs .lt. psadiprs(150) ) THEN |
|---|
| 385 | tonpsadiabat = thte*(prs/1000.)**GAMMA_RIP |
|---|
| 386 | RETURN |
|---|
| 387 | ENDIF |
|---|
| 388 | |
|---|
| 389 | |
|---|
| 390 | ! Otherwise, look for the given thte/prs point in the lookup table. |
|---|
| 391 | |
|---|
| 392 | DO jtch = 1,150-1 |
|---|
| 393 | IF ( thte.ge.psadithte(jtch) .AND. thte.lt.psadithte(jtch+1) ) THEN |
|---|
| 394 | jt = jtch |
|---|
| 395 | GOTO 213 |
|---|
| 396 | END IF |
|---|
| 397 | END DO |
|---|
| 398 | jt = -1 |
|---|
| 399 | 213 CONTINUE |
|---|
| 400 | |
|---|
| 401 | DO ipch = 1,150-1 |
|---|
| 402 | if ( prs.le.psadiprs(ipch) .AND. prs.gt.psadiprs(ipch+1) ) THEN |
|---|
| 403 | ip = ipch |
|---|
| 404 | GOTO 215 |
|---|
| 405 | END IF |
|---|
| 406 | ENDDO |
|---|
| 407 | ip = -1 |
|---|
| 408 | 215 CONTINUE |
|---|
| 409 | |
|---|
| 410 | |
|---|
| 411 | IF ( jt.eq.-1 .OR. ip.eq.-1 ) THEN |
|---|
| 412 | PRINT*, 'Outside of lookup table bounds. prs,thte=',prs,thte |
|---|
| 413 | STOP |
|---|
| 414 | ENDIF |
|---|
| 415 | |
|---|
| 416 | |
|---|
| 417 | fracjt = (thte-psadithte(jt))/(psadithte(jt+1)-psadithte(jt)) |
|---|
| 418 | fracjt2 = 1.-fracjt |
|---|
| 419 | fracip = (psadiprs(ip)-prs)/(psadiprs(ip)-psadiprs(ip+1)) |
|---|
| 420 | fracip2 = 1.-fracip |
|---|
| 421 | |
|---|
| 422 | IF ( psaditmk(ip,jt ).gt.1e9 .OR. psaditmk(ip+1,jt ).gt.1e9 .OR. & |
|---|
| 423 | psaditmk(ip,jt+1).gt.1e9 .OR. psaditmk(ip+1,jt+1).gt.1e9 ) THEN |
|---|
| 424 | PRINT*, 'Tried to access missing tmperature in lookup table.' |
|---|
| 425 | PRINT*, 'Prs and Thte probably unreasonable. prs,thte=',prs,thte |
|---|
| 426 | STOP |
|---|
| 427 | ENDIF |
|---|
| 428 | |
|---|
| 429 | tonpsadiabat = fracip2*fracjt2*psaditmk(ip ,jt )+ & |
|---|
| 430 | fracip *fracjt2*psaditmk(ip+1,jt )+ & |
|---|
| 431 | fracip2*fracjt *psaditmk(ip ,jt+1)+ & |
|---|
| 432 | fracip *fracjt *psaditmk(ip+1,jt+1) |
|---|
| 433 | |
|---|
| 434 | |
|---|
| 435 | END FUNCTION tonpsadiabat |
|---|
| 436 | |
|---|
| 437 | |
|---|
| 438 | END MODULE module_calc_cape |
|---|
