subroutine cflux(t000,z000d,r000,rlo,rmd,rhi,sst,vmax,pcore, + rlat,sspeed,ibiasc,iyear,ibasin,lulog,flux) c This routine calculates the flux parameter (cflux) in mks units. c The parameter is a difference between fluxes with RH equal to the c surface value and that with RH mixed between a deep tropospheric c layer. c c Input: t000 = 1000 hPa T (deg C) c z000d = 1000 hPa Z deviation (m) c r000 = 1000 hPa RH (%) c rlo = 850 to 700 hPa average RH (%) c rmd = 700 to 500 hPa average RH (%) c rhi = 500 to 300 hPa average RH (%) c sst = sea surface temperature (deg C) c vmax = maximum wind (kt) c pcore = minimum sea level pressure (hPa), c set to 9999. if missing, recalculated from vmax c rlat = latitude (deg N) c sspeed = storm translational speed (kt) c ibiasc = 1 to include RH, T bias corrections, =0 to skip c iyear = 4-digit year (only used for bias corrections) c lulog = Unit number for log file c common /pcons/ rd,rv,eps,rlv,ce,A,B,cctk c c Set idprint=1 to print intermediate variables for debugging idprint = 0 c c *** 1. Specify physical constants rd = 287.0 rv = 462.0 eps = rd/rv rlv = 2.4417e+6 c ce = 0.0012 ce = 0.00118 A = 17.27 B = 237.7 cctk = 273.15 cktm = 0.5144 gammad = 0.0098 c c *** 2. Estimate minimum sea level pressure (hPa) from pressure/wind relationship c if pcore is not provided. if (pcore .gt. 1050.0) then pcore = 1015.0 - (100.0/135.0)*(vmax-25.0) endif c c *** 3. Find standard height of 1000 hPa level and add it to z000d p = 1000.0 call stndz(p,zstd,tstd,thstd) z000 = zstd+z000d c c *** 4. Apply bias corrections c if (ibiasc .eq. 1) then c Apply those that depend on the year if (iyear .le. 1995) t000 = t000 - 2.2281 if (iyear .le. 1996) r000 = r000 - 5.0962 if (iyear .le. 1994) rlo = rlo + 9.6483 if (iyear .le. 1994) rmd = rmd +10.5621 endif c c Calculate sfc RH rsfc = r000 + 4.1022 c c Check to make sure RH never exceeds 100% if (rsfc .gt. 100.) r000 = 100.0 if (rlo .gt. 100.) rlo = 100.0 if (rmd .gt. 100.) rmd = 100.0 c c *** 5a,b. Find z=10m T and Td in the environment tsfce = t000 + gammad*(z000-10.0) call ttotd(tsfce,rsfc,tdec,tdek) c c *** 6. Calculate sfc T in the core assuming Td in the core c equals Td in the environment and RH=95% in the core rhcore = 95.0 call tdtot(tdec,rhcore,tcorec,tcorek) c c *** 7. Calculate specific humidity in the core in g/g (qcore) call esatcal(tcorec,escore) wcoresat = escore*eps/(pcore-escore) wcore = (rhcore/100.0)*wcoresat qcore = wcore/(1.0+wcore) c c *** 8. Adjust SST based on Cione SST cooling equation if (ibasin .eq. 1) then c = cktm*sspeed ssta = 1.12228*sst + 0.14256*c - 0.077859*rlat - 3.37056 else ssta = sst endif c c *** 9. Calculate saturated specific humidity at T=ssta call esatcal(ssta,esst) wsst = esst*eps/(pcore-esst) qsst = wsst/(1.0+wsst) c c *** 10. Vertically average RH c rhavg = (rhcore + rlo + rmd + rhi)/4.0 rhavg = (r000 + rlo + rmd)/3.0 c c *** 11. Repeat steps 5b through 7 with rhavg c c ++11.5b call ttotd(tsfce,rhavg,tdeca,tdeka) c c ++ 11.6. Calculate sfc T in the core assuming Td in the core c equals Td in the environment and RH=95% in the core rhcore = 95.0 call tdtot(tdeca,rhcore,tcoreca,tcoreka) c c ++ 11.7. Calculate specific humidity in the core in g/g (qcore) call esatcal(tcoreca,escorea) wcoresata = escorea*eps/(pcore-escorea) wcorea = (rhcore/100.0)*wcoresata qcorea = wcorea/(1.0+wcorea) c c *** 12. Calculate sfc density for rsfc and adjusted rsfc (=ravg) call rhocal(tcorec ,pcore,wcore ,rhocore) call rhocal(tcoreca,pcore,wcorea,rhocorea) c c *** 13. Calculate flux variable term1 = rlv*ce*vmax*cktm term2 = rhocore *(qsst-qcore) term3 = rhocorea*(qsst-qcorea) flux = term1*(term3-term2) c c *** 14. Debug print statements if (idprint .eq. 1) then write(lulog,400) zstd,z000 400 format('std z1000: ',f6.1,' z1000: ',f6.1) c write(lulog,401) pcore,rsfc 401 format('pcore: ',f6.1,/, + 'rhsfc: ',f6.1) write(lulog,402) tsfce,tdec,tcorec 402 format('tsfc: ',f6.2,' tdsfc: ',f6.2,' tcore: ',f6.2) c write(lulog,410) sst,ssta 410 format('SST: ',f6.2,' SSTa: ',f6.2) c write(lulog,415) qcore*1000.0,qsst*1000.0 415 format('qcore: ',f6.2,' qsst: ',f6.2) c write(lulog,421) rhavg 421 format('rhavg: ',f6.1) c write(lulog,422) tsfce,tdeca,tcoreca 422 format('tsfc: ',f6.2,' tdsfca:',f6.2,' tcoreca:',f6.2) c write(lulog,425) 1000.0*qcorea,1000.0*qsst 425 format('qcorea: ',f6.2,' qsst: 'f6.2) c write(lulog,430) rhocore,rhocorea,flux 430 format('rhocore: ',f6.3,' rhocorea: ',f6.3,/, + 'flux: ',e11.4) endif c return end subroutine ttotd(tc,rh,tdc,tdk) c This routine calculates the dewpoint temperature in deg c (tdc) c and deg K (tdk) given the temperature t in deg c (tc) and RH in % c common /pcons/ rd,rv,eps,rlv,ce,A,B,cctk c if (rh .lt. 1.0) rh = 1.0 c alpha = A*tc/(B+tc) + alog(rh/100.0) tdc = B*alpha/(A-alpha) tdk = tdc + 273.15 c return end subroutine tdtot(tdc,rh,tc,tk) c This routine calculates the temperature in deg c (tc) and c deg K (tk) given the dewpoint temperature in deg c (tdc) c and the RH in % c common /pcons/ rd,rv,eps,rlv,ce,A,B,cctk c if (rh .lt. 1.0) rh = 1.0 c alpha = A*tdc/(B+tdc) tc = B*(alpha - alog(rh/100.0)) tc = tc/(A-alpha+alog(rh/100.0)) tk = tc + cctk c return end subroutine esatcal(t,es) c This routine calculates the saturation vapor pressure in hPa (es) c given the temperature in deg C (t) c es = 6.112*exp(17.67*t/(t+243.5)) c return end subroutine rhocal(T,P,w,rho) c This routine calculates the density (rho) in kg/m3. The virtual temperature c correction is applied. c c Input: T (temperature in deg C) c P (total pressure in hPa) c w (mixing ratio g/g) c common /pcons/ rd,rv,eps,rlv,ce,A,B,cctk c c Convert T to Kelvin tk = T + cctk c c Convert P to Pa ppa = 100.0*p c c Virtual T tvk = tk*(1.0 + w/eps)/(1 + w) c c Density rho = ppa/(rd*tvk) c return end