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