#if defined(L08_1A) ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! Routine to calculate the spectral albedo of the land surface using ! the two stream approach of Sellers, 1995. !********************************************************************** SUBROUTINE albpft (p_field,land_field, & land_index,tile_index,tile_pts,ilayers, & getprofile, & albsoil,cosz,lai,alb_type, & fapar_dir,fapar_dif,fapar_dir2dif, & fapar_dif2dif,fapar_dir2dir,fsun) USE nstypes USE pftparm USE switches, ONLY : & ! imported scalars with intent(in) can_rad_mod USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook IMPLICIT NONE ! Subroutine arguments ! Scalar arguments with intent(in): INTEGER, INTENT(IN) :: & p_field &! Total number of grid points. ,land_field &! Number of land points. ,ilayers ! Number of layers over which ! the PAR absorption profile is to ! be calculated. LOGICAL, INTENT(IN) :: getprofile ! TRUE means profiles through the ! canopy are calculated ! ! FALSE means no profiles calculated ! Array arguments with intent(in): INTEGER, INTENT(IN) :: & land_index(land_field) & ! Index of land points. ,tile_pts(ntype) & ! Number of land points which ! ! include the surface type. ,tile_index(land_field,ntype) ! Indices of land points which ! ! include the surface type. REAL, INTENT(IN) :: & albsoil(land_field) & ! Soil albedo. ,cosz(p_field) & ! Cosine of the zenith angle. ,lai(land_field,npft) ! Leaf area index. ! Array arguments with intent(out): REAL, INTENT(OUT) :: & alb_type(land_field,ntype,4) & ! ! Albedos for surface types. ! ! (*,*,1) - Direct beam visible ! ! (*,*,2) - Diffuse visible ! ! (*,*,3) - Direct beam near-IR ! ! (*,*,4) - Diffuse near-IR ,fapar_dir(land_field,npft,ilayers) & ! ! Profile of absorbed PAR ! ! - Direct (fraction/LAI) ,fapar_dif(land_field,npft,ilayers) & ! ! Profile of absorbed PAR ! ! - Diffuse (fraction/LAI) ,fapar_dir2dif(land_field,npft,ilayers) & ! ! Profile of scattered PAR & ! ! -direct to diffuse ! ! (fraction/LAI) ,fapar_dif2dif(land_field,npft,ilayers) & ! ! Profile of absorbed PAR ! ! - Diffuse (fraction/LAI) ,fapar_dir2dir(land_field,npft,ilayers) & ! ! Profile of absorbed only direct PAR ! ! - (fraction/LAI) ,fsun(land_field,npft,ilayers) ! fraction of sunlit leaves ! Local arrays: REAL & albudif(land_field,2) & ! Diffuse albedo of the underlying ! ! surface. ,albudir(land_field,2) & ! Direct albedo of the underlying ! ! surface. ,alpl(2) & ! Leaf reflection coefficient. ,om(2) & ! Leaf scattering coefficient. ,rnet_dir(0:ilayers) & ! ! Net downward PAR at layer ! ! boundaries - incident direct ,rnet_dif(0:ilayers) & ! ! Net downward PAR at layer ! ! boundaries - incident diffuse ,taul(2) ! Leaf transmission coefficient. ! Local scalars: REAL & betadir & ! Upscatter parameter for direct beam. ,betadif & ! Upscatter parameter for diffuse beam ,coszm & ! Mean value of COSZ. ,k & ! Optical depth per unit leaf area. ,g & ! Relative projected leaf area in ! direction cosz. ,salb & ! Single scattering albedo. ,sqcost & ! Cosine squared of the mean leaf ! angle to the horizontal. ,b,c,ca,d,f,h,u1 & ! Miscellaneous variables from ,p1,p2,p3,p4,d1 & ! Sellers (1985). ,h1,h2,h3,h7,h8 & ! ,s1,s2,sig ! !----------------------------------------------------------------------- ! Additional work variables to calculate profiles of absorbed PAR !----------------------------------------------------------------------- REAL & dlai & ! LAI increment. ,la & ! Cumulative LAI through canopy. ,drdird_dlai,drdiru_dlai & ,drdifd_dlai,drdifu_dlai & ,u2,u3,d2,h4,h5,h6,h9,h10 & ! Rate of change of fluxes with LAI ! (W/m2/LAI). ,rdiru,rdird,rdifu,rdifd &! work variable ,dabeer_dla(0:ilayers) ! Attenuation of direct solar ! beam per unit LAI INTEGER & band,i,j,l,n ! Loop counters. INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 REAL(KIND=jprb) :: zhook_handle IF (lhook) CALL dr_hook('ALBPFT',zhook_in,zhook_handle) DO l=1,land_field albudif(l,1) = albsoil(l) albudif(l,2) = albsoil(l) albudir(l,1) = albsoil(l) albudir(l,2) = albsoil(l) END DO ! Initialise to zero at top. fapar_dir2dir(:,:,:) = 0.0 fapar_dir2dif(:,:,:) = 0.0 fapar_dif2dif(:,:,:) = 0.0 fsun(:,:,:) = 0.0 DO n=1,npft om(1) = omega(n) om(2) = omnir(n) alpl(1) = alpar(n) alpl(2) = alnir(n) DO band=1,2 ! Visible and near-IR bands taul(band) = om(band) - alpl(band) DO j=1,tile_pts(n) l = tile_index(j,n) i = land_index(l) IF (orient(n) == 0) THEN sqcost = 1./3. g = 0.5 coszm = 1.0 k = g / 0.01 salb = 0.5*om(band) IF (cosz(i) > 0.01) THEN k = g / cosz(i) salb = 0.5*om(band) * & ( 1. - cosz(i)*LOG((cosz(i)+1.)/cosz(i)) ) END IF ELSE IF (orient(n) == 1) THEN sqcost = 1. coszm = 1. ! K = G / COSZ(I) = 1.0 because G = COSZ(I) for horizontal leaves k = 1.0 salb = om(band)/4.0 END IF betadir = (1. + coszm*k)/(om(band)*coszm*k)*salb c = 0.5*( alpl(band) + taul(band) + & (alpl(band) - taul(band))*sqcost ) betadif = c / om(band) b = 1. - (1. - betadif)*om(band) d = om(band)*coszm*k*betadir f = om(band)*coszm*k*(1. - betadir) h = SQRT(b*b - c*c) / coszm sig = (coszm*k)**2 + c*c - b*b IF ( ABS(sig) < 1.0e-4 ) sig = SIGN( 1.0e-4,sig ) u1 = b - c/albudif(l,band) ca = c*albudir(l,band)/albudif(l,band) s1 = EXP(-h*lai(l,n)) s2 = EXP(-k*lai(l,n)) p1 = b + coszm*h p2 = b - coszm*h p3 = b + coszm*k p4 = b - coszm*k d1 = p1*(u1 - coszm*h)/s1 - p2*(u1 + coszm*h)*s1 h1 = -d*p4 - c*f h2 = ( (d - p3*h1/sig) * (u1 - coszm*h) / s1 - & (d - ca - (u1 + coszm*k)*h1/sig)*p2*s2 ) / d1 h3 = - ( (d - p3*h1/sig) * (u1 + coszm*h)*s1 - & (d - ca - (u1 + coszm*k)*h1/sig)*p1*s2 ) / d1 h7 = (c/d1)*(u1 - coszm*h) / s1 h8 = - (c/d1)*(u1 + coszm*h) * s1 alb_type(l,n,2*band-1) = h1/sig + h2 + h3 ! Direct beam alb_type(l,n,2*band) = h7 + h8 ! Diffuse !----------------------------------------------------------------------- ! If required calculate the profile of absorbed PAR through the canopy ! of direct and diffuse beams (BEWARE: assumes PAR is band 1) !----------------------------------------------------------------------- IF ( getprofile .AND. band==1 ) THEN u2 = b-c*albudif(l,band) u3 = f+c*albudif(l,band) d2 = (u2+coszm*h)/s1-(u2-coszm*h)*s1 h4 = -f*p3-c*d h5 = -1./d2*(h4/sig*(u2+coszm*h)/s1 & +(u3-h4/sig*(u2-coszm*k))*s2) h6 = 1./d2*(h4/sig*(u2-coszm*h)*s1 & +(u3-h4/sig*(u2-coszm*k))*s2) h9 = 1./d2*(u2+coszm*h)/s1 h10 = -1./d2*(u2-coszm*h)*s1 !----------------------------------------------------------------------- ! Two-stream equations for direct and diffuse upward and downward beams: ! RDIRU(I)=H1/SIG*EXP(-K*LA)+H2*EXP(-H*LA)+H3*EXP(H*LA) ! RDIRD(I)=(H4/SIG+1)*EXP(-K*LA)+H5*EXP(-H*LA)+H6*EXP(H*LA) ! RDIFU(I)=H7*EXP(-H*LA)+H8*EXP(H*LA) ! RDIFD(I)=H9*EXP(-H*LA)+H10*EXP(H*LA) !----------------------------------------------------------------------- dlai=lai(l,n)/REAL(ilayers) IF ( can_rad_mod /= 5 ) THEN !----------------------------------------------------------------------- ! Differentiate these equations to calculate PAR absorption per unit ! LAI down through the canopy. Centre derivatives in the centre of each ! LAI layer. !----------------------------------------------------------------------- la=0.5*dlai DO i=1,ilayers drdiru_dlai=-k*h1/sig*EXP(-k*la)-h*h2*EXP(-h*la) & +h*h3*EXP(h*la) drdird_dlai=-k*(h4/sig+1)*EXP(-k*la)-h*h5*EXP(-h*la) & +h*h6*EXP(h*la) drdifu_dlai=-h*h7*EXP(-h*la)+h*h8*EXP(h*la) drdifd_dlai=-h*h9*EXP(-h*la)+h*h10*EXP(h*la) fapar_dir(l,n,i)=-drdird_dlai+drdiru_dlai fapar_dif(l,n,i)=-drdifd_dlai+drdifu_dlai la=la+dlai END DO !layers ELSE !----------------------------------------------------------------------- ! can_rad_mod = 5 ! Don't use derivatives. !----------------------------------------------------------------------- la=dlai rnet_dir(0)=(h4-h1)/sig+(h5-h2)+(h6-h3) rnet_dif(0)=(h9-h7)+(h10-h8) DO i=1,ilayers rdiru=h1/sig*EXP(-k*la)+h2*EXP(-h*la)+h3*EXP(h*la) rdird=h4/sig*EXP(-k*la)+h5*EXP(-h*la)+h6*EXP(h*la) rdifu=h7*EXP(-h*la)+h8*EXP(h*la) rdifd=h9*EXP(-h*la)+h10*EXP(h*la) rnet_dir(i)=rdird-rdiru rnet_dif(i)=rdifd-rdifu fapar_dir(l,n,i)=(rnet_dir(i-1)-rnet_dir(i))/dlai fapar_dif(l,n,i)=(rnet_dif(i-1)-rnet_dif(i))/dlai dabeer_dla(i)=(EXP(-k*(la-dlai))-EXP(-k*la))/dlai fapar_dir2dir(l,n,i)=(1.- om(1))*dabeer_dla(i) fapar_dir2dif(l,n,i)= om(1)*dabeer_dla(i)+fapar_dir(l,n,i) fapar_dif2dif(l,n,i)=fapar_dif(l,n,i) fsun(l,n,i)=EXP(-k*la)*(EXP(k*dlai)-1.)/(k*dlai) la=la+dlai END DO !layers END IF ! can_rad_mod END IF ! abs. par (band=1+getProfile) END DO ! points END DO ! bands END DO ! pfts IF (lhook) CALL dr_hook('ALBPFT',zhook_out,zhook_handle) RETURN END SUBROUTINE albpft #endif