module StieglitzSnow ! reichle, 12 Aug 2014 - moved "*_calc_asnow()" to here from catchment.F90 ! - renamed "get_tf_nd()" to "StieglitzSnow_calc_tpsnow()" ! - removed "relayer()" (obsolete) ! - renamed "N_sm" ("soil moisture") to "N_zones" (could be LandIce...) ! - cleaned up MAPL constants ! - additional minor clean-up USE catch_constants, ONLY: & PIE => MAPL_PI, & ! - ALHE => MAPL_ALHL, & ! J/kg @15C ALHM => MAPL_ALHF, & ! J/kg TF => MAPL_TICE, & ! K RHOW => MAPL_RHOWTR ! kg/m^3 USE catch_constants, ONLY: MAPL_LANDICE public :: snowrt ! used by LandIce, Catchment public :: TRID ! used by LandIce public :: SNOW_ALBEDO ! used by LandIce, Catchment, and LDAS public :: StieglitzSnow_calc_asnow ! used by Catchment, and LDAS public :: StieglitzSnow_calc_tpsnow ! used by Catchment, and LDAS public :: StieglitzSnow_echo_constants ! used by LDAS public :: get_tf0d ! for now, to be unified w/ StieglitzSnow_calc_tpsnow, reichle, 12 Aug 2014 ! constants specific to StieglitzSnow real, private, parameter :: MINSWE = 0.013 ! kg/m^2 min SWE to avoid immediate melt real, public, parameter :: WEMIN = 26. ! kg/m^2 minimum SWE in areal fraction real, private, parameter :: cpw = 2065.22 ! @ 0 C [J/kg/K] real, private, parameter :: DZ1MAX = 0.08 ! m real, private, parameter :: RHOMA = 500. ! kg/m^3 maximum snow density real, private, parameter :: SNWALB_VISMIN = 0.5 real, private, parameter :: SNWALB_NIRMIN = 0.3 !+++ spectrally integrated values for VIS and NIR ++++++++++++++++++++ real, private, parameter :: RIV=0.018, RIN=0.017 !+++ so as to explain Abs. Co. for ice of 10 [m-1] +++++++++++++++++++ real, private, parameter :: AICEV=0.149, AICEN=19.851 real, private, parameter :: DENICE=917. contains subroutine snowrt(N_zones, N_snow, N_constit, tileType, & t1,area,tkgnd,precip,snowf,ts,dts, eturb,dedtc,hsturb,dhsdtc, & hlwtc,dhlwtc,desdtc,hlwout,raddn,zc1,totdepos,wss, & wesn,htsnn,sndz,fices,tpsn, rconstit, & areasc,areasc0,pre,fhgnd,evap,shflux,lhflux,hcorr,ghfluxsno & , sndzsc, wesnprec, sndzprec, sndz1perc & , wesnperc, wesndens, wesnrepar, mltwtr & , excs, drho0, wesnbot, tksno, dtss, maxsndepth, rhofs, targetthick) !********************************************************************* ! AUTHORS: M. Stieglitz, M. Suarez, R. Koster & S. Dery. ! VERSION: 2003b - This version last updated: 05/30/03. !********* ! INPUTS: !********* ! N_zones : number of zones in the horizontal dimension (eg, 3 for Catchment, 1 for LandIce) ! N_snow : number of snow layers ! N_constit : ! tileType : surface type of the tile ! t1 : Temperature of catchment zones [C] ! ts : Air temperature [K] ! area : Fraction of snow-free area in each catchment zone [0-1] ! precip : Precipitation (Rain+snowfall) [kg/m^2/s == mm/s] ! snowf : Snowfall per unit area of catchment [kg/m^2/s == mm/s] ! dts : Time step [s] ! eturb : Evaporation per unit area of snow [kg/m^2/s == mm/s] ! dedtc : d(eturb)/d(ts) [kg/m^2/s/K] ! hsturb : Sensible heat flux per unit area of snow [W/m^2] ! dhsdtc : d(hsturb)/d(ts) [W/m^2/K] ! hlwtc : Emitted IR per unit area of snow [W/m^2] ! dhlwtc : d(hlwtc)/d(ts) [W/m^2/K] ! raddn : Net solar + incident terrestrial per unit area of snow [W/m^2] ! tkgnd : Thermal diffusivity of soil in catchment zones [W/m/K] ! zc1 : Half-thickness of top soil layer [m] !*** Bin Zhao added ************************************************* ! maxsndepth : Maximum snow depth beyond which snow gets thrown away ! rhofs : fresh snow density ! targetthick : the target thickness distribution relayer redistribute mass ! and energy to; currently its value is surface type dependent ! for catchment, the 1st array element the target thickness ! the rest define a sigma distribution; ! for landice, it is an array with specified thicknesses !********* ! UPDATES: !********* ! wesn : Layer water contents per unit area of catchment [kg/m^2] ! htsnn : Layer heat contents relative to liquid water at 0 C [J/m^2] ! sndz : Layer depths [m] ! rconstit : !********* ! OUTPUTS: !********* ! tpsn : Layer temperatures [C] ! fices : Layer frozen fraction [0-1] ! areasc : Areal snow coverage at beginning of step [0-1] ! areasc0: Areal snow coverage at end of step [0-1] ! pre : Liquid water flow from snow base [kg/m^2/s] ! fhgnd : Heat flux at snow base at catchment zones [W/m^2] ! hlwout : Final emitted IR flux per unit area of snow [W/m^2] ! lhflux : Final latent heat flux per unit area of snow [W/m^2] ! shflux : Final sensible heat flux per unit area of snow [W/m^2] ! evap : Final evaporation per unit area of snow [kg/m^2/s] !*** Bin Zhao added ************************************************* ! sndzsc : top layer thickness change due to sublimation/condensation ! wesnprec : top layer water content change due to precip (different from precip itself) ! sndzprec : top layer thickness change due to precip ! sndz1perc : top layer thickness change due to percolation ! wesnperc : layer water content change due to percolation ! wesndens : layer water content change due to densification ! wesnrepar : layer water content change due to relayer ! mltwtr : total melt water production rate ! excs : frozen part of water content from densification excess ! drho0 : layer density change due to densification ! wesnbot : excessive water content due to thickness exceeding maximum depth ! tksno : layer conductivity ! dtss : top layer temperature change !********************************************************************* ! NOTA: By convention, wesn is representative for a catchment area ! equal to 1 whereas sndz is relative to the area covered by snow only. !********************************************************************* implicit none ! real, parameter :: lhv = 2.4548E6 ! 2.5008e6 ! @ 0 C [J/kg] ! real, parameter :: lhs = 2.8368E6 ! 2.8434e6 ! @ 0 C [J/kg] ! real, parameter :: lhf = (lhs-lhv) ! @ 0 C [J/kg] !rr real, parameter :: cpw_liquid = 4185. ! [J/kg/K] ! real, parameter :: tfrz = 273.16 ! @ 0 C [K] ! real, parameter :: rhofs = 150. ! [kg/m^3] ! real, parameter :: rhoma = 500. ! [kg/m^3] ! real, parameter :: rhow = 1000. ! [kg/m^3] ! real, parameter :: wemin = 13. ! [kg/m^2] real, parameter :: snfr = 0.01 ! holding capacity real, parameter :: small = 1.e-6 ! small number ! integer, parameter :: nlay = 3 ! number of layers ! integer, parameter :: N_zones = 3 ! number of zones ! real , parameter :: MIN_SNOW_MASS = .013 ! kg/M**2 equiv to 0.1% area integer, intent(in) :: N_zones, N_snow, N_constit, tileType real, intent(in ) :: t1(N_zones),area(N_zones),tkgnd(N_zones) real, intent(in) :: totdepos(N_constit) real, intent(in ) :: ts,precip,snowf,dts,dedtc,raddn,hlwtc real, intent(in ) :: dhsdtc,desdtc,dhlwtc,eturb,hsturb,zc1,wss real, intent(inout):: wesn(N_snow),htsnn(N_snow),sndz(N_snow) real, intent(inout):: rconstit(N_snow,N_constit) real, intent(out) :: tpsn(N_snow),fices(N_snow),fhgnd(N_zones) real, intent(out) :: hlwout,lhflux,shflux,areasc0,evap,areasc,pre real, intent(out) :: ghfluxsno real, intent(out) :: wesnprec !real, intent(out) :: wesnsc, wesnprec real, intent(out) :: sndzsc, sndzprec real, intent(out) :: sndz1perc real, intent(out) :: wesnperc(N_snow) real, intent(out) :: wesndens(N_snow) real, intent(out) :: wesnrepar(N_snow) real, intent(out) :: mltwtr real, intent(out) :: excs(N_snow) real, intent(out) :: drho0(N_snow) real, intent(out) :: wesnbot real, intent(out) :: tksno(N_snow) real, intent(out) :: dtss real, intent(in) :: maxsndepth real, intent(in) :: rhofs real, intent(in) :: targetthick(N_snow) !Locals real :: tsx, mass,snowd,rainf,denom,alhv,lhturb,dlhdtc,hcorr, & enew,eold,tdum,fnew,tnew,icedens,densfac,hnew,scale,t1ave, & flxnet,fdum,dw,waterin,waterout,snowin,snowout, mtwt, & waterbal,precision,flow,term,dz,w(0:N_snow),HTSPRIME real :: excsdz, excswe, sndzsum, melti, mtwt0, mtwt1 real, dimension(size(wesn) ) :: cmpc,dens real, dimension(size(wesn) ) :: tksn real, dimension(size(wesn) ) :: dtc,q,cl,cd,cr real, dimension(size(wesn)+1) :: fhsn,df real, dimension(size(wesn) ) :: htest,ttest,ftest logical, dimension(size(wesn) ) :: ice1,tzero, ice10,tzero0 real :: topthick real, dimension(size(wesn)-1) :: thickdist !real, dimension(size(wesn) ) :: wesn0, wesn1, wesn2 !real, dimension(size(wesn) ) :: sndz0, sndz1, sndz2 real, dimension(size(wesn) ) :: dens0 integer :: i,izone logical :: logdum,kflag snowd = sum(wesn) snowin = snowd ghfluxsno = 0. !rr correction for "cold" snow tsx = min(ts-tf,0.)*cpw !rr correction for heat content of rain !rr tsx_rain = max(ts-tf,0.)*cpw_liquid df = 0. dtc = 0. tpsn = 0. fices = 0. areasc = 0. areasc0= 0. pre = 0. fhgnd = 0. hlwout = 0. shflux = 0. lhflux = 0. evap = 0. excs = 0. hcorr = 0. dens = rhofs rainf = precip - snowf ! [kg/m^2/s] !wesnsc = 0. sndzsc = 0. wesnprec = 0. sndzprec = 0. sndz1perc = 0. wesnperc = 0. wesndens = 0. wesnrepar = 0. wesnbot = 0. !tksno = condice dtss = 0. excswe = 0. !wesn0 = wesn !sndz0 = sndz if(snowd <= MINSWE) then ! no snow ! Assume initial (very small) snow water melts; new area is based ! on new snowfall areasc = min(snowd/wemin,1.) areasc0 = 0. pre = snowd/dts + areasc*rainf wesn = 0. hcorr = hcorr + raddn*areasc + sum(htsnn)/dts htsnn = 0. sndz = 0. mltwtr = snowd/dts if(snowf > 0.) then ! only initialize with non-liquid part of precip ! liquid part runs off (above) wesn = snowf*dts/float(N_snow) htsnn = (tsx-alhm)*wesn areasc0 = min((snowf*dts)/wemin,1.) !sndz = wesn/rhofs !*** should have fractional snow cover taken into account sndz = wesn/(max(areasc0,small)*rhofs) ! hcorr = hcorr - (tsx-alhm)*snowf ! randy hcorr = hcorr - tsx*snowf ! randy !call relayer(N_snow, htsnn, wesn, sndz) select case (tileType) case (MAPL_LANDICE) call FindTargetThickDist(N_snow, sndz, targetthick, topthick, thickdist) case default topthick = targetthick(1) thickdist = targetthick(2:N_snow) end select call relayer2(N_snow, N_constit, topthick, thickdist, & htsnn, wesn, sndz, rconstit) call StieglitzSnow_calc_tpsnow(N_snow, htsnn, wesn, tpsn, fices) endif return ! if there was no snow at start of time step endif !**** Determine the fractional snow coverage areasc = min(snowd/wemin,1.) !**** Set the mean density & diffusivity of the layers do i=1,N_snow if(sndz(i) > 0) dens(i) = max(wesn(i)/(areasc*sndz(i)),rhofs) enddo tksn = 3.2217e-06*dens**2 tksno = tksn dens0 = dens !**** Determine temperature & frozen fraction of snow layers call StieglitzSnow_calc_tpsnow(N_snow, htsnn, wesn, tpsn, fices) mtwt = sum(wesn*(1.-fices)) !**** Calculate the ground-snow energy flux at 3 zones denom = 1./(sndz(N_snow)*0.5-zc1) fhgnd = -sqrt(tkgnd*tksn(N_snow))*area*denom*(tpsn(N_snow)-t1) fhsn(N_snow+1) = sum(fhgnd) do i=1,N_zones df(N_snow+1)=df(N_snow+1)-sqrt(tkgnd(i)*tksn(N_snow))*area(i)*denom enddo !**** Ensure against excessive heat flux between ground and snow: !**** if heat flux found to cause the lowest snow layer temperature !**** to "overshoot" (e.g. to become higher than the ground temperature !**** when it had been lower), reduce the heat flux. If the lowest !**** snow layer starts off at zero and the new temperature is greater !**** than zero, reduce the heat flux to melt only half of the lowest !**** layer snow. !**** t1ave=sum(t1*area)/sum(area) htest=htsnn htest(N_snow)=htest(N_snow)+fhsn(N_snow+1)*dts*areasc call StieglitzSnow_calc_tpsnow(N_snow, htest, wesn, ttest, ftest) scale=1. if((t1ave-tpsn(N_snow))*(t1ave-ttest(N_snow)) .lt. 0.) then scale=0.5*(tpsn(N_snow)-t1ave)/(tpsn(N_snow)-ttest(N_snow)) endif if(tpsn(N_snow) .eq. 0. .and. ttest(N_snow) .gt. 0. .and. & abs(fhsn(N_snow+1)) .gt. 1.e-10) then scale=(-0.5*htsnn(N_snow)/(dts*areasc))/fhsn(N_snow+1) endif fhsn(N_snow+1)=fhsn(N_snow+1)*scale df(N_snow+1)=df(N_snow+1)*scale fhgnd=fhgnd*scale !**** Calculate heat fluxes between snow layers. do i=2,N_snow df(i) = -sqrt(tksn(i-1)*tksn(i))/((sndz(i-1)+sndz(i))*0.5) fhsn(i)= df(i)*(tpsn(i-1)-tpsn(i)) enddo ghfluxsno = fhsn(2) !**** Effective heat of vaporization includes bringing snow to 0 C alhv = alhe + alhm !randy ! alhv = alhe + fices(1)*alhm + tpsn(1)*cpw !randy !**** Initial estimate of latent heat flux change with Tc lhturb = alhv*eturb dlhdtc = alhv*dedtc !**** Initial estimate of net surface flux & its change with Tc fhsn(1) = lhturb + hsturb + hlwtc - raddn df(1) = -(dlhdtc + dhsdtc + dhlwtc) !**** Prepare array elements for solution & coefficient matrices. !**** Terms are as follows: left (cl), central (cd) & right (cr) !**** diagonal terms in coefficient matrix & solution (q) terms. do i=1,N_snow call get_tf0d(htsnn(i),wesn(i),tdum,fdum, ice1(i),tzero(i)) if(ice1(i)) then cl(i) = df(i) cd(i) = cpw*wesn(i)/dts - df(i) - df(i+1) cr(i) = df(i+1) q(i) = fhsn(i+1)-fhsn(i) else cl(i) = 0. cd(i) = 1. cr(i) = 0. q(i) = 0. endif enddo cl(1) = 0. cr(N_snow) = 0. do i=1,N_snow-1 if(.not.ice1(i)) cl(i+1) = 0. enddo do i=2,N_snow if(.not.ice1(i)) cr(i-1) = 0. enddo !**** Solve the tri-diagonal matrix for implicit change in Tc. call TRID(dtc,cl,cd,cr,q,N_snow) !**** Check temperature changes for passages across critical points,i.e. !**** If implicit change has taken layer past melting/freezing, correct. do i=1,N_snow if(tpsn(i)+dtc(i) > 0. .or. htsnn(i)+wesn(i)*cpw*dtc(i) > 0.) then dtc(i)=-tpsn(i) endif if(.not.ice1(i)) dtc(i)=0. enddo !**** Further adjustments; compute new values of h associated with !**** all adjustments. eold=sum(htsnn) do i=1,N_snow !**** Quick check for "impossible" condition: if(.not.tzero(i) .and. .not.ice1(i)) then write(*,*) 'bad snow condition: fice,tpsn =',fices(i),tpsn(i) stop endif !**** Condition 1: layer starts fully frozen (temp < 0.) if(.not.tzero(i)) then tnew=tpsn(i)+dtc(i) fnew=1. endif !**** Condition 2: layer starts with temp = 0, fices < 1. ! Corrections for flxnet calculation: Koster, March 18, 2003. if(.not.ice1(i)) then tnew=0. if(i==1) flxnet= fhsn(i+1)+df(i+1)*(dtc(i)-dtc(i+1)) & -fhsn(i)-df(i)*dtc(i) if(i > 1 .and. i < N_snow) flxnet= & fhsn(i+1)+df(i+1)*(dtc(i)-dtc(i+1)) & -fhsn(i)-df(i)*(dtc(i-1)-dtc(i)) if(i==N_snow) flxnet=fhsn(i+1)+df(i+1)*dtc(i) & -fhsn(i)-df(i)*(dtc(i-1)-dtc(i)) HTSPRIME=HTSNN(I)+AREASC*FLXNET*DTS call get_tf0d(HTSPRIME,wesn(i), tdum,fnew,logdum,logdum) fnew=amax1(0., amin1(1., fnew)) endif !**** Condition 3: layer starts with temp = 0, fices = 1. ! Corrections for flxnet calculation: Koster, March 18, 2003. if(ice1(i) .and. tzero(i)) then if(dtc(i) < 0.) then tnew=tpsn(i)+dtc(i) fnew=1. endif if(dtc(i) >= 0.) then tnew=0. if(i==1) flxnet=fhsn(i+1)+df(i+1)*(dtc(i)-dtc(i+1)) & -fhsn(i)-df(i)*dtc(i) if(i > 1 .and. i < N_snow) flxnet= & fhsn(i+1)+df(i+1)*(dtc(i)-dtc(i+1)) & -fhsn(i)-df(i)*(dtc(i-1)-dtc(i)) if(i==N_snow) flxnet=fhsn(i+1)+df(i+1)*dtc(i) & -fhsn(i)-df(i)*(dtc(i-1)-dtc(i)) HTSPRIME=HTSNN(I)+AREASC*FLXNET*DTS call get_tf0d(HTSPRIME,wesn(i), tdum,fnew,logdum,logdum) fnew=amax1(0., amin1(1., fnew)) endif endif !**** Now update heat fluxes & compute sublimation or deposition. if(i == 1) then dtss = dtc(1) lhflux = lhturb + dlhdtc*dtc(1) shflux = hsturb + dhsdtc*dtc(1) hlwout = hlwtc + dhlwtc*dtc(1) evap = lhflux/alhv dw = -evap*dts*areasc if(-dw > wesn(1) ) then dw = -wesn(1) evap = -dw/(dts*areasc) ! shflux=shflux+(lhflux-evap*alhv) hcorr=hcorr+(lhflux-evap*alhv)*areasc lhflux=evap*alhv endif wesn(1) = wesn(1) + dw denom = 1./dens(1) if(dw > 0.) denom = 1./rhoma sndz(1) = sndz(1) + dw*denom !wesnsc = dw sndzsc = dw*denom endif if(i == N_snow) then do izone=1,N_zones fhgnd(izone)=fhgnd(izone)+area(izone)*df(N_snow+1)*dtc(N_snow) enddo endif !**** Now update thermodynamic quantities. htsnn(i)=(cpw*tnew-fnew*alhm)*wesn(i) tpsn(i) = tnew fices(i)= fnew enddo !**** Store excess heat in hcorr. enew=sum(htsnn) hcorr=hcorr-((enew-eold)/dts+areasc*(lhflux+shflux+hlwout-raddn) & -areasc*(fhsn(N_snow+1)+df(N_snow+1)*dtc(N_snow)) & ) call StieglitzSnow_calc_tpsnow(N_snow, htsnn, wesn, tpsn, fices) mltwtr = max(0., sum(wesn*(1.-fices)) - mtwt) mltwtr = mltwtr / dts !mtwt0 = sum(wesn*(1.-fices)) + areasc*rainf*dts mtwt0 = sum(wesn*fices) !rr!**** Add rainwater and snow at ts., bal. budget with shflux. !rr (tried and failed 19 Jun 2003, reichle) !rr !rr wesn (1) = wesn (1) + (rainf*areasc+snowf)*dts !rr htsnn(1) = htsnn(1) + (tsx -alhm)*(snowf*dts) + tsx_rain*rainf*dts !rr sndz (1) = sndz (1) + (snowf/rhofs)*dts !rr ! shflux = shflux + tsx*snowf ! randy !rr hcorr = hcorr - (tsx-alhm)*snowf - tsx_rain*rainf ! randy !**** Add rainwater at 0 C, snow at ts., bal. budget with shflux. wesn (1) = wesn (1) + (rainf*areasc+snowf)*dts htsnn(1) = htsnn(1) + (tsx -alhm)*(snowf*dts) sndz (1) = sndz (1) + (snowf/rhofs)*dts wesnprec = (rainf*areasc+snowf)*dts sndzprec = (snowf/rhofs)*dts ! shflux = shflux + tsx*snowf ! randy ! hcorr = hcorr - (tsx-alhm)*snowf ! randy hcorr = hcorr - tsx*snowf ! randy snowd=sum(wesn) call StieglitzSnow_calc_tpsnow(N_snow, htsnn, wesn, tpsn, fices) !**** Move meltwater through the pack. !**** Updated by Koster, August 27, 2002. pre = 0. flow = 0. wesnperc = wesn do i=1,N_snow if(flow > 0.) then wesn (i) = wesn(i) + flow call StieglitzSnow_calc_tpsnow(N_snow, htsnn, wesn, tpsn, fices) endif pre = max((1.-fices(i))*wesn(i), 0.) flow = 0. if(snowd > wemin) then icedens=wesn(i)*fices(i)/(sndz(i)+1.e-20) densfac=amax1(0., amin1(1., icedens/rhofs)) term=densfac*snfr*(sndz(i)*rhow-wesn(i)*fices(i)) if(pre > term) then pre = min(pre - term, wesn(i)) wesn(i) = wesn(i) - pre flow = pre endif else wesn(i) = wesn(i) - pre flow = pre endif !**** Adjust top layer snow depth to get proper density values !**** But limit this change for large throughflow (STEPH 06/19/03) if(i==1)then dz=min(flow/dens(i),0.5*sndz(i)) sndz(i)=sndz(i)-dz sndz1perc = -dz endif enddo wesnperc = wesn - wesnperc pre = flow/dts snowd=sum(wesn) !**** Update snow density by compaction (Pitman et al. 1991) excs = 0. mass = 0. w = 0. drho0 = 0. wesndens = wesn if(snowd > wemin) then ! Compaction only after full coverage. do i=1,N_snow dens(i) = rhofs if(sndz(i)>0.) dens(i) = max(wesn(i)/(sndz(i)),rhofs) enddo drho0 = dens cmpc = exp(14.643 - (4000./min(tpsn+tf,tf))-.02*dens) do i=1,N_snow w(i) = wesn(i) mass = mass + 0.5*(w(i)+w(i-1)) dens(i) = dens(i)*(1. + (dts*0.5e-7*9.81)*mass*cmpc(i)) !**** Clip densities below maximum value, adjust quantities accordingly !**** while conserving heat & mass (STEPH 06/21/03). if(dens(i) > rhoma) then excs(i) = (dens(i)-rhoma)*sndz(i) wesn(i) = wesn(i) - excs(i) hnew = (cpw*tpsn(i)-fices(i)*alhm)*wesn(i) hcorr= hcorr+(htsnn(i)-hnew)/dts htsnn(i)= hnew dens(i) = rhoma endif enddo drho0 = dens - drho0 endif wesndens = wesn - wesndens !pre = pre + sum(excs)/dts pre = pre + sum(excs*max(1.-fices,0.0))/dts excs = excs * fices / dts snowd=sum(wesn) areasc0 = max(small, min(snowd/wemin,1.) ) sndz = (wesn/areasc0)/dens sndzsum = sum(sndz) if(sndzsum > maxsndepth) then excsdz = sndzsum - maxsndepth excswe = dens(N_snow) * excsdz wesn(N_snow) = wesn(N_snow) - excswe hnew = (cpw*tpsn(N_snow)-fices(N_snow)*alhm)*wesn(N_snow) htsnn(N_snow)= hnew sndz(N_snow) = sndz(N_snow) - excsdz wesnbot = excswe endif !**** Restore layers to sigma values. wesnrepar = wesn do i=1,N_snow call get_tf0d(htsnn(i),wesn(i),tdum,fdum,ice10(i),tzero0(i)) enddo !sndz1 = sndz !call relayer(N_snow, htsnn, wesn, sndz) select case (tileType) case (MAPL_LANDICE) call FindTargetThickDist(N_snow, sndz, targetthick, topthick, thickdist) case default topthick = targetthick(1) thickdist = targetthick(2:N_snow) end select !wesn1 = wesn call relayer2(N_snow, N_constit, topthick, thickdist, & htsnn, wesn, sndz, rconstit) !wesn2 = wesn !sndz2 = sndz wesnrepar = wesn - wesnrepar call StieglitzSnow_calc_tpsnow(N_snow, htsnn, wesn, tpsn, fices) !**** Check that (ice10,tzero) conditions are conserved through !**** relayering process (or at least that (fices,tpsn) conditions don't !**** go through the (1,0) point); excess goes to hcorr. do i=1,N_snow kflag=.false. if(ice10(i).and.tzero0(i) .and. & (fices(i) .ne. 1. .or. tpsn(i) .ne. 0.) ) kflag=.true. if(.not.ice10(i).and.tzero0(i) .and. & (fices(i) .eq. 1. .and. tpsn(i) .lt. 0.) ) kflag=.true. if(ice10(i).and. .not.tzero0(i) .and. & (fices(i) .ne. 1. .and. tpsn(i) .eq. 0.) ) kflag=.true. if(kflag) then hnew=-alhm*wesn(i) hcorr=hcorr+(htsnn(i)-hnew)/dts htsnn(i)=hnew tpsn(i)=0. fices(i)=1. endif enddo !**** Reset fractional area coverage. areasc0 = min(sum(wesn)/wemin,1.) !**** Final check for water balance. waterin = (rainf*areasc+snowf)*dts + max(dw,0.) waterout = pre*dts - min(dw,0.) snowout = sum(wesn) + sum(excs) + excswe waterbal = snowin + waterin - waterout - snowout precision = snowout*small #if 0 if((waterbal > precision).and.(waterbal > small) .or. pre < -1.e-13 ) then write(*,*) 'Warning: Imbalance in snow water budget!', waterbal write(*,*) 'waterin = ', waterin write(*,*) 'snowin = ', snowin write(*,*) 'waterout = ', waterout write(*,*) 'snowout = ', snowout write(*,*) 'dw = ', dw write(*,*) 'excswe = ', excswe write(*,*) 'sum(excs) = ', sum(excs) write(*,*) 'snowf*dts = ', snowf*dts write(*,*) 'sum(wesn) = ', sum(wesn) write(*,*) (wesn(i), i=1,N_snow) write(*,*) 'sum(sndz) = ', sum(sndz) write(*,*) (sndz(i), i=1,N_snow) write(*,*) 'dens0 = ' write(*,*) (dens0(i), i=1,N_snow) !write(*,*) 'sum(wesn0) = ', sum(wesn0) !write(*,*) (wesn0(i), i=1,N_snow) write(*,*) 'sum(wesn1) = ', sum(wesn1) write(*,*) (wesn1(i), i=1,N_snow) write(*,*) 'sum(wesn2) = ', sum(wesn2) write(*,*) (wesn2(i), i=1,N_snow) !write(*,*) 'sum(sndz0) = ', sum(sndz0) !write(*,*) (sndz0(i), i=1,N_snow) write(*,*) 'sum(sndz1) = ', sum(sndz1) write(*,*) (sndz1(i), i=1,N_snow) write(*,*) 'sum(sndz2) = ', sum(sndz2) write(*,*) (sndz2(i), i=1,N_snow) !stop endif #endif return ! end snow end subroutine snowrt ! ********************************************************************** subroutine FindTargetThickDist(N_snow, sndz, dzmax, topthick, thickdist) integer, intent(in) :: N_snow real, intent(in) :: sndz(N_snow) real, intent(in) :: dzmax(N_snow) real, intent(out) :: topthick real, intent(out), dimension(N_snow-1) :: thickdist real, dimension(N_snow) :: sndzt real :: totald, dzdiff, restthick integer :: i integer, dimension(N_snow) :: mark logical :: lth_satisfy totald = sum(sndz) sndzt = totald/float(N_snow) if(sndzt(1) < dzmax(1)) then topthick = dzmax(1) do i=2,N_snow thickdist(i-1) = 1.0/real(N_snow-1,kind=4) enddo else mark = 0 do lth_satisfy = .true. do i=1,N_snow if(mark(i) == 0 .and. sndzt(i) > dzmax(i)) then sndzt(i) = dzmax(i) mark(i) = 1 lth_satisfy = .false. endif enddo if(lth_satisfy) exit dzdiff = 0.0 do i=1,N_snow if(mark(i) == 1) then dzdiff = dzdiff + sndzt(i) endif enddo restthick = (totald-dzdiff)/float(N_snow-sum(mark)) do i=1,N_snow if(mark(i) == 0) then sndzt(i) = restthick endif enddo enddo topthick = sndzt(1) totald = totald - topthick do i=2,N_snow thickdist(i-1) = sndzt(i)/totald enddo endif return end subroutine FindTargetThickDist ! ********************************************************************** subroutine relayer2(N_snow, N_constit, thick_toplayer, thickdist, & htsnn, wesn, sndz, rconstit) implicit none integer, intent(in) :: N_snow, N_constit real, intent(in) :: thick_toplayer real, intent(in), dimension(N_snow-1) :: thickdist real, intent(inout) :: htsnn(N_snow),wesn(N_snow),sndz(N_snow) real, intent(inout) :: rconstit(N_snow,N_constit) !integer, intent(in) :: pid, ktile real, dimension(size(sndz),2+N_Constit) :: h, s integer :: i, j, k, ilow, ihigh ! real, parameter :: dz1max = 0.08 ! [m] ! real, parameter :: wemin = 13.0 ! [kg/m2] !real, parameter :: small = 1.e-20 !real :: areasc,dz real :: areasc,dz real :: totalthick real, dimension(size(sndz)) :: tol_old, bol_old, tol_new, & bol_new real, dimension(size(sndz)) :: thickness !**** thick_toplayer: the assigned (final) thickness of the topmost layer (m) !**** thickdist: the assigned (final) distribution of thickness in layers !**** 2 through N_snow, in terms of fraction !**** h: array holding specific heat, water, and constituent contents !**** s: array holding the total and final heat, water, and constit. contents !**** ilow: first layer used in a particular relayering calculation !**** ihigh: final layer used in a particular relayering calculation !**** totalthick: total thickness of layers 2 through N_snow !**** thickness: array holding final thicknesses (m) of the snow layers !**** tol_old(i): depth (from surface) of the top of layer i, before & !**** relayering !**** bol_old(i): depth (from surface) of the bottom of layer i, before & !**** relayering !**** tol_old(i): depth (from surface) of the top of layer i, after & !**** relayering !**** bol_old(i): depth (from surface) of the bottom of layer i, after & !**** relayering !thickness(1)=thick_toplayer !totalthick=sum(sndz)-thick_toplayer !if(sum(sndz) .lt. thick_toplayer) & ! write(*,*) 'Snow layer thickness inconsistency' totalthick=sum(sndz) thickness(1)=amin1(totalthick*0.9, thick_toplayer) totalthick=totalthick-thickness(1) do i=1,N_snow-1 thickness(i+1)=thickdist(i)*totalthick enddo !**** Initialize some variables. h = 0. s = 0. dz = 0. !areasc = min(sum(wesn)/wemin,1.) !**** Compute specific heat & water contents of old layers. do i=1,N_snow if (sndz(i) > 0.) then h(i,1) = htsnn(i)/sndz(i) h(i,2) = wesn(i)/sndz(i) do k=1,N_Constit h(i,2+k)=rconstit(i,k)/sndz(i) enddo endif enddo !**** Determine old and new boundary depths (cumulative from top) !**** (tol refers to "top of layer", bol refers to "bottom of layer" tol_old(1)=0. bol_old(1)=sndz(1) tol_new(1)=0. bol_new(1)=thickness(1) do i=2,N_snow tol_old(i)=bol_old(i-1) bol_old(i)=bol_old(i-1)+sndz(i) tol_new(i)=bol_new(i-1) bol_new(i)=bol_new(i-1)+thickness(i) enddo !**** Redistribute quantities !**** Step 1: Do top layer ihigh=1 do k=1,N_snow if(bol_old(k) .lt. bol_new(1)) ihigh=k+1 enddo do k=1,ihigh if(k .lt. ihigh) dz=sndz(k) if(k .eq. ihigh) dz=bol_new(1)-tol_old(k) !s(1,:)=s(1,:)+h(1,:)*dz s(1,:)=s(1,:)+h(k,:)*dz enddo !**** Step 2: Do remaining layers do i=2,N_snow ilow=ihigh do k=ilow,N_snow if(bol_old(k) .lt. bol_new(i)) ihigh=k+1 enddo if(ihigh .eq. N_snow+1) ihigh=N_snow ! Account for potential truncation problem do k=ilow,ihigh if(k .eq. ilow .and. k .lt. ihigh) dz=bol_old(k)-tol_new(i) if(k .eq. ilow .and. k .eq. ihigh) dz=bol_new(i)-tol_new(i) if(k .gt. ilow .and. k .lt. ihigh) dz=bol_old(k)-tol_old(k) if(k .gt. ilow .and. k .eq. ihigh) dz=bol_new(i)-tol_old(k) s(i,:)=s(i,:)+h(k,:)*dz enddo enddo htsnn = s(:,1) wesn = s(:,2) do k=1,N_Constit rconstit(:,k)=s(:,2+k) enddo sndz=thickness ! if(sum(wesn) < wemin) sndz = sndz /(areasc + small) return end subroutine relayer2 ! ********************************************************************** subroutine get_tf0d(h,w,t,f,ice1,tzero) implicit none !RR real, parameter :: cpw = 2065.22 ! @ 0 C [J/kg/K] -- ALREADY DEFINED ABOVE ! real, parameter :: lhv = 2.4548E6 ! 2.5008e6 ! @ 0 C [J/kg] ! real, parameter :: lhs = 2.8368E6 ! 2.8434e6 ! @ 0 C [J/kg] !rr real, parameter :: lhv = 2.5008e6 ! @ 0 C [J/kg] !rr real, parameter :: lhs = 2.8434e6 ! @ 0 C [J/kg] ! real, parameter :: lhf = (lhs-lhv) ! @ 0 C [J/kg] real, parameter :: tfac=1./cpw real, parameter :: ffac=1./alhm real, intent(in ) :: w, h real, intent(out) :: t, f logical, intent(out) :: ice1,tzero real :: hbw hbw=0. if(w > 0.) hbw = h/w if(hbw < -1.00001*alhm) then t = (hbw+alhm)*tfac f = 1. ice1=.true. tzero=.false. elseif(hbw > -0.99999*alhm) then t = 0. f =-hbw*ffac ice1=.false. tzero=.true. else t = 0. f = 1. ice1=.true. tzero=.true. endif if(f < 0.) then t = hbw*tfac f = 0. endif if(w == 0.) then t = 0. f = 0. endif return end subroutine get_tf0d ! ********************************************************************** subroutine StieglitzSnow_calc_tpsnow(N,h,w,t,f) ! renamed for clarity: get_tf_nd() --> StieglitzSnow_calc_tpsnow() ! reichle, 12 Aug 2014 ! n-dimensional version of get_tf ! ! avoid slow "where" statements ! ! can be called for any number of layers or catchments, for example ! 1.) call StieglitzSnow_calc_tpsnow( ncatm, htsnn1(1:ncatm), wesn1(1:ncatm), ! tpsn(1:ncatm),f(1:ncatm) ) ! ! 2.) call StieglitzSnow_calc_tpsnow(N_snow, h, w, t, f) ! reichle, 22 Aug 2002 ! reichle, 29 Apr 2003 (updated parameter values) integer, intent(in) :: N real, dimension(n), intent(in) :: h, w real, dimension(n), intent(out) :: t, f ! local variables !RR real, parameter :: cpw = 2065.22 ! @ 0 C [J/kg/K] -- ALREADY DEFINED ABOVE ! real, parameter :: lhv = 2.4548E6 ! 2.5008e6 ! @ 0 C [J/kg] ! real, parameter :: lhs = 2.8368E6 ! 2.8434e6 ! @ 0 C [J/kg] !rr real, parameter :: lhv = 2.5008e6 ! @ 0 C [J/kg] !rr real, parameter :: lhs = 2.8434e6 ! @ 0 C [J/kg] ! real, parameter :: lhf = (lhs-lhv) ! @ 0 C [J/kg] real, parameter :: tfac=1./cpw real, parameter :: ffac=1./alhm integer :: i real :: hbw do i=1,N if(w(i) .gt. 0.0) then hbw = h(i)/w(i) else hbw = 0. endif if(hbw .lt. -alhm) then t(i) = (hbw+alhm)*tfac f(i) = 1. elseif(hbw .gt. -alhm) then t(i) = 0. f(i) = -hbw*ffac else t(i) = 0. f(i) = 1. endif if(f(i) .lt. 0.) then t(i) = hbw*tfac f(i) = 0. endif if(w(i) .eq. 0.) then t(i) = 0. f(i) = 0. endif end do return end subroutine StieglitzSnow_calc_tpsnow ! ******************************************************************** subroutine StieglitzSnow_calc_asnow( N_snow, NTILES, wesnn, asnow ) ! Calculate diagnostic snow area from prognostic SWE ! ! reichle, Nov 3, 2004 ! reichle, 2 Apr 2012 - revised for use without catch_types structures ! reichle, 12 Aug 2014 - moved to here from catchment.F90 ! [formerly known as catch_calc_asnow()] ! ! ---------------------------------------------------------------- implicit none integer, intent(in) :: N_snow, NTILES real, dimension(N_snow,NTILES), intent(in) :: wesnn real, dimension( NTILES), intent(out) :: asnow ! ----------------------------------------------------------- asnow = min( sum(wesnn,1)/wemin, 1. ) end subroutine StieglitzSnow_calc_asnow ! ********************************************************************** SUBROUTINE TRID(X,DD,D,RD,B,N) IMPLICIT NONE INTEGER,INTENT(IN) :: N REAL*4, INTENT(IN), DIMENSION(N) :: DD, RD REAL*4, INTENT(INOUT), DIMENSION(N) :: D, B REAL*4, INTENT(OUT),DIMENSION(N) :: X integer I,J real*4 RSF RSF=0. DO 10 I=2,N J=N+1-I if(D(J+1).ne.0.) RSF=RD(J)/D(J+1) D(J)=D(J)-DD(J+1)*RSF 10 B(J)=B(J)- B(J+1)*RSF if(D(1).ne.0.) X(1)=B(1)/D(1) DO 20 J=2,N 20 if(D(J).ne.0.) X(J)=(B(J)-DD(J)*X(J-1))/D(J) RETURN END SUBROUTINE TRID !======================================================================= ! Version 5.0.2 by Teppei J. Yasuanari on 02/14/2011 SUBROUTINE SNOW_ALBEDO (NCH, N_snow, N_constit, ITYP, VLAI, ZTH, & RHOFRESH, & SNWALB_VISMAX, SNWALB_NIRMAX, SLOPE, & WESN, HTSNN, SNDZ, & AVISDR, ANIRDR, AVISDF, ANIRDF, & ASNVDR, ASNNDR, ASNVDF, ASNNDF, & RCONSTIT, UM, RTS, PARDIR, PARDIF, & ABVIS, ABNIR & ) !USE sibalb_coeff !USE catch_constants IMPLICIT NONE INTEGER, INTENT(IN) :: NCH, N_snow, N_constit ! -------------------------------------- ! INPUTS: ! NCH: number of tiles considered ! N_snow: number of snow layers ! N_constit: number of constituents ! ITYP: vegetation type ! VLAI: the leaf area index. ! ZTH: The cosine of the solar zenith angle. ! RHOFRESH: density of fresh snow ! SNWALB_VISMAX: max of visible snow albedo ! SNWALB_NIRMAX: max of NIR snow albedo ! SLOPE: slope the albedo decreases higher density ! SLOPE > 0: it gets recomputed ! SLOPE < 0: it is directly used in the computation; ! for example,to recover the old formulation, ! set SLOPE=-0.0006, SNWALB_VISMAX=0.7, SNWALB_NIRMAX=0.5 ! WESN: snow water equivalent in each layer ! HTSNN: heat content of each layer ! SNDZ: depth of each layer ! UM: wind speed ! RTS: surface temperature ! PARDIR: photosynthetically active radiation, direct ! PARDIF: photosynthetically active radiation, diffuse ! RCONSTIT: array of constituent masses ! ABVIS: const array (not used when there is no constituents) ! ABNIR: const array (not used when there is no constituents) ! AVISDR: visible, direct albedo (snow-free). ! ANIRDR: near infra-red, direct albedo (snow-free). ! AVISDF: visible, diffuse albedo (snow-free). ! ANIRDF: near infra-red, diffuse albedo (snow-free). INTEGER, INTENT(IN), DIMENSION(NCH) :: ITYP REAL, INTENT(IN) :: RHOFRESH REAL, INTENT(IN) :: SNWALB_VISMAX, SNWALB_NIRMAX, SLOPE REAL, INTENT(IN), DIMENSION(NCH) :: AVISDR, ANIRDR, AVISDF, & ANIRDF, VLAI, ZTH REAL, INTENT(IN), DIMENSION(N_Snow,NCH) :: WESN, HTSNN, SNDZ REAL, INTENT(IN), DIMENSION(NCH), OPTIONAL :: UM, RTS, PARDIR, PARDIF REAL, INTENT(IN), DIMENSION(N_Constit), OPTIONAL :: ABVIS, ABNIR REAL, INTENT(IN),DIMENSION(N_snow,N_Constit,NCH), OPTIONAL :: RCONSTIT ! -------------------------------------- ! OUTPUTS: ! ASNVDR: snow albedo, visible direct ! ASNNDR: snow albedo, near-infrared direct ! ASNVDF: snow albedo, visible diffuse ! ASNNDF: snow albedo, near-infrared diffuse REAL, INTENT(OUT), DIMENSION(NCH) :: ASNVDR, ASNNDR, ASNVDF, ASNNDF ! -------------------------------------- ! Other variables as needed. Includes: ! SSA: snow specific surface area ! RHO_FS: fresh snow density ! EFFG: effective ice thickness (m) INTEGER :: I,M,J,K,K2 INTEGER, PARAMETER :: NTYPS_SIB=9 REAL :: rho_fs,DEGSZA,SD,SZASIN,COS50,SSALBV,SSALBN,AV,AN, & WSS, TS, FAC, FVEG, TOTDEP, SWE, DENS_EXC, AREASC, & DENSITY, ASNVDR_VEG, ASNNDR_VEG, ASNVDF_VEG, & ASNNDF_VEG, SUM1, SUM2, GK_B REAL, DIMENSION(NTYPS_SIB) :: SNWMID DATA SNWMID /50.,30.,45.,20.,30.,20.,2.,2.,2./ ! ********************************************************************* !FPP$ EXPAND (COEFFSIB) if(N_Constit > 0) then call ALB_WITH_IMPURITY (NCH, N_snow, N_constit, ZTH, & SNWALB_VISMAX, SNWALB_NIRMAX, & WESN, HTSNN, SNDZ, UM, RTS, PARDIR, PARDIF, & ABVIS, ABNIR, & ASNVDR, ASNNDR, ASNVDF, ASNNDF, & RCONSTIT & ) endif if(SLOPE < 0.0) then GK_B = SLOPE else GK_B = (0.85808-0.6)/(RHOFRESH-RHOMA) endif DO I=1,NCH SWE=SUM(WESN(:,I)) TOTDEP=SNDZ(1,I) AREASC = MIN(SWE/WEMIN,1.) !DENSITY=(SWE/(AREASC+1.e-20)) / (TOTDEP+1.e-20) !*** only use top layer density to dentermine albedo DENSITY=(WESN(1,I)/(AREASC+1.e-20)) / (TOTDEP+1.e-20) DENS_EXC=MAX(0., DENSITY-RHOFRESH) !********************************************************************* ! Use these when you use the original snow albedo model (comment out for alb) !ASNVDR(I) = MAX(SNWALB_VISMIN, SNWALB_VISMAX - DENS_EXC*.0006) !ASNNDR(I) = MAX(SNWALB_NIRMIN, SNWALB_NIRMAX - DENS_EXC*.0006) !ASNVDF(I) = MAX(SNWALB_VISMIN, SNWALB_VISMAX - DENS_EXC*.0006) !ASNNDF(I) = MAX(SNWALB_NIRMIN, SNWALB_NIRMAX - DENS_EXC*.0006) ASNVDR(I) = MAX(SNWALB_VISMIN, SNWALB_VISMAX + GK_B*DENS_EXC) ASNNDR(I) = MAX(SNWALB_NIRMIN, SNWALB_NIRMAX + GK_B*DENS_EXC) ASNVDF(I) = MAX(SNWALB_VISMIN, SNWALB_VISMAX + GK_B*DENS_EXC) ASNNDF(I) = MAX(SNWALB_NIRMIN, SNWALB_NIRMAX + GK_B*DENS_EXC) ! ACCOUNT FOR VEGETATION MASKING, FOR EACH COMPONENT ! A) FIRST DO MASKING IN VEGETATED FRACTION: FAC = SWE / (SWE + SNWMID(ITYP(I))) ASNVDR_VEG=AVISDR(I) + (ASNVDR(I)-AVISDR(I))*FAC ASNNDR_VEG=ANIRDR(I) + (ASNNDR(I)-ANIRDR(I))*FAC ASNVDF_VEG=AVISDF(I) + (ASNVDF(I)-AVISDF(I))*FAC ASNNDF_VEG=ANIRDF(I) + (ASNNDF(I)-ANIRDF(I))*FAC ! B) NOW ACCOUNT FOR SUBGRID VEGETATION FRACTION FVEG=AMIN1( 1., VLAI(I)/2. ) ASNVDR(I)=ASNVDR(I)*(1.-FVEG)+ASNVDR_VEG*FVEG ASNNDR(I)=ASNNDR(I)*(1.-FVEG)+ASNNDR_VEG*FVEG ASNVDF(I)=ASNVDF(I)*(1.-FVEG)+ASNVDF_VEG*FVEG ASNNDF(I)=ASNNDF(I)*(1.-FVEG)+ASNNDF_VEG*FVEG ENDDO RETURN END SUBROUTINE SNOW_ALBEDO SUBROUTINE ALB_WITH_IMPURITY (NCH, N_snow, N_constit, ZTH, & SNWALB_VISMAX, SNWALB_NIRMAX, & WESN, HTSNN, SNDZ, UM, RTS, PARDIR, PARDIF, & ABVIS, ABNIR, & ASNVDR, ASNNDR, ASNVDF, ASNNDF, & RCONSTIT & ) !USE sibalb_coeff !USE catch_constants IMPLICIT NONE INTEGER, INTENT(IN) :: NCH, N_snow, N_constit ! -------------------------------------- ! INPUTS: ! NCH: number of tiles considered ! ZTH: The cosine of the solar zenith angle. ! WESN: snow water equivalent in each layer ! HTSNN: heat content of each layer ! SNDZ: depth of each layer ! UM: wind speed ! RTS: surface temperature ! PARDIR: photosynthetically active radiation, direct ! PARDIF: photosynthetically active radiation, diffuse ! RCONSTIT: array of constituent masses ! AVISDR: visible, direct albedo (snow-free). ! ANIRDR: near infra-red, direct albedo (snow-free). ! AVISDF: visible, diffuse albedo (snow-free). ! ANIRDF: near infra-red, diffuse albedo (snow-free). REAL, INTENT(IN) :: SNWALB_VISMAX, SNWALB_NIRMAX REAL, INTENT(IN), DIMENSION(NCH) :: ZTH, UM, RTS, PARDIR, PARDIF REAL, INTENT(IN), DIMENSION(N_Constit) :: ABVIS, ABNIR REAL, INTENT(IN), DIMENSION(N_snow,N_Constit,NCH), OPTIONAL :: RCONSTIT REAL, INTENT(IN), DIMENSION(N_Snow,NCH) :: WESN, HTSNN, SNDZ ! -------------------------------------- ! OUTPUTS: ! ASNVDR: snow albedo, visible direct ! ASNNDR: snow albedo, near-infrared direct ! ASNVDF: snow albedo, visible diffuse ! ASNNDF: snow albedo, near-infrared diffuse REAL, INTENT(OUT), DIMENSION(NCH) :: ASNVDR, ASNNDR, ASNVDF, ASNNDF ! -------------------------------------- ! Other variables as needed. Includes: ! SSA: snow specific surface area ! RHO_FS: fresh snow density ! EFFG: effective ice thickness (m) INTEGER :: I,M,J,K,K2 INTEGER, PARAMETER :: NTYPS_SIB=9 REAL :: rho_fs,DEGSZA,SD,SZASIN,COS50,SSALBV,SSALBN,AV,AN, & WSS, TS, FAC, FVEG, TOTDEP, SWE, DENS_EXC, AREASC, & DENSITY, ASNVDR_VEG, ASNNDR_VEG, ASNVDF_VEG, & ASNNDF_VEG, SUM1, SUM2 REAL, DIMENSION(NCH) :: SZTH REAL, DIMENSION(N_snow,NCH) :: ABSCOV, ABSCON, EFFG, SSA, DENEL REAL, DIMENSION(N_snow,N_Constit,NCH) :: CONCENT REAL, DIMENSION(3) :: CWESN, CHTSNN, CSNDZ REAL, DIMENSION(3,NCH) :: TPSN, FICES REAL, DIMENSION(3) :: CTPSN, CFICES ! ********************************************************************* !FPP$ EXPAND (COEFFSIB) DO I=1,NCH SZTH(I)=ZTH(I) DEGSZA=ACOS(SZTH(I))*180./PIE SZASIN=SQRT(1.-(SZTH(I)**2.0)) COS50=COS(2.*PIE*50./360.) ! When it is cloud-covered, SZTA is set to 50 degree. THE VALUE ! USED HERE (0.1) CAN BE TUNED! IF(pardir(i)/(pardir(i)+pardif(i)+1.e-20) < 0.1) SZTH(I)=COS50 ! CTPSN: LAYER TEMPERATURE [degree C] DO M=1,N_Snow CWESN(M)=WESN(M,I) CHTSNN(M)=HTSNN(M,I) END DO CALL StieglitzSnow_calc_tpsnow(N_snow,CHTSNN,CWESN,CTPSN,CFICES) DO M=1,N_Snow TPSN(M,I)=CTPSN(M) FICES(M,I)=CFICES(M) END DO !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! SNOW ALBEDOES SWE=SUM(WESN(:,I)) TOTDEP=SUM(SNDZ(:,I)) AREASC = MIN(SWE/WEMIN,1.) DENSITY=(SWE/(AREASC+1.e-20)) / (TOTDEP+1.e-20) WSS=UM(I) TS=RTS(I) CALL WFSDEN(WSS,TS,RHO_FS) DENS_EXC=MAX(0., DENSITY-RHO_FS) !********************************************************************* IF(SWE > 0.01) THEN !========== SNOW CASE ========== DO M=1,N_SNOW ! Dry snow density in each snow layer [kg m-3] (N_Snow layers) DENEL(M,I)=(FICES(M,I)*WESN(M,I)/(AREASC+1.e-20))/(SNDZ(M,I) + 1.e-20) DENEL(M,I)=MIN(RHOMA,DENEL(M,I)) ! Mass concentrations of dust,BC, and OC [kg kg-1] (3 layers) IF(PRESENT(RCONSTIT)) THEN DO K=1,N_Constit CONCENT(M,K,I)=RCONSTIT(M,K,I)/(WESN(M,I) + 1.e-20) ENDDO ELSE DO K=1,N_Constit CONCENT(M,K,I)=0. ENDDO ENDIF ! Snow specific surface area (SSA) [m2 kg-1] ! Equation based on Yamazaki et al. [1991,1993] ! The equation is also based on the observed data ! by Narita [1971], Teionkagaku SSA(M,I)=10.0**(-15.32*((DENEL(M,I)/1000.)**3.0) & +16.65*((DENEL(M,I)/1000.)**2.0) & -7.3*(DENEL(M,I)/1000.)+2.23) ! Another Equation for SSA based on the equation (1) ! in Domine et al., JGR, vol.112, 2007 !! SSA(M,I)=( -308.2D0*LOG(DENEL(M,I)/1000.D0)-206.D0 )/10.D0 ! When each snow layer has excess water (snow gets wet) ! SSA is decreased to 60% values of the original value. IF(TPSN(M,I).GE.-0.001) SSA(M,I)=0.6*SSA(M,I) ! effective ice thickness ! comparable to effective snow grain radius ! if EFFG is multiplied by 3 EFFG(M,I)=2.D0/(DENICE*SSA(M,I) + 1.e-20) ! Calculation of Absorption coefficients [m-1] ! Equations for mass con. [kg kg-1] SUM2=0. do K=1,N_Constit SUM2=SUM2+CONCENT(M,K,I) ENDDO ! VIS range SUM1=0. do K=1,N_Constit SUM1=SUM1+CONCENT(M,K,I)*1.E03*ABVIS(K) ENDDO ABSCOV(M,I)=DENICE*SUM1+(1.-SUM2)*AICEV ! NIR range SUM1=0. do K=1,N_Constit SUM1=SUM1+CONCENT(M,K,I)*1.E03*ABNIR(K) ENDDO ABSCON(M,I)=DENICE*SUM1+(1.-SUM2)*AICEN ENDDO ! VIS & NIR CASES BY NEW SNOW ALBEDO MODEL ASNVDR(I)=ALB(NCH,I,N_Snow,RIV,ABSCOV,EFFG,DENEL,SNDZ,DENICE) ASNNDR(I)=ALB(NCH,I,N_Snow,RIN,ABSCON,EFFG,DENEL,SNDZ,DENICE) ! Diffuse components are calculated with the assumed SZTH of 50 degrees ASNVDF(I)=ASNVDR(I) ASNNDF(I)=ASNNDR(I) !+++++ Taking the effect of solar zenith angle (SZA) into account ++++++ ! The second terms for VIS and NIR in equations (6) & (7) ! [Marks and Dozier, Water Resources Research, Vol. 28, 3043-3054] ! For VIS [Note: unit for EFFG is um] ASNVDR(I)=ASNVDR(I) & -( (SQRT(1.5*(EFFG(1,I)*1.0E06))*1.375E-03)*(1.-COS50) ) & +( (SQRT(1.5*(EFFG(1,I)*1.0E06))*1.375E-03)*(1.-SZTH(I)) ) ! For NIR [Note: unit for EFFG is um] ASNNDR(I)=ASNNDR(I) & -( (( SQRT(1.5*(EFFG(1,I)*1.0E06))*2.0E-03)+0.1)*(1.-COS50) ) & +( (( SQRT(1.5*(EFFG(1,I)*1.0E06))*2.0E-03)+0.1)*(1.-SZTH(I)) ) ELSE !========== NO SNOW CASE ========== ! Use these when you use the original snow albedo model ASNVDR(I) = MAX(SNWALB_VISMIN, SNWALB_VISMAX - DENS_EXC*.0006) ASNNDR(I) = MAX(SNWALB_NIRMIN, SNWALB_NIRMAX - DENS_EXC*.0006) ASNVDF(I) = MAX(SNWALB_VISMIN, SNWALB_VISMAX - DENS_EXC*.0006) ASNNDF(I) = MAX(SNWALB_NIRMIN, SNWALB_NIRMAX - DENS_EXC*.0006) ENDIF ENDDO RETURN END SUBROUTINE ALB_WITH_IMPURITY ! ********************************************************************** !**** ------------------------------------------------------------------ !**** //////////////// Added by Teppei J. Yasunari ///////////////START1 !**** ------------------------------------------------------------------ !======================================================================= != To determine RHOFS = !======================================================================= ! Version 5.0.0 by Teppei J. Yasuanari on 09/23/2010 SUBROUTINE WFSDEN(UM,RTS,RHO_FS) REAL, INTENT(IN) :: UM,RTS ! REAL, INTENT(IN) :: UM,RTS,TSNOW,RPRES,RQST ! REAL, INTENT(IN),DIMENSION(3) :: TPSN REAL, INTENT(OUT) :: RHO_FS REAL :: ARHOFS,V1,T2,SR,ESAT,TW,E INTEGER :: IFG !--------------------------------------------------------------! ! YOU CAN CHOOSE ONE OF 10 TYPES OF "RHOFS" ! !--------------------------------------------------------------! ! ESAT : Saturated water vapor pressure [hPa] ! E : Water vapor pressure [hPa] ! TW : Wet bulb temperature [degree C] ! ! V1=UM ! T2=RTS-273.15 ! SR=TSNOW*3600. ! ! Saturated water vapor pressure !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Fresh snow density by Melloh et al. (2002) ! Hydrological Processes, 16, 3571-3584 ! ! RHOFS TYPE 1 ! ! IF(T2.GE.0.) THEN ! TW=T2-((1013.25/(0.667*RPRES))*(ESAT-E)) ! ELSE IF (T2.LT.0.) THEN ! TW=T2-((1013.25/(0.589*RPRES))*(ESAT-E)) ! END IF ! ! TWW=TW ! ARHOFS=(0.05+(0.0017*(((TW+273.15)-258.16)**1.5)))*1000. ! ARHOFS=MAX(ARHOFS,90.) ! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Multiple regression equation for new snow density ! based on Kajikawa et al. (2004), Sepyo, 66, p.561-565. ! ! RHOFS TYPE 2 ! Spatial dendrite (R**2=0.948; 10% confidence limit) ! ARHOFS=13.3+(53.9*SR)+(6.54*V1) ! ! RHOFS TYPE 3 ! Stellar crystal (R**2=0.817; 5% confidence limit) ! ARHOFS=23.4+(37.5*SR)+(7.32*V1)+(0.579*T2) ! ! RHOFS TYPE 4 ! Rimed stellar crystal (R**2=0.442; 5% confidence limit) ! ARHOFS=41.2+(8.26*SR)+(5.16*V1)+(0.422*T2) ! ! RHOFS TYPE 5 ! Rimed spatial dendrites (R**2=0.369; 5% confidence limit) ! ARHOFS=67.5+(23.4*SR)-(1.29*V1)+(3.65*T2) ! ! RHOFS TYPE 6 ! Rimed spatial dendrites 2 (R**2=0.321; 5% confidence limit) ! ARHOFS=43.1*EXP(0.106*V1) ! !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! 2001 J. Hydrometeorology (Boone and Etchevers) ! Fresh snow density based on wind speed and air temeperature ! ! RHOFS TYPE 7 ! ARHOFS=109.+(6.*(RTS-273.16))+(26.*(V1**0.5)) ! !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Equation 18 by Yamazaki et al. [1998] ! Seppyo, 60, 131-141. ! RHOFS TYPE 8 ARHOFS=67.+(13.*UM) IFG=2 IF(RTS-TF >= 0.) ARHOFS=200. !+++ USE Constant value of fresh snow density [kg m-3]+++ ! RHOFS TYPE 9 ! ARHOFS=110. ! Used in Yasunari et al. (2010) ! RHOFS TYPE 10 ! ARHOFS=150. ! Used in original snowpack model ! RHOFS does not exceed the maximum of 200. [kg m-3]. RHO_FS=MIN(ARHOFS,200.) ! write(*,*) 'FRESH SNOW DENSITY',RHOFS ! (comment out for den) rho_fs=150. RETURN END SUBROUTINE WFSDEN !======================================================================= != New snow albedo functions = !======================================================================= ! Version 5.0.0 by Teppei J. Yasuanari on 09/23/2010 FUNCTION ALB(NCH,I,L,RI,ABSCO,EFFG,DENEL,SNDZ,DENICE) IMPLICIT NONE INTEGER :: J,K INTEGER, INTENT(IN) :: L,I,NCH REAL :: ALB,DZ REAL, INTENT(IN) :: RI,DENICE REAL, INTENT(IN), DIMENSION(L,NCH) :: EFFG,ABSCO,SNDZ,DENEL REAL, DIMENSION(L) :: REF,TRANS,A,B,AMU, & TJALPHA,TJBETA,DI DO J=1,L ! Reflectance in one ice layer REF(J)=RI+( ((1.D0-RI)**2.D0)*RI & *EXP(-2.D0*ABSCO(J,I)*EFFG(J,I)) & /(1.D0-( (RI*EXP(-ABSCO(J,I)*EFFG(J,I)))**2.D0 )) ) ! Transparency in one ice layer TRANS(J)=( ((1.D0-RI)**2.D0)*EXP(-ABSCO(J,I)*EFFG(J,I)) & /(1.D0-( (RI*EXP(-ABSCO(J,I)*EFFG(J,I)))**2.D0 )) ) A(J)=( (1.D0-TRANS(J))/EFFG(J,I) )*( DENEL(J,I)/DENICE ) B(J)=( REF(J)/EFFG(J,I) )*( DENEL(J,I)/DENICE ) ! AMU(J)=SQRT( (A(J)**2.D0) - (B(J)**2.D0) ) AMU(J)=SQRT( MAX(0.D0,(A(J)**2.D0) - (B(J)**2.D0)) ) TJALPHA(J)=( A(J)-AMU(J) )/B(J) TJBETA(J)=( A(J)+AMU(J) )/B(J) END DO DI(L)=0. DO K=L-1,1,-1 IF(K==2) THEN DZ=SNDZ(1,I)+SNDZ(2,I) ELSE IF(K==1) THEN DZ=SNDZ(1,I) END IF DI(K)=( (TJALPHA(K)-TJBETA(K+1))*DI(K+1) & *EXP(2.D0*DZ*AMU(K+1))+TJALPHA(K)-TJALPHA(K+1) )& /( (TJBETA(K+1)-TJBETA(K))*DI(K+1) & *EXP(2.D0*DZ*AMU(K+1))+TJALPHA(K+1)-TJBETA(K) ) & *( EXP(-2.D0*DZ*AMU(K)) ) END DO ! Surface Snow Albedo ALB=RI+( ((1.D0-RI)**2.D0)*(TJBETA(1)*DI(1)+TJALPHA(1)) & /(1.D0+DI(1)-(RI*(TJBETA(1)*DI(1)+TJALPHA(1)))) ) RETURN END FUNCTION ALB ! ******************************************************************** subroutine StieglitzSnow_echo_constants(logunit) ! reichle, 12 Aug 2014 implicit none integer, intent(in) :: logunit write (logunit,*) write (logunit,*) '-----------------------------------------------------------' write (logunit,*) write (logunit,*) 'StieglitzSnow_echo_constants():' write (logunit,*) write (logunit,*) 'PIE = ', PIE write (logunit,*) 'ALHE = ', ALHE write (logunit,*) 'ALHM = ', ALHM write (logunit,*) 'TF = ', TF write (logunit,*) 'RHOW = ', RHOW write (logunit,*) write (logunit,*) 'MINSWE = ', MINSWE write (logunit,*) 'WEMIN = ', WEMIN write (logunit,*) 'CPW = ', CPW write (logunit,*) 'RHOMA = ', RHOMA write (logunit,*) 'DZ1MAX = ', DZ1MAX write (logunit,*) write (logunit,*) 'SNWALB_VISMIN = ', SNWALB_VISMIN write (logunit,*) 'SNWALB_NIRMIN = ', SNWALB_NIRMIN write (logunit,*) write (logunit,*) 'end StieglitzSnow_echo_constants()' write (logunit,*) write (logunit,*) '-----------------------------------------------------------' write (logunit,*) end subroutine StieglitzSnow_echo_constants end module StieglitzSnow