cnv21.F90

Go to the documentation of this file.

 
 
subroutine cnv21(ifl1,ifl2)
 
  use grib_mod
  use params
  integer,intent(in) :: ifl1,ifl2
 
  CHARACTER(len=1),allocatable,dimension(:) :: cgrib
  type(gribfield) :: gfld
  integer,dimension(200) :: jids,jpdt,jgdt
  integer :: kpds(200),kgds(200),kens(200),kprob(2)
  integer :: kclust(16),kmembr(80)
  integer :: currlen=0
  integer :: igds(5)=(/0,0,0,0,0/)
  real :: xprob(2)
  logical :: unpack=.true.
  !
  !
  ifli1=0
  jdisc=-1
  jids=-9999
  jpdt=-9999
  jgdt=-9999
  jpdtn=-1
  jgdtn=-1
  !
  icount=0
  jskp=0
  do
     call getgb2(ifl1,ifli1,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt, &
          unpack,jskp,gfld,iret)
     if (iret.ne.0) then
        if (iret.eq.99) exit
        print *,' getgb2 error = ',iret
        cycle
        !call errexit(17)
     endif
     icount=icount+1
     !
     !  Ensure that cgrib array is large enough
     !
     newlen=4*gfld%ngrdpts
     if (newlen.gt.currlen) then
        if (allocated(cgrib)) deallocate(cgrib)
        allocate(cgrib(newlen),stat=is)
        currlen=newlen
     endif
     !
     !   Construct GDS
     !
     igds(1)=gfld%griddef
     igds(2)=gfld%ngrdpts
     igds(3)=gfld%numoct_opt
     igds(4)=gfld%interp_opt
     igds(5)=gfld%igdtnum
     if (.NOT. associated(gfld%list_opt)) &
          allocate(gfld%list_opt(1))
     call gdt2gds(igds,gfld%igdtmpl,gfld%num_opt,gfld%list_opt, &
          kgds,igrid,iret)
     if (iret.ne.0) then
        print *,'cnv21: could not create gds'
        cycle
     endif
     !print *,' SAGT: NCEP GRID: ',igrid
     !
     !   Construct PDS
     !
     call makepds(gfld%discipline,gfld%idsect,gfld%ipdtnum, &
          gfld%ipdtmpl,gfld%ibmap,gfld%idrtnum, &
          gfld%idrtmpl,kpds,iret)
     if (iret.ne.0) then
        print *,'cnv21: could not create pds in GRIB1'
        cycle
     endif
     kpds(3)=igrid
     !
     !  Check for Coastal Ocean circulation and UKMET grib grid id.
     !  ON 388 defined grid id 238 same as grid 244
     !  If the process model is 45, and UK Met(74), the grid id is 2 or 45
     !  If the process model is 121, the grid id is 238
     !  If the process model is 123, the grid id is 244
     !
     if (kpds(1).eq.7.AND.kpds(2).eq.121) kpds(3)=238
     if (kpds(1).eq.7.AND.kpds(2).eq.123) kpds(3)=244
     if (kpds(1).eq.74) then
        if (kpds(2).eq.45.AND.kpds(3).eq.2)  kpds(3)=2
        if (kpds(2).eq.15.AND.kpds(3).eq.45) kpds(3)=45
        if (kpds(2).eq.45.AND.kpds(3).eq.45) kpds(3)=45
     end if
     !
     !   Construct Ensemble info, if necessary
     !
     if ((gfld%ipdtnum.ge.1.AND.gfld%ipdtnum.le.6).OR. &
          (gfld%ipdtnum.ge.9.AND.gfld%ipdtnum.le.14)) then
        call makepdsens(gfld%ipdtnum,gfld%ipdtmpl,kpds,kens,kprob, &
             xprob,kclust,kmembr,iret)
     endif
     !
     !   If not using bit-map, must assign dummy bit-map
     !
     if (gfld%ibmap.ne.0 .AND. gfld%ibmap.ne.254) then
        !gfld%bmap => dummy
        if ((gfld%idrtnum.eq.2 .OR. gfld%idrtnum.eq.3) .AND. &
             gfld%idrtmpl(7).ne.0) then       ! convert missings to bitmap
           allocate(gfld%bmap(gfld%ngrdpts))
           kpds(4)=ior(kpds(4),64)
           if (gfld%idrtmpl(7).eq.1) then
              call rdieee(gfld%idrtmpl(8),rmiss1,1)
              do i=1,gfld%ngrdpts
                 if (gfld%fld(i) .eq. rmiss1) then
                    gfld%bmap(i)=.false.
                 else
                    gfld%bmap(i)=.true.
                 endif
              enddo
           endif
           if (gfld%idrtmpl(7).eq.2) then
              call rdieee(gfld%idrtmpl(8),rmiss1,1)
              call rdieee(gfld%idrtmpl(9),rmiss2,1)
              do i=1,gfld%ngrdpts
                 if (gfld%fld(i).eq.rmiss1 .OR. &
                      gfld%fld(i).eq.rmiss2) then
                    gfld%bmap(i)=.false.
                 else
                    gfld%bmap(i)=.true.
                 endif
              enddo
           endif
        endif
        if ((gfld%idrtnum.eq.2 .OR. gfld%idrtnum.eq.3) .AND. &
             gfld%idrtmpl(7).eq.0) then       ! convert missings to bitmap
           allocate(gfld%bmap(gfld%ngrdpts))
           kpds(4)=ior(kpds(4),64)
           call rdieee(gfld%idrtmpl(8),rmiss1,1)
           if (rmiss1 .lt. -9999.0) then
              rmiss1=rmiss1*10.0
           else
              rmiss1=-9999.0
           endif
           do i=1,gfld%ngrdpts
              if (gfld%fld(i) .eq. rmiss1) then
                 gfld%bmap(i)=.false.
              else
                 gfld%bmap(i)=.true.
              endif
           enddo
        endif
     endif
     !
     !   Pack and write GRIB 1 field
     !
     ibs=gfld%idrtmpl(2)
     !print *,'SAGT:before putgbexn'
     if (.NOT. associated(gfld%bmap)) allocate(gfld%bmap(1))
     imug=0
     call putgbexn(ifl2,gfld%ngrdpts,kpds,kgds,kens,kprob, &
          xprob,kclust,kmembr,ibs,imug,gfld%bmap, &
          gfld%fld,iret)
     !print *,'SAGT:after putgbexn'
     if (iret.ne.0) then
        print *,' putgbexn error = ',iret
        cycle
        !call errexit(17)
     endif
 
     call gf_free(gfld)
 
  enddo
 
  if (allocated(cgrib)) deallocate(cgrib)
 
  return
end subroutine cnv21
 
subroutine makepds(idisc,idsect,ipdsnum,ipdstmpl,ibmap, &
    idrsnum,idrstmpl,kpds,iret)
 
  use params
 
  integer,intent(in) :: idsect(*),ipdstmpl(*),idrstmpl(*)
  integer,intent(in) :: ipdsnum,idisc,idrsnum,ibmap
  integer,intent(out) :: kpds(*)
  integer,intent(out) :: iret
 
  iret=0
  ipos=0
  kpds(1:24)=0
  if ( (ipdsnum.lt.0).OR.(ipdsnum.gt.15) ) then
    print *,'makepds: Don:t know GRIB2 PDT 4.',ipdsnum
    iret=2
    return
  endif
 
  kpds(1)=idsect(1)
  kpds(2)=ipdstmpl(5)
  kpds(3)=255
  kpds(4)=128
  if ( ibmap.ne.255 ) kpds(4)=kpds(4)+64
  if ( ibmap.ge.1.AND.ibmap.le.253 ) then
    print *,'makepds: Don:t know about predefined bit-map ',ibmap
    iret=1
    return
  endif
  call param_g2_to_g1(idisc,ipdstmpl(1),ipdstmpl(2),kpds(5), &
      kpds(19))
  !
  !  Special parameters for ICAO WAFS (Max Icing, TP and CAT)
  !
  If (ipdstmpl(16).eq.2.and.ipdstmpl(1).eq.19.and. &
      ipdstmpl(2).eq.20) kpds(5) = 169
  If (ipdstmpl(16).eq.2.and.ipdstmpl(1).eq.19.and. &
      ipdstmpl(2).eq.21) kpds(5) = 171
  If (ipdstmpl(16).eq.2.and.ipdstmpl(1).eq.19.and. &
      ipdstmpl(2).eq.22) kpds(5) = 173
  !
  !  Special parameters for NAM (NMMB)
  !
  If (idisc.eq.0.and.ipdstmpl(1).eq.2) then
    if (ipdstmpl(2).eq.220) then
        kpds(5) = 237
        kpds(19) = 129
    end if
    if (ipdstmpl(2).eq.221) then
        kpds(5) = 238
        kpds(19) = 129
    end if
    if (ipdstmpl(2).eq.222) then
        kpds(5) = 253
        kpds(19) = 129
    end if
    if (ipdstmpl(2).eq.223) then
        kpds(5) = 254
        kpds(19) = 129
    end if
  endif
  !
  If (idisc.eq.0.and.ipdstmpl(2).eq.16 &
      .and.ipdstmpl(3).eq.198) then
    kpds(5) = 235
    kpds(19) = 129
  endif
  !
  If (idisc.eq.0.and.ipdstmpl(2).eq.7 &
      .and.ipdstmpl(3).eq.199) then
    kpds(5) = 236
    kpds(19) = 129
  endif
  !
  !  Special parameters for ICAO Height at CB Base and Top
  !  in GRIB1 Table 140
  !
  If (ipdstmpl(1).eq.3.and.ipdstmpl(2).eq.3) then
    If (ipdstmpl(10).eq.11) then
        kpds(19) = 140
        kpds(5)  = 179
    end if
    If (ipdstmpl(10).eq.12) then
        kpds(19) = 140
        kpds(5)  = 180
    end if
  end if
  !
  call levelcnv(ipdstmpl,kpds(6),kpds(7))      ! level
  kpds(8)=mod(idsect(6),100)
  if ( kpds(8).eq.0 ) kpds(8)=100
  kpds(9)=idsect(7)                            ! Year
  kpds(10)=idsect(8)                           ! Month
  kpds(11)=idsect(9)                           ! Day
  kpds(12)=idsect(10)                          ! Hour
  if ( ipdstmpl(8).ne.13 ) then
    kpds(13)=ipdstmpl(8)                      ! Time Unit
  else
    kpds(13)=254
  endif
  kpds(14)=ipdstmpl(9)                         ! P1
  if ( ipdsnum.le.7 ) then                     ! P2
    kpds(15)=0
    kpds(16)=0
    kpds(20)=0
    if ( kpds(14).eq.0 ) kpds(16)=1
    if ( kpds(14).gt.255 ) kpds(16)=10
    if ( ipdstmpl(5).eq.77.OR.ipdstmpl(5).eq.81.OR. &
          ipdstmpl(5).eq.96.OR.ipdstmpl(5).eq.80.OR. &
          ipdstmpl(5).eq.82.OR.ipdstmpl(5).eq.120.OR. &
          ipdstmpl(5).eq.47.OR.ipdstmpl(5).eq.11 ) then
        kpds(16)=10
    end if
    if (ipdstmpl(5).eq.84.AND.kpds(5).eq.154)kpds(16) = 10
    !
    !          NOAA Wave Watch III and Coastal Ocean Circulation
    !          and Alaska Waters Regional Wave Model
    !
    if ( ipdstmpl(5).eq.88.OR.ipdstmpl(5).eq.121 &
          .OR.ipdstmpl(5).eq.122.OR.ipdstmpl(5).eq.123 &
          .OR.ipdstmpl(5).eq.124.OR.ipdstmpl(5).eq.125 &
          .OR.ipdstmpl(5).eq.131.OR.ipdstmpl(5).eq.45 &
          .OR.ipdstmpl(5).eq.11 ) then
        kpds(16) = 0
        !
        ! Level Surface set to 1
        !
        if (kpds(5).eq.80.OR.kpds(5).eq.82.OR. &
            kpds(5).eq.88.OR.kpds(5).eq.49.OR. &
            kpds(5).eq.50) kpds(7)=1  ! Level Surface
        if (ipdstmpl(5).eq.122.OR.ipdstmpl(5).eq.124.OR. &
            ipdstmpl(5).eq.131.OR.ipdstmpl(5).eq.123.OR. &
            ipdstmpl(5).eq.125.OR.ipdstmpl(5).eq.88.OR. &
            ipdstmpl(5).eq.121) kpds(7)=1
        if (idsect(1).eq.54.AND.ipdstmpl(5).eq.45) kpds(16) = 10
    endif
  else
    selectcase (ipdsnum)
    case(8)
        ipos=24
    case(9)
        ipos=31
    case(10)
        ipos=25
    case(11)
        ipos=27
    case(12)
        ipos=26
    case(13)
        ipos=40
    case(14)
        ipos=39
    end select
    kpds(15)=ipdstmpl(ipos+3)+kpds(14)  ! P2
    selectcase (ipdstmpl(ipos))
    case (255)
        kpds(16)=2
    case (0)
        kpds(16)=3
    case (1)
        kpds(16)=4
    case (2)
        kpds(16)=2
    case (3)
        kpds(16)=2
    case (4)
        kpds(16)=5
    case (51)
        kpds(16)=51
    end select
    kpds(20)=ipdstmpl(ipos-1)
  endif
  if (ipdstmpl(9) .ge. 252) then
    if (ipdstmpl(ipos+3).eq.3) then
        kpds(13)= 10                          ! Forecast time unit is 3-hour
        kpds(14)=ipdstmpl(9)/3                ! Time range P1
        kpds(15)=ipdstmpl(ipos+3)/3+kpds(14)  ! Time range P2
    else if (ipdstmpl(ipos+3).eq.6) then
        kpds(13)= 11                          ! Forecast time unit is 6-hour
        kpds(14)=ipdstmpl(9)/6                ! Time range P1
        kpds(15)=ipdstmpl(ipos+3)/6+kpds(14)  ! Time range P2
    else if (ipdstmpl(ipos+3).eq.12) then
        kpds(13)= 12                          ! Forecast time unit is 12-hour
        kpds(14)=ipdstmpl(9)/12               ! Time range P1
        kpds(15)=ipdstmpl(ipos+3)/12+kpds(14) ! Time range P2
    end if
  end if
  if (ipdsnum .eq. 8 .AND. ipdstmpl(9) .eq. 0) then
    if (ipdstmpl(ipos+3).ge.252) then
        kpds(13)= 10                          ! Forecast time unit is hour
        kpds(14)=ipdstmpl(9)/3                ! Time range P1
        kpds(15)=ipdstmpl(ipos+3)/3+kpds(14)  ! Time range P2
    end if
  end if
  !
  !  Checking total preciptation for 15-hr or 18-hr or 21-hr or 24-hr accumulation
  !  after forecast hour F240
  !
  if (ipdstmpl(9) .ge. 240 )then
    if ( ipdstmpl(ipos+3).eq.15 .OR. ipdstmpl(ipos+3).eq.18 &
          .OR. ipdstmpl(ipos+3).eq.21 .OR. &
          ipdstmpl(ipos+3).eq.24 ) then
        kpds(13)= 10                          ! Forecast time unit is 3-hour
        kpds(14)=ipdstmpl(9)/3                ! Time range P1
        kpds(15)=ipdstmpl(ipos+3)/3+kpds(14)  ! Time range P2
    end if
  end if
  !
  !   Checking Unit of Time Range for FNMOC (APCP)
  !
  if (ipdstmpl(4).eq.58 .AND. ipdsnum.eq.11 .AND. &
      (ipdstmpl(1).eq.1 .AND.ipdstmpl(2).eq.8) &
      .AND. (ipdstmpl(10).eq.1)) then
    if (ipdstmpl(9) .ge. 252) then
        kpds(13)= 11      !  Forecast time unit is 6-hour
        kpds(14)=ipdstmpl(9)/6      ! Time range P1
        kpds(15)=ipdstmpl(ipos+3)/6+kpds(14)  ! Time range P2
    else
        kpds(13)= 1       !  Forecast time unit is 1-hour
        kpds(14)=ipdstmpl(9)  ! Time range P1
    end if
  end if
  !
  !   Special case for FNMOC (TMAX and TMIN)
  !
  if (ipdstmpl(4).eq.58 .AND. ipdsnum.eq.11 .AND. &
      (ipdstmpl(1).eq.0 &
      .AND.ipdstmpl(2).eq.0).AND.(ipdstmpl(10).eq.103)) then
    kpds(16) = 2
    !   For Maximum Temperature
    If (ipdstmpl(27).eq.2 .AND. ipdstmpl(1).eq.0 .AND. &
          ipdstmpl(2).eq.0) kpds(5) = 15
    !   For Minimum Temperature
    If (ipdstmpl(27).eq.3 .AND. ipdstmpl(1).eq.0 .AND. &
          ipdstmpl(2).eq.0) kpds(5) = 16
  end if
  !
  !   Special case for WAFS (Mean/MAx IP,CTP and CAT)
  !
  if (ipdstmpl(5).eq.96.AND.((ipdstmpl(1).eq.19) &
      .AND.(ipdstmpl(2).eq.20.or.ipdstmpl(2).eq.21.or. &
      ipdstmpl(2).eq.22)).AND.(ipdstmpl(10).eq.100)) then
    kpds(16) = 10
  end if
  !
  kpds(17)=0
  kpds(18)=1                                   ! GRIB edition
  kpds(21)=(idsect(6)/100)+1                   ! Century
  if ( kpds(8).eq.100 ) kpds(21)=idsect(6)/100
  kpds(22)=idrstmpl(3)                         ! Decimal scale factor
  kpds(23)=idsect(2)                           ! Sub-center
  return
end subroutine makepds
 
subroutine levelcnv(ipdstmpl,ltype,lval)
  integer,intent(in) :: ipdstmpl(*)
  integer,intent(out) :: ltype,lval
 
  ltype=255
  lval=0
  ltype1=ipdstmpl(10)
  ltype2=ipdstmpl(13)
 
  if ( ltype1.eq.10.AND.ltype2.eq.255 ) then
     ltype=200
     lval=0
  elseif ( ltype1.eq.11.AND.ltype2.eq.255 ) then
     ltype=216
     lval=0
  elseif ( ltype1.eq.12.AND.ltype2.eq.255 ) then
     ltype=217
     lval=0
  elseif ( ltype1.lt.100.AND.ltype2.eq.255 ) then
     ltype=ltype1
     lval=0
  elseif ( ltype1.eq.1.AND.ltype2.eq.8 ) then
     ltype=ltype1
     lval=0
  elseif ( ltype1.eq.10.AND.ltype2.eq.255 ) then
     ltype=200
     lval=0
  elseif ( ltype1.eq.235.AND.ltype2.eq.255 ) then
     ltype=235
     rscal1=10.**(-ipdstmpl(11))
     lval=nint(real(ipdstmpl(12))*rscal1)
  elseif ( ltype1.ge.200.AND.ltype2.eq.255 ) then
     ltype=ltype1
     lval=0
  elseif (ltype1.eq.100.AND.ltype2.eq.255 ) then
     ltype=100
     rscal1=10.**(-ipdstmpl(11))
     lval=nint(real(ipdstmpl(12))*rscal1/100.)
  elseif (ltype1.eq.100.AND.ltype2.eq.100 ) then
     ltype=101
     rscal1=10.**(-ipdstmpl(11))
     lval1=nint(real(ipdstmpl(12))*rscal1/1000.)
     rscal2=10.**(-ipdstmpl(14))
     lval2=nint(real(ipdstmpl(15))*rscal2/1000.)
     lval=(lval1*256)+lval2
  elseif (ltype1.eq.101.AND.ltype2.eq.255 ) then
     ltype=102
     lval=0
  elseif (ltype1.eq.102.AND.ltype2.eq.255 ) then
     ltype=103
     rscal1=10.**(-ipdstmpl(11))
     lval=nint(real(ipdstmpl(12))*rscal1)
  elseif (ltype1.eq.102.AND.ltype2.eq.102 ) then
     ltype=104
     rscal1=10.**(-ipdstmpl(11))
     lval1=nint(real(ipdstmpl(12))*rscal1)
     rscal2=10.**(-ipdstmpl(14))
     lval2=nint(real(ipdstmpl(15))*rscal2)
     lval=(lval1*256)+lval2
  elseif (ltype1.eq.103.AND.ltype2.eq.255 ) then
     ltype=105
     rscal1=10.**(-ipdstmpl(11))
     lval=nint(real(ipdstmpl(12))*rscal1)
  elseif (ltype1.eq.103.AND.ltype2.eq.103 ) then
     ltype=106
     rscal1=10.**(-ipdstmpl(11))
     lval1=nint(real(ipdstmpl(12))*rscal1/100.)
     rscal2=10.**(-ipdstmpl(14))
     lval2=nint(real(ipdstmpl(15))*rscal2/100.)
     lval=(lval1*256)+lval2
  elseif (ltype1.eq.104.AND.ltype2.eq.255 ) then
     ltype=107
     rscal1=10.**(-ipdstmpl(11))
     lval=nint(real(ipdstmpl(12))*rscal1*10000.)
  elseif (ltype1.eq.104.AND.ltype2.eq.104 ) then
     ltype=108
     rscal1=10.**(-ipdstmpl(11))
     lval1=nint(real(ipdstmpl(12))*rscal1*100.)
     rscal2=10.**(-ipdstmpl(14))
     lval2=nint(real(ipdstmpl(15))*rscal2*100.)
     lval=(lval1*256)+lval2
  elseif (ltype1.eq.105.AND.ltype2.eq.255 ) then
     ltype=109
     lval=ipdstmpl(12)
  elseif (ltype1.eq.105.AND.ltype2.eq.105 ) then
     ltype=110
     rscal1=10.**(-ipdstmpl(11))
     lval1=nint(real(ipdstmpl(12))*rscal1)
     rscal2=10.**(-ipdstmpl(14))
     lval2=nint(real(ipdstmpl(15))*rscal2)
     lval=(lval1*256)+lval2
  elseif (ltype1.eq.106.AND.ltype2.eq.255 ) then
     ltype=111
     rscal1=10.**(-ipdstmpl(11))
     lval=nint(real(ipdstmpl(12))*rscal1*100.)
  elseif (ltype1.eq.106.AND.ltype2.eq.106 ) then
     ltype=112
     rscal1=10.**(-ipdstmpl(11))
     lval1=nint(real(ipdstmpl(12))*rscal1*100.)
     rscal2=10.**(-ipdstmpl(14))
     lval2=nint(real(ipdstmpl(15))*rscal2*100.)
     lval=(lval1*256)+lval2
  elseif (ltype1.eq.107.AND.ltype2.eq.255 ) then
     ltype=113
     rscal1=10.**(-ipdstmpl(11))
     lval=nint(real(ipdstmpl(12))*rscal1)
  elseif (ltype1.eq.107.AND.ltype2.eq.107 ) then
     ltype=114
     rscal1=10.**(-ipdstmpl(11))
     lval1=475-nint(real(ipdstmpl(12))*rscal1)
     rscal2=10.**(-ipdstmpl(14))
     lval2=475-nint(real(ipdstmpl(15))*rscal2)
     lval=(lval1*256)+lval2
  elseif (ltype1.eq.108.AND.ltype2.eq.255 ) then
     ltype=115
     rscal1=10.**(-ipdstmpl(11))
     lval=nint(real(ipdstmpl(12))*rscal1/100.)
  elseif (ltype1.eq.108.AND.ltype2.eq.108 ) then
     ltype=116
     rscal1=10.**(-ipdstmpl(11))
     lval1=nint(real(ipdstmpl(12))*rscal1/100.)
     rscal2=10.**(-ipdstmpl(14))
     lval2=nint(real(ipdstmpl(15))*rscal2/100.)
     lval=(lval1*256)+lval2
  elseif (ltype1.eq.109.AND.ltype2.eq.255 ) then
     ltype=117
     rscal1=10.**(-ipdstmpl(11))
     lval=nint(real(ipdstmpl(12))*rscal1*1000000000.)
  elseif (ltype1.eq.111.AND.ltype2.eq.255 ) then
     ltype=119
     rscal1=10.**(-ipdstmpl(11))
     lval=nint(real(ipdstmpl(12))*rscal1*10000.)
  elseif (ltype1.eq.111.AND.ltype2.eq.111 ) then
     ltype=120
     rscal1=10.**(-ipdstmpl(11))
     lval1=nint(real(ipdstmpl(12))*rscal1*100.)
     rscal2=10.**(-ipdstmpl(14))
     lval2=nint(real(ipdstmpl(15))*rscal2*100.)
     lval=(lval1*256)+lval2
  elseif (ltype1.eq.160.AND.ltype2.eq.255 ) then
     ltype=160
     rscal1=10.**(-ipdstmpl(11))
     lval=nint(real(ipdstmpl(12))*rscal1)
  elseif ((ltype1.ge.236.AND.ltype1.le.239).AND. &
       (ltype2.ge.236.AND.ltype2.le.239)) then
     ltype=ltype1
     rscal1=10.**(-ipdstmpl(11))
     lval1=nint(real(ipdstmpl(12))*rscal1)
     rscal2=10.**(-ipdstmpl(14))
     lval2=nint(real(ipdstmpl(15))*rscal2)
     lval=(lval1*256)+lval2
  else
     print *,'levelcnv: GRIB2 Levels ',ltype1,ltype2, &
          ' not recognized.'
     ltype=255
  endif
 
  !  High resolution stuff
  !        elseif (ltype.eq.121) then
  !           ipdstmpl(10)=100
  !           ipdstmpl(12)=(1100+(lval/256))*100
  !           ipdstmpl(13)=100
  !           ipdstmpl(15)=(1100+mod(lval,256))*100
  !        elseif (ltype.eq.125) then
  !           ipdstmpl(10)=103
  !           ipdstmpl(11)=-2
  !           ipdstmpl(12)=lval
  !        elseif (ltype.eq.128) then
  !           ipdstmpl(10)=104
  !           ipdstmpl(11)=-3
  !           ipdstmpl(12)=1100+(lval/256)
  !           ipdstmpl(13)=104
  !           ipdstmpl(14)=-3
  !           ipdstmpl(15)=1100+mod(lval,256)
  !        elseif (ltype.eq.141) then
  !           ipdstmpl(10)=100
  !           ipdstmpl(12)=(lval/256)*100
  !           ipdstmpl(13)=100
  !           ipdstmpl(15)=(1100+mod(lval,256))*100
 
  return
end subroutine levelcnv
 
subroutine makepdsens(ipdsnum,ipdstmpl,kpds,kens,kprob, &
    xprob,kclust,kmembr,iret)
  use params
 
  integer,intent(in) :: ipdstmpl(*)
  integer,intent(in) :: ipdsnum
  integer,intent(inout) :: kpds(*)
  integer,intent(out) :: kens(5),kprob(2)
  integer,intent(out) :: kclust(16),kmembr(80)
  real,intent(out) :: xprob(2)
  integer,intent(out) :: iret
 
  iret=0
  kpds(23)=2          !  subcenter = ensemble
 
  kens(1:5)=0
  kprob(1:2)=0
  xprob(1:2)=0.
  kclust(1:16)=0
  kmembr(1:80)=0
  !
  !  Individual Ensemble Fcst
  !
  if (ipdsnum.eq.1.OR.ipdsnum.eq.11) then
    kens(1)=1
    selectcase (ipdstmpl(16))
    case(0)
      kens(2)=1
      kens(3)=1
    case(1)
      kens(2)=1
      kens(3)=2
    case(2)
      kens(2)=2
      kens(3)=ipdstmpl(17)
    case(3)
      kens(2)=3
      kens(3)=ipdstmpl(17)
    case(4)
      kens(2)=3
      kens(3)=ipdstmpl(17)
    case(192)
      kens(2)=3
      kens(3)=ipdstmpl(17)
    end select
    kens(4)=1
    kens(5)=255
    !
    !  Probability Fcst
    !
  elseif (ipdsnum.eq.5.OR.ipdsnum.eq.9) then
    kens(1)=1
    kens(2)=5
    kens(3)=0
    kens(4)=0
    kens(5)=255
    kprob(1)=kpds(5)
    kpds(5)=191
    kprob(2)=ipdstmpl(18)+1
    if (kprob(2).eq.1) then
      rscale=10.**ipdstmpl(19)
      xprob(1)=real(ipdstmpl(20))/rscale
      xprob(2)=0.0
    elseif (kprob(2).eq.2) then
      xprob(1)=0.0
      rscale=10.**ipdstmpl(21)
      xprob(2)=real(ipdstmpl(22))/rscale
    elseif (kprob(2).eq.3) then
      rscale=10.**ipdstmpl(19)
      xprob(1)=real(ipdstmpl(20))/rscale
      rscale=10.**ipdstmpl(21)
      xprob(2)=real(ipdstmpl(22))/rscale
    endif
    kclust(1)=ipdstmpl(17)
    !
    !  Derived Ensemble Fcst
    !
  elseif (ipdsnum.eq.2.OR.ipdsnum.eq.12) then
    kens(1)=1
    kens(2)=5
    kens(3)=0
    selectcase (ipdstmpl(16))
    case(0)
      kens(4)=1
    case(1)
      kens(4)=2
    case(2)
      kens(4)=11
    case(3)
      kens(4)=12
    end select
    !kens(5)=89
    kens(5)=0
    kclust(1)=ipdstmpl(17)
  else
    print *,'makepdsens: Don:t know GRIB2 PDT 4.',ipdsnum
    iret=2
  endif
 
  return
end subroutine makepdsens
 
subroutine gdt2gds(igds,igdstmpl,idefnum,ideflist,kgds, &
    igrid,iret)
 
  !
  integer,intent(in) :: idefnum
  integer,intent(in) :: igds(*),igdstmpl(*),ideflist(*)
  integer,intent(out) :: kgds(*),igrid,iret
 
  integer :: kgds72(200),kgds71(200),idum(200),jdum(200)
 
  iret=0
  if (igds(5).eq.0) then       !  Lat/Lon grid
    kgds(1)=0
    kgds(2)=igdstmpl(8)            ! Ni
    kgds(3)=igdstmpl(9)            ! Nj
    kgds(4)=igdstmpl(12)/1000      ! Lat of 1st grid point
    kgds(5)=igdstmpl(13)/1000      ! Long of 1st grid point
    kgds(6)=0                      ! resolution and component flags
    if (igdstmpl(1)==2) kgds(6)=64
    if (btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5))  &
        kgds(6)=kgds(6)+128
    if (btest(igdstmpl(14),3)) kgds(6)=kgds(6)+8
    kgds(7)=igdstmpl(15)/1000      ! Lat of last grid point
    kgds(8)=igdstmpl(16)/1000      ! Long of last grid point
    kgds(9)=igdstmpl(17)/1000      ! Di
    kgds(10)=igdstmpl(18)/1000     ! Dj
    kgds(11)=igdstmpl(19)          ! Scanning mode
    kgds(12)=0
    kgds(13)=0
    kgds(14)=0
    kgds(15)=0
    kgds(16)=0
    kgds(17)=0
    kgds(18)=0
    kgds(19)=0
    kgds(20)=255
    kgds(21)=0
    kgds(22)=0
    !
    !  Process irreg grid stuff, if necessary
    !
    if (idefnum.ne.0) then
      if (igdstmpl(8).eq.-1) then
          kgds(2)=65535
          kgds(9)=65535
      endif
      if (igdstmpl(9).eq.-1) then
          kgds(3)=65535
          kgds(10)=65535
      endif
      kgds(19)=0
      kgds(20)=33
      if (kgds(1).eq.1.OR.kgds(1).eq.3) kgds(20)=43
      kgds(21)=igds(2)                   ! num of grid points
      do j=1,idefnum
          kgds(21+j)=ideflist(j)
      enddo
    endif
  elseif (igds(5).eq.10) then       !  Mercator grid
    kgds(1)=1                 ! Grid Definition Template number
    kgds(2)=igdstmpl(8)            ! Ni
    kgds(3)=igdstmpl(9)            ! Nj
    kgds(4)=igdstmpl(10)/1000      ! Lat of 1st grid point
    kgds(5)=igdstmpl(11)/1000      ! Long of 1st grid point
    kgds(6)=0                      ! resolution and component flags
    if (igdstmpl(1)==2) kgds(6)=64
    if (btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5))  &
        kgds(6)=kgds(6)+128
    if (btest(igdstmpl(12),3)) kgds(6)=kgds(6)+8
    kgds(7)=igdstmpl(14)/1000      ! Lat of last grid point
    kgds(8)=igdstmpl(15)/1000      ! Long of last grid point
    kgds(9)=igdstmpl(13)/1000      ! Lat intersects earth
    kgds(10)=0
    kgds(11)=igdstmpl(16)          ! Scanning mode
    kgds(12)=igdstmpl(18)/1000     ! Di
    kgds(13)=igdstmpl(19)/1000     ! Dj
    kgds(14)=0
    kgds(15)=0
    kgds(16)=0
    kgds(17)=0
    kgds(18)=0
    kgds(19)=0
    kgds(20)=255
    kgds(21)=0
    kgds(22)=0
  elseif (igds(5).eq.30) then       ! Lambert Conformal Grid
    kgds(1)=3
    kgds(2)=igdstmpl(8)            ! Nx
    kgds(3)=igdstmpl(9)            ! Ny
    kgds(4)=igdstmpl(10)/1000      ! Lat of 1st grid point
    kgds(5)=igdstmpl(11)/1000      ! Long of 1st grid point
    kgds(6)=0                      ! resolution and component flags
    if (igdstmpl(1)==2) kgds(6)=64
    if (btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5))  &
        kgds(6)=kgds(6)+128
    if (btest(igdstmpl(12),3)) kgds(6)=kgds(6)+8
    kgds(7)=igdstmpl(14)/1000      ! Lon of orientation
    kgds(8)=igdstmpl(15)/1000      ! Dx
    kgds(9)=igdstmpl(16)/1000      ! Dy
    kgds(10)=igdstmpl(17)          ! Projection Center Flag
    kgds(11)=igdstmpl(18)          ! Scanning mode
    kgds(12)=igdstmpl(19)/1000     ! Lat in 1
    kgds(13)=igdstmpl(20)/1000     ! Lat in 2
    kgds(14)=igdstmpl(21)/1000     ! Lat of S. Pole of projection
    kgds(15)=igdstmpl(22)/1000     ! Lon of S. Pole of projection
    kgds(16)=0
    kgds(17)=0
    kgds(18)=0
    kgds(19)=0
    kgds(20)=255
    kgds(21)=0
    kgds(22)=0
  elseif (igds(5).eq.40) then       !  Gaussian Lat/Lon grid
    kgds(1)=4
    kgds(2)=igdstmpl(8)            ! Ni
    kgds(3)=igdstmpl(9)            ! Nj
    kgds(4)=igdstmpl(12)/1000      ! Lat of 1st grid point
    kgds(5)=igdstmpl(13)/1000      ! Long of 1st grid point
    kgds(6)=0                      ! resolution and component flags
    if (igdstmpl(1)==2) kgds(6)=64
    if (btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5))  &
        kgds(6)=kgds(6)+128
    if (btest(igdstmpl(14),3)) kgds(6)=kgds(6)+8
    kgds(7)=igdstmpl(15)/1000      ! Lat of last grid point
    kgds(8)=igdstmpl(16)/1000      ! Long of last grid point
    kgds(9)=igdstmpl(17)/1000      ! Di
    kgds(10)=igdstmpl(18)          ! N - Number of parallels
    kgds(11)=igdstmpl(19)          ! Scanning mode
    kgds(12)=0
    kgds(13)=0
    kgds(14)=0
    kgds(15)=0
    kgds(16)=0
    kgds(17)=0
    kgds(18)=0
    kgds(19)=0
    kgds(20)=255
    kgds(21)=0
    kgds(22)=0
  elseif (igds(5).eq.20) then       ! Polar Stereographic Grid
    kgds(1)=5
    kgds(2)=igdstmpl(8)            ! Nx
    kgds(3)=igdstmpl(9)            ! Ny
    kgds(4)=igdstmpl(10)/1000      ! Lat of 1st grid point
    kgds(5)=igdstmpl(11)/1000      ! Long of 1st grid point
    kgds(6)=0                      ! resolution and component flags
    if (igdstmpl(1)==2) kgds(6)=64
    if (btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5))  &
        kgds(6)=kgds(6)+128
    if (btest(igdstmpl(12),3)) kgds(6)=kgds(6)+8
    kgds(7)=igdstmpl(14)/1000      ! Lon of orientation
    kgds(8)=igdstmpl(15)/1000      ! Dx
    kgds(9)=igdstmpl(16)/1000      ! Dy
    kgds(10)=igdstmpl(17)          ! Projection Center Flag
    kgds(11)=igdstmpl(18)          ! Scanning mode
    kgds(12)=0
    kgds(13)=0
    kgds(14)=0
    kgds(15)=0
    kgds(16)=0
    kgds(17)=0
    kgds(18)=0
    kgds(19)=0
    kgds(20)=255
    kgds(21)=0
    kgds(22)=0
  elseif (igds(5).eq.204) then      ! Curvilinear Orthogonal
    kgds(1)=204
    kgds(2)=igdstmpl(8)            ! Ni
    kgds(3)=igdstmpl(9)            ! Nj
    kgds(4)=0
    kgds(5)=0
    kgds(6)=0                      ! resolution and component flags
    if (igdstmpl(1)==2) kgds(6)=64
    if (btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5)) &
        kgds(6)=kgds(6)+128
    if (btest(igdstmpl(14),3)) kgds(6)=kgds(6)+8
    kgds(7)=0
    kgds(8)=0
    kgds(9)=0
    kgds(10)=0
    kgds(11)=igdstmpl(19)          ! Scanning mode
    kgds(12)=0
    kgds(13)=0
    kgds(14)=0
    kgds(15)=0
    kgds(16)=0
    kgds(17)=0
    kgds(18)=0
    kgds(19)=0
    kgds(20)=255
    kgds(21)=0
    kgds(22)=0
    !
    !  Process irreg grid stuff, if necessary
    !
    if (idefnum.ne.0) then
      if (igdstmpl(8).eq.-1) then
          kgds(2)=65535
          kgds(9)=65535
      endif
      if (igdstmpl(9).eq.-1) then
          kgds(3)=65535
          kgds(10)=65535
      endif
      kgds(19)=0
      kgds(20)=33
      if (kgds(1).eq.1.OR.kgds(1).eq.3) kgds(20)=43
      kgds(21)=igds(2)                   ! num of grid points
      do j=1,idefnum
          kgds(21+j)=ideflist(j)
      enddo
    endif
  elseif (igds(5).eq.32768) then    ! Rotate Lat/Lon grid
    kgds(1)=203                      ! Arakawa Staggerred E/B grid
    kgds(2)=igdstmpl(8)            ! Ni
    kgds(3)=igdstmpl(9)            ! Nj
    kgds(4)=igdstmpl(12)/1000      ! Lat of 1st grid point
    kgds(5)=igdstmpl(13)/1000      ! Lon of 1st grid point
    kgds(6)=0                      ! resolution and component flags
    if (igdstmpl(1)==2) kgds(6)=64
    if (btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5)) &
        kgds(6)=kgds(6)+128
    if (btest(igdstmpl(14),3)) kgds(6)=kgds(6)+8
    kgds(7)=igdstmpl(15)/1000      ! Lat of last grid point
    kgds(8)=igdstmpl(16)/1000      ! Lon of last grid point
    kgds(9)=igdstmpl(17)/1000      ! Di
    kgds(10)=igdstmpl(18)/1000     ! Dj
    kgds(11)=igdstmpl(19)          ! Scanning mode
    kgds(12)=0
    kgds(13)=0
    kgds(14)=0
    kgds(15)=0
    kgds(16)=0
    kgds(17)=0
    kgds(18)=0
    kgds(19)=0
    kgds(20)=255
    kgds(21)=0
    kgds(22)=0
    !
    !  Process irreg grid stuff, if necessary
    !
    if (idefnum.ne.0) then
      if (igdstmpl(8).eq.-1) then
          kgds(2)=65535
          kgds(9)=65535
      endif
      if (igdstmpl(9).eq.-1) then
          kgds(3)=65535
          kgds(10)=65535
      endif
      kgds(19)=0
      kgds(20)=33
      if (kgds(1).eq.1.OR.kgds(1).eq.3) kgds(20)=43
      kgds(21)=igds(2)                   ! num of grid points
      do j=1,idefnum
          kgds(21+j)=ideflist(j)
      enddo
    endif
  elseif (igds(5).eq.32769) then    ! Rotate Lat/Lon grid
    kgds(1)=205                    ! Arakawa Staggerred for Non-E Stagger grid
    kgds(2)=igdstmpl(8)            ! Ni
    kgds(3)=igdstmpl(9)            ! Nj
    kgds(4)=igdstmpl(12)/1000      ! Lat of 1st grid point
    kgds(5)=igdstmpl(13)/1000      ! Lon of 1st grid point
    kgds(6)=0                      ! resolution and component flags
    if (igdstmpl(1)==2) kgds(6)=64
    if (btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5)) &
        kgds(6)=kgds(6)+128
    if (btest(igdstmpl(14),3)) kgds(6)=kgds(6)+8
    kgds(7)=igdstmpl(15)/1000      ! Lat of last grid point
    kgds(8)=igdstmpl(16)/1000      ! Lon of last grid point
    kgds(9)=igdstmpl(17)/1000      ! Di
    kgds(10)=igdstmpl(18)/1000     ! Dj
    kgds(11)=igdstmpl(19)          ! Scanning mode
    kgds(12)=igdstmpl(20)/1000
    kgds(13)=igdstmpl(21)/1000
    kgds(14)=0
    kgds(15)=0
    kgds(16)=0
    kgds(17)=0
    kgds(18)=0
    kgds(19)=0
    kgds(20)=255
    kgds(21)=0
    kgds(22)=0
  else
    print *,'gdt2gds: Unrecognized GRIB2 GDT = 3.',igds(5)
    iret=1
    kgds(1:22)=0
    return
  endif
  !
  !   Can we determine NCEP grid number ?
  !
  igrid=255
  do j=254,1,-1
    !do j=225,225
    kgds71=0
    kgds72=0
    call w3fi71(j,kgds71,ierr)
    if (ierr.ne.0) cycle
    ! convert W to E for longitudes
    if (kgds71(3).eq.0) then    ! lat/lon
      if (kgds71(7).lt.0) kgds71(7)=360000+kgds71(7)
      if (kgds71(10).lt.0) kgds71(10)=360000+kgds71(10)
    elseif (kgds71(3).eq.1) then    ! mercator
      if (kgds71(7).lt.0) kgds71(7)=360000+kgds71(7)
      if (kgds71(10).lt.0) kgds71(10)=360000+kgds71(10)
    elseif (kgds71(3).eq.3) then     ! lambert conformal
      if (kgds71(7).lt.0) kgds71(7)=360000+kgds71(7)
      if (kgds71(9).lt.0) kgds71(9)=360000+kgds71(9)
      if (kgds71(18).lt.0) kgds71(18)=360000+kgds71(18)
    elseif (kgds71(3).eq.4) then     ! Guassian lat/lon
      if (kgds71(7).lt.0) kgds71(7)=360000+kgds71(7)
      if (kgds71(10).lt.0) kgds71(10)=360000+kgds71(10)
    elseif (kgds71(3).eq.5) then     ! polar stereographic
      if (kgds71(7).lt.0) kgds71(7)=360000+kgds71(7)
      if (kgds71(9).lt.0) kgds71(9)=360000+kgds71(9)
    endif
    call r63w72(idum,kgds,jdum,kgds72)
    if (kgds72(3).eq.3) kgds72(14)=0    ! lambert conformal fix
    if (kgds72(3).eq.1) kgds72(15:18)=0    ! mercator fix
    if (kgds72(3).eq.5) kgds72(14:18)=0    ! polar str fix
    !           print *,' kgds71(',j,')= ', kgds71(1:30)
    !           print *,' kgds72       = ', kgds72(1:30)
    if (all(kgds71.eq.kgds72)) then
      igrid=j
      return
    endif
  enddo
 
  return
end subroutine gdt2gds
 
SUBROUTINE putgbexn(LUGB,KF,KPDS,KGDS,KENS, &
    KPROB,XPROB,KCLUST,KMEMBR,IBS,NBITS,LB,F,IRET)
 
  INTEGER KPDS(200),KGDS(200),KENS(200)
  INTEGER KPROB(2),KCLUST(16),KMEMBR(80)
  REAL XPROB(2)
  LOGICAL*1 LB(KF)
  REAL F(KF)
  !      PARAMETER(MAXBIT=16)
  parameter(maxbit=24)
  INTEGER IBM(KF),IPDS(200),IGDS(200),IBDS(200)
  CHARACTER PDS(400),GRIB(1000+KF*(MAXBIT+1)/8)
 
  !  GET W3FI72 PARAMETERS
  !print *,'SAGT: start putgbexn'
  CALL r63w72(kpds,kgds,ipds,igds)
  ibds=0
 
  !  COUNT VALID DATA
  kbm=kf
  IF(ipds(7).NE.0) THEN
    kbm=0
    DO i=1,kf
      IF(lb(i)) THEN
          ibm(i)=1
          kbm=kbm+1
      ELSE
          ibm(i)=0
      ENDIF
    ENDDO
    IF(kbm.EQ.kf) ipds(7)=0
  ENDIF
 
  !  GET NUMBER OF BITS AND ROUND DATA
  IF(nbits.GT.0) THEN
    nbit=nbits
  ELSE
    IF(kbm.EQ.0) THEN
      DO i=1,kf
          f(i)=0.
      ENDDO
      nbit=0
    ELSE
      !print *,'SAGT:',IPDS(7),IBS,IPDS(25),KF
      !print *,'SAGT:',count(ibm.eq.0),count(ibm.eq.1)
      CALL setbit(ipds(7),-ibs,ipds(25),kf,ibm,f,fmin,fmax,nbit)
      nbit=min(nbit,maxbit)
    ENDIF
  ENDIF
 
  !  CREATE PRODUCT DEFINITION SECTION
  CALL w3fi68(ipds,pds)
  IF(ipds(24).EQ.2) THEN
    ilast=45
    IF (ipds(8).EQ.191.OR.ipds(8).EQ.192) ilast=55
    IF (kens(2).EQ.5) ilast=76
    IF (kens(2).EQ.5) ilast=86
    IF (kens(2).EQ.4) ilast=86
    CALL pdsens(kens,kprob,xprob,kclust,kmembr,ilast,pds)
  ENDIF
 
  !  PACK AND WRITE GRIB DATA
  igflag=1
  igrid=kpds(3)
  if (igrid.ne.255) igflag=0
  !print *,minval(f(1:kf)),maxval(f(1:kf))
  !print *,nbit,kf
  !print *,(ipds(j),j=1,28)
  !write(6,fmt='(28z2)') (pds(j),j=1,28)
  !print *,(kgds(j),j=1,28)
  !print *,(igds(j),j=1,28)
  icomp=0
  CALL w3fi72(0,f,0,nbit,1,ipds,pds, &
      igflag,igrid,igds,icomp,0,ibm,kf,ibds, &
      kfo,grib,lgrib,iret)
  IF(iret.EQ.0) CALL wryte(lugb,lgrib,grib)
 
  RETURN
END SUBROUTINE putgbexn
 
SUBROUTINE setbit(IBM,IBS,IDS,LEN,MG,G,GMIN,GMAX,NBIT)
 
  dimension mg(len),g(len)
 
  !  ROUND FIELD AND DETERMINE EXTREMES WHERE BITMAP IS ON
  s=2.**ibs*10.**ids
  IF(ibm.EQ.0) THEN
     gmax=g(1)
     gmin=g(1)
     DO i=2,len
        gmax=max(gmax,g(i))
        gmin=min(gmin,g(i))
     ENDDO
  ELSE
     i1=1
     DO WHILE(i1.LE.len.AND.mg(i1).EQ.0)
        i1=i1+1
     ENDDO
     IF(i1.LE.len) THEN
        DO i=1,i1-1
           g(i)=0.
        ENDDO
        gmax=g(i1)
        gmin=g(i1)
        DO i=i1+1,len
           IF(mg(i).NE.0) THEN
              gmax=max(gmax,g(i))
              gmin=min(gmin,g(i))
           ELSE
              g(i)=0.
           ENDIF
        ENDDO
     ELSE
        DO i=1,len
           g(i)=0.
        ENDDO
        gmax=0.
        gmin=0.
     ENDIF
  ENDIF
 
  !  COMPUTE NUMBER OF BITS
  nbit=log((gmax-gmin)*s+0.9)/log(2.)+1.
 
  RETURN
END SUBROUTINE setbit