program vint c c ABSTRACT: This program interpolates from various pressure levels c onto regularly-spaced, 50-mb vertical levels. The intent is that c we can use data with relatively coarse vertical resolution to c get data on the necessary 50-mb intervals that we need for Bob c Hart's cyclone phase space. For each model, we will need to read c in a control file that contains the levels that we are c interpolating from. c c Written by Tim Marchok USE params USE grib_mod implicit none type(gribfield) :: holdgfld integer, parameter :: lugb=11,lulv=16,lugi=31,lout=51,maxlev=200 integer kpds(200),kgds(200) integer nlevsin,iriret,iogret,kf,iggret,igdret,iidret,ixo,k,n integer iha,iho,iva,irfa,iodret,ifcsthour,iia,iparm,nlevsout integer gribver,g2_jpdtn integer ilevs(maxlev) real, allocatable :: xinpdat(:,:),xoutdat(:,:),xoutlevs_p(:) logical(1), allocatable :: valid_pt(:),readflag(:) namelist/timein/ifcsthour,iparm,gribver,g2_jpdtn c read (5,NML=timein,END=201) 201 continue print *,' ' print *,'*----------------------------------------------------*' print *,' ' print *,' +++ Top of vint +++' print *,' ' print *,'After namelist read, input forecast hour = ',ifcsthour print *,' input grib parm = ',iparm print *,' GRIB version= ',gribver print *,' GRIB2 JPDTN= g2_jpdtn= ' & ,g2_jpdtn if (iparm == 7 .or. iparm == 156) then nlevsout = 13 ! dealing with height else nlevsout = 5 ! dealing with temperature endif allocate (xoutlevs_p(nlevsout),stat=ixo) if (ixo /= 0) then print *,' ' print *,'!!! ERROR in vint allocating the xoutlevs_p array.' print *,'!!! ixo= ',ixo print *,' ' goto 899 endif do k = 1,nlevsout xoutlevs_p(k) = 300. + float((k-1)*50) enddo ilevs = -999 call read_input_levels (lulv,maxlev,nlevsin,ilevs,iriret) if (iriret /= 0) then print *,' ' print *,'!!! ERROR in vint. ' print *,'!!! RETURN CODE FROM read_input_levels /= 0' print *,'!!! RETURN CODE = iriret = ',iriret print *,'!!! EXITING....' print *,' ' goto 899 endif call open_grib_files (lugb,lugi,lout,gribver,iogret) if (iogret /= 0) then print '(/,a45,i4,/)','!!! ERROR: in vint open_grib_files, rc= ' & ,iogret goto 899 endif call getgridinfo (lugb,lugi,kf,kpds,kgds,holdgfld,ifcsthour,iparm & ,gribver,g2_jpdtn,iggret) allocate (xinpdat(kf,nlevsin),stat=iha) allocate (xoutdat(kf,nlevsout),stat=iho) allocate (valid_pt(kf),stat=iva) allocate (readflag(nlevsin),stat=irfa) if (iha /= 0 .or. iho /= 0 .or. iva /= 0 .or. irfa /= 0) then print *,' ' print *,'!!! ERROR in vint.' print *,'!!! ERROR allocating the xinpdat, readflag, or the' print *,'!!! valid_pt array, iha= ',iha,' iva= ',iva print *,'!!! irfa= ',irfa,' iho= ',iho print *,' ' goto 899 endif print *,'hold check, holdgfld%ipdtlen = ',holdgfld%ipdtlen do n = 1,holdgfld%ipdtlen print *,'hold check, n= ',n,' holdgfld%ipdtmpl= ' & ,holdgfld%ipdtmpl(n) enddo call getdata (lugb,lugi,kf,valid_pt,nlevsin,ilevs,maxlev & ,readflag,xinpdat,ifcsthour,iparm,gribver,g2_jpdtn & ,igdret) call interp_data (kf,valid_pt,nlevsin,ilevs,maxlev,readflag & ,xinpdat,xoutdat,xoutlevs_p,nlevsout,iidret) call output_data (lout,kf,kpds,kgds,holdgfld,xoutdat,valid_pt & ,xoutlevs_p,nlevsout,gribver,iodret) deallocate (xinpdat) deallocate (xoutdat) deallocate (valid_pt) deallocate (readflag) deallocate (xoutlevs_p) 899 continue c stop end c c--------------------------------------------------------------------- c c--------------------------------------------------------------------- subroutine read_input_levels (lulv,maxlev,nlevsin,ilevs,iriret) c c ABSTRACT: This subroutine reads in a text file that contains c the number of input pressure levels for a given model. The c format of the file goes like this, from upper levels to c lower, for example: c c 1 200 c 2 400 c 3 500 c 4 700 c 5 850 c 6 925 c 7 1000 c c implicit none integer lulv,nlevsin,maxlev,iriret,inplev,ict,lvix integer ilevs(maxlev) c iriret=0 ict = 0 do while (.true.) print *,'Top of while loop in vint read_input_levels' read (lulv,85,end=130) lvix,inplev if (inplev > 0 .and. inplev <= 1000) then ict = ict + 1 ilevs(ict) = inplev else print *,' ' print *,'!!! ERROR: Input level not between 0 and 1000' print *,'!!! in vint. inplev= ',inplev print *,'!!! STOPPING EXECUTION' STOP 91 endif print *,'vint readloop, ict= ',ict,' inplev= ',inplev enddo 85 format (i4,1x,i4) 130 continue nlevsin = ict print *,' ' print *,'Total number of vint levels read in = ',nlevsin c return end c--------------------------------------------------------------------- c c--------------------------------------------------------------------- subroutine getgridinfo (lugb,lugi,kf,kpds,kgds,holdgfld,ifcsthour & ,iparm,gribver,g2_jpdtn,iggret) c c ABSTRACT: The purpose of this subroutine is just to get the max c values of i and j and the dx and dy grid spacing intervals for the c grid to be used in the rest of the program. So just read the c grib file to get the lon and lat data. Also, get the info for c the data grid's boundaries. This boundary information will be c used later in the tracking algorithm, and is accessed via Module c grid_bounds. c C INPUT: C lugb The Fortran unit number for the GRIB data file C lugi The Fortran unit number for the GRIB index file c ifcsthour input forecast hour to search for c iparm input grib parm to search for c gribver integer (1 or 2) to indicate if using GRIB1 / GRIB2 c g2_jpdtn If GRIB2 data being read, this is the value for JPDTN c that is input to getgb2. C C OUTPUT: c kf Number of gridpoints on the grid c kpds pds array for a GRIB1 record c kgds gds array for a GRIB1 record c holdgfld info for a GRIB2 record c C iggret The return code from this subroutine c USE params USE grib_mod implicit none c type(gribfield) :: gfld,prevfld,holdgfld integer,dimension(200) :: jids,jpdt,jgdt logical(1), allocatable :: lb(:) integer, parameter :: jf=4000000 integer jpds(200),jgds(200) integer kpds(200),kgds(200) integer :: listsec1(13) integer ila,ifa,iret,ifcsthour,imax,jmax,jskp,jdisc integer lugb,lugi,kf,j,k,iggret,iparm,gribver,g2_jpdtn integer jpdtn,jgdtn,npoints,icount,ipack,krec integer :: listsec0(2)=(/0,2/) integer :: igds(5)=(/0,0,0,0,0/),previgds(5) integer :: idrstmpl(200) integer :: currlen=1000000 logical :: unpack=.true. logical :: open_grb=.false. real, allocatable :: f(:) real dx,dy c iggret = 0 allocate (lb(jf),stat=ila) allocate (f(jf),stat=ifa) if (ila /= 0 .or. ifa /= 0) then print *,' ' print *,'!!! ERROR in vint.' print *,'!!! ERROR in getgridinfo allocating either lb or f' print *,'!!! ila = ',ila,' ifa= ',ifa iggret = 97 return endif if (gribver == 2) then ! Search for a record from a GRIB2 file ! ! --- Initialize Variables --- ! gfld%idsect => NULL() gfld%local => NULL() gfld%list_opt => NULL() gfld%igdtmpl => NULL() gfld%ipdtmpl => NULL() gfld%coord_list => NULL() gfld%idrtmpl => NULL() gfld%bmap => NULL() gfld%fld => NULL() jdisc=0 ! meteorological products jids=-9999 jpdtn=g2_jpdtn ! 0 = analysis or forecast; 1 = ens fcst jgdtn=0 ! lat/lon grid jgdt=-9999 jpdt=-9999 npoints=0 icount=0 jskp=0 c Search for Temperature or GP Height by production template.... JPDT(1:15)=(/-9999,-9999,-9999,-9999,-9999,-9999,-9999,-9999 & ,-9999,-9999,-9999,-9999,-9999,-9999,-9999/) if (iparm == 7) then ! GP Height jpdt(1) = 3 ! Param category from Table 4.1 jpdt(2) = 5 ! Param number from Table 4.2-0-3 elseif (iparm == 11) then ! Temperature jpdt(1) = 0 ! Param category from Table 4.1 jpdt(2) = 0 ! Param category from Table 4.2 endif jpdt(9) = ifcsthour call getgb2(lugb,lugi,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt & ,unpack,krec,gfld,iret) if ( iret.ne.0) then print *,' ' print *,' ERROR: getgb2 error in getgridinfo = ',iret endif c Determine packing information from GRIB2 file c The default packing is 40 JPEG 2000 ipack = 40 print *,' gfld%idrtnum = ', gfld%idrtnum ! Set DRT info ( packing info ) if ( gfld%idrtnum.eq.0 ) then ! Simple packing ipack = 0 elseif ( gfld%idrtnum.eq.2 ) then ! Complex packing ipack = 2 elseif ( gfld%idrtnum.eq.3 ) then ! Complex & spatial packing ipack = 31 elseif ( gfld%idrtnum.eq.40.or.gfld%idrtnum.eq.15 ) then ! JPEG 2000 packing ipack = 40 elseif ( gfld%idrtnum.eq.41 ) then ! PNG packing ipack = 41 endif print *,'After check of idrtnum, ipack= ',ipack print *,'Number of gridpts= gfld%ngrdpts= ',gfld%ngrdpts print *,'Number of elements= gfld%igdtlen= ',gfld%igdtlen print *,'PDT num= gfld%ipdtnum= ',gfld%ipdtnum print *,'GDT num= gfld%igdtnum= ',gfld%igdtnum imax = gfld%igdtmpl(8) print *,'at A' jmax = gfld%igdtmpl(9) print *,'at B' dx = float(gfld%igdtmpl(17))/1.e6 print *,'at C' dy = float(gfld%igdtmpl(17))/1.e6 print *,'at D' kf = gfld%ngrdpts print *,'at E' holdgfld = gfld else ! Search for a record from a GRIB1 file jpds = -1 jgds = -1 j=0 jpds(5) = iparm ! Get a record for the input parm selected jpds(6) = 100 ! Get a record on a standard pressure level jpds(14) = ifcsthour call getgb(lugb,lugi,jf,j,jpds,jgds, & kf,k,kpds,kgds,lb,f,iret) if (iret.ne.0) then print *,' ' print *,'!!! ERROR in vint getgridinfo calling getgb' print *,'!!! Return code from getgb = iret = ',iret iggret = iret return else iggret=0 imax = kgds(2) jmax = kgds(3) dx = float(kgds(9))/1000. dy = float(kgds(10))/1000. endif endif print *,' ' print *,'In vint getgridinfo, grid dimensions follow:' print *,'imax= ',imax,' jmax= ',jmax print *,' dx= ',dx,' dy= ',dy print *,'number of gridpoints = ',kf deallocate (lb); deallocate(f) return end c--------------------------------------------------------------------- c c--------------------------------------------------------------------- subroutine getdata (lugb,lugi,kf,valid_pt,nlevsin,ilevs,maxlev & ,readflag,xinpdat,ifcsthour,iparm,gribver,g2_jpdtn & ,igdret) c c ABSTRACT: This subroutine reads the input GRIB file for the c tracked parameters. USE params USE grib_mod implicit none c type(gribfield) :: gfld,prevfld CHARACTER(len=8) :: pabbrev integer,dimension(200) :: jids,jpdt,jgdt logical(1) valid_pt(kf),lb(kf),readflag(nlevsin) integer, parameter :: jf=4000000 integer ilevs(maxlev) integer jpds(200),jgds(200),kpds(200),kgds(200) integer lugb,lugi,kf,nlevsin,maxlev,igdret,jskp,jdisc integer i,j,k,ict,np,lev,ifcsthour,iret,iparm,gribver,g2_jpdtn integer jpdtn,jgdtn,npoints,icount,ipack,krec integer pdt_4p0_vert_level,pdt_4p0_vtime,mm integer :: listsec0(2)=(/0,2/) integer :: listsec1(13) integer :: igds(5)=(/0,0,0,0,0/),previgds(5) integer :: idrstmpl(200) integer :: currlen=1000000 logical :: unpack=.true. logical :: open_grb=.false. real f(kf),xinpdat(kf,nlevsin),xtemp(kf) real dmin,dmax,firstval,lastval c igdret=0 ict = 0 level_loop: do lev = 1,nlevsin print *,' ' print *,'In vint getdata read loop, lev= ',lev,' level= ' & ,ilevs(lev) if (gribver == 2) then ! ! --- Initialize Variables --- ! gfld%idsect => NULL() gfld%local => NULL() gfld%list_opt => NULL() gfld%igdtmpl => NULL() gfld%ipdtmpl => NULL() gfld%coord_list => NULL() gfld%idrtmpl => NULL() gfld%bmap => NULL() gfld%fld => NULL() jdisc=0 ! meteorological products jids=-9999 jpdtn=g2_jpdtn ! 0 = analysis or forecast; 1 = ens fcst jgdtn=0 ! lat/lon grid jgdt=-9999 jpdt=-9999 npoints=0 icount=0 jskp=0 c Search for input parameter by production template 4.0. This c vint program is used primarily for temperature, but still we c will leave that as a variable and not-hard wire it in case we c choose to average something else in the future. ! We are looking for Temperature or GP Height here. This ! block of code, or even the smaller subset block of code that ! contains the JPDT(1) and JPDT(2) assignments, can of course ! be modified if this program is to be used for interpolating ! other variables.... ! Set defaults for JPDT, then override in array ! assignments below... JPDT(1:15)=(/-9999,-9999,-9999,-9999,-9999,-9999,-9999,-9999 & ,-9999,-9999,-9999,-9999,-9999,-9999,-9999/) print *,' ' print *,'In getdata vint, iparm= ',iparm if (iparm == 7) then ! GP Height jpdt(1) = 3 ! Param category from Table 4.1 jpdt(2) = 5 ! Param number from Table 4.2-0-3 elseif (iparm == 11) then ! Temperature jpdt(1) = 0 ! Param category from Table 4.1 jpdt(2) = 0 ! Param category from Table 4.2 endif JPDT(9) = ifcsthour JPDT(10) = 100 ! Isobaric surface requested (Table 4.5) JPDT(12) = ilevs(lev) * 100 ! value of specific level print *,'before getgb2 call, value of unpack = ',unpack do mm = 1,15 print *,'VINT getdata mm= ',mm,' JPDT(mm)= ',JPDT(mm) enddo call getgb2(lugb,lugi,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt & ,unpack,krec,gfld,iret) print *,'iret from getgb2 in getdata = ',iret print *,'after getgb2 call, value of unpacked = ' & ,gfld%unpacked print *,'after getgb2 call, gfld%ndpts = ',gfld%ndpts print *,'after getgb2 call, gfld%ibmap = ',gfld%ibmap if ( iret == 0) then c Determine packing information from GRIB2 file c The default packing is 40 JPEG 2000 ipack = 40 print *,' gfld%idrtnum = ', gfld%idrtnum ! Set DRT info ( packing info ) if ( gfld%idrtnum.eq.0 ) then ! Simple packing ipack = 0 elseif ( gfld%idrtnum.eq.2 ) then ! Complex packing ipack = 2 elseif ( gfld%idrtnum.eq.3 ) then ! Complex & spatial & ! packing ipack = 31 elseif ( gfld%idrtnum.eq.40.or.gfld%idrtnum.eq.15 ) then ! JPEG 2000 packing ipack = 40 elseif ( gfld%idrtnum.eq.41 ) then ! PNG packing ipack = 41 endif print *,'After check of idrtnum, ipack= ',ipack print *,'Number of gridpts= gfld%ngrdpts= ',gfld%ngrdpts print *,'Number of elements= gfld%igdtlen= ',gfld%igdtlen print *,'GDT num= gfld%igdtnum= ',gfld%igdtnum kf = gfld%ndpts ! Number of gridpoints returned from read do np = 1,kf xinpdat(np,lev) = gfld%fld(np) xtemp(np) = gfld%fld(np) if (gfld%ibmap == 0) then valid_pt(np) = gfld%bmap(np) else valid_pt(np) = .true. endif enddo readflag(lev) = .TRUE. c call bitmapchk(kf,gfld%bmap,gfld%fld,dmin,dmax) call bitmapchk(kf,valid_pt,xtemp,dmin,dmax) if (ict == 0) then c do np = 1,kf c valid_pt(np) = gfld%bmap(np) c enddo ict = ict + 1 endif firstval=gfld%fld(1) lastval=gfld%fld(kf) print *,' ' print *,' SECTION 0: discipl= ',gfld%discipline & ,' gribver= ',gfld%version print *,' ' print *,' SECTION 1: ' do j = 1,gfld%idsectlen print *,' sect1, j= ',j,' gfld%idsect(j)= ' & ,gfld%idsect(j) enddo if ( associated(gfld%local).AND.gfld%locallen.gt.0) then print *,' ' print *,' SECTION 2: ',gfld%locallen,' bytes' else print *,' ' print *,' SECTION 2 DOES NOT EXIST IN THIS RECORD' endif print *,' ' print *,' SECTION 3: griddef= ',gfld%griddef print *,' ngrdpts= ',gfld%ngrdpts print *,' numoct_opt= ',gfld%numoct_opt print *,' interp_opt= ',gfld%interp_opt print *,' igdtnum= ',gfld%igdtnum print *,' igdtlen= ',gfld%igdtlen print *,' ' print '(a17,i3,a2)',' GRID TEMPLATE 3.',gfld%igdtnum,': ' do j=1,gfld%igdtlen print *,' j= ',j,' gfld%igdtmpl(j)= ',gfld%igdtmpl(j) enddo print *,' ' print *,' PDT num (gfld%ipdtnum) = ',gfld%ipdtnum print *,' ' print '(a20,i3,a2)',' PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ' do j=1,gfld%ipdtlen print *,' sect 4 j= ',j,' gfld%ipdtmpl(j)= ' & ,gfld%ipdtmpl(j) enddo c Print out values for data representation type print *,' ' print '(a21,i3,a2)',' DATA REP TEMPLATE 5.',gfld%idrtnum & ,': ' do j=1,gfld%idrtlen print *,' sect 5 j= ',j,' gfld%idrtmpl(j)= ' & ,gfld%idrtmpl(j) enddo c Get parameter abbrev for record that was retrieved pdt_4p0_vtime = gfld%ipdtmpl(9) pdt_4p0_vert_level = gfld%ipdtmpl(12) pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1) & ,gfld%ipdtmpl(2)) print *,' ' write (6,131) 131 format (' rec# param level byy bmm bdd bhh ' & ,'fhr npts firstval lastval minval ' & ,' maxval') print '(i5,3x,a8,2x,6i5,2x,i8,4g12.4)' & ,krec,pabbrev,pdt_4p0_vert_level/100,gfld%idsect(6) & ,gfld%idsect(7),gfld%idsect(8),gfld%idsect(9) & ,pdt_4p0_vtime,gfld%ndpts,firstval,lastval,dmin,dmax do np = 1,kf xinpdat(np,lev) = gfld%fld(np) enddo else print *,' ' print *,'!!! ERROR: GRIB2 VINT READ IN GETDATA FAILED FOR ' & ,'LEVEL LEV= ',LEV print *,' ' readflag(lev) = .FALSE. do np = 1,kf xinpdat(np,lev) = -99999.0 enddo endif else ! Reading a GRIB1 file.... jpds = -1 jgds = -1 j=0 jpds(5) = iparm ! grib parameter id to read in jpds(6) = 100 ! level id to indicate a pressure level jpds(7) = ilevs(lev) ! actual level of the layer jpds(14) = ifcsthour ! lead time to search for call getgb (lugb,lugi,jf,j,jpds,jgds, & kf,k,kpds,kgds,lb,f,iret) print *,' ' print *,'After vint getgb call, j= ',j,' k= ',k,' level= ' & ,ilevs(lev),' iret= ',iret if (iret == 0) then readflag(lev) = .TRUE. call bitmapchk(kf,lb,f,dmin,dmax) if (ict == 0) then do np = 1,kf valid_pt(np) = lb(np) enddo ict = ict + 1 endif write (6,31) 31 format (' rec# parm# levt lev byy bmm bdd bhh fhr ' & ,'npts minval maxval') print '(i4,2x,8i5,i8,2g12.4)', & k,(kpds(i),i=5,11),kpds(14),kf,dmin,dmax do np = 1,kf xinpdat(np,lev) = f(np) enddo else print *,' ' print *,'!!! ERROR: VINT READ FAILED FOR LEVEL LEV= ',LEV print *,' ' readflag(lev) = .FALSE. do np = 1,kf xinpdat(np,lev) = -99999.0 enddo endif endif enddo level_loop c return end c c----------------------------------------------------------------------- c c----------------------------------------------------------------------- subroutine interp_data (kf,valid_pt,nlevsin,ilevs,maxlev,readflag & ,xinpdat,xoutdat,xoutlevs_p,nlevsout,iidret) c c ABSTRACT: This routine interpolates data in between available c pressure levels to get data resolution at the 50-mb c resolution that we need for the cyclone phase space c diagnostics. implicit none logical(1) valid_pt(kf),readflag(nlevsin) integer ilevs(maxlev) integer nlevsin,nlevsout,maxlev,kf,kout,kin,k,n,kup,klo integer iidret real xinpdat(kf,nlevsin),xoutdat(kf,nlevsout) real xoutlevs_p(nlevsout),xoutlevs_lnp(nlevsout) real xinlevs_p(nlevsin),xinlevs_lnp(nlevsin) real pdiff,pdiffmin,xu,xo,xl,yu,yl c iidret=0 print *,' ' print *,'*----------------------------------------------*' print *,' Listing of standard output levels follows....' print *,'*----------------------------------------------*' print *,' ' do k = 1,nlevsout xoutlevs_lnp(k) = log(xoutlevs_p(k)) write (6,81) k,xoutlevs_p(k),xoutlevs_lnp(k) enddo 81 format (1x,'k= ',i3,' p= ',f6.1,' ln(p)= ',f9.6) do k = 1,nlevsin xinlevs_p(k) = float(ilevs(k)) xinlevs_lnp(k) = log(xinlevs_p(k)) enddo c ----------------------------------------------------------------- c We want to loop through for all the *output* levels that we need. c We may have some input levels that match perfectly, often at c least the standard levels like 500, 700, 850. For these levels, c just take the data directly from the input file. For other c output levels that fall between the input levels, we need to c find the nearest upper and lower levels. output_loop: do kout = 1,nlevsout print *,' ' print *,'+------------------------------------------------+' print *,'Top of vint output_loop, kout= ',kout,' pressure= ' & ,xoutlevs_p(kout) ! Loop through all of the input levels and find the level ! that is closest to the output level from the *upper* side. ! And again, in this upper loop, if we hit a level that ! exactly matches a needed output level, just copy that data ! and then cycle back to the top of output_loop. kup = -999 klo = -999 pdiffmin = 9999.0 inp_loop_up: do kin = 1,nlevsin if (xinlevs_p(kin) == xoutlevs_p(kout)) then print *,' ' print *,'+++ Exact level found. kout= ',kout print *,'+++ level= ',xoutlevs_p(kout) print *,'+++ Data copied. No interpolation needed.' if (readflag(kin)) then do n = 1,kf xoutdat(n,kout) = xinpdat(n,kin) enddo cycle output_loop else print *,' ' print *,'!!! ERROR: readflag is FALSE in interp_data for' print *,'!!! level kin= ',kin,', which is a level that ' print *,'!!! exactly matches a required output level, and' print *,'!!! the user has identified as being an input ' print *,'!!! level with valid data for this model. We ' print *,'!!! will get the data from a different level.' endif else pdiff = xoutlevs_p(kout) - xinlevs_p(kin) if (pdiff > 0.) then ! We have a level higher than outlev if (pdiff < pdiffmin) then pdiffmin = pdiff kup = kin endif endif endif enddo inp_loop_up pdiffmin = 9999.0 inp_loop_lo: do kin = 1,nlevsin pdiff = xinlevs_p(kin) - xoutlevs_p(kout) if (pdiff > 0.) then ! We have a level lower than outlev if (pdiff < pdiffmin) then pdiffmin = pdiff klo = kin endif endif enddo inp_loop_lo if (kup == -999 .or. klo == -999) then print *,' ' print *,'!!! ERROR: While interpolating, could not find ' print *,'!!! either an upper or lower input level to use' print *,'!!! for interpolating *from*.' print *,'!!! kup= ',kup,' klo= ',klo print *,' ' print *,'!!! STOPPING....' stop 91 endif if (.not. readflag(kup) .or. .not. readflag(klo)) then print *,' ' print *,'!!! ERROR: In interp_data, either the upper or the' print *,'!!! lower input level closest to the target output' print *,'!!! level did not have valid data read in.' print *,'!!! ' write (6,91) ' upper level k= ',kup,xinlevs_p(kup) & ,xinlevs_lnp(kup) write (6,101) xoutlevs_p(kout),xoutlevs_lnp(kout) write (6,91) ' lower level k= ',klo,xinlevs_p(klo) & ,xinlevs_lnp(klo) print *,'!!! readflag upper = ',readflag(kup) print *,'!!! readflag lower = ',readflag(klo) print *,'!!! EXITING....' stop 92 endif print *,' ' write (6,91) ' upper level k= ',kup,xinlevs_p(kup) & ,xinlevs_lnp(kup) write (6,101) xoutlevs_p(kout),xoutlevs_lnp(kout) write (6,91) ' lower level k= ',klo,xinlevs_p(klo) & ,xinlevs_lnp(klo) 91 format (1x,a17,1x,i3,' pressure= ',f6.1,' ln(p)= ',f9.6) 101 format (13x,'Target output pressure= ',f6.1,' ln(p)= ',f9.6) !-------------------------------------------------------------- ! Now perform the linear interpolation. Here is the notation ! used in the interpolation: ! ! xu = ln of pressure at upper level ! xo = ln of pressure at output level ! xl = ln of pressure at lower level ! yu = data value at upper level ! yl = data value at lower level !-------------------------------------------------------------- xu = xinlevs_lnp(kup) xo = xoutlevs_lnp(kout) xl = xinlevs_lnp(klo) do n = 1,kf yu = xinpdat(n,kup) yl = xinpdat(n,klo) xoutdat(n,kout) = ((yl * (xo - xu)) - (yu * (xo - xl))) & / (xl - xu) enddo enddo output_loop c return end c c---------------------------------------------------------------------- c c---------------------------------------------------------------------- subroutine output_data (lout,kf,kpds,kgds,holdgfld,xoutdat & ,valid_pt,xoutlevs_p,nlevsout,gribver,iodret) c c ABSTRACT: This routine writes out the output data on the c specified output pressure levels. USE params USE grib_mod implicit none CHARACTER(len=1),pointer,dimension(:) :: cgrib type(gribfield) :: holdgfld logical(1) valid_pt(kf),bmap(kf) integer lout,kf,lugb,lugi,iodret,nlevsout,igoret,ipret,lev integer gribver,ierr,ipack,lengrib,npoints,newlen,idrsnum integer numcoord,ica,n,j integer :: idrstmpl(200) integer :: currlen=1000000 integer :: listsec0(2)=(/0,2/) integer :: igds(5)=(/0,0,0,0,0/),previgds(5) integer kpds(200),kgds(200) integer(4), parameter::idefnum=1 integer(4) ideflist(idefnum),ibmap real coordlist real xoutdat(kf,nlevsout),xoutlevs_p(nlevsout) c iodret=0 call baopenw (lout,"fort.51",igoret) print *,'baopenw: igoret= ',igoret if (igoret /= 0) then print *,' ' print *,'!!! ERROR in vint in sub output_data opening' print *,'!!! **OUTPUT** grib file. baopenw return codes:' print *,'!!! grib file 1 return code = igoret = ',igoret STOP 95 return endif levloop: do lev = 1,nlevsout if (gribver == 2) then ! Write data out as a GRIB2 message.... allocate(cgrib(currlen),stat=ica) if (ica /= 0) then print *,' ' print *,'ERROR in output_data allocating cgrib' print *,'ica= ',ica iodret=95 return endif ! Ensure that cgrib array is large enough if (holdgfld%ifldnum == 1 ) then ! start new GRIB2 message npoints=holdgfld%ngrdpts else npoints=npoints+holdgfld%ngrdpts endif newlen=npoints*4 if ( newlen.gt.currlen ) then ccc if (allocated(cgrib)) deallocate(cgrib) if (associated(cgrib)) deallocate(cgrib) allocate(cgrib(newlen),stat=ierr) c call realloc (cgrib,currlen,newlen,ierr) if (ierr == 0) then print *,' ' print *,'re-allocate for large grib msg: ' print *,' currlen= ',currlen print *,' newlen= ',newlen currlen=newlen else print *,'ERROR returned from 2nd allocate cgrib = ',ierr stop 95 endif endif ! Create new GRIB Message listsec0(1)=holdgfld%discipline listsec0(2)=holdgfld%version print *,'output, holdgfld%idsectlen= ',holdgfld%idsectlen do j = 1,holdgfld%idsectlen print *,' sect1, j= ',j,' holdgfld%idsect(j)= ' & ,holdgfld%idsect(j) enddo call gribcreate(cgrib,currlen,listsec0,holdgfld%idsect,ierr) if (ierr.ne.0) then write(6,*) ' ERROR creating new GRIB2 field (gribcreate)= ' & ,ierr stop 95 endif previgds=igds igds(1)=holdgfld%griddef igds(2)=holdgfld%ngrdpts igds(3)=holdgfld%numoct_opt igds(4)=holdgfld%interp_opt igds(5)=holdgfld%igdtnum if (igds(3) == 0) then ideflist = 0 endif call addgrid (cgrib,currlen,igds,holdgfld%igdtmpl & ,holdgfld%igdtlen,ideflist,idefnum,ierr) if (ierr.ne.0) then write(6,*) ' ERROR from addgrid adding GRIB2 grid = ',ierr stop 95 endif holdgfld%ipdtmpl(12) = int(xoutlevs_p(lev)) * 100 ipack = 40 idrsnum = ipack idrstmpl = 0 idrstmpl(2)= holdgfld%idrtmpl(2) idrstmpl(3)= holdgfld%idrtmpl(3) idrstmpl(6)= 0 idrstmpl(7)= 255 numcoord=0 coordlist=0.0 ! Only needed for hybrid vertical coordinate, ! not here, so set it to 0.0 ! 0 - A bit map applies to this product and is specified in ! this section ! 255 - A bit map does not apply to this product ibmap=255 ! Bitmap indicator (see Code Table 6.0) print *,' ' print *,'output, holdgfld%ipdtlen= ',holdgfld%ipdtlen do n = 1,holdgfld%ipdtlen print *,'output, n= ',n,' holdgfld%ipdtmpl= ' & ,holdgfld%ipdtmpl(n) enddo print *,'output, kf= ',kf c do n = 1,kf c print *,'output, n= ',n,' xoutdat(n)= ',xoutdat(n) c enddo call addfield (cgrib,currlen,holdgfld%ipdtnum,holdgfld%ipdtmpl & ,holdgfld%ipdtlen,coordlist & ,numcoord & ,idrsnum,idrstmpl,200 & ,xoutdat(1,lev),kf,ibmap,bmap,ierr) if (ierr /= 0) then write(6,*) ' ERROR from addfield adding GRIB2 data = ',ierr stop 95 endif ! Finalize GRIB message after all grids ! and fields have been added. It adds the End Section ( "7777" ) call gribend(cgrib,currlen,lengrib,ierr) call wryte(lout,lengrib,cgrib) if (ierr == 0) then print *,' ' print *,'+++ GRIB2 write successful. ' print *,' Len of message = currlen= ',currlen print *,' Len of entire GRIB2 message = lengrib= ' & ,lengrib else print *,' ERROR from gribend writing GRIB2 msg = ',ierr stop 95 endif else ! Write data out as a GRIB1 message.... kpds(7) = int(xoutlevs_p(lev)) print *,'In vint, just before call to putgb, kf= ',kf call putgb (lout,kf,kpds,kgds,valid_pt,xoutdat(1,lev),ipret) print *,'In vint, just after call to putgb, kf= ',kf if (ipret == 0) then print *,' ' print *,'+++ IPRET = 0 after call to putgb in vint' print *,' ' else print *,' ' print *,'!!!!!! ERROR in vint.' print *,'!!!!!! ERROR: IPRET NE 0 AFTER CALL TO PUTGB !!!' print *,'!!!!!! Level index= ',lev print *,'!!!!!! pressure= ',xoutlevs_p(lev) print *,' ' endif write(*,980) kpds(1),kpds(2) write(*,981) kpds(3),kpds(4) write(*,982) kpds(5),kpds(6) write(*,983) kpds(7),kpds(8) write(*,984) kpds(9),kpds(10) write(*,985) kpds(11),kpds(12) write(*,986) kpds(13),kpds(14) write(*,987) kpds(15),kpds(16) write(*,988) kpds(17),kpds(18) write(*,989) kpds(19),kpds(20) write(*,990) kpds(21),kpds(22) write(*,991) kpds(23),kpds(24) write(*,992) kpds(25) write(*,880) kgds(1),kgds(2) write(*,881) kgds(3),kgds(4) write(*,882) kgds(5),kgds(6) write(*,883) kgds(7),kgds(8) write(*,884) kgds(9),kgds(10) write(*,885) kgds(11),kgds(12) write(*,886) kgds(13),kgds(14) write(*,887) kgds(15),kgds(16) write(*,888) kgds(17),kgds(18) write(*,889) kgds(19),kgds(20) write(*,890) kgds(21),kgds(22) 980 format(' kpds(1) = ',i7,' kpds(2) = ',i7) 981 format(' kpds(3) = ',i7,' kpds(4) = ',i7) 982 format(' kpds(5) = ',i7,' kpds(6) = ',i7) 983 format(' kpds(7) = ',i7,' kpds(8) = ',i7) 984 format(' kpds(9) = ',i7,' kpds(10) = ',i7) 985 format(' kpds(11) = ',i7,' kpds(12) = ',i7) 986 format(' kpds(13) = ',i7,' kpds(14) = ',i7) 987 format(' kpds(15) = ',i7,' kpds(16) = ',i7) 988 format(' kpds(17) = ',i7,' kpds(18) = ',i7) 989 format(' kpds(19) = ',i7,' kpds(20) = ',i7) 990 format(' kpds(21) = ',i7,' kpds(22) = ',i7) 991 format(' kpds(23) = ',i7,' kpds(24) = ',i7) 992 format(' kpds(25) = ',i7) 880 format(' kgds(1) = ',i7,' kgds(2) = ',i7) 881 format(' kgds(3) = ',i7,' kgds(4) = ',i7) 882 format(' kgds(5) = ',i7,' kgds(6) = ',i7) 883 format(' kgds(7) = ',i7,' kgds(8) = ',i7) 884 format(' kgds(9) = ',i7,' kgds(10) = ',i7) 885 format(' kgds(11) = ',i7,' kgds(12) = ',i7) 886 format(' kgds(13) = ',i7,' kgds(14) = ',i7) 887 format(' kgds(15) = ',i7,' kgds(16) = ',i7) 888 format(' kgds(17) = ',i7,' kgds(18) = ',i7) 889 format(' kgds(19) = ',i7,' kgds(20) = ',i7) 890 format(' kgds(20) = ',i7,' kgds(22) = ',i7) endif enddo levloop c return end c c----------------------------------------------------------------------- c c----------------------------------------------------------------------- subroutine open_grib_files (lugb,lugi,lout,gribver,iret) C ABSTRACT: This subroutine must be called before any attempt is C made to read from the input GRIB files. The GRIB and index files C are opened with a call to baopenr. This call to baopenr was not C needed in the cray version of this program (the files could be C opened with a simple Cray assign statement), but the GRIB-reading C utilities on the SP do require calls to this subroutine (it has C something to do with the GRIB I/O being done in C on the SP, and C the C I/O package needs an explicit open statement). C C INPUT: C lugb The Fortran unit number for the GRIB data file C lugi The Fortran unit number for the GRIB index file C lout The Fortran unit number for the output grib file c gribver integer (1 or 2) to indicate if using GRIB1 / GRIB2 C C OUTPUT: C iret The return code from this subroutine implicit none character fnameg*7,fnamei*7,fnameo*7 integer iret,gribver,lugb,lugi,lout,igoret,iioret,iooret iret=0 fnameg(1:5) = "fort." fnamei(1:5) = "fort." fnameo(1:5) = "fort." write(fnameg(6:7),'(I2)') lugb write(fnamei(6:7),'(I2)') lugi write(fnameo(6:7),'(I2)') lout call baopenr (lugb,fnameg,igoret) call baopenr (lugi,fnamei,iioret) call baopenw (lout,fnameo,iooret) print *,' ' print *,'vint: baopen: igoret= ',igoret,' iioret= ',iioret & ,' iooret= ',iooret if (igoret /= 0 .or. iioret /= 0 .or. iooret /= 0) then print *,' ' print *,'!!! ERROR in vint.' print *,'!!! ERROR in sub open_grib_files opening grib file' print *,'!!! or grib index file. baopen return codes:' print *,'!!! grib file return code = igoret = ',igoret print *,'!!! index file return code = iioret = ',iioret print *,'!!! output file return code = iooret = ',iooret iret = 93 return endif return end c c------------------------------------------------------------------- c c------------------------------------------------------------------- subroutine bitmapchk (n,ld,d,dmin,dmax) c c This subroutine checks the bitmap for non-existent data values. c Since the data from the regional models have been interpolated c from either a polar stereographic or lambert conformal grid c onto a lat/lon grid, there will be some gridpoints around the c edges of this lat/lon grid that have no data; these grid c points have been bitmapped out by Mark Iredell's interpolater. c To provide another means of checking for invalid data points c later in the program, set these bitmapped data values to a c value of -999.0. The min and max of this array are also c returned if a user wants to check for reasonable values. c logical(1) ld dimension ld(n),d(n) c dmin=1.E15 dmax=-1.E15 c do i=1,n if (ld(i)) then dmin=min(dmin,d(i)) dmax=max(dmax,d(i)) else d(i) = -999.0 endif enddo c return end