!----------------------------------- subroutine sfc_land & !................................... ! --- inputs: & ( im, km, ps, u1, v1, t1, q1, smc, soiltyp, & & sigmaf, vegtype, sfcemis, dlwflx, swnet, delt, & & zorl, tg3, cm, ch, prsl1, prslki, slimsk, & & ddvel, flag_iter, flag_guess, & ! --- input/outputs: & weasd, tskin, tprcp, srflag, stc, canopy, tsurf, & ! --- outputs: & qsurf, snowmt, gflux, zsoil, rhscnpy, rhsmc, & & aim, bim, cim, drain, evap, hflx, ep, cmm, chh, & & evbs, evcw, trans, sbsno, snowc, stm, snohf, & & twilt, tref & ) ! ===================================================================== ! ! description: osu land surface model ! ! ! ! usage: ! ! ! ! call sfc_land ! ! inputs: ! ! ( im, km, ps, u1, v1, t1, q1, smc, soiltyp, ! ! sigmaf, vegtype, sfcemis, dlwflx, swnet, delt, ! ! zorl, tg3, cm, ch, prsl1, prslki, slimsk, ! ! ddvel, flag_iter, flag_guess, ! ! input/outputs: ! ! weasd, tskin, tprcp, srflag, stc, canopy, tsurf, ! ! outputs: ! ! qsurf, snowmt, gflux, zsoil, rhscnpy, rhsmc, ! ! aim, bim, cim, drain, evap, hflx, ep, cmm, chh, ! ! evbs, evcw, trans, sbsno, snowc, stm, snohf, ! ! twilt, tref ) ! ! ! ! subprograms called: none ! ! ! ! ! ! program history log: ! ! xxxx -- original version created from Hula Lu's progtm ! ! ! 200x -- sarah lu modified (need description) ! ! oct 2006 -- h. wei modified (need description) ! ! apr 2009 -- y.-t. hou modified to include surface emissivity ! ! effect on lw radiation. also replaced slrad (a ! ! confussing term) with sfc net sw flux swnet ! ! that is redefined as (du-up). rewrite the code ! ! and add program documentation block. ! ! sep 2009 -- s. moorthi some additional modification ! ! ! ! ! ! ==================== defination of variables ==================== ! ! ! ! inputs: size ! ! im, km - integer, horiz dimension and num of soil layers 1 ! ! ps - real, surface pressure im ! ! u1, v1 - real, u/v component of surface layer wind im ! ! t1 - real, surface layer mean temperature ( k ) im ! ! q1 - real, surface layer mean specific humidity im ! ! smc - real, soil moisture content (fractional) im,km ! ! soiltyp - integer, soil type (integer index) im ! ! sigmaf - real, areal fractional cover of green vegetation im ! ! vegtype - integer, vegetation type (integer index) im ! ! sfcemis - real, sfc lw emissivity ( fraction ) im ! ! dlwflx - real, total sky sfc downward lw flux ( w/m**2 ) im ! ! swnet - real, total sky sfc netsw flx into ground(w/m**2) im ! ! delt - real, time interval (second) 1 ! ! zorl - real, surface roughness im ! ! tg3 - real, deep soil temperature im ! ! cm - real, surface exchange coeff for momentum (m/s) im ! ! ch - real, surface exchange coeff heat & moisture(m/s) im ! ! prsl1 - real, surface layer mean pressure im ! ! prslki - real, im ! ! slimsk - real, sea/land/ice mask (=0/1/2) im ! ! ddvel - real, im ! ! flag_iter- logical, im ! ! flag_guess-logical, im ! ! ! ! input/outputs: ! ! weasd - real, water equivalent accumulated snow depth (mm)im ! ! tskin - real, ground surface skin temperature ( k ) im ! ! tprcp - real, total precipitation im ! ! srflag - real, snow/rain flag for precipitation im ! ! stc - real, soil temp (k) im,km ! ! canopy - real, canopy moisture content (m) im ! ! tsurf - real, surface skin temperature (after iteration) im ! ! ! ! outputs: ! ! qsurf - real, specific humidity at sfc im ! ! snowmt - real, snow melt (m) im ! ! gflux - real, soil heat flux (w/m**2) im ! ! zsoil - real, soil depth im,km ! ! rhscnpy - real, im ! ! rhsmc - real, im,km ! ! aim - real, tridiagonal matrix coeff for soil moist im,km ! ! bim - real, tridiagonal matrix coeff for soil moist im,km ! ! cim - real, tridiagonal matrix coeff for soil moist im,km ! ! drain - real, subsurface runoff im ! ! evap - real, evaperation from latent heat flux im ! ! hflx - real, sensible heat flux im ! ! ep - real, potential evaporation im ! ! cmm - real, im ! ! chh - real, im ! ! evbs - real, direct soil evaporation (m/s) im ! ! evcw - real, canopy water evaporation (m/s) im ! ! trans - real, im ! ! sbsno - real, sublimation/deposit from snopack (m/s) im ! ! snowc - real, fractional snow cover im ! ! stm - real, total soil column moisture content (m) im ! ! snohf - real, snow/freezing-rain latent heat flux (w/m**2)im ! ! twilt - real, dry soil moisture threshold im ! ! tref - real, soil moisture threshold im ! ! ! ! ==================== end of description ===================== ! ! use machine , only : kind_phys use funcphys, only : fpvs use physcons, only : grav => con_g, sbc => con_sbc, cp => con_cp, & & hvap => con_hvap, hfus => con_hfus, & & eps => con_eps, epsm1 => con_epsm1, & & t0c => con_t0c, rvrdm1 => con_fvirt, & & rd => con_rd ! implicit none ! ! --- constant parameters: real (kind=kind_phys), parameter :: cpinv = 1.0/cp real (kind=kind_phys), parameter :: hvapi = 1.0/hvap real (kind=kind_phys), parameter :: elocp = hvap/cp real (kind=kind_phys), parameter :: dfsnow = 0.31 real (kind=kind_phys), parameter :: ch2o = 4.2e6 real (kind=kind_phys), parameter :: csoil = 1.26e6 real (kind=kind_phys), parameter :: scanop = 0.5 real (kind=kind_phys), parameter :: cfactr = 0.5 real (kind=kind_phys), parameter :: zbot =-3.0 real (kind=kind_phys), parameter :: topt = 298.0 real (kind=kind_phys), parameter :: rhoh2o = 1000.0 real (kind=kind_phys), parameter :: ctfil1 = 0.5 real (kind=kind_phys), parameter :: ctfil2 = 1.0-ctfil1 real (kind=kind_phys), parameter :: snomin = 1.0e-9 ! --- input: integer, intent(in) :: im, km, soiltyp(im), vegtype(im) real (kind=kind_phys), dimension(im), intent(in) :: ps, u1, v1, & & t1, q1, sigmaf, sfcemis, dlwflx, swnet, zorl, tg3, cm, ch, & & prsl1, prslki, slimsk, ddvel real (kind=kind_phys), dimension(im,km), intent(in) :: smc real (kind=kind_phys), intent(in) :: delt logical, intent(in) :: flag_iter(im), flag_guess(im) ! --- in/out: real (kind=kind_phys), dimension(im), intent(inout) :: weasd, & & tskin, tprcp, srflag, canopy, tsurf real (kind=kind_phys), dimension(im,km), intent(inout) :: stc ! --- output: real (kind=kind_phys), dimension(im), intent(out) :: qsurf, & & snowmt, gflux, rhscnpy, drain, evap, hflx, ep, chh, cmm, & & evbs, evcw, trans, sbsno, snowc, stm, snohf, twilt, & & tref real (kind=kind_phys), dimension(im,km), intent(out) :: zsoil, & & rhsmc, aim, bim, cim ! --- external functions: real (kind=kind_phys) :: funcdf, funckt, ktsoil ! --- locals: real (kind=kind_phys), dimension(im) :: weasd_old, tprcp_old, & & srflag_old, tskin_old, canopy_old, wind, canfac, ddz, & & ddz2, delta, dew, df1, dft0, dft1, dft2, dmdz, dmdz2, & & dtdz1, dtdz2, ec, edir, etpfac, factsnw, fx, gx, hcpct, & & partlnd, q0, qs1, qss, rcap, rch, rho, rs, rsmall, & & slwd, smcz, snoevp, snowd, term1, term2, theta1, tv1, & & tsea, xx, yy, zz, kt1, kt2 real (kind=kind_phys), dimension(im,km) :: stc_old, et, stsoil, & & ai, bi, ci, rhstc real (kind=kind_phys) :: bfact, cc, delt2, df2, eth, ff, g, rcq, & & rcs, rct, rsi, rss, smcdry, t12, t14, tflx, tem integer :: i, k logical :: flag(im), flagsnw(im) ! --- local data arrays: ! the 13 vegetation types are: ! 1 ... broadleave-evergreen trees (tropical forest) ! 2 ... broadleave-deciduous trees ! 3 ... broadleave and needle leave trees (mixed forest) ! 4 ... needleleave-evergreen trees ! 5 ... needleleave-deciduous trees (larch) ! 6 ... broadleave trees with groundcover (savanna) ! 7 ... groundcover only (perenial) ! 8 ... broadleave shrubs with perenial groundcover ! 9 ... broadleave shrubs with bare soil ! 10 ... dwarf trees and shrubs with ground cover (trunda) ! 11 ... bare soil ! 12 ... cultivations (use parameters from type 7) ! 13 ... glacial real(kind=kind_phys),dimension(13), save :: rsmax, rsmin, rgl, hs data rsmax / 13*5000.0 / data rsmin / 150., 100., 125., 150., 100., 70., 40., & & 300., 400., 150., 999., 040., 999. / data rgl / 5*30., 65., 4*100., 999., 100., 999. / data hs / 41.69, 54.53, 51.93, 47.35, 47.35, 54.53, 36.35, & & 3*42.00, 999.0, 36.35, 999.0 / real(kind=kind_phys), dimension(9), save :: smdry, smref, smwlt data smdry / .07, .14, .22, .08, .18, .16, .12, .10, .07 / data smref / .283,.387,.412,.312,.338,.382,.315,.329,.283 / data smwlt / .029,.119,.139,.047,.010,.103,.069,.066,.029 / ! !===> ... begin here ! delt2 = delt * 2.0 ! --- ... set default flag for land do i = 1, im flag(i) = ( slimsk(i) == 1.0 ) enddo ! --- ... save land-related prognostic fields for guess run do i = 1, im if (flag(i) .and. flag_guess(i)) then weasd_old(i) = weasd(i) tskin_old (i) = tskin(i) canopy_old(i) = canopy(i) tprcp_old (i) = tprcp(i) srflag_old(i) = srflag(i) do k = 1, km stc_old(i,k) = stc(i,k) enddo endif enddo ! --- ... initialize variables. all units are supposedly m.k.s. unless ! specifie ps is in pascals ! wind is wind speed, theta1 is adiabatic surface temp from ! level 1, rho is density, qs1 is sat. hum. at level1 and qss ! is sat. hum. at surface ! surface roughness length is converted to m from cm ! net sw flux swnet is dn-up, and dlw is positive dnwd ! qs1 = fpvs(t1) ! qss = fpvs(tskin) do i = 1, im if (flag_iter(i) .and. flag(i)) then slwd(i) = swnet(i) + dlwflx(i) wind(i) = sqrt(u1(i)*u1(i) + v1(i)*v1(i)) & & + max(0.0, min(ddvel(i), 30.0)) wind(i) = max(wind(i), 1.0) q0(i) = max(q1(i), 1.e-8) ! tsurf(i) = tskin(i) tsea(i) = tsurf(i) theta1(i) = t1(i) * prslki(i) tv1(i) = t1(i) * (1.0 + rvrdm1*q0(i)) ! thv1(i) = theta1(i) * (1.0 + rvrdm1*q0(i)) ! tvs(i) = tsea(i) * (1.0 + rvrdm1*q0(i)) rho(i) = prsl1(i) / (rd * tv1(i)) qs1(i) = fpvs(t1(i)) qs1(i) = eps*qs1(i) / (prsl1(i) + epsm1*qs1(i)) qs1(i) = max(qs1(i), 1.e-8) q0(i) = min(qs1(i), q0(i)) qss(i) = fpvs(tskin(i)) !!! change to tsurf? qss(i) = eps*qss(i) / (ps(i) + epsm1*qss(i)) ! rs = plantr rs(i) = 0. if (vegtype(i) > 0.0) rs(i) = rsmin(vegtype(i)) canopy(i) = max(canopy(i), 0.0) factsnw(i) = 10.0 ! --- ... snow depth in water equivalent is converted from mm to m unit snowd(i) = weasd(i) * 0.001 flagsnw(i) = .false. ! --- ... when snow depth is less than 1 mm, a patchy snow is assumed ! and soil is allowed to interact with the atmosphere. ! we should eventually move to a linear combination of soil and ! snow under the condition of patchy snow. if (snowd(i)>0.001 .or. slimsk(i)==2.0) rs(i) = 0.0 if (snowd(i)>0.001) flagsnw(i) = .true. endif ! end if_flag_iter_block enddo ! end do_i_loop do i = 1, im if (flag_iter(i) .and. flag(i)) then zsoil(i,1) = -0.10 do k = 2, km zsoil(i,k) = zsoil(i,k-1) + (-2.0 - zsoil(i,1)) / (km - 1) enddo ! --- ... wei: use the same soil layer structure as noah if running with 4-layer if (km > 0.2)then zsoil(i,2) = -0.4 zsoil(i,3) = -1.0 zsoil(i,4) = -2.0 endif endif enddo do i = 1, im if (flag_iter(i) .and. flag(i)) then drain(i) = 0.0 endif enddo do k = 1, km do i = 1, im if (flag_iter(i) .and. flag(i)) then et (i,k) = 0.0 rhsmc(i,k) = 0.0 aim (i,k) = 0.0 bim (i,k) = 1.0 cim (i,k) = 0.0 stsoil(i,k) = stc(i,k) endif enddo enddo do i = 1, im if (flag_iter(i) .and. flag(i)) then edir(i) = 0.0 ec (i) = 0.0 evap(i) = 0.0 hflx(i) = 0.0 ep (i) = 0.0 fx (i) = 0.0 snowmt(i) = 0.0 gflux (i) = 0.0 rhscnpy(i)= 0.0 etpfac(i) = 0.0 canfac(i) = 0.0 evbs (i) = 0.0 evcw (i) = 0.0 trans(i) = 0.0 sbsno(i) = 0.0 snowc(i) = 0.0 snohf(i) = 0.0 endif enddo ! --- ... rcp = rho cp ch v do i = 1, im if (flag_iter(i) .and. flag(i)) then rch(i) = rho(i) * cp * ch(i) * wind(i) cmm(i) = cm(i) * wind(i) chh(i) = rho(i) * ch(i) * wind(i) endif enddo ! --- ... compute soil/snow/ice heat flux in preparation for surface ! energy balance calculation do i = 1, im if (flag_iter(i) .and. flag(i)) then smcz(i) = 0.5 * (smc(i,1) + 0.20) dft0(i) = ktsoil(smcz(i), soiltyp(i)) endif enddo do i = 1, im if (flag_iter(i) .and. flag(i)) then if (flagsnw(i)) then ! --- ... when snow covered, ground heat flux comes from snow tflx = min(t1(i), tsea(i)) gflux(i) = -dfsnow * (tflx - stsoil(i,1)) & & / (factsnw(i) * max(snowd(i), 0.001)) else gflux(i) = dft0(i) * (stsoil(i,1) - tsea(i)) & & / (-0.5 * zsoil(i,1)) endif gflux(i) = max(gflux(i), -200.0) gflux(i) = min(gflux(i), 200.0) endif enddo do i = 1, im if (flag_iter(i) .and. flag(i)) then partlnd(i) = 1.0 if (snowd(i)>0.0 .and. snowd(i)<=0.001) then partlnd(i) = 1.0 - snowd(i) / 0.001 endif endif enddo do i = 1, im if (flag_iter(i) .and. flag(i)) then snoevp(i) = 0.0 if (snowd(i) > 0.001) partlnd(i) = 0.0 endif enddo ! --- ... compute potential evaporation for land and sea ice do i = 1, im if (flag_iter(i) .and. flag(i)) then t12 = t1(i) * t1(i) t14 = t12 * t12 ! --- ... rcap = fnet - sigma t**4 + gflx - rho cp ch v (t1-theta1) rcap(i) = slwd(i) - sfcemis(i)*sbc*t14 + gflux(i) & & - rch(i)*(t1(i) - theta1(i)) ! --- ... rsmall = 4 sigma t**3 / rch(i) + 1 rsmall(i) = 4.0*sfcemis(i)*sbc*t1(i)*t12 / rch(i) + 1.0 ! --- ... delta = l / cp * dqs/dt delta(i) = elocp*eps*hvap*qs1(i) / (rd*t12) ! --- ... potential evapotranspiration ( watts / m**2 ) and ! potential evaporation term1(i) = elocp*rsmall(i)*rch(i) * (qs1(i) - q0(i)) term2(i) = rcap(i) * delta(i) ep(i) = (elocp*rsmall(i)*rch(i) * (qs1(i) - q0(i)) & & + rcap(i)*delta(i)) ep(i) = ep(i) / (rsmall(i) + delta(i)) endif enddo ! --- ... actual evaporation over land in three parts : edir, et, and ec ! direct evaporation from soil, the unit goes from m s-1 to kg m-2 s-1 do i = 1, im flag(i) = (slimsk(i)==1.0) .and. (ep(i)>0.0) enddo do i = 1, im if (flag_iter(i))then if (flag(i)) then df1(i) = funcdf(smc(i,1),soiltyp(i)) kt1(i) = funckt(smc(i,1),soiltyp(i)) endif if (flag(i) .and. stc(i,1)0.0) enddo ! --- ... change of snow depth due to evaporation or sublimation ! convert evap from kg m-2 s-1 to m s-1 to determine the reduction of s do i = 1, im if (flag_iter(i) .and. flag(i)) then bfact = snowd(i) / (delt*ep(i) / (hvap*rhoh2o)) bfact = min(bfact, 1.0) ! --- ... the evaporation of snow if (ep(i) <= 0.0) bfact = 1.0 if (snowd(i) <= 0.001) then ! evap = (snowd(i)/0.001) * bfact*ep(i) + evap ! snoevp(i) = bfact*ep(i) * (1.0 - partlnd(i)) ! evap = evap + snoevp(i) snoevp(i) = bfact * ep(i) ! evap = evap + snoevp(i) * (1.0 - partlnd(i)) evap(i) = evap(i) + snoevp(i) * (1.0 - partlnd(i)) else ! evap(i) = bfact * ep(i) snoevp(i) = bfact * ep(i) evap(i) = snoevp(i) endif tsea(i) = t1(i) & & + (rcap(i) - gflux(i) - dfsnow*(t1(i) - stsoil(i,1)) & & / (factsnw(i) * max(snowd(i), 0.001)) & ! & + theta1 - t1 & ! & - bfact * ep(i)) / (rsmall(i) * rch(i) & & - snoevp(i)) / (rsmall(i) * rch(i) & & + dfsnow / (factsnw(i) * max(snowd(i), 0.001))) ! snowd(i) = snowd(i) - bfact*ep(i)*delt / (rhoh2o*hvap) snowd(i) = snowd(i) - snoevp(i)*delt / (rhoh2o*hvap) snowd(i) = max(snowd(i), 0.0) endif enddo ! --- ... snow melt (m) do i = 1, im flag(i) = (slimsk(i)==1.0) .and. (snowd(i)>0.0) enddo do i = 1, im if (flag_iter(i)) then if (flag(i) .and. tsea(i)>t0c) then snowmt(i) = rch(i)*rsmall(i)*delt * (tsea(i) - t0c) & & / (rhoh2o*hfus) snowmt(i) = min(snowmt(i), snowd(i)) snowd(i) = snowd(i) - snowmt(i) snowd(i) = max(snowd(i), 0.0) tsea (i) = max(t0c, tsea(i) - hfus*snowmt(i)*rhoh2o & & / (rch(i)*rsmall(i)*delt)) endif endif enddo ! --- ... we need to re-evaluate evaporation because of snow melt ! the skin temperature is now bounded to 0 deg c ! qss = fpvs(tsea) do i = 1, im flag(i) = (slimsk(i) == 1.0) if (flag_iter(i) .and. flag(i))then ! if (snowd(i) > 0.0) then if (snowd(i) > snomin) then !jfe qss(i) = 1000.0 * fpvs(tsea(i)) qss(i) = fpvs(tsea(i)) qss(i) = eps*qss(i) / (ps(i) + epsm1*qss(i)) evap(i) = elocp*rch(i) * (qss(i) - q0(i)) endif endif enddo ! --- ... prepare tendency terms for the soil moisture field without ! precipitat. the unit of moisture flux needs to become m s-1 ! for soil moisture. hence the factor of rhoh2o do i = 1, im if (flag_iter(i)) then if (flag(i)) then df1(i) = funcdf(smcz(i),soiltyp(i)) kt1(i) = funckt(smcz(i),soiltyp(i)) endif if (flag(i) .and. stc(i,1) 0.001) then if (flag(i) .and. flagsnw(i)) then yy(i) = stsoil(i,1) ! --- ... heat flux from snow is explicit in time zz(i) = 1.0 xx(i) = dfsnow * (stsoil(i,1) - tsea(i)) & & / (-factsnw(i) * max(snowd(i), 0.001)) endif endif enddo ! --- ... compute the forcing and the implicit matrix elements for update ! ch2o is the heat capacity of water and csoil is the heat capacity ! of soil do i = 1, im if (flag_iter(i) .and. flag(i)) then smcz(i) = max(smc(i,1), smc(i,2)) dtdz1(i) = (stsoil(i,1) - stsoil(i,2)) / (-0.5*zsoil(i,2)) dft1(i) = ktsoil(smcz(i),soiltyp(i)) hcpct(i) = smc(i,1)*ch2o + (1.0 - smc(i,1)) * csoil dft2(i) = dft1(i) ddz(i) = 1.0 / (-0.5*zsoil(i,2)) ! --- ... ai, bi, and ci are the elements of the tridiagonal matrix for the ! implicit update of the soil temperature ai(i,1) = 0.0 bi(i,1) = dft1(i)*ddz(i) / (-zsoil(i,1)*hcpct(i)) ci(i,1) = -bi(i,1) bi(i,1) = bi(i,1) & & + dft0(i) / (0.5*zsoil(i,1)**2 * hcpct(i)*zz(i)) rhstc(i,1) = (dft1(i)*dtdz1(i) - xx(i))/(zsoil(i,1)*hcpct(i)) endif enddo do k = 2, km do i = 1, im if (flag_iter(i) .and. flag(i)) then hcpct(i) = smc(i,k)*ch2o + (1.0 - smc(i,k)) * csoil endif enddo if (k < km) then do i = 1, im if (flag_iter(i) .and. flag(i)) then dtdz2(i) = (stsoil(i,k) - stsoil(i,k+1)) & & / (0.5 * (zsoil(i,k-1) - zsoil(i,k+1))) smcz(i) = max(smc(i,k), smc(i,k+1)) dft2(i) = ktsoil(smcz(i),soiltyp(i)) ddz2(i) = 2.0 / (zsoil(i,k-1) - zsoil(i,k+1)) ci(i,k) = -dft2(i) * ddz2(i) & & / ((zsoil(i,k-1) - zsoil(i,k)) * hcpct(i)) endif enddo else ! if_k_block ! --- ... at the bottom, climatology is assumed at 2m depth for land and ! freezing temperature is assumed for sea ice at z(km) do i = 1, im if (flag_iter(i) .and. flag(i)) then dtdz2(i) = (stsoil(i,k) - tg3(i)) & & / (0.5 * (zsoil(i,k-1) + zsoil(i,k)) - zbot) dft2(i) = ktsoil(smc(i,k),soiltyp(i)) ci(i,k) = 0.0 endif enddo endif ! end if_k_block do i = 1, im if (flag_iter(i) .and. flag(i)) then rhstc(i,k) = (dft2(i)*dtdz2(i) - dft1(i)*dtdz1(i)) & & / ((zsoil(i,k) - zsoil(i,k-1)) * hcpct(i)) ai(i,k) = -dft1(i) * ddz(i) & & / ((zsoil(i,k-1) - zsoil(i,k)) * hcpct(i)) bi(i,k) = -(ai(i,k) + ci(i,k)) dft1(i) = dft2(i) dtdz1(i) = dtdz2(i) ddz(i) = ddz2(i) endif enddo enddo ! end do_k_loop ! --- ... solve the tri-diagonal matrix do k = 1, km do i = 1, im if (flag_iter(i) .and. flag(i)) then rhstc(i,k) = rhstc(i,k) * delt2 ai(i,k) = ai(i,k) * delt2 bi(i,k) = 1.0 + bi(i,k)*delt2 ci(i,k) = ci(i,k) * delt2 endif enddo enddo ! --- ... forward elimination do i = 1, im if (flag_iter(i) .and. flag(i)) then ci(i,1) = -ci(i,1) / bi(i,1) rhstc(i,1) = rhstc(i,1) / bi(i,1) endif enddo do k = 2, km do i = 1, im if (flag_iter(i) .and. flag(i)) then cc = 1.0 / (bi(i,k) + ai(i,k)*ci(i,k-1)) ci(i,k) = -ci(i,k) * cc rhstc(i,k) = (rhstc(i,k) - ai(i,k)*rhstc(i,k-1)) * cc endif enddo enddo ! --- ... backward substituttion do i = 1, im if (flag_iter(i) .and. flag(i)) then ci(i,km) = rhstc(i,km) endif enddo do k = km-1, 1 do i = 1, im if (flag_iter(i) .and. flag(i)) then ci(i,k) = ci(i,k)*ci(i,k+1) + rhstc(i,k) endif enddo enddo ! --- ... update soil and ice temperature do k = 1, km do i = 1, im if (flag_iter(i) .and. flag(i)) then stsoil(i,k) = stsoil(i,k) + ci(i,k) endif enddo enddo ! --- ... update surface temperature for snow free surfaces do i = 1, im if (flag_iter(i)) then if (flag(i) .and. .not.flagsnw(i)) then tsea(i) = (yy(i) + (zz(i) - 1.0) * stsoil(i,1)) / zz(i) endif endif enddo ! --- ... time filter for soil and skin temperature do i = 1, im if (flag_iter(i) .and. flag(i)) then tsurf(i) = ctfil1*tsea(i) + ctfil2*tsurf(i) endif enddo do k = 1, km do i = 1, im if (flag_iter(i) .and. flag(i)) then stc(i,k) = ctfil1*stsoil(i,k) + ctfil2*stc(i,k) endif enddo enddo ! --- ... gflux calculation do i = 1, im if (flag_iter(i) .and. flag(i)) then if (flagsnw(i)) then gflux(i) = -dfsnow * (tsurf(i) - stc(i,1)) & & / (factsnw(i) * max(snowd(i), 0.001)) else gflux(i) = dft0(i) * (stc(i,1) - tsurf(i)) & & / (-0.5*zsoil(i,1)) endif endif enddo ! --- ... calculate sensible heat flux do i = 1, im if (flag_iter(i) .and. flag(i)) then hflx(i) = rch(i) * (tsurf(i) - theta1(i)) endif enddo ! --- ... the rest of the output do i = 1, im if (flag_iter(i) .and. flag(i)) then qsurf(i) = q1(i) + evap(i) / (elocp*rch(i)) evbs (i) = edir(i) evcw (i) = ec(i) sbsno(i) = snoevp(i) snowc(i) = 1.0 - partlnd(i) stm (i) = -smc(i,1) * zsoil(i,1) snohf(i) = dfsnow * (t1(i) - stsoil(i,1)) trans(i) = et(i,1) do k = 2, km stm(i) = stm(i) + smc(i,k) * (zsoil(i,k-1) - zsoil(i,k)) trans(i) = trans(i) + et(i,k) enddo ! --- ... convert snow depth back to mm of water equivalent weasd(i) = snowd(i) * 1000.0 endif enddo do i = 1, im if (flag_iter(i) .and. flag(i)) then tem = 1.0 / rho(i) hflx(i) = hflx(i) * tem * cpinv evap(i) = evap(i) * tem * hvapi endif enddo ! --- ... restore land-related prognostic fields for guess run do i = 1, im if (flag(i)) then if (flag_guess(i)) then weasd(i) = weasd_old(i) tskin(i) = tskin_old(i) canopy(i) = canopy_old(i) tprcp(i) = tprcp_old(i) srflag(i) = srflag_old(i) do k = 1, km stc(i,k) = stc_old(i,k) enddo else tskin(i) = tsurf(i) endif endif enddo return end