12 use mapl_constantsmod, kp => mapl_r8
16 use machine,
only : kp => kind_phys
52 public hlatv, tmin, hlatf, rgasv, pcf, cp, epsqs, ttrice
57 parameter(plenest=250)
63 real(kp) estbl(plenest)
76 integer,
parameter :: iulog=6
80 real(kp) function estblf( td )
84 real(kp),
intent(in) :: td
90 e = max(min(td,tmax),tmin)
93 estblf = (tmin+ai-e+1._kp)* estbl(i)-(tmin+ai-e)* estbl(i+1)
106 subroutine gestbl(tmn ,tmx ,trice ,ip ,epsil , latvap ,latice , &
107 & rh2o ,cpair ,tmeltx )
112 real(kp),
intent(in) :: tmn
113 real(kp),
intent(in) :: tmx
114 real(kp),
intent(in) :: epsil
115 real(kp),
intent(in) :: trice
116 real(kp),
intent(in) :: latvap
117 real(kp),
intent(in) :: latice
118 real(kp),
intent(in) :: rh2o
119 real(kp),
intent(in) :: cpair
120 real(kp),
intent(in) :: tmeltx
150 lentbl = int(tmax-tmin+2.000001_kp)
151 if (lentbl .gt. plenest)
then
155 write(*,*)
"AHHH wv_sat"
164 if (ttrice /= 0.0_kp)
then
176 call gffgch(t,estbl(n),tmelt,itype)
179 do n=lentbl+1,plenest
180 estbl(n) = -99999.0_kp
190 pcf(1) = 5.04469588506e-01_kp
191 pcf(2) = -5.47288442819e+00_kp
192 pcf(3) = -3.67471858735e-01_kp
193 pcf(4) = -8.95963532403e-03_kp
194 pcf(5) = -7.78053686625e-05_kp
2149000
format(
'GESTBL: FATAL ERROR ********************************* &
215 &',/,
' TMAX AND TMIN REQUIRE A LARGER DIMENSION ON THE LENGTH', &
216 &
' OF THE SATURATION VAPOR PRESSURE TABLE ESTBL(PLENEST)',/, &
217 &
' TMAX, TMIN, AND PLENEST => ', 2f7.2, i3)
227 subroutine aqsat(t ,p ,es ,qs ,ii , ilen ,kk ,kstart ,kend )
232 integer,
intent(in) :: ii
233 integer,
intent(in) :: kk
234 integer,
intent(in) :: ilen
235 integer,
intent(in) :: kstart
236 integer,
intent(in) :: kend
237 real(kp),
intent(in) :: t(ii,kk)
238 real(kp),
intent(in) :: p(ii,kk)
242 real(kp),
intent(out) :: es(ii,kk)
243 real(kp),
intent(out) :: qs(ii,kk)
252 omeps = 1.0_kp - epsqs
255 es(i,k) = min(estblf(t(i,k)),p(i,k))
259 qs(i,k) = min(1.0_kp, epsqs*es(i,k)/(p(i,k)-omeps*es(i,k)))
288 integer,
intent(in) :: ii
289 integer,
intent(in) :: kk
290 integer,
intent(in) :: ilen
291 integer,
intent(in) :: kstart
292 integer,
intent(in) :: kend
293 real(kp),
intent(in) :: t(ii,kk)
294 real(kp),
intent(in) :: p(ii,kk)
298 real(kp),
intent(out) :: es(ii,kk)
299 real(kp),
intent(out) :: qs(ii,kk)
308 omeps = 1.0_kp - epsqs
313 es(i,k) = min(
polysvp(t(i,k),0), p(i,k))
316 es(i,k) = min(fpvsl(t(i,k)), p(i,k))
321 qs(i,k) = min(1.0_kp, epsqs*es(i,k)/(p(i,k)-omeps*es(i,k)))
348 subroutine aqsatd(t, p, es, qs, gam, ii, ilen, kk, kstart, kend)
353 integer,
intent(in) :: ii
354 integer,
intent(in) :: ilen
355 integer,
intent(in) :: kk
356 integer,
intent(in) :: kstart
357 integer,
intent(in) :: kend
359 real(kp),
intent(in) :: t(ii,kk)
360 real(kp),
intent(in) :: p(ii,kk)
365 real(kp),
intent(out) :: es(ii,kk)
366 real(kp),
intent(out) :: qs(ii,kk)
367 real(kp),
intent(out) :: gam(ii,kk)
386 omeps = 1.0_kp - epsqs
389 es(i,k) = min(p(i,k), estblf(t(i,k)))
393 qs(i,k) = min(1.0_kp, epsqs*es(i,k)/(p(i,k)-omeps*es(i,k)))
410 if ((.not. icephs) .or. (ttrice == 0.0_kp))
go to 10
411 trinv = 1.0_kp/ttrice
424 lflg = (tc >= -ttrice .and. tc < 0.0_kp)
425 weight = min(-tc*trinv,1.0_kp)
426 hlatsb = hlatv + weight*hlatf
427 hlatvp = hlatv - 2369.0_kp*tc
428 if (t(i,k) < tmelt)
then
434 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) &
439 desdt = hltalt*es(i,k)/(rgasv*t(i,k)*t(i,k)) + tterm*trinv
440 gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/(cp*es(i,k)*(p(i,k) &
442 if (qs(i,k) == 1.0_kp) gam(i,k) = 0.0_kp
456 hlatvp = hlatv - 2369.0_kp*(t(i,k)-tmelt)
458 hlatsb = hlatv + hlatf
462 if (t(i,k) < tmelt)
then
467 desdt = hltalt*es(i,k)/(rgasv*t(i,k)*t(i,k))
468 gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/(cp*es(i,k)*(p(i,k) &
470 if (qs(i,k) == 1.0_kp) gam(i,k) = 0.0_kp
484 subroutine vqsatd(t ,p ,es ,qs ,gam , len )
489 integer,
intent(in) :: len
490 real(kp),
intent(in) :: t(len)
491 real(kp),
intent(in) :: p(len)
495 real(kp),
intent(out) :: es(len)
496 real(kp),
intent(out) :: qs(len)
497 real(kp),
intent(out) :: gam(len)
518 omeps = 1.0_kp - epsqs
520 es(i) = min(estblf(t(i)), p(i))
524 qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
529 qs(i) = min(1.0_kp,qs(i))
542 if ((.not. icephs) .or. (ttrice.eq.0.0_kp))
go to 10
543 trinv = 1.0_kp/ttrice
554 lflg = (tc >= -ttrice .and. tc < 0.0_kp)
555 weight = min(-tc*trinv,1.0_kp)
556 hlatsb = hlatv + weight*hlatf
557 hlatvp = hlatv - 2369.0_kp*tc
558 if (t(i) < tmelt)
then
564 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) &
569 desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv
570 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
571 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
582 hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt)
584 hlatsb = hlatv + hlatf
588 if (t(i) < tmelt)
then
593 desdt = hltalt*es(i)/(rgasv*t(i)*t(i))
594 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
595 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
611 integer,
intent(in) :: len
612 real(kp),
intent(in) :: t(len)
613 real(kp),
intent(in) :: p(len)
618 real(kp),
intent(out) :: es(len)
619 real(kp),
intent(out) :: qs(len)
620 real(kp),
intent(out) :: gam(len)
637 omeps = 1.0_kp - epsqs
640 es(i) = min(fpvsl(t(i)), p(i))
642 es(i) = min(
polysvp(t(i),0), p(i))
647 qs(i) = min(1.0_kp, epsqs*es(i) / (p(i)-omeps*es(i)))
668 hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt)
670 if (t(i) < tmelt)
then
675 desdt = hltalt*es(i)/(rgasv*t(i)*t(i))
676 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
677 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
705 polysvp = murphykoop_svp_ice(t)
718 polysvp = 10._kp**(-9.09718_kp*(273.16_kp/t-1._kp)-3.56654_kp* &
719 & log10(273.16_kp/t)+0.876793_kp*(1._kp-t/273.16_kp)+ &
720 & log10(6.1071_kp))*100._kp
727 polysvp = 10._kp**(-7.90298_kp*(373.16_kp/t-1._kp)+ 5.02808_kp* &
728 &log10(373.16_kp/t)- 1.3816e-7_kp*(10._kp**(11.344_kp*(1._kp-t/ &
729 &373.16_kp))-1._kp)+ 8.1328e-3_kp*(10._kp**(-3.49149_kp*(373.16_kp/ &
730 &t-1._kp))-1._kp)+ log10(1013.246_kp))*100._kp
740 integer function fqsatd(t ,p ,es ,qs ,gam , len )
746 integer,
intent(in) :: len
747 real(kp),
intent(in) :: t(len)
748 real(kp),
intent(in) :: p(len)
750 real(kp),
intent(out) :: es(len)
751 real(kp),
intent(out) :: qs(len)
752 real(kp),
intent(out) :: gam(len)
754 call vqsatd(t ,p ,es ,qs ,gam , len )
759 real(kp) function qsat_water(t,p)
764 real(kp) ps, ts, e1, e2, f1, f2, f3, f4, f5, f
772 e1 = 11.344_kp*(1.0_kp - t/ts)
773 e2 = -3.49149_kp*(ts/t - 1.0_kp)
774 f1 = -7.90298_kp*(ts/t - 1.0_kp)
775 f2 = 5.02808_kp*log10(ts/t)
776 f3 = -1.3816_kp*(10.0_kp**e1 - 1.0_kp)/10000000.0_kp
777 f4 = 8.1328_kp*(10.0_kp**e2 - 1.0_kp)/1000.0_kp
779 f = f1 + f2 + f3 + f4 + f5
780 es = (10.0_kp**f)*100.0_kp
782 qsat_water = epsqs*es/(p-(1.-epsqs)*es)
783 if(qsat_water < 0.) qsat_water = 1.
785 end function qsat_water
791 integer,
intent(in) :: len
794 real(kp) qsat_water(len)
796 real(kp),
parameter :: t0inv = 1._kp/273._kp
802 es = 611._kp*exp(coef*(t0inv-1./t(i)))
803 qsat_water(i) = epsqs*es/(p(i)-(1.-epsqs)*es)
804 if(qsat_water(i) < 0.) qsat_water(i) = 1.
817 real(kp),
parameter :: t0inv = 1._kp/273._kp
818 es = 611.*exp((hlatv+hlatf)/rgasv*(t0inv-1./t))
819 qsat_ice = epsqs*es/(p-(1.-epsqs)*es)
827 integer,
intent(in) :: len
832 real(kp),
parameter :: t0inv = 1._kp/273._kp
836 coef = (hlatv+hlatf)/rgasv
838 es = 611.*exp(coef*(t0inv-1./t(i)))
839 qsat_ice(i) = epsqs*es/(p(i)-(1.-epsqs)*es)
859 integer,
intent(in) :: len
860 real(kp),
intent(in) :: t(len)
861 real(kp),
intent(in) :: p(len)
866 real(kp),
intent(out) :: es(len)
867 real(kp),
intent(out) :: qs(len)
870 real(kp),
intent(out) :: dqsdt(len)
893 omeps = 1.0_kp - epsqs
896 es(i) = min(
polysvp(t(i),0), p(i))
899 es(i) = min(fpvsl(t(i)), p(i))
904 qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
909 qs(i) = min(1.0_kp,qs(i))
925 hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt)
927 if (t(i) < tmelt)
then
932 desdt = hltalt*es(i)/(rgasv*t(i)*t(i))
933 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
934 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
936 dqsdt(i) = (cp/hltalt)*gam(i)
952 real(kp),
intent(in) :: t, p
957 real(kp),
intent(out) :: es, qs, dqsdt
963 real(kp) omeps, hltalt, hlatsb, hlatvp, desdt, gam
967 omeps = 1.0_kp - epsqs
973 es = min(p, fpvsl(t))
978 qs = min(1.0_kp, epsqs*es/(p-omeps*es))
996 hlatvp = hlatv - 2369.0_kp*(t-tmelt)
1003 desdt = hltalt*es/(rgasv*t*t)
1004 gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
1005 if (qs >= 1.0_kp) gam = 0.0_kp
1007 dqsdt = (cp/hltalt)*gam
1036 integer,
intent(in) :: len
1037 real(kp),
intent(in) :: t(len)
1038 real(kp),
intent(in) :: p(len)
1042 real(kp),
intent(out) :: es(len)
1043 real(kp),
intent(out) :: qs(len)
1046 real(kp),
intent(out) :: dqsdt(len)
1073 omeps = 1.0_kp - epsqs
1076 es(i) = min(p(i), estblf(t(i)))
1079 es(i) = min(p(i), fpvsi(t(i)))
1084 qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
1089 qs(i) = min(1.0_kp,qs(i))
1101 if ((.not. icephs) .or. (ttrice == 0.0_kp))
go to 10
1102 trinv = 1.0_kp/ttrice
1113 lflg = (tc >= -ttrice .and. tc < 0.0_kp)
1114 weight = min(-tc*trinv,1.0_kp)
1115 hlatsb = hlatv + weight*hlatf
1116 hlatvp = hlatv - 2369.0_kp*tc
1117 if (t(i) < tmelt)
then
1123 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) &
1128 desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv
1129 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
1130 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
1132 dqsdt(i) = (cp/hltalt)*gam(i)
1144 hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt)
1146 hlatsb = hlatv + hlatf
1150 if (t(i) < tmelt)
then
1155 desdt = hltalt*es(i)/(rgasv*t(i)*t(i))
1156 gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
1157 if (qs(i) == 1.0_kp) gam(i) = 0.0_kp
1159 dqsdt(i) = (cp/hltalt)*gam(i)
1191 real(kp),
intent(in) :: t, p
1195 real(kp),
intent(out) :: es, qs, dqsdt
1219 omeps = 1.0_kp - epsqs
1227 es = min(fpvs(t), p)
1232 qs = epsqs*es/(p - omeps*es)
1250 if ((.not. icephs) .or. (ttrice == 0.0_kp))
go to 10
1251 trinv = 1.0_kp/ttrice
1263 lflg = (tc >= -ttrice .and. tc < 0.0_kp)
1264 weight = min(-tc*trinv,1.0_kp)
1265 hlatsb = hlatv + weight*hlatf
1266 hlatvp = hlatv - 2369.0_kp*tc
1273 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) &
1278 desdt = hltalt*es/(rgasv*t*t) + tterm*trinv
1279 gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
1280 if (qs == 1.0_kp) gam = 0.0_kp
1282 dqsdt = (cp/hltalt)*gam
1297 hlatvp = hlatv - 2369.0_kp*(t-tmelt)
1299 hlatsb = hlatv + hlatf
1308 desdt = hltalt*es/(rgasv*t*t)
1309 gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
1310 if (qs == 1.0_kp) gam = 0.0_kp
1312 dqsdt = (cp/hltalt)*gam
1362 real(kp),
intent(in) :: t ,tmelt
1366 integer,
intent(inout) :: itype
1368 real(kp),
intent(out) :: es
1396 tr = abs(real(itype,kp))
1403 if (tr > 40.0_kp)
then
1408 if(t < (tmelt - tr) .and. itype == 1)
go to 10
1414 e1 = 11.344_kp*(1.0_kp - t/ts)
1415 e2 = -3.49149_kp*(ts/t - 1.0_kp)
1416 f1 = -7.90298_kp*(ts/t - 1.0_kp)
1417 f2 = 5.02808_kp*log10(ts/t)
1418 f3 = -1.3816_kp*(10.0_kp**e1 - 1.0_kp)/10000000.0_kp
1419 f4 = 8.1328_kp*(10.0_kp**e2 - 1.0_kp)/1000.0_kp
1421 f = f1 + f2 + f3 + f4 + f5
1422 es = (10.0_kp**f)*100.0_kp
1425 if(t >= tmelt .or. itype == 0)
go to 20
1431 term1 = 2.01889049_kp/(t0/t)
1432 term2 = 3.56654_kp*log(t0/t)
1433 term3 = 20.947031_kp*(t0/t)
1434 es = 575.185606e10_kp*exp(-(term1 + term2 + term3))
1436 if (t < (tmelt - tr))
go to 20
1440 weight = min((tmelt - t)/tr,1.0_kp)
1441 es = weight*es + (1.0_kp - weight)*eswtr
1447900
format(
'GFFGCH: FATAL ERROR ******************************',/, &
1448 &
'TRANSITION RANGE FOR WATER TO ICE SATURATION VAPOR',
' PRESSURE, &
1449 & TR, EXCEEDS MAXIMUM ALLOWABLE VALUE OF',
' 40.0 DEGREES C',/, &
1457 real(kp),
intent(in) :: tx
1464 es = exp(54.842763_kp - (6763.22_kp / t) - (4.210_kp * log(t)) + &
1465 & (0.000367_kp * t) + (tanh(0.0415_kp * (t - 218.8_kp)) * &
1466 & (53.878_kp - (1331.22_kp / t) - (9.44523_kp * log(t)) + &
1467 & 0.014025_kp * t)))
1471 function murphykoop_svp_ice(tx)
result(es)
1472 real(kp),
intent(in) :: tx
1480 es = exp(9.550426_kp - (5723.265_kp / t) + (3.53068_kp * &
1481 & log(t)) - (0.00728332_kp * t))
1483 end function murphykoop_svp_ice
1492 real(kp),
intent(in) :: t, p
1496 real(kp),
intent(out) :: es, qs, dqsdt
1503 real(kp) omeps, hltalt, hlatsb, hlatvp, desdt, gam
1507 omeps = 1.0_kp - epsqs
1513 es = min(fpvsi(t),p)
1518 qs = min(1.0_kp, epsqs*es/(p-omeps*es))
1536 hltalt = hlatv + hlatf
1537 desdt = hltalt*es/(rgasv*t*t)
1538 if (qs < 1.0_kp)
then
1539 gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
1544 dqsdt = (cp/hltalt)*gam
subroutine, public aqsat(t, p, es, qs, ii, ilen, kk, kstart, kend)
This subroutine is the utility procedure to look up and returen saturation vapor pressure from precom...
subroutine, public aqsat_water(t, p, es, qs, ii, ilen, kk, kstart, kend)
This subroutine includes the utility procedure to look up and return saturation vapor pressure from p...
subroutine, public aqsatd(t, p, es, qs, gam, ii, ilen, kk, kstart, kend)
This subroutine include the utility procedure to look up and returen saturation vapor pressure from p...
subroutine, public vqsatd(t, p, es, qs, gam, len)
This subroutine is the utility procedure to look up and return saturation vapor pressure from precomp...
subroutine, public vqsatd_water(t, p, es, qs, gam, len)
subroutine, public vqsat_water(t, p, qsat_water, len)
This subroutine.
subroutine, public vqsatd2_water_single(t, p, es, qs, dqsdt)
real(kp) function, public qsat_ice(t, p)
real(kp) function murphykoop_svp_water(tx)
DONIF USe Murphy and Koop (2005) (Written by Andrew Gettelman)
subroutine, public vqsatd2(t, p, es, qs, dqsdt, len)
real(kp) function, public polysvp(t, typ)
subroutine gffgch(t, es, tmelt, itype)
subroutine, public vqsatd2_ice_single(t, p, es, qs, dqsdt)
subroutine, public gestbl(tmn, tmx, trice, ip, epsil, latvap, latice, rh2o, cpair, tmeltx)
This subroutine builds saturation vapor pressure table for later procedure.
subroutine, public vqsatd2_water(t, p, es, qs, dqsdt, len)
subroutine, public vqsatd2_single(t, p, es, qs, dqsdt)
subroutine, public vqsat_ice(t, p, qsat_ice, len)
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...