      subroutine albsfc(sm,sice,sno,vegfrc,hphi,hlam,albedo,albase  &
                       ,mxsnal,idat)

!********************************************************************!
!                                                                    !
!     calculat surface albedo                                        !
!                                                                    !
!********************************************************************!
      implicit none
      include "param_o.h"
      include "const.h"

      real,parameter::dtr=3.141592654/180.
      real,parameter::H1=1.0,D5=5.E-1,D01=1.00E-2,HM1=-1.0,    &
         H90=90.0,H360=360.0,D00=0.0,radfac=57.295777951
      real,parameter::snup=0.040,salp=2.6
      REAL*4:: ALBC1(361,180),ALBC2 (361,180),ALBC3(361,180),  &
         ALBC4(361,180),ALBC(361,180),SNALMX(361,180),s1,s2,   &
         wgt1,wgt2
      real::elat,elon,dif,w1,w2,ar1,ar2,ar3,ar4,arx,rsnow,     &
            snofac
      real,dimension(0:im+1,0:jm+1,nm)::albase,mxsnal,sm,sice, &
                hphi,hlam,vegfrc,albedo,sno
      integer*4::juld,julm(13)
      integer::i,j,n,ilon1,ilon2,ilat1,ilat2,ilonx,ilatx
      DATA JULM/0,31,59,90,120,151,181,212,243,273,304,334,365/

!----------------------------------------
!c
!C READ MAX SNOW ALBEDO FILE
!C 90N-20N VIA DAVE ROBINSON, JCAM, 1985, VOL. 24, 402-411
!C 20N-90S VIA LAND-SFC TYPE 'CORRELATED' TO ROBINSON DATABASE
!C SNALMX = GLOBAL 1-DEG x 1-DEG MAXIMUM SNOW ALBEDO
!C Units in percent, later converted to fraction
!C values over sea=0.0 (non-land mass)
!C values over land between 21.0 and 80.0
        write(6,*) 'snalmx values read in:'
      READ(20)SNALMX

!-----------------------------------------------------------------------
! READ ALBEDO FILES
! GLOBAL 1-DEG x 1-DEG (4) SEASONAL SNOWFREE ALBEDO
! 90N-90S VIA MATTHEWS ***get reference
! ALBC1 = 3-MONTH AVERAGE CENTERED ON 31 Jan
! ALBC2 = 3-MONTH AVERAGE CENTERED ON 30 Apr
! ALBC3 = 3-MONTH AVERAGE CENTERED ON 31 Jul
! ALBC4 = 3-MONTH AVERAGE CENTERED ON 31 Oct
! Units in percent, later converted to fraction
! values over sea=6.0 (non-land mass)
! values over land between 11.0 and 32.0 (exception, glacier=0.75)
! *NOTE: in future, replace w/high-res 0.144-deg albedo NESDIS database
!   ...and follow similar spatial/temporal averaging for greenness frac
!

      READ(21)ALBC1
      READ(22)ALBC2
      READ(23)ALBC3
      READ(24)ALBC4

! FIND JULIAN DAY OF YEAR TO DO THE TEMPORAL INTERPOLATION
      JULD = JULM(IDAT(1)) + IDAT(2)
      IF(JULD.LE.32) THEN
        S1 = 32 - JULD
        S2 = JULD + 30
        WGT1 = S2/(S1+S2)
        WGT2 = S1/(S1+S2)
        DO J = 1,180
        DO I = 1,361
        ALBC(I,J) = WGT1 * ALBC1(I,J) + WGT2 * ALBC4(I,J)
        END DO
        END DO
      ELSE IF(JULD.LE.121.AND.JULD.GT.32) THEN
        S1 = 121 - JULD
        S2 = JULD - 32
        WGT1 = S2/(S1+S2)
        WGT2 = S1/(S1+S2)
        DO J = 1,180
        DO I = 1,361
        ALBC(I,J) = WGT1 * ALBC2(I,J) + WGT2 * ALBC1(I,J)
        END DO
        END DO
      ELSE IF(JULD.LE.213.AND.JULD.GT.121) THEN
        S1 = 213 - JULD
        S2 = JULD - 121
        WGT1 = S2/(S1+S2)
        WGT2 = S1/(S1+S2)
        DO J = 1,180
        DO I = 1,361
        ALBC(I,J) = WGT1 * ALBC3(I,J) + WGT2 * ALBC2(I,J)
        END DO
        END DO
      ELSE IF(JULD.LE.305.AND.JULD.GT.213) then
        S1 = 305 - JULD
        S2 = JULD - 213
        WGT1 = S2/(S1+S2)
        WGT2 = S1/(S1+S2)
        DO J = 1,180
        DO I = 1,361
        ALBC(I,J) = WGT1 * ALBC4(I,J) + WGT2 * ALBC3(I,J)
        END DO
        END DO
      ELSE
        S1 = 365 - JULD + 32
        S2 = JULD - 305
        WGT1 = S2/(S1+S2)
        WGT2 = S1/(S1+S2)
        DO J = 1,180
        DO I = 1,361
        ALBC(I,J) = WGT1 * ALBC1(I,J) + WGT2 * ALBC4(I,J)
        END DO
        END DO
      END IF
!
! BEGIN SPATIAL INTERPOLATION FOR BASELINE SNOWFREE ALBEDO AND MAX SNOW
! ALBEDOS
!
      do n=1,nm
      DO J=1,JM
      DO I=1,IM
!new
          ALBASE(I,J,n)=20.
          MXSNAL(I,J,n)=55.
!new
          IF (SM(I,J,n) .LT. 0.9.or.sice(i,j,n).eq.1.) THEN
!
!-----------------------------------------------------------------------
! IF LANDMASS, DETERMINE LAT,LON AND WEIGHTS
!
      ELAT=90.+hphi(I,J,n)/DTR
      ELON=hlam(I,J,n)/DTR
      IF(ELON.GT.360.)ELON=ELON-360.
      ILON1=INT(ELON)
      DIF=ELON-ILON1
      IF(DIF.GT.D5)ILON1=ILON1+1
      IF(ILON1.EQ.D00)ILON1=360
      ILON2=ILON1+1
      ILAT1=INT(ELAT)
      DIF=ELAT-ILAT1
      IF(DIF.GT.D5)ILAT1=MIN(ILAT1+1,179)
      ILAT2=ILAT1+1
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
      ILAT1=MAX(1,MIN(180,ILAT1))
      ILAT2=MAX(1,MIN(180,ILAT2))
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME
! FIX ME

      W1=ELON-ILON1+D5
      IF(W1.LT.D00)W1=W1+360.
      W2=ELAT-ILAT1+D5
!
      AR1=W1*W2
      AR2=W1*(H1-W2)
      AR3=(H1-W1)*(H1-W2)
      AR4=(H1-W1)*W2
!-----------------------------------------------------------------------
! INTERPOLATE BASELINE SNOWFREE ALBEDO TO E GRID
! INPUT:  ALBC (GLOBAL 1-DEG x 1-DEG)
! OUTPUT: ALBASE (SPECIFIED ETA GRID)
!
! INTERPOLATE MAX SNOW ALBEDO TO E GRID
! INPUT:  SNALMX (GLOBAL 1-DEG x 1-DEG)
! OUTPUT: MXSNAL (SPECIFIED ETA GRID)
!
! DOING BOTH BASELINE SNOWFREE ALBEDO AND MAX SNOW ALBEDO INTERPOLATIONS
! IN THE SAME BLOCK IS POSSIBLE SINCE THEY HAVE IDENTICAL LAND-SEA MASKS
!            IF ( (ALBC(ILON2,ILAT2) .NE. 0.) .AND.
!     .           (ALBC(ILON2,ILAT1) .NE. 0.) .AND.
!     .           (ALBC(ILON1,ILAT1) .NE. 0.) .AND.
!     .           (ALBC(ILON1,ILAT2) .NE. 0.) ) THEN
! the quarterly mathews albedo data base sea values = 6.0
! beginning lat/long=-90,0 (southpole, prime meredian)
! lowest land value = 11.0
! max non-glacial value about 32.0 +/-
! glacial value = 75.0
!
            IF ( (ALBC(ILON2,ILAT2) .GT. 7.) .AND.  &
                 (ALBC(ILON2,ILAT1) .GT. 7.) .AND.  &
                 (ALBC(ILON1,ILAT1) .GT. 7.) .AND.  &
                 (ALBC(ILON1,ILAT2) .GT. 7.) ) THEN
!-----------------------------------------------------------------------
! ALL 4 SURROUNDING POINTS ARE LAND POINTS
              ALBASE(I,J,n)=AR1*ALBC(ILON2,ILAT2)+  &
                            AR2*ALBC(ILON2,ILAT1)+  &
                            AR3*ALBC(ILON1,ILAT1)+  &
                            AR4*ALBC(ILON1,ILAT2)
              MXSNAL(I,J,n)=AR1*SNALMX(ILON2,ILAT2)+  &
                            AR2*SNALMX(ILON2,ILAT1)+  &
                            AR3*SNALMX(ILON1,ILAT1)+  &
                            AR4*SNALMX(ILON1,ILAT2)
!
            ELSE
!-----------------------------------------------------------------------
! ONE OR MORE OF THE 4 SURROUNDING POINT ARE NOT LAND POINTS
! TAKE NEAREST NEIGHBOR LAND POINT IN THE FOLLOWING ORDER:
! (ILON2,ILAT2),(ILON2,ILAT1),(ILON1,ILAT1),(ILON1,ILAT2)
!
              ARX=-999.
              IF (ALBC(ILON2,ILAT2) .GT. 7.) THEN
                IF (AR1 .GT. ARX) THEN
                  ARX=AR1
                  ILONX=ILON2
                  ILATX=ILAT2
                ENDIF
              ENDIF
              IF (ALBC(ILON2,ILAT1) .GT. 7.) THEN
                IF (AR2 .GT. ARX) THEN
                  ARX=AR2
                  ILONX=ILON2
                  ILATX=ILAT1
                ENDIF
              ENDIF
              IF (ALBC(ILON1,ILAT1) .GT. 7.) THEN
                IF (AR3 .GT. ARX) THEN
                  ARX=AR3
                  ILONX=ILON1
                  ILATX=ILAT1
                ENDIF
              ENDIF
              IF (ALBC(ILON1,ILAT2) .GT. 7.) THEN
                IF (AR4 .GT. ARX) THEN
                  ARX=AR4
                  ILONX=ILON1
                  ILATX=ILAT2
                ENDIF
              ENDIF
!-----------------------------------------------------------------------
! Use nearest land neighbor:
              IF (ARX .GT. 0.) THEN
                ALBASE(I,J,n)=ALBC(ILONX,ILATX)
                MXSNAL(I,J,n)=SNALMX(ILONX,ILATX)
              ELSE
!-----------------------------------------------------------------------
! NO SURROUNDING POINTS ARE LAND (ARX=-999):
! SET DEFAULT SNOWFREE ALBEDO=20 PERCENT, DEFAULT SNOWALB=55 PERCENT
                ALBASE(I,J,n)=20.
                MXSNAL(I,J,n)=55.
!                PRINT *,'AT: LAT,LON',ELAT,ELON
!                PRINT *,'SNOWFREE ALBEDO SET TO DEFAULT VALUE OF 20%'
!                PRINT *,'MAX SNOW ALBEDO SET TO DEFAULT VALUE OF 55%'
!-----------------------------------------------------------------------
! end of ALBASE,MXSNAL land points<4 block
              ENDIF
!-----------------------------------------------------------------------
! end of ALBASE,MXSNAL interpolation block
            ENDIF
!-----------------------------------------------------------------------
! end of land (SM=0) block
          ENDIF
!-----------------------------------------------------------------------
! CONVERT ALBEDO UNITS: PERCENT TO FRACTION
          ALBASE(I,J,n)=ALBASE(I,J,n)*D01
          MXSNAL(I,J,n)=MXSNAL(I,J,n)*D01
      ENDDO
      ENDDO
      ENDDO

      do n=1,nm
      do i=1,im
      do j=1,jm
        if(sm(i,j,n).gt.0.9)then
	  albedo(i,j,n)=0.06
	  albase(i,j,n)=0.06
        endif
	  if(sice(i,j,n).eq.1.)then
	    albedo(i,j,n)=0.6
	    albase(i,j,n)=0.6
          endif

      enddo
      enddo
      enddo

!-----------------------------------------------------------------------
! DETERMINE ALBEDO OVER LAND
      do n=1,nm
      DO J=1,JM
      DO I=1,IM
          IF(SM(I,J,n).LT.0.9.AND.SICE(I,J,n).LT.0.9) THEN
! SNOWFREE ALBEDO
            IF ( (SNO(I,J,n) .EQ. 0.0) .OR.  &
                 (ALBASE(I,J,n) .GE. MXSNAL(I,J,n) ) ) THEN
              ALBEDO(I,J,n) = ALBASE(I,J,n)
           ELSE
! MODIFY ALBEDO IF SNOWCOVER:
! BELOW SNOWDEPTH THRESHOLD...
              IF (SNO(I,J,n) .LT. SNUP) THEN
                RSNOW = SNO(I,J,n)/SNUP
                SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP))
! ABOVE SNOWDEPTH THRESHOLD...
              ELSE
                SNOFAC = 1.0
              ENDIF
! CALCULATE ALBEDO ACCOUNTING FOR SNOWDEPTH AND VEGFRC
              ALBEDO(I,J,n) = ALBASE(I,J,n)  &
                + (1.0-VEGFRC(I,J,n))*SNOFAC*(MXSNAL(I,J,n)-ALBASE(I,J,n))
            ENDIF
          END IF
      ENDDO
      ENDDO
      ENDDO

      END SUBROUTINE ALBSFC 
