MODULE ProVegMap
 IMPLICIT  NONE
 PRIVATE
 TYPE GridDef
    CHARACTER(LEN=20) :: GRIDDEF  
    INTEGER           :: ncols  
    INTEGER           :: nlins  
    REAL              :: X1  
    REAL              :: Y2  
    REAL              :: resX  
    REAL              :: resY  
    REAL              :: Lati  
    REAL              :: Loni  
    REAL              :: nodatavalue
 END TYPE GridDef  
 TYPE(GridDef), ALLOCATABLE :: GridProv(:)
 TYPE(GridDef), ALLOCATABLE :: GridBamg(:)

 REAL         , ALLOCATABLE :: MapVegProVeg(:,:)
 REAL         , ALLOCATABLE :: MapVegBamLSMK(:,:)
 REAL         , ALLOCATABLE :: MapVegBamIbis(:,:)
 REAL         , ALLOCATABLE :: MapVegBamIbisNew(:,:)

 PUBLIC :: InitProVeg
 PUBLIC :: RunProVeg
 PUBLIC :: FinalizeProVeg
CONTAINS 
 SUBROUTINE InitProVeg()
  IMPLICIT NONE
  CHARACTER(LEN=255) :: FileName
  CHARACTER(LEN=255) :: FileName2
  CHARACTER(LEN=255) :: FileName3
  CHARACTER(LEN=255) :: FileName4
  CHARACTER(LEN=50 ) :: Line=''
  CHARACTER(LEN=255) :: Line1
  CHARACTER(LEN=255) :: Line2
  CHARACTER(LEN=255) :: Line3
  REAL               :: lat
  INTEGER            :: ic 
  INTEGER            :: lrec
  INTEGER            :: iline
  INTEGER            :: ivar
  INTEGER            :: j
  FileName='./inputdata/proveg_map/ProVeg_Brasil_GRR.spr'
  FileName2='./inputdata/bam_map/ibismsk.form'
  FileName3='./inputdata/proveg_map/ProVeg_Brasil_GRR.bin'
  FileName4='./inputdata/bam_map/LandSeaMaskNavy.bin'

  ALLOCATE(GridProv(1))
  ALLOCATE(GridBamg(1))
  OPEN(1,FILE=TRIM(FileName), ACCESS='SEQUENTIAL',FORM='FORMATTED',STATUS='OLD', ACTION='READ')
  READ(1,'(12/1x)')
  READ(1,'(a)')Line1
  READ(1,'(a)')Line2
  READ(1,'(a)')Line3
  ivar=0
  iline=0
  DO ic=1,LEN(TRIM(Line2))+10
      IF(Line2(ic:ic) /= '')THEN
         iline=iline+1
         Line(iline:iline)=Line2(ic:ic)
      ELSE 
         IF(LEN(TRIM(Line)) /= 0 .and. ivar==0)THEN
            ivar=ivar+1
            READ(Line,*)GridProv(1)%GRIDDEF  
         ELSE IF(LEN(TRIM(Line)) /= 0 .and. ivar==1)THEN
             ivar=ivar+1
             READ(Line,*)GridProv(1)%ncols  
         ELSE IF(LEN(TRIM(Line)) /= 0 .and. ivar==2)THEN
            ivar=ivar+1
            READ(Line,*)GridProv(1)%nlins  
         ELSE IF(LEN(TRIM(Line)) /= 0 .and. ivar==3)THEN
            ivar=ivar+1
            READ(Line,*)GridProv(1)%X1  
         ELSE IF(LEN(TRIM(Line)) /= 0 .and. ivar==4)THEN
            ivar=ivar+1
            READ(Line,*)GridProv(1)%Y2  
         ELSE IF(LEN(TRIM(Line)) /= 0 .and. ivar==5)THEN
            ivar=ivar+1
            READ(Line,*)GridProv(1)%resX  
         ELSE IF(LEN(TRIM(Line)) /= 0 .and. ivar==6)THEN
            ivar=ivar+1
            READ(Line,*)GridProv(1)%resY  
         ELSE IF(LEN(TRIM(Line)) /= 0 .and. ivar==7)THEN
            ivar=ivar+1
            READ(Line,*)GridProv(1)%nodatavalue
         END IF
         Line=''
         iline=0
      ENDIF
  END DO

  GridProv(1)%Loni=GridProv(1)%X1 
  lat=GridProv(1)%Y2
  DO j=2,GridProv(1)%nlins
   lat=lat-GridProv(1)%resY
  END DO
  GridProv(1)%Lati=lat

  PRINT*,GridProv(1)%ncols
  PRINT*,GridProv(1)%nlins
  PRINT*,GridProv(1)%X1 
  PRINT*,GridProv(1)%Y2 
  PRINT*,GridProv(1)%Loni
  PRINT*,GridProv(1)%Lati
  PRINT*,GridProv(1)%nodatavalue
  ALLOCATE(MapVegProVeg(GridProv(1)%ncols,GridProv(1)%nlins))
  !DO j=1,GridProv(1)%nlins
  !      READ(1,*)(MapVegProVeg(i,GridProv(1)%nlins-j+1),i=1,GridProv(1)%ncols)
  !END DO
  CLOSE(1,STATUS='KEEP')
  INQUIRE(IOLENGTH=lrec)MapVegProVeg
  OPEN(1,FILE=TRIM(FileName3),ACCESS='DIRECT',FORM='UNFORMATTED',RECL=lrec,STATUS='OLD',ACTION='READ') 
  READ(1,rec=1)MapVegProVeg
  CLOSE(1,STATUS='KEEP')

  GridBamg(1)%ncols=720
  GridBamg(1)%nlins=360
  GridBamg(1)%X1=-179.75
  GridBamg(1)%Y2=89.7500
  GridBamg(1)%Loni=GridBamg(1)%X1
  GridBamg(1)%Lati=-89.7500
  GridBamg(1)%resX=0.5000000000
  GridBamg(1)%resY=0.5000000000
  GridBamg(1)%nodatavalue=-999.0
  PRINT*,GridBamg(1)%ncols
  PRINT*,GridBamg(1)%nlins
  PRINT*,GridBamg(1)%X1 
  PRINT*,GridBamg(1)%Y2 
  PRINT*,GridBamg(1)%Loni
  PRINT*,GridBamg(1)%Lati
  PRINT*,GridBamg(1)%nodatavalue

  ALLOCATE(MapVegBamLSMK(GridBamg(1)%ncols,GridBamg(1)%nlins))
  INQUIRE(IOLENGTH=lrec)MapVegBamLSMK
  OPEN(1,FILE=TRIM(FileName4), ACCESS='DIRECT',FORM='UNFORMATTED',RECL=lrec,STATUS='OLD',ACTION='READ') 
  READ(1,rec=1)MapVegBamLSMK
  MapVegBamLSMK = cshift(MapVegBamLSMK,GridBamg(1)%ncols/2,DIM=1)
  CLOSE(1,STATUS='KEEP')

  ALLOCATE(MapVegBamIbis(GridBamg(1)%ncols,GridBamg(1)%nlins))
  ALLOCATE(MapVegBamIbisNew(GridBamg(1)%ncols,GridBamg(1)%nlins))
  INQUIRE(IOLENGTH=lrec)MapVegBamIbis
  OPEN(1,FILE=TRIM(FileName2), ACCESS='DIRECT',FORM='UNFORMATTED',RECL=lrec,STATUS='OLD',ACTION='READ') 
  READ(1,rec=1)MapVegBamIbis
  MapVegBamIbisNew=MapVegBamIbis
  CLOSE(1,STATUS='KEEP')

 END  SUBROUTINE InitProVeg


 SUBROUTINE RunProVeg()
 IMPLICIT NONE
   CHARACTER(LEN=255) :: FileName
   CHARACTER(LEN=255) :: FileName2
   CHARACTER(LEN=255) :: FileName3
   CHARACTER(LEN=255) :: FileName4
   REAL               :: CoorProVegLat(GridProv(1)%nlins)
   REAL               :: CoorProVegLon(GridProv(1)%ncols)
   REAL               :: CoorProVegDiffLat(GridProv(1)%nlins)
   REAL               :: CoorProVegDiffLon(GridProv(1)%ncols)
   REAL               :: CoorBamGblLat(GridBamg(1)%nlins)
   REAL               :: CoorBamGblLon(GridBamg(1)%ncols)
   REAL               :: CountSSIBVeg(0:13)

   INTEGER            :: nLatFirstPt,nLatLastPt
   INTEGER            :: nLonFirstPt,nLonLastPt
   REAL               :: nlat
   REAL               :: nlon
   INTEGER            :: lrec,ii,jj,i,j,iveg
   FileName ='./outputdata/ProVeg_Brasil_GRR.bin'
   FileName2='./outputdata/BamVeg_Global_GRR.bin'
   FileName3='./outputdata/BamVeg_Global_MERGE_GRR.bin'
   FileName4='./outputdata/ibismsk_new.form'

   INQUIRE(IOLENGTH=lrec)MapVegProVeg
   OPEN(11,FILE=TRIM(FileName),ACCESS='DIRECT',FORM='UNFORMATTED',RECL=lrec,STATUS='UNKNOWN',ACTION='WRITE') 
    DO j=1,GridProv(1)%nlins
       DO i=1,GridProv(1)%ncols
          IF(MapVegProVeg(i,j) < 0.0 .or. MapVegProVeg(i,j) > 13.0 )THEN
             MapVegProVeg(i,j)=GridProv(1)%nodatavalue
          END IF
      END DO
   END DO
   WRITE(11,rec=1)MapVegProVeg
   
   CoorProVegLat(1) = GridProv(1)%Lati
   DO j=2,GridProv(1)%nlins
      CoorProVegLat(j) = CoorProVegLat(j-1) + GridProv(1)%resY
   END DO

   CoorProVegLon(1) = GridProv(1)%Loni
   DO i=2,GridProv(1)%ncols
      CoorProVegLon(i) = CoorProVegLon(i-1) + GridProv(1)%resX
   END DO
   PRINT*,MAXVAL(CoorProVegLat),MINVAL(CoorProVegLat)
   PRINT*,MAXVAL(CoorProVegLon),MINVAL(CoorProVegLon)

   CoorBamGblLat(1) = GridBamg(1)%Lati
   DO j=2,GridBamg(1)%nlins
      CoorBamGblLat(j) = CoorBamGblLat(j-1) + GridBamg(1)%resY
   END DO

   CoorBamGblLon(1) = GridBamg(1)%Loni
   DO i=2,GridBamg(1)%ncols
      CoorBamGblLon(i) = CoorBamGblLon(i-1) + GridBamg(1)%resX
   END DO
   PRINT*,MAXVAL(CoorBamGblLat),MINVAL(CoorBamGblLat)
   PRINT*,MAXVAL(CoorBamGblLon),MINVAL(CoorBamGblLon)


   DO jj=1,GridBamg(1)%nlins
      nlat=CoorBamGblLat(jj) - (GridBamg(1)%resY/2)
      CoorProVegDiffLat=ABS(CoorProVegLat-nlat)
      nLatFirstPt=MINLOC(CoorProVegDiffLat,dim=1)
      nlat=CoorBamGblLat(jj) + (GridBamg(1)%resY/2)
      CoorProVegDiffLat=ABS(CoorProVegLat-nlat)
      nLatLastPt=MINLOC(CoorProVegDiffLat,dim=1)
      DO ii=1,GridBamg(1)%ncols
         nlon=CoorBamGblLon(ii) - (GridBamg(1)%resX/2)
         CoorProVegDiffLon=ABS(CoorProVegLon-nlon)
         nLonFirstPt=MINLOC(CoorProVegDiffLon,dim=1)
         nlon=CoorBamGblLon(ii) + (GridBamg(1)%resX/2)
         CoorProVegDiffLon=ABS(CoorProVegLon-nlon)
         nLonLastPt=MINLOC(CoorProVegDiffLon,dim=1)
         CountSSIBVeg=0.0
         IF(nLatFirstPt /= nLatLastPt .AND. nLonFirstPt/=nLonLastPt)THEN
            DO j=nLatFirstPt,nLatLastPt
               DO i=nLonFirstPt,nLonLastPt
                  IF(MapVegProVeg(i,j) /= GridProv(1)%nodatavalue)THEN
                     iveg=INT(MapVegProVeg(i,j))
                     CountSSIBVeg(iveg) = CountSSIBVeg(iveg)  + 1.0
                  END IF
               END DO
            END DO
             iveg=MAXLOC(CountSSIBVeg(0:13) ,dim=1)-1
             IF(CountSSIBVeg(iveg) /= 0.0)THEN
                IF     (iveg ==1)THEN
                   MapVegBamIbisNew(ii,jj) = 1
                ELSE IF(iveg ==2)THEN
                   MapVegBamIbisNew(ii,jj) = 2
                ELSE IF(iveg ==3)THEN
                   MapVegBamIbisNew(ii,jj) = 3
                ELSE IF(iveg ==4)THEN
                   MapVegBamIbisNew(ii,jj) = 4
                ELSE IF(iveg ==5)THEN
                   MapVegBamIbisNew(ii,jj) = 7
                ELSE IF(iveg ==6)THEN
                   MapVegBamIbisNew(ii,jj) = 9
                ELSE IF(iveg ==7)THEN
                   MapVegBamIbisNew(ii,jj) = 10
                ELSE IF(iveg ==8)THEN
                   MapVegBamIbisNew(ii,jj) = 11
                ELSE IF(iveg ==9)THEN
                   MapVegBamIbisNew(ii,jj) = 12
                ELSE IF(iveg ==10)THEN
                   MapVegBamIbisNew(ii,jj) = 13
                ELSE IF(iveg ==11)THEN
                   MapVegBamIbisNew(ii,jj) = 14
                ELSE IF(iveg ==12)THEN
                   MapVegBamIbisNew(ii,jj) = 10
                ELSE IF(iveg ==13)THEN
                   MapVegBamIbisNew(ii,jj) = 15
                END IF
                PRINT*,'inner',iveg,CountSSIBVeg(iveg) 
             END IF 
         ELSE
             iveg=MAXLOC(CountSSIBVeg(0:13) ,dim=1)-1
         END IF
      END DO
   END DO


   INQUIRE(IOLENGTH=lrec)MapVegBamIbisNew
   OPEN(1,FILE=TRIM(FileName4), ACCESS='DIRECT',FORM='UNFORMATTED',RECL=lrec,STATUS='UNKNOWN',ACTION='WRITE') 
   WRITE(1,rec=1)MapVegBamIbisNew
   CLOSE(1,STATUS='KEEP')




   INQUIRE(IOLENGTH=lrec)MapVegBamIbis
   OPEN(12,FILE=TRIM(FileName2),ACCESS='DIRECT',FORM='UNFORMATTED',RECL=lrec,STATUS='UNKNOWN',ACTION='WRITE') 
   DO jj=1,GridBamg(1)%nlins
      DO ii=1,GridBamg(1)%ncols
         IF(MapVegBamLSMK(ii,jj) == 1.0 )THEN
            MapVegBamIbis(ii,jj)=MapVegBamIbis(ii,jj)
         ELSE
            MapVegBamIbis(ii,jj)=0.0
         END IF
      END DO
   END DO
   WRITE(12,rec=1)MapVegBamIbis
   CLOSE(12,STATUS='KEEP')

   INQUIRE(IOLENGTH=lrec)MapVegBamIbisNew
   OPEN(12,FILE=TRIM(FileName3),ACCESS='DIRECT',FORM='UNFORMATTED',RECL=lrec,STATUS='UNKNOWN',ACTION='WRITE') 
   DO jj=1,GridBamg(1)%nlins
      DO ii=1,GridBamg(1)%ncols
         IF(MapVegBamLSMK(ii,jj) == 1.0 )THEN
            MapVegBamIbisNew(ii,jj)=MapVegBamIbisNew(ii,jj)
         ELSE
            MapVegBamIbisNew(ii,jj)=0.0
         END IF
      END DO
   END DO
   WRITE(12,rec=1)MapVegBamIbisNew
   CLOSE(12,STATUS='KEEP')


 END SUBROUTINE RunProVeg
 
 
 SUBROUTINE  FinalizeProVeg()
 IMPLICIT NONE
   DEALLOCATE(GridProv)
   DEALLOCATE(GridBamg)
   DEALLOCATE(MapVegProVeg)
   DEALLOCATE(MapVegBamIbis)

 END SUBROUTINE  FinalizeProVeg

END MODULE ProVegMap



PROGRAM Main
 USE ProVegMap, Only : InitProVeg,RunProVeg,FinalizeProVeg
 IMPLICIT NONE

 CALL Init()
 CALL Run()
 CALL Finalize()
CONTAINS

 SUBROUTINE Init()
  IMPLICIT NONE
  CALL InitProVeg()
 END  SUBROUTINE Init

 SUBROUTINE Run()
  IMPLICIT NONE
  CALL RunProVeg() 

 END  SUBROUTINE Run

 SUBROUTINE Finalize()
  IMPLICIT NONE
  CALL FinalizeProVeg()
 END  SUBROUTINE Finalize

END PROGRAM Main
