program SIMFLX c c This programs estimates the surface fluxes of heat c moisture and momentum using the similsrity theory. c surface fluxes calculations require too many surface c parameters which need to be defined at each grid point. c because of the voluminous data required ,this program c is designed to run over one grid point only .It is a c column model for surface flux calculations. c Althought the surface fluxes are calculated using c only the anemometer level and surface variables, the c code presented here is set closely compatible to that c which would normally run in a Numerical Weather Prediction c model, where data comes into a vertical array on sigma c levels for each grid point (sigma = p/ps) where p is the c pressure at the given level and ps is the surface pressure. c parameter ( nk = 12 ) real t(nk),q(nk),sh(nk),s(nk) common/wet/eps,ae(2),be(2),ht(2),tc common/sfc/rib,ustr,ol,cd,ch,cq,z2,z1,sigz2,zv,zt c c Parameters definition . c c parsol : Ice absorption coefficient (over ocean) c cfonte : Snow melting coefficient (over land ) c humsol : Soil moisture availability c albedo : Surface albedo c flxiri : Outgoing longwave radiation at surface c flxvis : Incoming Shortwave radiation at surface c z0 : Roughness length c vent : Wind speed at lowest model level c ts : Surface temperature (sst over oceans) c ps : Surface pressure c c Land Sea matrix c c ntypes = -1 .... Land c ntypes = 0 .... Ocean c ntypes = 1 .... Land with snow c ntypes = 2 .... Ocean with ice c c Define the sigma levels c data s /.1,.2,.3,.4,.5,.6,.7,.8,.85,.9,.95,.99/ c c Define temperature and specific humidity profile for c the given grid point. c data t/220.16,218.06,223.71,235.78,247.31,256.73, & 264.38,269.16,271.98,274.71,277.10,278.41/ c data q/ 4.0412000e-05,1.7426601e-05,2.6705300e-05 & ,9.5105803e-05,2.6963701e-04,6.8341498e-04 & ,1.6110500e-03,2.1570900e-03,2.3708199e-03 & ,2.5780201e-03,2.9228700e-03,3.2207901e-03/ c c Define surface parameters. c data parsol,z0,ntypes,humsol/ 0.45,1.26,-1,0.15/ data ts,ps,albedo /278.96,88230.5,8.05e-02/ data vent,flxvis,flxiri /9.43,567.93,178.22/ c c Compute intermediate sigma levels . c do 8100 k = 1,nk-1 8100 sh(k) = sqrt ( s(k)*s(k+1) ) sh(nk) = sqrt ( s(nk) ) c c Define moisture constants. c call WETCNS c c Call surface fluxes routine. c call flxsrf(flxqdm,flxcha,flxvap,ts,parsol,t(nk) & ,flxiri,flxvis,vent,z0,ntypes,albedo & ,humsol,q(nk),ps,s(nk),sh(nk)) c c Display Output surface fluxes . c print *,flxqdm,flxcha,flxvap stop end subroutine WETCNS c c This subroutine computes latent heat c of evaporation with respect to water c and ice . c common /wet/eps,ae(2),be(2),ht(2),tc data cp,elv,eli/1004.5,2.500e+6,2.834e+6/ data eps,e0,t0,rgocp/.622,610.71,273.,.286/ tc = 263. ht(1) = elv/cp ht(2) = eli / cp do 8110 i = 1, 2 be(i) = eps * ht(i) / rgocp 8110 ae(i) = be(i) / t0 + alog(e0) return end subroutine FLXSRF (flxqdm,flxcha,flxvap,tss,parsol,ta, & flxiri,flxvis,vent,z0,ntypes,albedo, & humsol,qa,ps,sigmav,sigmat) c c This subroutine computes the surface fluxes of c latent and sensible heat . c parameter ( nk = 12 ) common /sfc/ rib,ustr,ol,cd,ch,cq,z2,z1,sigz2,zv,zt data cp,capa,grav,stebol /1004.5,.286,9.81,5.6678e-8/ data cfonte /15./ c c Calculations are done using 2 levels : c Level 1 refers to the ground surface c Level 2 refers to the lowest model level c ( or measurment level) c rocp = ps/(0.5*(ta + tss)*capa) rho = rocp/cp prg = ps/(rho*grav) zt = (1. - sigmat)*prg zv = (1. - sigmav)*prg u1 = 0.0 u2 = vent t2 = ta * sigmat ** (-capa) q2 = qa c c Ocean case. c if (ntypes .ne. 0) goto 200 q1 = humspc(tss,tss,ps) call SFLX (u1,u2,tss,t2,q1,q2,z1,z2,ntypes,rib, & flxcha,flxvap,flxqdm,rho,humsol) z0 = z1 goto 1000 c c Treatment of land points including c frozen land and sea. Compute net c radiative flux through ice-covered c oceans and melting of ice over c frozen land or sea. c c Set the first guess temperature . c 200 tss = t2 c c Compute the net radiation absorbed at c the ground surface. c raynet = (1.- albedo)*flxvis + & flxiri - stebol*tss**4 c c Test the type of grid point : c Land,land with ice,ocean with ice. c if (ntypes .ne. 2) goto 210 c c Ice covered ocean. c raynet = raynet + parsol * (273.15 - tss) parsol = 0. 210 if (ntypes.lt.1 .or. tss.lt.273.15) goto 220 c c Land with snow c raynet = raynet + cfonte*(273.15-tss) 220 continue flxsrt = (1.-albedo)*flxvis c c Compute the ground temperature by solving the c surface energy balance. c call TG (flxsrt,flxiri,humsol,u1,u2,tss,t2 & ,q1,q2,z0,z2,ps,raynet,rho,parsol) c c Compute surface fluxes from similarity theory. c call SFLX (u1,u2,tss,t2,q1,q2,z0,z2,ntypes,rib, & flxcha,flxvap,flxqdm,rho,humsol) 1000 sigz2 = 1. - (z2/prg) return end subroutine SFLX (u1,u2,t1,t2,q1,q2,z1,z2, & ntypes,rib,fsx,flt,fm,rho,gw) parameter ( nk = 12) c c This subroutine returns surface sensible (fsx),and c latent (flt) heat ,and momentum fluxes via stability c dependent bulk aerodynamic formulation. c c Input parameters : c c ntypes : type of surface c u1,u2 : winds in m/s at level 1 & 2 c t1,t2 : potential temp. at level 1 & 2 c q1,q2 : specific humidity at level 1 & 2 c common /sfc/ rb,ustr,ol,cd,ch,cq,x7,x8,sigz2,zv,zt data g,vkc,alpha,cp/9.81,.35,.04,1004.5/ indx = 0 chp = 0. cqp = 0. if (ntypes .ne. 0) goto 15 c c Ocean case . c z2 = zt c c Set initial value for cd. c cd = 0.0011 c c Compute first guess for the friction c velocity using u2 and then z0. c ustrsq = cd*(u2**2) z0 = alpha*ustrsq/g if (z0 .lt. 1.0e-4) z0 = 1.0e-4 z1 = z0 if (z1 .ge. z2) z2 = z2+z1 c c Complete itterations to obtain the c friction velocity ,and the roughness c lenght from Charnock formula. c do 8130 i = 1, 3 u2dum = u2*alog(z2/z1)/alog(zv/z1) call SFXPAR (u1,u2dum,t1,t2,z1,z2,ustr,tstr & ,rib,ol,cd,ch,cq,chp,cqp,indx) z0 = alpha*(ustr**2)/g if (z0 .lt. 1.0e-4) z0 = 1.0e-4 z1 = z0 if (z1 .ge. z2) z2 = z2+z1 8130 continue chmax = -2.0e-03 cdmax = 2.0e-03 goto 25 c c Land case. c 15 u2dum = u2*alog(z2/z1)/alog(zv/z1) chmax = -5.0e-03 cdmax = 5.0e-03 c c Compute surface parameters and turbulent c exchange coefficients. c 25 call SFXPAR (u1,u2dum,t1,t2,z1,z2,ustr,tstr & ,rib,ol,cd,ch,cq,chp,cqp,indx) cd = amin1(cd,cdmax) ch = amax1(ch,chmax) cq = ch u22 = u2dum u22 = amax1(5.0,u22) c c Compute surface fluxes. c fsx : Sensible heat flux c flt : Latent heat flux c fm : Momentum flux c fsx = rho*ch*cp*u22*(t2-t1) flt = rho*cq*u22*gw*(q2-q1)*2.500e+6 fm = rho*cd*u22*u22 return end subroutine TG (fsg,flg,gw,u1,u2,ts,t2,qs,q2,z1 & ,z2,ps,raynet,rho,parsol) c c This subroutine uses the coupled surface energy c balance and similarity theory to obtain the c surface temperature. c parameter ( nk = 12) common/sfc/rib,ustr,ol,cd,ch,cq,y1,y2,sigz2,zv,zt data stebol /5.6678e-08/ data cp,z0max/1004.5,5./ c ts : computed surface temperature . c qs : computed specific humidity . c c First guess ts is t2 afterwards(during model integration), c the ts at pevious time step is used as a first guess . c c A minimum speed of 1 m/sec is assumed for computing c surface flux . c c t2 and q2 are extrapolated to the lowest sigma level c assuming constant rh and potential temperature . c z2 = zt hl = htbycp(ts)*cp velt = u2*alog(z2/z1)/alog(zv/z1) velt = amax1(5.0,velt) indx = 1 x2 = fsg + flg x3 = 0.622*hl/287. c c Initially, set rnet = raynet ,to account for the c depletion of energy due to ice/snow melting and c absorbtion. Subsequently, compute rnet as , c rnet = x2 - stebol * ts**4 c rnet = raynet c c Aplly Newton-Raphson to solve the surface energy c balance.iteration for surface temperature are coupled c to surface parameters (Loop 100).For each new value of c the ground temperature ,the surface characteristics are c updated until convergence is obtained. c do 8120 il= 1, 7 call SFXPAR (u1,velt,ts,t2,z1,z2,ustr,tstr & ,rib,ol,cd,ch,cq,chp,cqp,indx) c ch = amax1(ch,-5.0e-03) cq = ch swch = 1. c c Compute heat flux (xs) and moisture flux (xl) terms c xs = rho*ch*cp*velt xl = rho*cq*hl*gw*velt qg = humspc(ts,ts,ps) tx = ts*ts c c fts and dfdts are the functions f(ts) and its c derivative (d/dts( f(ts) ) terms required for c the Newton-Raphson iteration method . c dfdts = -4.*stebol*tx*ts+xs+xl*qg*swch*x3/tx & -rho*cp*velt*chp*(t2-ts)-rho*hl*velt*cqp*swch & *gw*(q2-qg) fts = rnet-xs*(t2-ts)-xl*swch*(q2-qg) ts1 = ts - (fts/dfdts) if (abs(ts1-ts) .lt. 0.1) goto 101 ts = ts1 rnet = x2 - stebol*ts*ts*ts*ts 8120 continue ts = t2 goto 201 101 ts = ts1 c c Set a limit on ts in case of non convergence. c 201 if (ts .ge. 323.) then print *,'surface temperature greater than 50 ' ts = 323. endif c c Compute saturation specific humidity at c surface temperature ts. c qs = humspc(ts,ts,ps) return end function HUMSPC (td,t,p) c c This subroutine computes the specific c humidity as a function of dewpoint td c and pressure p depending on phase. c common/wet/eps,ae(2),be(2),ht(2),tc iph = 1 if (t.le.tc) iph = 2 e = exp(ae(iph) - be(iph)/td ) humspc = eps * e / p return end subroutine SFXPAR (u1,u2,t1,t2,z1,z2,ustr,tstr & ,rib,ol,cd,ch,cq,chp,cqp,indx) c c This subroutine returns the following parameters : c c ustr : friction velocity c tstr : characteristic potential temp. c ol : Monin-Obukhov's length c cd : bulk-transfer stress coef. c ch : bulk-transfer sensible heat flux coef. c cq : bulk-transfer latent heat flux coef. c c Constants definitions : c c g : gravitational acceleration c vkc : Von karman's constant c c Derived parametric constants c c rib : bulk richardson's number c b : buoyancy parameter c c Input parameters : c c u1,u2 : winds at heights z1,z2 in m/sec c t1,t2 : potential temps at z1,z2 in deg k c z1,z2 : heights in m c indx : controls computation of second order c terms in the surface energy balance. c data g,vkc/9.81,0.35/ c c Compute dz, du and beta c if du < 1.0 m/s set to 1.0 m/s c dz = z2-z1 du = u2-u1 b = g/((t1+t2)*0.5) if (du .lt. 1.) du = 1.0 dt2t1 = t2-t1 rib = dt2t1*b*dz/(du**2) c c Cap Rib since large Rib produces c small transfer coefficients . c if (rib .gt. 0.212) rib = 0.212 eta = alog(z2/z1) c c Consider only two regimes : c c neutral/stable rib .ge. 0 c unstable rib .lt. 0 c if (rib .ge. 0.) then c c Function fm is computed using eq.(8.39) c fm = 1. / ((1. + 4.7*rib)**2) fh = fm end if if (rib .lt. 0.) then c c cmm and chh are computed using eqs.(8.36) c with C* = 7.4 and C* = 5.3 ,respectively. c fm and fh are computed from eqs.(8.40) and c (8.41). c cmm = 7.4*vkc*vkc*9.4*sqrt(z2/z1)/(eta*eta) chh = 5.3*vkc*vkc*9.4*sqrt(z2/z1)/(eta*eta) fm = 1. - 9.4*rib/(1.+cmm*sqrt(abs(rib))) fh = 1. - 9.4*rib/(1.+chh*sqrt(abs(rib))) end if c c Compute remaining quantities . c ustrsq = vkc*vkc*du*du *fm/( eta*eta) tstrus = vkc*vkc*du*dt2t1*fh/(0.74*eta*eta) ustr = sqrt(ustrsq) tstr = tstrus/ustr cd = vkc*vkc*fm/( eta*eta) ch = -vkc*vkc*fh/(0.74*eta*eta) cq = ch ol = 999. if (tstrus .ne. 0.) then ol = ustrsq/(b*vkc*tstrus/sqrt(ustrsq)) end if chp = 0. cqp = 0. return end function HTBYCP (t) c c This funciotn computes latent heat c over cp depending on phase. c common/wet/eps,ae(2),be(2),ht(2),tc if (t.gt.tc) htbycp = ht(1) if (t.le.tc) htbycp = ht(2) return end