!CCC 1. split one file to different forecast hour file according toNCEP standard file naming format !CCC 2. change kpds5=164 in hysplit to kpds5=89 in nesdes grib data format. PROGRAM CONVERT_HYSPLIT IMPLICIT NONE !jp INTEGER, PARAMETER :: JMAX = 534, IMAX = 801, MAXDATA=IMAX*JMAX INTEGER, PARAMETER :: JPDS_DIM = 200 ! default SIZE for JPDS array in libw3.a INTEGER, PARAMETER :: NHOUR = 48 REAL,ALLOCATABLE,DIMENSION(:,:) :: DATA REAL,ALLOCATABLE,DIMENSION(:) :: VAR LOGICAL*1, ALLOCATABLE, DIMENSION(:) :: LB !jf CHARACTER (LEN=120) :: INFILE, OUTFILE CHARACTER (LEN=120), DIMENSION(JPDS_DIM) :: CKPDS INTEGER, DIMENSION(JPDS_DIM) :: KPDS, KGDS INTEGER, DIMENSION( 5) :: KENS INTEGER, DIMENSION( 2) :: KPROB REAL, DIMENSION( 2) :: XPROB INTEGER, DIMENSION(16) :: KCLUST INTEGER, DIMENSION(80) :: KMEMBR INTEGER :: ICNT, ISEEK, LLGRIB, LLSKIP, IERR INTEGER :: I, J, K, L, M, N, JF, IIN, IOUT INTEGER :: IMAX, JMAX, MAXDATA INTEGER :: NLEN REAL :: BMAX, BMIN CHARACTER (LEN=10):: ENVVAR CHARACTER (LEN=16):: FHEAD CHARACTER (LEN=4):: CYCLE !CHARACTER :: FHEAD READ(5,*) FHEAD !READ(5,80) FHEAD !80 FORMAT(A16) READ(5,81) CYCLE 81 FORMAT(A4) READ(5,82) JMAX READ(5,82) IMAX 82 FORMAT(I4) JF = IMAX*JMAX ALLOCATE(VAR(JF)) ALLOCATE(DATA(JF,NHOUR)) ALLOCATE(LB(JF)) !INFILE='grib_pbl.1hr_227' !ZZ CYCLE='t06z' !WRITE(*,*) TRIM(INFILE),' ', TRIM(CYCLE) IIN=10 ENVVAR='XLFUNIT_ ' WRITE(ENVVAR(9:10),FMT='(I2)') IIN print*,"ENVVAR=",ENVVAR CALL GETENV(ENVVAR,INFILE) print*,"INFILE=",INFILE CALL BAOPENR(IIN,INFILE,IERR) IF ( IERR /= 0 ) STOP 'GB OPEN FAILED' ISEEK=0 CALL SKGB(IIN,ISEEK,LLGRIB,LLSKIP) DO WHILE ( LLGRIB > 0 ) ! llgrib CALL RDGB(IIN,LLGRIB,LLSKIP,KPDS,KGDS,JF,LB,VAR) CALL SKGB(IIN,ISEEK,LLGRIB,LLSKIP) DO I = 1, NHOUR IF ( KPDS(5) == 164 .AND. KPDS(15) == I ) THEN DO J = 1, 25 WRITE(CKPDS(J),*)KPDS(J) END DO WRITE(*,'(''RD '',I2.2,'' KPDS = '', 25(A))') I,(TRIM(CKPDS(J)),J=1,25) DATA(:,I)=VAR(:) EXIT END IF END DO ENDDO ! end of while loop for llgrib CALL BACLOSE(IIN,IERR) IF ( IERR /= 0 ) STOP 'GB CLOSE FAILED' DO I = 1, NHOUR IOUT = I + 50 WRITE(OUTFILE,'(A,".",A,''.f'',I2.2)')trim(FHEAD), CYCLE, I WRITE(*,'(''PACKING '',A)') TRIM(OUTFILE) CALL BAOPENW(IOUT,OUTFILE,IERR) KPDS(5) = 89 KPDS(6) = 105 KPDS(7) = 5000 KPDS(14) = I KPDS(22) = 5 DO J = 1, JF IF ( DATA(J,I) < -90. ) THEN !jp VAR(J)= -99.0 ! conditional cerification VAR(J)= 0.0 ! un-conditional verification ELSE VAR(J)=10.0**(DATA(J,I)) ! convert from log10 C -> concentration C ! unit: micro-gram/m3 END IF END DO KPROB(1) = -1 !: OCT 46 KPROB(2) = -1 !: OCT 47 XPROB(1) = -1. !: OCT 48-51 XPROB(2) = -1. !: OCT 52-55 KENS(1) = -1 !: OCT 41 KENS(2) = -1 !: OCT 42 KENS(3) = -1 !: OCT 43 KENS(4) = -1 KENS(5) = 255 KCLUST(1)= 1 DO J = 1, 25 WRITE(CKPDS(J),*)KPDS(J) END DO WRITE(*,'(''WRT KPDS = '', 25(A))') (TRIM(CKPDS(J)),J=1,25) CALL PUTGBEX(IOUT,JF,KPDS,KGDS,KENS,KPROB,XPROB,KCLUST,KMEMBR,LB,VAR,IERR) CALL BACLOSE(IOUT,IERR) END DO CONTAINS !********************************************************************** SUBROUTINE RDGB(LUGB,LGRIB,LSKIP,KPDS,KGDS,NDATA,LBMS,DATA) IMPLICIT NONE ! ! READ GRIB FILE ! INPUT ! LUGB - LOGICAL UNIT TO READ ! LGRIB - LENGTH OF GRIB RECORD ! LSKIP - BYTES TO SKIP FOR GRIB RECORD ! OUTPUT ! KPDS(22) - UNPACKED PRODUCT DEFINITION SECTION ! KGDS(20) - UNPACKED GRID DEFINITION SECTION ! NDATA - NUMBER OF DATA POINTS ! LBMS(NDATA) - LOGICAL BIT MAP ! DATA(NDATA) - DATA UNPACKED ! INTEGER, DIMENSION(:), INTENT(OUT) :: KPDS, KGDS LOGICAL*1, DIMENSION(:), INTENT(OUT) :: LBMS REAL, DIMENSION(:), INTENT(OUT) :: DATA INTEGER, INTENT(IN) :: LUGB, LSKIP, LGRIB INTEGER, INTENT(OUT) :: NDATA CHARACTER (LEN=1), DIMENSION(LGRIB) :: GRIB INTEGER, DIMENSION(20) :: KPTR INTEGER :: LREAD, IRET NDATA=0 CALL BAREAD(LUGB,LSKIP,LGRIB,LREAD,GRIB) IF ( LREAD < LGRIB ) RETURN CALL W3FI63(GRIB,KPDS,KGDS,LBMS,DATA,KPTR,IRET) IF(IRET /= 0) RETURN NDATA=KPTR(10) RETURN END SUBROUTINE RDGB !******************************************** !******************************************** SUBROUTINE SKGB(LUGB,ISEEK,LGRIB,LSKIP) IMPLICIT NONE ! ! SEEK FOR NEXT GRIB1 RECORD WITHIN THE NEXT LSEEK=4096 BYTES ! INPUT ! LUGB - LOGICAL UNIT TO READ ! ISEEK - BYTES TO SKIP BEFORE SEARCH (SET TO 0 AT START) ! OUTPUT ! ISEEK - NUMBER OF BYTES READ SO FAR ! LGRIB - LENGTH OF GRIB RECORD (0 IF NOT FOUND) ! LSKIP - BYTES TO SKIP FOR GRIB RECORD ! INTEGER, PARAMETER :: LSEEK=4096 CHARACTER (LEN=LSEEK) :: C INTEGER, INTENT(IN) :: LUGB INTEGER, INTENT(INOUT) :: ISEEK INTEGER, INTENT(OUT):: LSKIP, LGRIB INTEGER :: I, LREAD, IRET, MOVA2I CALL BAREAD(LUGB,ISEEK,LSEEK,LREAD,C) DO I=0,LREAD-8 IF(C(I+1:I+4) == 'GRIB'.AND.MOVA2I(C(I+8:I+8)) == 1) THEN LGRIB=MOVA2I(C(I+5:I+5))*65536 & +MOVA2I(C(I+6:I+6))*256 & +MOVA2I(C(I+7:I+7)) LSKIP=ISEEK+I ISEEK=LSKIP+LGRIB RETURN ENDIF ENDDO LGRIB=0 LSKIP=0 ISEEK=ISEEK+LSEEK ! CALL W3TAGE('SREF_COM_GRIB') RETURN END SUBROUTINE SKGB END PROGRAM CONVERT_HYSPLIT