Kd.zip_LANDSAT8 MATLAB_landsat matlab_landsat8_landsat8__matlab

Program IOP C INVERSION of IOP BASED on Rrs and REMOTELY RETRIEVED Kd C LOISEL AND POTEAU, IOCCG REPORT N°5, 2006 C LOISEL AND STRAMSKI, APPLIED OPTICS, Vol 39, N°18, 2000 C LOISEL et al, APPLIED OPTICS, Vol 40, N°15, 2001 DIMENSION RRS(5) DIMENSION CTERRS(2,3,5), CTEKD(3,3,5), AW(5), BBW(5) DIMENSION A(5), BB(5), XKD(5) CHARACTER*120 FILENAME C LUTs OPEN(10,FILE='./LUT_RRS') OPEN(11,FILE='./LUT_KD') OPEN(12,FILE='./LUT_AW') DO J=1,3 DO K=1,5 READ(10,*) (CTERRS(I,J,K), I=1,2) ENDDO ENDDO DO J=1,3 DO K=1,5 READ(11,*) (CTEKD(I,J,K),I=1,3) ENDDO ENDDO DO K=1,5 READ(12,*) AW(K), BBW(K) ENDDO CLOSE(10) CLOSE(11) CLOSE(12) OPEN(1,FILE='Result_KD') OPEN(2,FILE='Result_A') OPEN(3,FILE='Result_BB') C INPUTs: Rrs et ASOL WRITE(6,*) '--------------------------------------------------' WRITE(6,*) ' INPUTS ' WRITE(6,*) '--------------------------------------------------' WRITE(6,*) WRITE(6,*) 'Rrs = (Lw(0+)/Ed(0+)) at 410, 440, 490, 510, 550 nm' WRITE(6,*) 'Solar Zenith Angle (degree) if SZA>60 then SZA=60' WRITE(6,*) '' WRITE(6,*) '--------------------------------------------------' WRITE(6,*) ' INPUT FILE FORMAT ' WRITE(6,*) '--------------------------------------------------' WRITE(6,*) '' WRITE(6,*) 'Rrs410, Rrs440, Rrs490, Rrs510, Rrs550, SZA degree' WRITE(6,*) '' WRITE(6,*) '--------------------------------------------------' WRITE(6,*) ' OUTPUT FILE FORMAT ' WRITE(6,*) '--------------------------------------------------' WRITE(6,*) '' WRITE(6,*) 'Result_KD: Kd at 410, 440, 490, 510, 550 nm' WRITE(6,*) 'Result_A : a_total at 410, 440, 490, 510, 550 nm' WRITE(6,*) 'Result_BB: bb_total at 410, 440, 490, 510, 550 nm' WRITE(6,*) '--------------------------------------------------' WRITE(6,*) '' WRITE(6,*) 'PLEASE INPUT FILE NAME' READ(5,*) FILENAME OPEN(48,FILE=FILENAME,STATUS='OLD') DO I=1,10000000 READ(48,*,END=111) (RRS(K), K=1,5), ASOL CALL IOPIOP (CTERRS,CTEKD,RRS,ASOL,XKD,AW,BBW,A,BB) WRITE(1,991) (XKD(K), K=1,5) WRITE(2,991) (A(K), K=1,5) WRITE(3,991) (BB(K), K=1,5) 991 FORMAT(f7.5,2x,f7.5,2x,f7.5,2x,f7.5,2x,f7.5) 111 ENDDO CLOSE(1) CLOSE(2) CLOSE(3) END SUBROUTINE IOPIOP (CTERRS,CTEKD,RRS,ASOL,XKD,AW,BBW,A,BB) DIMENSION RRS(5) DIMENSION CTERRS(2,3,5), CTEKD(3,3,5), AW(5),BBW(5) DIMENSION R(5), XKD(5), R0(3), XKD0(3) DIMENSION A(5),BB(5) REAL EC50,HILLSLOPE,h,eta,MINH,MAXH PI=ACOS(-1.0) XMU=COS(ASIN(SIN(ASOL*PI/180.)/1.34)) RAP=RRS(2)/RRS(5) DO K=1,5 ! LAMBDA LOOP DO J=1,3 ! MUS LOOP R0(J) = CTERRS(1,J,K)*RRS(K)**CTERRS(2,J,K) XKD0(J) =10**((CTEKD(1,J,K)*log10(RAP)+CTEKD(2,J,K))/ + (CTEKD(3,J,K)+log10(RAP))) ENDDO ! END OF MUS LOOP IF(ASOL.LT.60) THEN THETAS=ASOL ELSE THETAS=60 ENDIF CALL INTERPOL(K,THETAS,R0,XKD0,R,XKD) BW=2*BBW(K) ETA=0.05 BALPHA = -0.0042 *ASOL+0.3594 AALPHA = -0.0007*ASOL+0.1024 ALPHA=AALPHA*log10(ETA)+BALPHA DELTA=0.871+0.40*ETA-1.83*ETA**2. BB(K) = XKD(K)*10**(ALPHA)*R(K)**DELTA DO LL=1,2 ETA = BW/(BW+((BB(K)-BBW(K))/0.0183)) IF(BB(k).LT.BBW(K)) THEN ETA = BW/(BW+(((BBW(K)+BBW(K)/10)-BBW(K))/0.0183)) ENDIF ALPHA=(AALPHA)*log10(ETA)+(BALPHA) DELTA=0.871+0.40*ETA-1.83*ETA**2. BB(K) = XKD(K)*10**(ALPHA)*R(K)**DELTA ENDDO ETA = BW/(BW+((BB(K)-BBW(K))/0.0183)) IF(BB(k).LT.BBW(K)) THEN ETA = BW/(BW+(((BBW(K)+BBW(K)/10)-BBW(K))/0.0183)) ENDIF MINH = -0.0034*ASOL+1.248 MAXH = -0.0084*ASOL+1.7223 EC50 = 0.0033*ASOL-1.9837 HILLSLOPE = 0.0034*ASOL**2-0.0674*ASOL-9.34 H = MINH+(MAXH-MINH)/(1+10**(((EC50)-log10(ETA))*HILLSLOPE)) H=10**H A(K) = XMU*XKD(K)/(1+H*(R(K)/(1-R(K))))**0.5 ENDDO ! END LAMBDA LOOP END SUBROUTINE INTERPOL(K,THETAS,R0,XKD0,R,XKD) DIMENSION R(5), XKD(5), R0(3), XKD0(3) IF(THETAS.GE.0.AND.THETAS.LE.30) THEN R(K)=THETAS*(R0(1)-R0(2))/(0-30)+(-30*R0(1))/(0-30) XKD(K)=THETAS*(XKD0(1)-XKD0(2))/(0-30)+(-30*XKD0(1))/(0-30) ELSEIF(THETAS.GT.30.AND.THETAS.LE.60) THEN R(K)=THETAS*(R0(2)-R0(3))/(30-60)+(30*R0(3)-60*R0(2))/(30-60) XKD(K)=THETAS*(XKD0(2)-XKD0(3))/(30-60)+ + (30*XKD0(3)-60*XKD0(2))/(30-60) ENDIF END