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