c....................................................................... c PROGRAM: bulk_v2.5b c COARE bulk flux version 2.5 c c 1. Based on Liu et al. 1979 and Liu's original computer code. c c 2. Adopt blending between Kansas and Free-convection forms by Fairall, c Rogers and Bradley October 20 1993 (ref. Fairall et al. JGR, 101, 3747-3764, 1996. c c 3. Modified by Fairall and Bradley June 22 1994 to include c a) cool skin effect (needs radiative fluxes) c equations quoted are from Fairall et al. 1995b (JGR, 1995, to appear) c b) sensible heat flux (ocean cooling) due to rain at wet bulb temp. c formalism by Gosnel, Webster and Fairall (JGR, 1995, to appear) c c) a simplified version of Price, Weller and Pinkel (PWP) model for solar c warm layer c This option requires a long data file in time order, with resolution at c least 1 hour, including radiative fluxes c c 4. Version 2f includes minor modifications and corrections to warm layer and c cool skin effect. (2/95) c c 5. Version 2.5 further code modifications (changes to Tau to output in kinematic units c which is required by the warm layer calculation), 6% decrease of scalar transfer c coefficients by adusting roughness Reynolds number coefficients, c additional outputs: (8/95) c c The Gustiness factor, wg, the Webb correction to the latent heat flux c and the rain momentum flux are now output. c c The form of the transfer coefficients has been changed c to correspond directly with the form of the equations summarized below. c c Latitude and Longitude is now a required input. c c Note the wind speed should be input relative to the sea surface. c c.................................................................... c Summary and form of the basic flux equations (Fairall et al. 1995, to be published c in JGR) c c Sensible heat flux: c c QS=-rhoa*Cpa*Ustar*Tstar = rhoa*Cpa*Ch*S*(Ts-Pt) c c Latent heat flux: c c QE=-rhoa*Le*Ustar*Qstar = rhoa*Le*Ce*S*(qs-q) c c Wind stress c c Tau_i=-rhoa*Ustar**2 = rhoa*Cd*S*(us_i-u_i) c c where c Ch, Ce, and Cd are the transfer coefficients for sensible c heat, latent heat and stress, Pt is the potential temperature, c q is the water vapor mixing ratio, and u_i is one of the horizontal c wind components relative to the fixed earth measured at some c atmospheric reference height, z. S is the average value of the wind speed c relative to the sea surface, which includes a gustiness factor, wg, such c that c c S=SQR(u**2 + v**2 + wg**2). c c Ts is sea surface interface temperature (skin temperature), qs is the c interfacial value of the saturation mixing ratio, us_i is the surface current. c Note that the code requires the wind components relative to the sea surface (i.e., us_i=0). c c c.................................................................... c input: c intime (time, either yymnddhhmmss.ssss or dd.dddd) c hUm (wind measurement height) m c hTm (T&rh measurement height) m c hUs (wind standard height) m - use this to change your mean data to different ref. Height. c hTs (T&rh standard height) m c ts_depth (depth of sst instrument - a positive value) m - for cool skin/warm layer correction c ws (wind speed relative to the sea surface) m/s c sst (sea surface temp.) deg. C c atb (air temperature) deg. C c qq (specific humidity) g/kg or qq (RH as decimal) - Note, code works internally in kg/kg!! c rs (shortwave radiation) W/m2 c rl (downwelling longwave) W/m2 c rain (average rainrate) mm/hour c pp (pressure) mb c zi (boundary-layer depth; use 600m if data unavailable) m c Jcool (=1 for cool skin calculation; =0 if SST measured by IR radiometer) c Jwarm (=1 for warm layer calculation; =0 if SST measured by IR radiometer) c intime (time variable) togatime in yymnddhhmmss.ssss - GMT c glat (latitude) degrees [latitude north +ve, latitude south -ve] c glon (longitude) degrees [Longitude east +ve, longitude west -ve] c.................................................................... c c output: c QF W/m**2 (Sensible heat flux) - turbulent part only c QE W/m**2 (Latent heat flux) - turbulent part only c RF W/m**2 (Rainfall heat flux) c TAU N /m**2 (wind stress) - turbulent part only c Ustar m/s (velocity scaling parameter - friction velocity) c Qstar kg/kg (humidity scaling parameter) c Tstar C (temperature scaling parameter) c CD - drag coefficient c CH - transfer coefficient for heat c CE - transfer coefficient for moisture c RR - Roughness Reynolds number c RT - Roughness Reynolds number for temperature c RQ - Roughness Reynolds number for moisture c ZL - ht/L where L is the Obukhov length c Z0 - roughness length c zot - roughness length for temperature c zoq - roughness length for humidity c dt_wrm - total warm layer temperature difference C c dter - cool skin temperature difference C c T0 - skin temperature C (T0 = sst - dter + dt_wrm*ts_depth/zwarm) c wg - gustiness factor m/s c Taur - momentum flux due to rain N/m**2 c Ef_webb - Webb correction to the latent heat flux W/m**2 c-------------------------------------------------------------- c c subroutine bulk_mat_v25b(m,intime, & hUm,hTm,hUs,hTs,ts_depth, & ws,sst,atb,qq,ws_h,Ta_h,qq_h, & rs,rl,rain,pp,zi,Jcool,Jwarm,glat,glon,QH,QE,RF,TAU, & Ustar,Tstar,Qstar,CD,CH,CE,RR,RT,RQ, & ZL,ZO,zot,zoq,dt_wrm,dter,T0,wg,Tau_r,Ef_webb) c-------------------------------------------------------------- real*8 hUm(m),hTm(m),hUs(m),hTs(m),ws_h(m),Ta_h(m),qq_h(m) real*8 ws(m),sst(m),atb(m),qq(m),pp(m),zi(m),rain(m) real*8 QH(m),QE(m),TAU(m),Ustar(m),Qstar(m),Tstar(m) real*8 rl(m),rs(m),RF(m),T0(m),ts_depth(m),wg(m) real*8 CD(m),CE(m),CH(m),RR(m),RT(m),RQ(m),Zl(m),ZO(m) real*8 zot(m),zoq(m),dt_wrm(m),dter(m) real*8 glat(m),glon(m),Ef_webb(m),Tau_r(m) real*8 intime(m),sol_time,tim c real*8 time_old,qcol_ac,tau_ac,tau_old,rf_old,hf_old,ef_old real*8 tk_pwp,fxp integer hh,mn,yy,mm,dd,jamset,jump,m character*17 chtime COMMON /old/time_old,qcol_ac,tau_ac,tau_old,rf_old,hf_old, & ef_old,jamset,jump,fxp,tk_pwp c initialize variables qcol_ac=0. tau_ac=0. time_old=0. jamset=0. tau_old=0. hf_old=0. ef_old=0. rf_old=0. jamset=0 jump=0 c c initial guesses to the warm layer parameters, these values simulate c what is expected in an early morning situation: fxp=0.5 implies a c shallow heating layer to start the integration; tk_pwp=19.0 implies c the thickness is a maximum from the day before and is not meant to c match this timestep's fxp. c fxp=0.5 tk_pwp=19.0 c c main loop of the flux routine c c c loop through data do 30 i=1,m c--------------------------------------------------------------------- c check warm layer, cool skin switches c if(Jwarm.gt.0) then if(i.eq.1) then Jwarm=2 else Jwarm=1 endif endif c c convert time to fractions of day (e.g. yymnddhhmmss.ssss -> dd.dddd) c skip this part if intime is in the right format c if(intime(i) .gt. 1e9) then write(chtime(1:17),'(f17.4)') intime(i) read(chtime,12) yy,mn,dd,hh,mm,ss 12 format(5i2,f7.4) tim=(float(hh)+float(mm)/60.+ss/3600.) endif c c and convert to local solar time in seconds c sol_time=mod(glon(i)/15.+tim+24.,24.)*3600 c c c call bulk flux routine c call bulk_flux(sol_time,hUm(i),hTm(i),hUs(i),hTs(i), & ws(i),sst(i),atb(i),qq(i),ws_h(i),Ta_h(i),qq_h(i), & rs(i),rl(i),rain(i),pp(i),zi(i),Jcool,Jwarm,QH(i), & QE(i),RF(i),TAU(i),Ustar(i),Tstar(i),Qstar(i), & CD(i),CH(i),CE(i),RR(i),RT(i),RQ(i),ZL(i),ZO(i), & zot(i),zoq(i),dt_wrm(i),dter(i), & T0(i),ts_depth(i),wg(i),TAU_r(i),EF_webb(i)) c c 30 continue c c end of loop c return end c c***************************************************************** subroutine bulk_flux(sol_time,hUm,hTm,hUs,hTs, & ws,sst,atb,qq,ws_h,Ta_h,qq_h,rs,rl,rainx,pp,zix,Jcoolx,Jwarmx, & HF,EF,RF,TAU,Ustar,Tstar,Qstar, & CD,CH,CE,RRx,RTx,RQx,ZLx,ZOx,zotx,zoqx, & dt_wrmx,dterx,T0,ts_depthx,wgx,TAU_r,Hl_webb) c real*8 hUm,hTm,hUs,hTs,ws_h,Ta_h,qq_h,ts_depthx real*8 ws,sst,atb,qq,pp,zix,rainx real*8 HF,EF,TAU,Ustar,Qstar,Tstar,TAU_r,Hl_webb real*8 rl,rs,RF,T0,sol_time,wgx,ee real*8 CD,CE,CH,RRx,RTx,RQx,Zlx,ZOx real*8 zotx,zoqx,dt_wrmx,dterx,qjoule,qr_out,dtime real*8 ZUs,ZTs,ZQs,U_hs,T_hs,Q_hs,ctd1,ctd2,rich real*8 dqs_dt,dwat,dtmp c real*8 U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,QA real*8 USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt real*8 al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,visa real*8 visw,von,fdg,DU_Wg,Wg,newtime real*8 time_old,qcol_ac,tau_ac,tau_old,rf_old,hf_old,ef_old real*8 fxp,tk_pwp,dsea integer ID,jump,jamset,ihumid COMMON /old/time_old,qcol_ac,tau_ac,tau_old,rf_old,hf_old, & ef_old,jamset,jump,fxp,tk_pwp COMMON/PIN/U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,ID COMMON/POUT/USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt COMMON/const/al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK, & visa,visw,von,fdg COMMON/wgust/DU_Wg,Wg c c initialize variables that appear in COMMON c Jcool=Jcoolx Jwarm=Jwarmx ZU=hUm !height of wind measurement ZT=hTm !height of temperature measurement ZQ=hTm !height of water vapor measurement ZUs=hUs !standard height of wind measurement ZTs=hTs !standard height of temperature measurement ZQs=hTs !standard height of water vapor measurement U=ws !wind speed m/s TS=sst !surface temp. Celsius T=atb !air temp. Celsius P=pp !pressure mb zi=zix !atmospheric boundary layer depth toK=273.16 Rnl= 0.97*(5.67e-8*(TS+toK)**4-rl) !Net longwave (up = +) Rns=0.945*rs !Net shortwave (into water) rain=rainx !rainfall ts_depth=ts_depthx !depth of sst measurement c c Warm Layer c if(Jwarm.ne.0) then newtime=sol_time !run time in secs if(Jwarm.eq.2) then jump=1 goto 16 !set time_old and pass thru' ASL elseif(newtime.gt.21600.and.jump.eq.1) then goto 16 !6 am too late to start elseif(newtime.lt.time_old) then !reset all var. at midnight jump=0 !test threshold q morning only jamset=0 fxp=0.5 tk_pwp=19.0 tau_ac=0.0 qcol_ac=0.0 dt_wrm=0.0 goto 16 else rich=.65 !critical Rich. No. ctd1=sqrt(2*rich*cpw/(al*grav*rhow)) !u*^2 integrated so ctd2=sqrt(2*al*grav/(rich*rhow))/(cpw**1.5) !has /rhow in both dtime=newtime-time_old !delta time qr_out=rnl+hf_old+ef_old+rf !flux out from previous pass q_pwp=fxp*rns-qr_out !effective net warming if(q_pwp.lt.50.and.jamset.eq.0) go to 16 !integration threshold jamset=1 tau_ac=tau_ac+tau_old*dtime !tau from previous pass c if(qcol_ac+q_pwp*dtime.gt.0) then do 10 index=1,5 !iterate for warm layer thickness fxp=1.-(0.28*0.014*(1-dexp(-tk_pwp/0.014)) ! & +0.27*0.357*(1-dexp(-tk_pwp/0.357)) ! & +.45*12.82*(1-dexp(-tk_pwp/12.82)))/tk_pwp !solar absorb. prof qjoule=(fxp*rns-qr_out)*dtime if((qcol_ac+qjoule.gt.0.0)) & tk_pwp=min(19.,ctd1*tau_ac/sqrt(qcol_ac+qjoule)) 10 continue else fxp=.76 tk_pwp=19 qjoule=(fxp*rns-qr_out)*dtime endif c qcol_ac=qcol_ac+qjoule !integrate heat input if(qcol_ac.gt.0) then dt_wrm=ctd2*(qcol_ac)**1.5/tau_ac !pwp model warming else dt_wrm=0. endif endif if(tk_pwp.lt.ts_depth) then !pwp layer deeper than sensor dsea=dt_wrm else dsea=dt_wrm*ts_depth/tk_pwp !linear temperature profile endif ts=ts+dsea 16 time_old=newtime endif c c dt_wrm=dsea c c end of warm layer c 15 call humidity(T,P,QA) !Teten's formula returns sat. air in mb if(qq.lt.2.) then !checks whether humidity in g/Kg or RH R=qq ee=QA*R !convert from RH using vapour pressure Q=.62197*(ee/(P-0.378*ee)) !Spec. humidity kg/kg else Q=qq/1000. !g/kg to kg/kg endif QA=.62197*(QA/(P-0.378*QA)) !convert from mb to spec. humidity kg/kg call humidity(TS,P,QS) !sea QS returned in mb QS=QS*0.98 !reduced for salinity Kraus 1972 p. 46 QS=.62197*(QS/(P-0.378*QS)) !convert from mb to spec. humidity kg/kg c c calculate atmospheric surface layer c call ASL(Jcool,IER) if(IER.ge.0) then c c compute surface stress (TAU), sensible heat flux (HF), c latent heat flux (EF) & other parameters c S=sqrt(u*u + wg*wg) !velocity incl. gustiness param. TAU=rhoa*USR*usr*u/S !kinematic units HF=-cpa*rhoa*USR*TSR EF=-xlv*rhoa*USR*QSR c c compute heat flux due to rainfall c dwat=2.11e-5*((T+toK)/toK)**1.94 !water vapour diffusivity dtmp=(1.+3.309e-3*T-1.44e-6*T*T)*0.02411/(rhoa*cpa) !heat diffusivity dqs_dt=QA*xlv/(rgas*(T+toK)**2) !Clausius-Clapeyron alfac= 1/(1+0.622*(dqs_dt*xlv*dwat)/(cpa*dtmp)) !wet bulb factor W=-1.61*USR*QSR-(1+1.61*Q)*USR*TSR/(T+toK) rF= rain*alfac*cpw*((TS-T)+(QS-Q)*xlv/cpa)/3600. c c Webb correction to latent heat flux c Hl_webb=rhoa*xlv*W*Q c c compute momentum flux due to rainfall c TAU_r=0.85*rain/3600*u c dterx=dter !cool skin parameters tktx=tkt ! T0=ts-dter tau_old=tau ef_old=ef hf_old=hf dt_wrmx=dt_wrm !warm layer parameter c c compute transfer coefficients c CD=(USR/S)**2 CH=USR*TSR/(S*(T-TS+.0098*zt+dter)) !revised 2e to 2f !included '+dter' CE=USR*QSR/(S*(Q-QS+dqer)) Ustar=USR Tstar=TSR Qstar=QSR RRx=RR RTx=RT RQx=RQ ZLx=ZL ZOx=ZO zotx=zot zoqx=zoq wgx=wg ihumid=0 if(qq .lt. 2) ihumid=1 call h_adjust(ZUs,ZTs,ZQs,U_hs,T_hs,Q_hs,ihumid) ws_h=U_hs Ta_h=T_hs qq_h=Q_hs else !input parameters out of range EF=-999. HF=-999. TAU=-999. TAUr=-999. EF_webb=-999. Ustar=-999. Tstar=-999. Qstar=-999. RRx=-999. RTx=-999. RQx=-999. ZLx=-999. ZOx=-999. ws_h=-999. Ta_h=-999. qq_h=-999. wg=-999. endif return !return to main program end c c ------------------------------------------------------------------ subroutine ASL(Jcool,IER) c c TO EVALUATE SURFACE FLUXES, SURFACE ROUGHNESS AND STABILITY OF c THE ATMOSPHERIC SURFACE LAYER FROM BULK PARAMETERS BASED ON c LIU ET AL. (79) JAS 36 1722-1735 c real*8 U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,D,S real*8 USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt real*8 al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,visa real*8 visw,von,fdg,DU_Wg,Wg real*8 PTZ,PQZ,PUZ,ZTl,ZLN,ZQl,wetc,bigc,cpv,dq,dqq,dt,dtt integer ID COMMON/PIN/U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,ID COMMON/POUT/USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt COMMON/const/al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK, & visa,visw,von,fdg COMMON/wgust/DU_Wg,Wg c Factors Beta=1.2 ! evaluated from Fairall's low windspeed turbulence data Von=0.4 ! von Karman's "constant" c fdg=1.00 ! Fairall's LKB rr to von karman adjustment fdg=1.00 !based on results from Flux workshop August 1995 toK=273.16 ! Celsius to Kelvin grav=9.72 ! gravity equatorial value (ref. IGPP-SIO) c Air constants and coefficients Rgas=287.1 !J/kg/K gas const. dry air xlv=(2.501-0.00237*TS)*1e+6 !J/kg latent heat of vaporization at TS Cpa=1004.67 !J/kg/K specific heat of dry air (Businger 1982) Cpv=Cpa*(1+0.84*Q) !Moist air - currently not used (Businger 1982) rhoa=P*100./(Rgas*(T+toK)*(1.+.61*Q)) !kg/m3 Moist air density ( " ) visa=1.326e-5*(1+6.542e-3*T+8.301e-6*T*T-4.84e-9*T*T*T) !m2/s !Kinematic viscosity of dry air - Andreas (1989) CRREL Rep. 89-11 c Cool skin constants al=2.1e-5*(ts+3.2)**0.79 !water thermal expansion coefft. be=0.026 !salinity expansion coefft. cpw=4000. !J/kg/K specific heat water rhow=1022. !kg/m3 density water visw=1.e-6 !m2/s kinematic viscosity water tcw=0.6 !W/m/K Thermal conductivity water bigc=16.*grav*cpw*(rhow*visw)**3/(tcw*tcw*rhoa*rhoa) wetc=.622*xlv*QS/(rgas*(TS+toK)**2) !correction for dq;slope of sat. vap. c Initialise everything IER=0 ZL=0. Dter=0. !cool skin Dt Dqer=0. !cool skin Dq c Initial guesses c next line originally required for surface current, input now requires relative wind c US=0. !surface current = 0. c Wg=0.5 !Gustiness factor initial guess ZO=.0005 !roughness initial guess tkt=.001 !guess sublayer thickness DU=U !assumes U is measured rel. to current DU_Wg=(DU**2.+Wg**2.)**.5 !include gustiness in wind spd. difference !equivalent to S in definition of fluxes DT=T-TS+.0098*zt !potential temperature diff DQ=Q-QS USR=.04*DU_Wg TSR=.04*DT QSR=.04*DQ if(DU_Wg.ne.0.)then TA=T+toK RI=grav*ZU*(DT+0.61*TA*DQ)/(TA*DU_Wg**2) else IER=-1 RI=-999. endif if(RI.gt.0.25) IER=-1 do 200 index=1,20 !iterate 20 times call ZETA(T,Q,USR,TSR,QSR,ZU,ZLN) ZL=ZLN PUZ=PSI(1,ZL) ZTL=ZL*ZT/ZU ZQL=ZL*ZQ/ZU PTZ=PSI(2,ZTL) PQZ=PSI(2,ZQL) ZO=0.011*USR*USR/grav + 0.11*visa/USR !after Smith 1988 USR=DU_Wg*von/(dlog(ZU/ZO)-PUZ) !Gustiness effect incl. RR=ZO*USR/VISA call LKB(RR,RT,1) if(RT.eq.-999.) then IER = -2 !error - return return endif 21 call LKB(RR,RQ,2) if(RQ.eq.-999.) then IER = -2 !error - return return endif 22 zot=rt*visa/usr zoq=rq*visa/usr S = (dlog(ZT/zot)-PTZ)/(von*fdg) !coeff fdg=1.04 included following D = (dlog(ZQ/zoq)-PQZ)/(von*fdg) !Fairall observations during COARE. !NOTE coeff changed to 1. !This has been adjusted to 0.94 !following discussion during FLux !workshop August 1995. dtt=(dt+dter) dqq=(dq+dqer) ! or we could calculate new sat. hum. tsr=dtt/S !! modify qsr=dqq/D !! fluxes TVSR=TSR*(1.+0.61*Q)+(0.61*TA*QSR) Bf=-grav/TA*USR*TVSR if(Bf.gt.0) then Wg=Beta*(Bf*zi)**0.333 else Wg=0. endif DU_Wg=(DU**2.+Wg**2.)**.5 !include gustiness in wind spd. if(Jcool.ne.0) then !Cool skin hsb=-rhoa*cpa*usr*tsr hlb=-rhoa*xlv*usr*qsr qout=rnl+hsb+hlb dels=rns*(.137+11*tkt-6.6e-5/tkt*(1-dexp(-tkt/8.0e-4))) ! Eq.16 ! Shortwave, correction ! from version 2e to 2f qcol=qout-dels if(qcol.gt.0.) then alq=Al*qcol+be*hlb*cpw/xlv ! Eq. 7 Buoy flux water xlamx=6/(1+(bigc*alq/usr**4)**.75)**.333 ! Eq 13 Saunders coeff. tkt=xlamx*visw/(sqrt(rhoa/rhow)*usr) ! Eq.11 Sublayer thickness dter=qcol*tkt/tcw ! Eq.12 Cool skin else dter=0. endif dqer=wetc*dter endif !end cool skin 200 continue !end iterations return !to main subroutine, bulk_flux end c c------------------------------------------------------------------ subroutine humidity(T,P,Qsat) c c Tetens' formula for saturation vp Buck(1981) JAM 20, 1527-1532 c real*8 T,P,Qsat c Qsat = (1.0007+3.46e-6*P)*6.1121*dexp(17.502*T/(240.97+T)) return end c c------------------------------------------------------------------ subroutine LKB(RR,RT,IFLAG) c c TO DETERMINE THE LOWER BOUNDARY VALUE RT OF THE c LOGARITHMIC PROFILES OF TEMPERATURE (IFLAG=1) c OR HUMIDITY (IFLAG=2) IN THE ATMOSPHERE FROM ROUGHNESS c REYNOLD NUMBER RR BETWEEN 0 AND 1000. OUT OF RANGE c RR INDICATED BY RT=-999. BASED ON LIU ET AL.(1979) c JAS 36 1722-1723 c------------------------------------------------------------------- c Scalar RR relation from Moana Wave data. c c real*8 A(9,2),B(9,2),RAN(9),RR,RT c integer iflag c DATA A/0.177,2.7e3,1.03,1.026,1.625,4.661,34.904,1667.19,5.88E5, c & 0.292,3.7e3,1.4,1.393,1.956,4.994,30.709,1448.68,2.98E5/ c DATA B/0.,4.28,0,-0.599,-1.018,-1.475,-2.067,-2.907,-3.935, c & 0.,4.28,0,-0.528,-0.870,-1.297,-1.845,-2.682,-3.616/ c DATA RAN/0.11,.16,1.00,3.0,10.0,30.0,100.,300.,1000./ c------------------------------------------------------------------- c c Scalar RR relation from Liu et al. c real*8 a(8,2),b(8,2),ran(8),RR,RT integer iflag c data a/0.177,1.376,1.026,1.625,4.661,34.904,1667.19,5.88e5, & 0.292,1.808,1.393,1.956,4.994,30.709,1448.68,2.98e5/ data b/0.,0.929,-0.599,-1.018,-1.475,-2.067,-2.907,-3.935, & 0.,0.826,-0.528,-0.870,-1.297,-1.845,-2.682,-3.616/ data ran/0.11,0.825,3.0,10.0,30.0,100.,300.,1000./ c I=1 if((RR.le.0.).or.(RR.ge.1000.)) goto 90 10 continue if(RR.le.RAN(I)) goto 20 I=I+1 goto 10 20 RT=A(I,IFLAG)*RR**B(I,IFLAG) goto 99 90 RT=-999. 99 return end c c----------------------------------------------------------------- function PSI(ID,ZL) C C TO EVALUATE THE STABILITY FUNCTION PSI FOR WIND SPEED (IFLAG=1) C OR FOR TEMPERATURE AND HUMIDITY PROFILES FROM STABILITY C PARAMETER ZL. SEE LIU ET AL (1979). c Modified to include convective form following Fairall (Unpublished) C real*8 Zl,chik,psik,f,chic,psic integer ID c if(ZL)10,20,30 10 F=1./(1+zl*zl) CHIK=(1.-16.*ZL)**0.25 if(ID.eq.1) goto 11 PSIK=2.*dlog((1.+CHIK*CHIK)/2.) goto 12 11 PSIK=2.*dlog((1.+CHIK)/2.)+dlog((1.+CHIK*CHIK)/2.) & -2.*ATAN(CHIK)+2.*ATAN(1.) 12 CHIC=(1.-12.87*ZL)**.333 !for very unstable conditions PSIC=1.5*dlog((CHIC*CHIC+CHIC+1.)/3.) & -(3.**.5)*ATAN((2*CHIC+1.)/(3.**.5)) & +4.*ATAN(1.)/(3.**0.5) c c match Kansas and free-conv. forms with weighting F c PSI= F*PSIK+(1-F)*PSIC goto 99 20 PSI=0. goto 99 30 continue PSI=-4.7*ZL 99 return end c c------------------------------------------------------------ subroutine ZETA(T,Q,USR,TSR,QSR,Z,ZL) C C TO EVALUATE OBUKHOVS STABILITY PARAMETER Z/L FROM AVERAGE C TEMP T IN DEG C, AVERAGE HUMIDITY Q IN GM/GM, HEIGHT IN M, C AND FRICTIONAL VEL,TEMP.,HUM. IN MKS UNITS C SEE LIU ET AL. (1979) C real*8 T,Q,OB,TVSR,TV,TA,sgn real*8 USR,TSR,QSR,Z,ZL real*8 al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,visa real*8 visw,von,fdg COMMON/const/al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK, & visa,visw,von,fdg c TA=T+toK TV=TA*(1.+0.61*Q) TVSR=TSR*(1.+0.61*Q)+0.61*TA*QSR sgn=sign(1.,tvsr) !added this to avoid program if(abs(tvsr) .lt. 1.e-3) then !failure when TVSR is very small tvsr=sgn*tvsr endif OB=TV*USR*USR/(grav*VON*TVSR) ZL=Z/OB c if(ZL .gt. 1000) ZL=1000. goto 99 c 10 ZL=0. 99 return end c------------------------------------------------------------------------- subroutine H_ADJUST(ZUs,ZTs,ZQs,U_hs,T_hs,Q_hs,IHUMID) c c This subroutine adjusts the U,T,Q variables to the specified c standard height (ZUs,ZTs,ZQs) using the loglayer profiles. c The DELTA correction (adjustment) is relative to the surface c measurement. Cronin 4/13/94 c real*8 ZUs,ZTs,ZQs,U_hs,T_hs,Q_hs,ZUsL,ZTsL,ZQsL real*8 U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,PUZs,PTZs,PQZs real*8 U_wg_hs,Rho_hs,Rho_avg,QA,Rho,P_hs,ee real*8 USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt real*8 al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,visa real*8 visw,von,fdg,DU_Wg,Wg,S,D integer ID,ihumid c COMMON/PIN/U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,ID COMMON/POUT/USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt COMMON/const/al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK, & visa,visw,von,fdg COMMON/wgust/DU_Wg,Wg call ZETA(T,Q,USR,TSR,QSR,ZUs,ZUsL) call ZETA(T,Q,USR,TSR,QSR,ZTs,ZTsL) call ZETA(T,Q,USR,TSR,QSR,ZQs,ZQsL) PUZs= PSI(1,ZUsL) PTZs= PSI(2,ZTsL) PQZs= PSI(2,ZQsL) S = (dlog(ZTs*USR/VISA/RT)-PTZs)/(von*fdg) D = (dlog(ZQs*USR/VISA/RQ)-PQZs)/(von*fdg) T_hs =TSR*S +TS - dter -.0098*ZTs Q_hs =(QSR*D + QS - dqer)*1000 U_wg_hs = USR*(dlog(ZUs/ZO) - PUZs)/0.4 if(U_wg_hs.ge.Wg) then U_hs = SQRT(U_wg_hs**2 - Wg**2) else U_hs = U_wg_hs endif c if(IHUMID.eq.1) then ! then need to convert sp hum into rh Q_hs = Q_hs/1000 ! sh kg/kg RHO=1./(287.*(T+273.16)*(1.+.61*Q))*P*100. P_hs = P - (RHO*grav*(ZTs - ZT))/100 !Approx hydrost.Pressure mb RHO_hs=1./(287.*(T_hs+273.16)*(1.+.61*Q_hs))*P_hs*100 RHO_avg = (RHO + RHO_hs)/2 P_hs = P -(RHO_avg*grav*(ZTs - ZT))/100 !hydrostatic Pressure call humidity(T_hs,P_hs,QA) !Teten's formula for Pvap,sat ee=Q_hs*P_hs/(.62197 + .378*Q_hs) !to get vapor pressure Q_hs = ee/QA !to get relative humidity endif return end