SUBROUTINE KMVEGT(TLM0D,TPH0D,DLMD,DPHD + ,WBD,SBD,SM,sice,vegmap) INCLUDE "parmeta" PARAMETER (IMT=IM*2-1,JMT=JM/2+1) PARAMETER (IMJM=IM*JM-JM/2) PARAMETER (IMA=7200,JMA=9000,IMJMA=IMA*JMA) PARAMETER (NTYP = 12) DIMENSION NSIBV(NTYP),NEGRD(IMJM,NTYP),NMAX(IMJM) DIMENSION NSIBG(6,3) ,NEGRDG(IMJM, 6) DIMENSION SM(IM,JM),sice(IM,JM),SMK(IMJM) DIMENSION vegmap(im,jm) REAL ALAT,ALON INTEGER NTY(IMA,JMA) C DATA NSIBV/1,2,3,4,5,6,7,8,9,10,11,12/ DATA (NSIBG(1,I), I=1,3) / 1, 2,0/ DATA (NSIBG(2,I), I=1,3) / 3, 4,0/ DATA (NSIBG(3,I), I=1,3) / 5, 6,0/ DATA (NSIBG(4,I), I=1,3) / 7, 8,0/ DATA (NSIBG(5,I), I=1,3) / 9,10,0/ DATA (NSIBG(6,I), I=1,3) /11,12,0/ C c DATA MISS/-99/,NZERO/0/ DATA MISS/0/,NZERO/0/ C C *** ALAT AND ALON ARE THE GEODETIC LATITUE C *** AND LONGITUDE OF THE 1KM LAT/LON GRID. C *** FOR EACH PAIR OF ALAT/ALON, FIND A K-INDEX ON ETA C *** E GRID H POINT. THEN AT EVERY H POINT, C *** WE HAVE THE NUMBERS OF OB FOR THE C *** 20 VEG TYPES AND THE NUMBERS OF OB FOR THE 6 VEG C *** TYPE GROUPS. AT LAST, DETERMINE THE DOMINANT VEG C *** TYPE IN THE DOMAINT VEG GROUP ON EVERY H GRID POINT C DO K = 1,IMJM NMAX (K)=MISS DO N = 1,NTYP NEGRD(K,N) =NZERO NEGRDG(K,N)=NZERO ENDDO ENDDO DO K=1,IMJM JX=(K-1)/IMT+1 IX=K-(JX-1)*IMT IF(IX.LE.IM)THEN I=IX J=2*JX-1 ELSE I=IX-IM J=2*JX ENDIF SMK(K)=SM(I,J) ENDDO write(65) SMK C*** C READ SEA-LAND MASK FOR 40KM ETA-GRID (1d) C READ(10) SM C C READ MASK in 2-D WITH 1KM RESOLUTION C READ VEGETATION TYPE MASK in 1-D WITH 1KM RESOLUTION C (INTEGER VARIABLE) C NNN=1 NNOB=1 c dxy=0.008333333 !resolucao do mapa de veg. (1km) c CJLG read(44,'(i2)')((NTY(i,j),i=1,ima),j=1,jma) CJLG write(63)NTY n=0 do 210 nj=1,jma alat = -59.99583453 + (nj-1)*dxy do 210 ni=1,ima alon = -89.15898363 + (ni-1)*dxy n=n+1 read(44,'(i2)')NTY(ni,nj) C c IF (NTY(ni,nj) .LE. 0 .OR. NTY(ni,nj) .GT. NTYP )then NNN=NNN+1 GO TO 210 ENDIF C C BE WARE THE WEST IS POSITIVE. C C The comment above is assuming the given longitudes C were positive. But since we have given the western C boundary as negative, we comment the line below C in which the longitude was set to negative C FLAT = ALAT FLON = -ALON C C FIND KINDEX OF OBSERVATION C ALN = FLON APH = FLAT C CJLG CALL TOM( APH, ALN, KINDEX) CALL TOM(DLMD,DPHD,WBD,SBD,APH,ALN,KINDEX) IF(KINDEX .LE. 0 .OR. KINDEX .GT. IMJM) THEN C WRITE(6,6002) KINDEX, APH, ALN, ALAT NNOB=NNOB+1 GOTO 210 ENDIF 6002 FORMAT(1x,'REJECTED: KINDEX, APH, ALN, ALAT',I5,1x,3F11.6) K= KINDEX IF(SMK(K) .gt. 0.5) GO TO 210 C C COUNT THE NUMBERS OF EVERY VEG TYPE IN ONE GRID BOX C DO 500 NV = 1, NTYP IF(NTY(ni,nj) .eq. NSIBV(NV) ) THEN NEGRD(K,NV)=NEGRD(K,NV) + 1 ENDIF 500 CONTINUE C C GROUP THE 20 VEG TYPES INTO 6 AND COUNT THE NUMBER OF EACH GROUP C DO 620 NG = 1, 6 DO 600 NGG = 1, 3 IF (NTY(ni,nj) .EQ. NSIBG(NG,NGG)) THEN NEGRDG(K,NG)= NEGRDG(K,NG)+1 ENDIF 600 CONTINUE 620 CONTINUE 210 CONTINUE c c fim do primeiro loop c DO 220 K=1,IMJM IF(SMK(K) .EQ. 1.0) GO TO 220 NTYY=0 C DETERMINE THE DOMINANT GROUP IN THE GRID BOX(exclude water) C MMAX = -99999 DO 650 NG=1,6 IF (NEGRDG(K,NG) .LE.MMAX .OR. NEGRDG(K,NG) .EQ. NZERO) THEN GO TO 650 ENDIF MMAX=NEGRDG(K,NG) NTYY=NG 650 CONTINUE IF(NTYY .EQ. 0) THEN C write(47,*) 'No available data at the grid point ',K GO TO 220 ENDIF C C NOW DETERMINE THE DOMINANT VEG TYPE IN THE DOMINANT GROUP C MMAX = -99999 NGG=1 DO 302 NV = 1, NTYP IF (NSIBV(NV) .NE. NSIBG(NTYY,NGG)) GO TO 302 IF (NGG .GT. 3) THEN GO TO 220 ELSEIF (NEGRD(K,NV) .LT. MMAX .OR. NEGRD(K,NV) .EQ. NZERO) THEN NGG=NGG+1 GO TO 302 ELSE NGGG=NGG NVV=NV NMAX(K)=NSIBG(NTYY,NGGG) NGG=NGG+1 MMAX=NEGRD(K,NV) GO TO 302 ENDIF 302 CONTINUE c write(48,*) 'K=',K,' NTYY= ',NTYY,' NVV= ',NVV,' NGG= ',NGGG 220 CONTINUE C CJLG WRITE(30) NMAX write(62) NMAX C c c Convert veg_map from k-index into 2-dim index im,jm c DO J=1,JM DO I=1,IM vegmap(I,J)=0. ENDDO ENDDO C DO K=1,IMJM JX=(K-1)/IMT+1 IX=K-(JX-1)*IMT IF(IX.LE.IM)THEN I=IX J=2*JX-1 ELSE I=IX-IM J=2*JX ENDIF vegmap(I,J)=NMAX(K) IF(IX.EQ.IMT)vegmap(I+1,J)=vegmap(I,J) ENDDO write(66)vegmap WRITE(61,1061) (NMAX(K),K=1,IMJM) 1061 FORMAT(1x,16(1x,I4)) C DO 720 NN = 1,IMJM if (NMAX(NN) .eq. MISS) GO TO 720 c WRITE(63,1062) NN, (NEGRDG(NN,N1),N1=1,6), NMAX(NN) c WRITE(62,1063) NN, (NEGRD(NN,MM),MM=1,20) 720 CONTINUE 1062 FORMAT(1x,I5,6(1x,I4), 6x,I4) 1063 FORMAT(1x,I5,20I4) C 9999 WRITE(6,*) 'Data length; ',N-1 WRITE(6,*) 'Missing data number: ', NNN-1 WRITE(6,*) 'Rejected data number: ',NNOB-1 C RETURN END C SUBROUTINE TOM(DLMD,DPHD,WBD,SBD,APHD,ALMD,KMIN) C ****************************************************************** C * * C * ROUTINE FOR CALCULATING THE K INDEX OF THE NEAREST * C * HEIGHT POINT ON THE E GRID TO THAT OF A POINT * C * GIVEN IN GEODETIC COORDINATES * C * LAT POS NORTH LON POS WEST * C ****************************************************************** INCLUDE "parmeta" INTEGER IM REAL DLMD,DPHD,WBD,SBD PARAMETER 1 (IM1=IM-1,NINC=2*IM-1,DTR=1.745329E-2,NMAX=4) DIMENSION AH(NMAX),KH(NMAX) WB=WBD*DTR SB=SBD*DTR DLM=DLMD*DTR DPH=DPHD*DTR C---------------------------------------------------------------------- AIM1=REAL(IM1) RDLM=1./DLM RDPH=1./DPH C C*** CONVERT LOCATION FROM GEODETIC TO TRANSFORMED LAT/LON C CJLG CALL TLLRAD(APHD,ALMD,TLM,TPH) CALL TLLRAD(TLM0D,TPH0D,APHD,ALMD,TLM,TPH) C C*** X1,Y1 ARE THE COORDINATES WITH RESPECT TO THE X,Y TRANSFORMED C*** SYSTEM C X1=(TLM-WB)*RDLM Y1=(TPH-SB)*RDPH C C*** X2,Y2 ARE THE COORDINATES WITH RESPECT TO THE X',Y' TRANSFORMED C*** SYSTEM (TRANSLATED BY AIM1 ALONG Y') C X2=.5*( X1+Y1) Y2=.5*(-X1+Y1)+AIM1 C C*** I2,J2 ARE THE COORDINATES OF THE SOUTHERN HEIGHT POINT IN C*** A CLUSTER OF FOUR C I2=INT(X2) J2=INT(Y2) C P=X2-I2 Q=Y2-J2 PQ=P*Q C JR=J2-IM1 I3=I2-JR J3=I2+JR C*** C*** INDEX K CORRESPONDING TO POINT (I2,J2) (SOUTHERNMOST OF THE C*** FOUR SURROUNDING HEIGHT POINTS) C*** K=J3*IM-J3/2+(I3+2)/2 C*** C*** FIND WHICH OF THE FOUR H POINTS IS NEAREST C*** AH(1)=P*Q AH(2)=(1.-P)*Q AH(3)=P*(1.-Q) AH(4)=(1.-P)*(1.-Q) KH(1)=K KH(2)=K+IM KH(3)=K+IM1 KH(4)=K+NINC C KMIN=KH(1) NM=1 100 NL=NMAX1-NM DO 200 NN=1,NL NX=NM+NN IF(AH(NM).GT.AH(NX))THEN KMIN=KH(NX) NM=NX IF(NX.LT.NMAX1)GO TO 100 ENDIF 200 CONTINUE C---------------------------------------------------------------------- RETURN END C C---------------------------------------------------------------------- SUBROUTINE TLLRAD(TLM0D,TPH0D,APHD,ALMD,TLM,TPH) C ****************************************************************** C * * C * ROUTINE FOR CONVERTING FROM GEODETIC COORDINATES (RADIANS) * C * TO TRANSFORMED COORDINATES (RADIANS) * C * * C ****************************************************************** P A R A M E T E R C 1 (DTR=1.745329E-2) TLM0=TLM0D*DTR TPH0=TPH0D*DTR C--------------------------------------------------------------------- STPH0=SIN(TPH0) CTPH0=COS(TPH0) ALM=ALMD*DTR APH=APHD*DTR RELM=ALM-TLM0 SRLM=SIN(RELM) CRLM=COS(RELM) C SPH=SIN(APH) CPH=COS(APH) CC=CPH*CRLM C ANUM=CPH*SRLM DENOM=CTPH0*CC+STPH0*SPH C---------------------------------------------------------------------- TLM=-ATAN2(ANUM,DENOM) TPH=ASIN(CTPH0*SPH-STPH0*CC) C---------------------------------------------------------------------- RETURN END