program SIMFLX_EXP c c This program shows the calculation of the surface fluxes c for a limited domain. Calculations are presented as if c they were in a numerical model except that there is no time c integration.This case would correspond to fluxes computation c for one time step for all grid points. c c Parameters definition : c c l : east-west dimension c m : south-north dimension c nk : number of vertical levels c slat : southernmost latitude c wlon : westernmost longitude c hz : time of the day in hours c day : day of the month c month : month of the year c dphi : grid distance in the north-south direction c dlambda : grid distance in the east - west direction c kount : counter of the number of time steps c pas : the time step in secondes. c stebol : stefan-boltzman constant c c c Variables definitions c c Input : c c wt : tridimensional temperature field c wq : tridimensional specific humidity field c t : temperature profile for the considered point c q : specific humdity // c s : sigma levels c sh : intermediate sigma levels c declsc : sun declination angle c alb : surface albedo c sfps : surface pressure c tss : surface temperature c ven : surface wind c rlon : longitude c alat : latitude c c Output : c c warm : heating rate c cool : cooling rate c rtrter : vertical profile of terrestrial radiative flux c rtrsol : vertical profile of solar radiative flux c flxvis : flux of shortwave radiation at surface c flxiri : flux of infra-red radiation at surface c irtop : infra-red radiation at top of atmosphere c c common blocks wetcns and rrr define the moisture c constants and other variables used in the radiative c transfer calculations. c parameter (l = 17,m = 10,nk = 12,nk2 = nk+2) real t(nk),q(nk),sh(nk),s(nk),declsc(2) real rtrter(nk),rtrsol(nk),irtop(l,m) real alb(l,m),sfps(l,m),tss(l,m),rlon(l,m) real ven(l,m),alat(l,m),cool(l,m,nk) real wt(l,m,nk),wq(l,m,nk),warm(l,m,nk) common/wetc/eps,ae(2),be(2),ht(2),tc common/rrr /fup(nk2,8),fdw(nk2,8),rw(nk2) &,fx(nk2),tl(nk2),aso(nk2,8),rso(nk2,8),albc(8),rsa,secz, &coe(8),dtc(nk2,8),fnet(nk2,8),cld(nk2,8),rhx(nk),dp(nk), &qz(nk),p(nk2),tz(nk),dpz(nk2),qs(nk),sso(8),dtw(nk2,8) dimension ntyp(l,m),pars(l,m),rlen(l,m),hums(l,m) real flxq(l,m),flxc(l,m),flxv(l,m),topo(l,m) common /rhh /rh(nk),dqq(nk) common /ipos /ihem,nbad,nbad1,nshal common /loc /ic,jc common /extra/rib,ustr,ol,cd,ch,cq,z2,z1,sigz2,zv,zt common /density/rho c c Input : c c ntyp : land-sea matrix c pars : soil parameter c rlen : roughness lenth c hums : soil humdity c topo : topography c c Output : c c flxq : momentum flux c flxc : heat flux c flxv : moisture flux c swc and swr are logical switchs for shortwaves c and longwaves radiation calculations,respectively c logical swc,swr data swr ,swc /.true.,.true./ data pi,dphi,dlambda /3.14159 ,5.5376 ,5.625 / data s /.1,.2,.3,.4,.5,.6,.7,.8,.85,.9,.95,.99/ data pas,kount/1.,1./ data stebol /5.6678e-8/ data slat,wlong,hz,day,month/2.77,-28.13,12.,20,12/ picon = pi/180. open (11,file='ts17.dat',status='old') open (12,file='ps17.dat',status='old') open (13,file='alb17.dat',status='old') open (14,file='vent17.dat',status='old') open (15,file='t17.dat',status='old') open (16,file='q17.dat',status='old') open (17,file='ntyp17.dat',status='old') open (18,file='pars17.dat',status='old') open (19,file='z017.dat',status='old') open (20,file='hums17.dat',status='old') open (21,file='xht17.dat ',status='old') open (22,file='warm.dat ',status='unknown ') open (23,file='cool.dat',status='unknown') open (24,file='surflx.dat',status='unknown') c 1000 format(6e13.6) 1001 format(10(1x,17i2,/)) 100 format(6e13.6) 101 format(10(1x,17i2,/)) c c In normal model operations variables in files 11 to 21 c are read in the initialization procedure and carried c by the model. In this example these variables are read in. c c Read temperature and specific humdity for all levels c do 8200 k = 1, nk read (15,1000) ((wt(i,j,k),i=1,l),j=1,m) read (16,1001) ((wq(i,j,k),i=1,l),j=1,m) 8200 continue c c Read surface variables . c read (11,100)((tss (i,j),i=1,l),j=1,m) read (12,100)((sfps (i,j),i=1,l),j=1,m) read (13,100)((alb (i,j),i=1,l),j=1,m) read (14,100)((ven (i,j),i=1,l),j=1,m) read (17,101)((ntyp (i,j),i=1,l),j=1,m) read (18,100)((pars (i,j),i=1,l),j=1,m) read (19,100)((rlen (i,j),i=1,l),j=1,m) read (20,100)((hums (i,j),i=1,l),j=1,m) read (21,100)((topo (i,j),i=1,l),j=1,m) c c Compute declination angle from day and month . c call CONRAD (declsc,day,month) c c Compute intermediate sigma levels (sh) c do 8202 k = 1, nk-1 sh(k) = sqrt ( s(k)*s(k+1) ) 8202 continue sh(nk) = sqrt ( s(nk) ) c c Main loop for all grid points . c do 8204 j = 1, m do 8204 i = 1, l c c Compute the position of the given grid point c alat(i,j) = (slat + (j-1)*dphi)*picon rlon(i,j) = (wlon + (i-1)*dlambda)*picon c c Prepare colone arrays for the given grid point c do 8206 k = 1, nk t(k) = wt(i,j,k) q(k) = wq(i,j,k) 8206 continue c c Prepare surface parameters c parsol = pars (i,j) z0 = rlen (i,j) ntypes = ntyp (i,j) humsol = hums (i,j) xht1 = topo (i,j) ts = tss (i,j) ps = sfps (i,j) rlong = rlon (i,j) sinlat = sin (alat(i,j)) coslat = cos (alat(i,j)) albedo = alb (i,j) vent = ven (i,j) f = 2.*7.29e-05*sinlat ic = i jc = j c c Prepare moisture constants c call WETCNS c c Call radiation routine to obtain the infra-red and shortwave c radiative fluxes at the surface . these are required inputs c for surface fluxes calculation. c call RAD (rtrsol,rtrter,flxvis,flxiri,ts,t,q,ps,sh,nk, 1 rlong,sinlat,coslat,s,hz,albedo,vent,declsc, 2 swr,swc,irtop,ic,jc) c c Save output of warm and cool for the grid point c do 8208 k = 1, nk warm(i,j,k) = rtrsol(k) cool(i,j,k) = rtrter(k) 8208 continue c c Call surface fluxes routine for the given grid point . c call flxsrf(flxqdm,flxcha,flxvap,ts,parsol,t(nk),flxiri, & flxvis,vent,z0,ntypes,albedo,humsol,f,ch1,ch2, & ch3,q(nk),ps,s(nk),pas,sh(nk),xht1,kount) c c Save surface fluxes output for the given grid point c flxq(i,j) = flxqdm flxc(i,j) = flxcha flxv(i,j) = flxvap 8204 continue 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 subroutine CONRAD (declsc,day,month) c c This subroutine computes the declination angle . c The Julian month is assumed to have length of c 365/12=30.4 c real declsc(2) data pi /3.14159/ ajour = 30.4 * (month-1) + day rdecl = .412 *cos((ajour+10.)*2.*pi /365.25 - pi) declsc(1) = sin (rdecl) declsc(2) = cos (rdecl) return end