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