module blending ! ABSTRACT: This program reads GRIB2 file and makes inventory ! of GRIB2 file ! ! PROGRAM HISTORY LOG: ! 2020-01-21 Y Mao ! ! SUBPROGRAMS CALLED: (LIST ALL CALLED FROM ANYWHERE IN CODES) ! LIBRARY: ! G2LIB - GB_INFO, GT_GETFLD, PRLEVEL, PRVTIME ! W3LIB - GBYTE, SKGB ! BACIO - BAOPENR, BAREAD, BACLOSE ! SYSTEM - IARGC FUNCTION RETURNS NUMBER OF ARGUMENT ON ! COMMAND LINE ! - GETARG ROUTINE RETURNS COMMAND LINE ARGUMENT ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 ! ! No bit-map will be used for blended output ! ! Turbulence blending: ! Max of UK US ! ! Icing severity blending: ! 0. Select pressures ! 1. Max of UK US ! 2. Blend by Gaussian Kernel Filter (sigma=1) ! 3. Keep the original matching data ! 4. Re-categorize by different thresholds ! 5. Not greater than max, and not smaller than min of US and UK ! ! CB blending: ! 1. extent: average ! 2. base: min ! 3. top: max contains subroutine process(gfile1,gfile2,gfile3) use grib_mod use params implicit none character(*), intent(in) :: gfile1,gfile2,gfile3 INTEGER :: NARG integer :: ifl1,ifl2,ifl3 integer, parameter :: msk1=32000 real, parameter :: EPSILON=0.000001 integer :: currlen=0,icount,itot integer :: iseek,lskip,lgrib,lengrib CHARACTER(len=1),allocatable,dimension(:) :: cgrib integer :: listsec0(3),listsec1(13) integer :: numfields,numlocal,maxlocal integer :: ierr,n,k,im,jm,i,j,jj,ij logical :: unpack,expand character(len=8) :: pabbrev integer :: jids(20),jpdt(20),jgdt(20) integer :: jpdtn,jgdtn type(gribfield) :: gfld,gfld2 real, allocatable :: usdata(:,:),ukdata(:,:),blddata(:) real :: missing character(3) :: whichblnd logical :: exclusive real :: fldmax,fldmin,sum unpack=.true. expand=.true. call getlun90(ifl1,1) call getlun90(ifl2,1) call getlun90(ifl3,1) print *, "blended file handles=",ifl1,ifl2,ifl3 CALL BAOPENR(ifl1,gfile1,ierr) if(ierr/=0)print*,'cant open ',trim(gfile1) CALL BAOPENR(ifl2,gfile2,ierr) if(ierr/=0)print*,'cant open ',trim(gfile2) call baopenw(ifl3,gfile3,ierr) print*,'Opened ',ifl3,'for grib2 data ', & trim(gfile3), 'return code is ',ierr itot=0 icount=0 iseek=0 do call skgb(ifl1,iseek,msk1,lskip,lgrib) if (lgrib==0) exit ! end loop at EOF or problem if (lgrib>currlen) then if (allocated(cgrib)) deallocate(cgrib) allocate(cgrib(lgrib),stat=ierr) currlen=lgrib endif call baread(ifl1,lskip,lgrib,lengrib,cgrib) if (lgrib/=lengrib) then print *,' degrib2: IO Error.' call errexit(9) endif iseek=lskip+lgrib icount=icount+1 PRINT * PRINT *,'GRIB MESSAGE ',icount,' starts at',lskip+1 PRINT * ! Unpack GRIB2 field call gb_info(cgrib,lengrib,listsec0,listsec1,& numfields,numlocal,maxlocal,ierr) if (ierr/=0) then write(6,*) ' ERROR querying GRIB2 message = ',ierr stop 10 endif itot=itot+numfields print *,' SECTION 0: ',(listsec0(j),j=1,3) print *,' SECTION 1: ',(listsec1(j),j=1,13) print *,' Contains ',numlocal,' Local Sections ', & ' and ',numfields,' data fields.' do n=1,numfields call gf_getfld(cgrib,lengrib,n,unpack,expand,gfld,ierr) if (ierr/=0) then write(6,*) ' ERROR extracting field = ',ierr cycle endif ! specify missing data values for different fields ! don't process fields other than icing, turublence, and CB ! 1. For icing severity and GTG, max of UK US ! 2. For CB top, max of US UK ! 3. For CB base, min of US UK ! 4. For CB extent, average of US UK if( gfld%ipdtmpl(1)==19 .and. gfld%ipdtmpl(2)==37)then ! ICESEV missing=-1. whichblnd='max' exclusive=.false. else if(gfld%ipdtmpl(1)==19 .and. gfld%ipdtmpl(2)==30)then ! EDPARM missing=-0.5 whichblnd='max' exclusive=.false. else if(gfld%ipdtmpl(1)==19 .and. gfld%ipdtmpl(2)==29)then ! CATEDR missing=-0.5 whichblnd='max' exclusive=.false. else if(gfld%ipdtmpl(1)==19 .and. gfld%ipdtmpl(2)==28)then ! MWTURB missing=-0.5 whichblnd='max' exclusive=.false. else if(gfld%ipdtmpl(1)== 6 .and. gfld%ipdtmpl(2)==25)then ! CB extent missing=-0.1 whichblnd='max' exclusive=.false. else if(gfld%ipdtmpl(2)== 3 .and. gfld%ipdtmpl(10)==11)then ! CB base missing=-1. whichblnd='min' exclusive=.false. else if(gfld%ipdtmpl(2)== 3 .and. gfld%ipdtmpl(10)==12)then ! CB top missing=-1. whichblnd='max' exclusive=.false. else cycle end if print * print *,' FIELD ',n if (n==1) then print *,' SECTION 0: ',gfld%discipline,gfld%version print *,' SECTION 1: ',(gfld%idsect(j),j=1,gfld%idsectlen) endif if ( associated(gfld%local).AND.gfld%locallen>0 ) then print *,' SECTION 2: ',gfld%locallen,' bytes' endif print *,' SECTION 3: ',gfld%griddef,gfld%ngrdpts, & gfld%numoct_opt,gfld%interp_opt,gfld%igdtnum print *,' GRID TEMPLATE 3.',gfld%igdtnum,': ', & (gfld%igdtmpl(j),j=1,gfld%igdtlen) if ( gfld%num_opt == 0 ) then print *,' NO Optional List Defining Number of Data Points.' else print *,' Section 3 Optional List: ', & (gfld%list_opt(j),j=1,gfld%num_opt) endif print *,' PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ', & (gfld%ipdtmpl(j),j=1,gfld%ipdtlen) pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),& gfld%ipdtmpl(2)) if ( gfld%num_coord == 0 ) then print *,' NO Optional Vertical Coordinate List.' else print *,' Section 4 Optional Coordinates: ',& (gfld%coord_list(j),j=1,gfld%num_coord) endif if ( gfld%ibmap /= 255 ) then print *,' Num. of Data Points = ',gfld%ndpts, & ' with BIT-MAP ',gfld%ibmap else print *,' Num. of Data Points = ',gfld%ndpts, & ' NO BIT-MAP ' endif print *,' DRS TEMPLATE 5.',gfld%idrtnum,': ', & (gfld%idrtmpl(j),j=1,gfld%idrtlen) fldmax=gfld%fld(1) fldmin=gfld%fld(1) sum=gfld%fld(1) do j=2,gfld%ndpts if (gfld%fld(j)>fldmax) fldmax=gfld%fld(j) if (gfld%fld(j)<fldmin) fldmin=gfld%fld(j) sum=sum+gfld%fld(j) enddo print *,' Data Values:' write(6,fmt='(" MIN=",f21.8," AVE=",f21.8, & " MAX=",f21.8)') fldmin,sum/gfld%ndpts,fldmax ! read UK's matching products print*,'reading UK WAFS' jids= -9999 ! jids(5)=listsec1(5) jids(6)=listsec1(6) ! year jids(7)=listsec1(7) ! mon jids(8)=listsec1(8) ! day jids(9)=listsec1(9) ! hr print*,'disipline,jpd= ',listsec0(1),jids(6:9) jpdtn=gfld%ipdtnum jpdt=-9999 jpdt(1)=gfld%ipdtmpl(1) ! cat number jpdt(2)=gfld%ipdtmpl(2) ! parm number jpdt(3)=gfld%ipdtmpl(3) ! (0-analysis, 1-forecast, or 2-analysis error) jpdt(9)=gfld%ipdtmpl(9) ! forecast hour jpdt(10)=gfld%ipdtmpl(10) ! level ID jpdt(12)=gfld%ipdtmpl(12) ! level value if(gfld%ipdtlen>=16) jpdt(16)=gfld%ipdtmpl(16) ! spatial statistical processing print*,'jpdtn,jpdt= ',jpdtn,jpdt(1:10) jgdtn=gfld%igdtnum jgdt=-9999 ! jgdt(1)=gfld%igdtmpl(1) jgdt(8)=gfld%igdtmpl(8) jgdt(9)=gfld%igdtmpl(9) print*,'jgdtn,jgdt= ',jgdtn,jgdt(1:9) call getgb2(ifl2,0,0,listsec0(1),jids,jpdtn,jpdt, & gfld%igdtnum,jgdt,.TRUE.,k,gfld2,ierr) print*,'US and UK dimensions= ',gfld%ndpts,gfld2%ndpts,"ierr=",ierr im=gfld%igdtmpl(8) jm=gfld%igdtmpl(9) allocate(blddata(im*jm)) if_ierr: if(ierr==0)then allocate(usdata(im,jm)) allocate(ukdata(im,jm)) do j=1,jm if(gfld%igdtmpl(12) == gfld2%igdtmpl(12)) then ! for template 3.0 jj=j else jj=jm-j+1 ! UK data is from south to north while US data from north to south endif do i=1,im usdata(i,j)=gfld%fld((j-1)*im+i) ukdata(i,j)=gfld2%fld((jj-1)*im+i) ij=(j-1)*im+i if(gfld%ibmap /= 255) then !If US with BIT-MAP if(.not. gfld%bmap((j-1)*im+i)) usdata(i,j)=missing endif if(gfld2%ibmap /= 255) then !If UK with BIT-MAP if(.not. gfld2%bmap((jj-1)*im+i)) ukdata(i,j)=missing end if blddata(ij)=generalblending(whichblnd,missing,exclusive,usdata(i,j),ukdata(i,j)) end do end do j=jm/2 print*,'sample 2D data ',(i,j,usdata(i,j),ukdata(i,j),i=210,220) ! Icing severity needs more processes. if(gfld%ipdtmpl(1)==19 .and. gfld%ipdtmpl(2)==37)then ! ICESEV ! 1. Blend by Gaussian Kernel Filter call gaussian_smooth(0,im,jm,1,1,missing,blddata) do j=1,jm do i=1,im ij=(j-1)*im+i ! 2. Keep the original matching data if(usdata(i,j) == ukdata(i,j)) blddata(ij)=usdata(i,j) ! 3. Re-categorize by different thresholds if(abs(blddata(ij)-missing)<=EPSILON) cycle if(blddata(ij) <= 0.8) then blddata(ij) = 0. elseif(blddata(ij) <= 1.5) then blddata(ij) = 1. elseif(blddata(ij) <= 2.4) then blddata(ij) = 2. elseif(blddata(ij) <= 3.2) then blddata(ij) = 3. else blddata(ij) = 4. end if ! 4. Not greater than max, and not smaller than min of US and UK blddata(ij)=min(blddata(ij),max(usdata(i,j),ukdata(i,j))) blddata(ij)=max(blddata(ij),min(usdata(i,j),ukdata(i,j))) end do end do end if call gf_free(gfld2) deallocate(usdata) deallocate(ukdata) else print*,'error code= ',ierr print*, pabbrev,' not found, writting US data as blended' blddata = gfld%fld if(gfld%ibmap /= 255) then where(.not. gfld%bmap) blddata = missing end if end if if_ierr ! MIN and MAX value after first step of blending. ! Used to set templete 5 elements do ij = 1, im*jm if(blddata(ij) /= missing) then i = ij fldmin = blddata(ij) fldmax = blddata(ij) exit end if end do do ij = i, im*jm if(blddata(ij) /= missing) then if(blddata(ij) > fldmax) fldmax = blddata(ij) if(blddata(ij) < fldmin) fldmin = blddata(ij) end if end do ! ngrdpts>=ndpts when bitmap is used (for underground gridpoints) call write_grib2(fldmin,fldmax,gfld%ngrdpts,blddata,ifl3,listsec1,& gfld%igdtnum,gfld%igdtlen,gfld%igdtmpl,& gfld%ipdtnum,gfld%ipdtlen,gfld%ipdtmpl,& gfld%idrtnum,gfld%idrtlen,gfld%idrtmpl) call gf_free(gfld) deallocate(blddata) enddo enddo print *," " print *, ' Total Number of Fields Found = ',itot call BACLOSE(ifl1, ierr) call BACLOSE(ifl2, ierr) call BACLOSE(ifl3, ierr) end subroutine process ! General blending of min, max or average of two values, ! when they are not missing values function generalblending(whichblnd,missing,exclusive,a,b) real :: generalblending character(3), intent(in) :: whichblnd real, intent(in) :: missing logical,intent(in) :: exclusive real, intent(in) :: a, b generalblending = missing if(a == missing .or. b == missing) then if(exclusive) then generalblending = missing elseif(a == missing) then generalblending = b else generalblending = a endif else if(whichblnd == 'max') then generalblending = max(a,b) elseif(whichblnd == 'min') then generalblending = min(a,b) elseif(whichblnd == 'avg') then generalblending = (a+b)/2. end if end if end function generalblending ! Abstract: smoothing a gridded field using Gaussian Kernel smoothing technique ! ! Reference: Amy Harless et al: "A report and feature-based verification study ! of the CAPS 2008 storm-scale ensemble forecast for severe convection ! weather", AMS Conference 2012 ! Input: ! iregion: 0 - global, 1 - regional ! im,jm: X and Y dimension of field A ! nbr: range of smoothing (in grids) ! s: sigma value of Gaussian smoothing ! A: gridded field to be smoothed ! Output: ! A: Smoothed field ! Programmer: ! 2015-12-02: Binbin Zhou, NCEP/EMC ! 2020-01-23: Y Mao, NCEP/EMC ! subroutine gaussian_smooth (iregion,im,jm,nbr,s,missing,A) implicit none integer, intent(in) :: iregion ! 0 - global, 1-regional integer, intent(in) :: im,jm integer, intent(in) :: nbr, s real, intent(in) :: missing real, intent(inout) :: A(im*jm) real :: B(im*jm) real :: f1,f2,G,Gsum,AxG integer :: i,j,ij,ii,ip,jp,ijp,i1,i2,j1,j2 f1=1./(3.14*s*s) f2=-0.5/(s*s) B=A do jp = 1,jm do ip = 1,im ijp=(jp-1)*im + ip i1 = ip - nbr i2 = ip + nbr j1 = jp - nbr j2 = jp + nbr if ( j1<=1 ) j1=1 if ( j2>=jm ) j2=jm Gsum=0. AxG=0. do j = j1,j2 do i = i1,i2 ii = i if( i<1 ) then if(iregion == 0 ) then ii=i+IM else ii=1 end if end if if ( i>im ) then if(iregion == 0 ) then ii=i-IM else ii=im end if end if ij = (j-1)*im + ii if(A(ij) /= missing) then G=f1*exp(f2*((ip-i)*(ip-i)+(jp-j)*(jp-j))) Gsum=Gsum+G AxG=AxG+G*A(ij) end if end do end do if(Gsum>0.0) B(ijp)=AxG/Gsum end do end do A=B return end subroutine gaussian_smooth subroutine write_grib2(min,max,npt,fld,lunout,listsec1in,& igdtnum,igdtlen,jgdt,ipdtnum,ipdtlen,jpdt,idrtnum,idrtlen,idrtmpl) !****************************************************************** ! prgmmr: pondeca org: np20 date: 2006-03-03 * ! * ! abstract: * ! use steve gilbert's gribcreate, addgrid, addfield, and gribend * ! subroutines to write data in Grib2 format * ! * ! program history log: * ! 2006-03-03 pondeca * ! 2020-01-24 Y Mao (remove nflds and its dependances) * ! * ! input argument list: * ! * ! 1. min: min value of array fld(npt) * ! * ! 2. max: max value of array fld(npt) * ! * ! 3. npt: size of data array to be written * ! * ! 4. fld: data array to be written * ! * ! 5. lunout: file handler/unit of output Grib2 file * ! * ! 6. listsec1in: array of reference time (year, month, day, hour, * ! minutes and seconds) * ! * ! 7. igdtnum: Grid Definition Template Number (Code Table 3.0) * ! * ! 8. igdtlen: length of grid template array * ! * ! 9. jgdt: array values of Grid Definition Template * ! * ! 10. ipdtnum: Product Definition Template Number (Code Table 4.0)* ! * ! 11. ipdtlen: length of product template array * ! * ! 12. jpdt: array values of Product Definition Template * ! * ! attributes: * ! language: f90 * !****************************************************************** implicit none real, intent(in) :: min,max integer, intent(in) :: npt real(4), intent(in) :: fld(npt) integer, intent(in) :: lunout integer, intent(in) :: listsec1in(:) integer, intent(in) :: igdtnum integer, intent(in) :: igdtlen integer, intent(in) :: jgdt(igdtlen) integer, intent(in) :: ipdtnum integer, intent(in) :: ipdtlen integer, intent(in) :: jpdt(ipdtlen) integer, intent(in) :: idrtnum integer, intent(in) :: idrtlen integer, intent(in) :: idrtmpl(idrtlen) integer(4), parameter :: idefnum=1 integer(4) :: max_bytes integer(4) :: listsec0(2),listsec1(13) integer(4) :: ierr,i,j,ij integer(4) :: lengrib integer(4) :: igds(5) integer(4) :: igdstmpl(igdtlen) integer(4) :: ideflist(idefnum) integer(4) :: ipdstmpl(ipdtlen) integer(4) :: numcoord integer(4) :: ibmap logical*1 :: bmap(npt) character*1,allocatable :: cgrib(:) real(4) :: coordlist max_bytes = npt*4 allocate(cgrib(max_bytes)) !==>initialize new GRIB2 message and pack ! GRIB2 sections 0 (Indicator Section) and 1 (Identification Section) listsec0(1)=0 ! Discipline-GRIB Master Table Number (see Code Table 0.0) listsec0(2)=2 ! GRIB Edition Number (currently 2) listsec1(:) = listsec1in(:) listsec1(1)=7 ! Id of orginating centre (Common Code Table C-1) listsec1(2)=4 !"EMC"! Id of orginating sub-centre (local table)/Table C of ON388 ! Yali Mao, GFS master table is 25 for WAFS at 0.25 deg (US unblended data is already set to 25) ! listsec1(3)=25 ! GRIB Master Tables Version Number (Code Table 1.0) listsec1(4)=1 ! per Brent! GRIB Local Tables Version Number (Code Table 1.1) listsec1(5)=1 ! Significance of Reference Time (Code Table 1.2) ! listsec1(6) ! Reference Time - Year (4 digits) ! listsec1(7) ! Reference Time - Month ! listsec1(8) ! Reference Time - Day ! listsec1(9) ! Reference Time - Hour listsec1(10) = 0 ! Reference Time - Minute listsec1(11) = 0 ! Reference Time - Second listsec1(12) = 0 ! Production status of data (Code Table 1.3) listsec1(13) = 1 ! Type of processed data (Code Table 1.4) ! 0 for analysis products and 1 for forecast products call gribcreate(cgrib,max_bytes,listsec0,listsec1,ierr) print*,'gribcreate status=',ierr !==> Pack up Grid Definition Section (Section 3) add to GRIB2 message. call apply_template_300(jgdt,igdstmpl) igds(1)=0 !Source of grid definition (see Code Table 3.0) igds(2)=npt !Number of grid points in the defined grid. igds(3)=0 !Number of octets needed for each additional grid points definition igds(4)=0 !Interpretation of list for optional points definition (Code Table 3.11) igds(5)=igdtnum !Grid Definition Template Number (Code Table 3.1) ideflist=0 !Used if igds(3) /= 0. Dummy array otherwise call addgrid(cgrib,max_bytes,igds,igdstmpl,igdtlen,ideflist,idefnum,ierr) print*,'addgrid status=',ierr !==> pack up sections 4 through 7 for a given field and add them to a GRIB2 message. ! They are Product Definition Section, Data Representation Section, Bit-Map Section ! and Data Section, respectively. call apply_template_40(jpdt,ipdstmpl) print*,'product template in new Grib file= ',ipdstmpl ! Use US unblended template 5 information (5.40) !!! idrtnum=40 !Data Representation Template Number ( see Code Table 5.0 ) !!! call apply_template_50(min,max,jpdt(1),jpdt(2),idrtmpl) numcoord=0 coordlist=0. !needed for hybrid vertical coordinate ibmap=255 ! Bitmap indicator ( see Code Table 6.0 ), WAFS Blended products do not use bit-map print *, "npt=",npt print *, "ipdtnum=",ipdtnum,ipdtlen,ipdstmpl print *, "coordlist=",coordlist,numcoord print *, "idrtnum=",idrtnum,idrtlen,idrtmpl call addfield(cgrib,max_bytes,ipdtnum,ipdstmpl,ipdtlen, & coordlist,numcoord,idrtnum,idrtmpl, & idrtlen,fld,npt,ibmap,bmap,ierr) print*,'addfield status=',ierr !==> finalize GRIB message after all grids ! and fields have been added. It adds the End Section ( "7777" ) call gribend(cgrib,max_bytes,lengrib,ierr) print*,'gribend status=',ierr print*,'length of the final GRIB2 message in octets =',lengrib call wryte(lunout, lengrib, cgrib) deallocate(cgrib) return end subroutine write_grib2 !=========================================================================== subroutine apply_template_300(jgdt,ifield3) implicit none integer,intent(in)::jgdt(:) integer(4),intent(out) :: ifield3(:) ifield3 = jgdt ifield3(11) = 0 return end subroutine apply_template_300 !=========================================================================== subroutine apply_template_40(jpdt,ifield4) implicit none integer,intent(in) :: jpdt(:) integer(4),intent(out) :: ifield4(:) ifield4 = jpdt !==> ifield4(1):parameter category (see Code Table 4.1) !==> ifield4(2):parameter number (see Code Table 4.2) !==> ifield4(3):type of generating process (see Code Table 4.3) ! 0 - analysis ! 2 - forecast ! 7 - analysis error !==>ifield4(4):background generating process identifier ! (defined by originating Center) ifield4(4) = 0 !hasn't been defined yet !==>ifield4(5):analysis or forecast generating process identifier ! (defined by originating Center) ifield4(5) = 96 !==>ifield4(6):hours of observational data cutoff after reference time !==>ifield4(7):minutes of observational data cutoff after reference time ifield4(6) = 0 ! per steve ifield4(7) = 0 !==>ifield4(8):indicator of unit of time range (see Code Table 4.4) ifield4(8) = 1 !==>ifield4(9):forecast time in units defined by ifield4(8) !==>ifield4(10):type of first fixed surface (see Code Table 4.5) !==>ifield4(11):scale factor of first fixed surface ifield4(11) = 0 !Because not saving any precision !==>ifield4(12):scaled value of first fixed surface !==>ifield4(13):type of second fixed surface(See Code Table 4.5) !==>ifield4(14):scale factor of second fixed surface !==>ifield4(15):scaled value of second fixed surface ifield4(13) = 255 ifield4(14) = 0 ifield4(15) = 0 if(size(jpdt)>=16) then !==> ifield4(16):Statistical process used within the spatial area (see Code Table 4.10) ifield4(17) = 3 ! Type of spatial processing ifield4(18) = 1 ! Number of data points used in spatial processing end if return end subroutine apply_template_40 !=========================================================================== subroutine apply_template_50(min,max,ncat,nparm,ifield5) implicit none real,intent(in) :: min,max integer,intent(in) :: ncat,nparm integer(4),intent(out) :: ifield5(:) ! reference value(R) (IEEE 32-bit floating-point value) ifield5(1)=0 ! Any value. Will be overwritten ifield5(2)=0! binary scale factor (E) ! decimal scale factor (D) if (ncat==19 .and. nparm==37) then ! ICESEV ifield5(3) = 0 else if (ncat==19 .and. nparm==30) then ! EDR ifield5(3) = 2 else if(ncat==3 .and. nparm==3) then ! CB base/top ifield5(3) = 0 else if (ncat==6 .and. nparm==25) then ! Cb ext ifield5(3) = 1 else ifield5(3) = 2 endif ! number of bits used for each packed value for simple packing ! or for each group reference value for complex packing or ! spatial differencing ifield5(4) = 0 ! Must reset to 0 ifield5(5) = 0 ! type of original field values (See Code Table 5.1) ! Rarely happens for WAFS data, just in case if(min == max) then ifield5(2)=0 ifield5(3)=0 end if return end subroutine apply_template_50 !=========================================================================== !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: GETLUN GET UNIQUE LOGICAL UNIT NUMBERS ! PRGMMR: SMITH, TRACY ORG: FSL/PROFS DATE: 90-06-15 ! ! ABSTRACT: THIS PROGRAM GETS UNIQUE LOGICAL UNIT NUMBERS FOR OPFILE ! OR RETURNS THEM TO THE POOL FOR CLFILE. ! ! PROGRAM HISTORY LOG: ! FORTRAN 90 VERSION IS GETLUN90: PONDECA, DATE: 2006-03-08 ! ! USAGE: CALL GETLUN(LUN,OPTN) ! INPUT ARGUMENT LIST: ! LUN - INTEGER LOGICAL UNIT NUMBER ! OPTN - INTEGER CNCT=1, DSCT=2. ! IF CONNECTING A FILE(CNCT) SET THE NUMBER TO ! NEGATIVE SO IT WON'T BE USED UNTIL AFTER ! DSCT SETS IT POSITIVE. ! ! OUTPUT ARGUMENT LIST: ! LUN - INTEGER LOGICAL UNIT NUMBER ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN-90 ! MACHINE: NAS-9000 !$$$ SUBROUTINE GETLUN90(LUN,OPTN) !* THIS PROGRAM GETS UNIQUE LOGICAL UNIT NUMBERS FOR OPFILE !* OR RETURNS THEM TO THE POOL FOR CLFILE IMPLICIT NONE INTEGER, PARAMETER :: CNCT=1,DSCT=2 INTEGER :: LUN,OPTN,I INTEGER :: NUM(80)=(/ & 99,98,97,96,95,94,93,92,91,90, & 89,88,87,86,85,84,83,82,81,80, & 79,78,77,76,75,74,73,72,71,70, & 69,68,67,66,65,64,63,62,61,60, & 59,58,57,56,55,54,53,52,51,50, & 49,48,47,46,45,44,43,42,41,40, & 39,38,37,36,35,34,33,32,31,30, & 29,28,27,26,25,24,23,22,21,20 /) !* START IF(OPTN == CNCT) THEN DO I=1,80 IF(NUM(I)>0) THEN LUN=NUM(I) NUM(I)=-NUM(I) return ENDIF END DO PRINT*, 'NEED MORE THAN 80 UNIT NUMBERS' ELSE IF(OPTN == DSCT) THEN !* MAKE THE NUMBER AVAILABLE BY SETTING POSITIVE DO I=1,80 IF(LUN == -NUM(I)) NUM(I)=ABS(NUM(I)) ENDDO END IF RETURN END SUBROUTINE GETLUN90 end module blending program main use blending implicit none character(60) :: gfile1,gfile2,gfile3 INTEGER :: NARG call start() ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! GET ARGUMENTS NARG=IARGC() IF(NARG /= 3) THEN CALL ERRMSG('blending: Incorrect usage') CALL ERRMSG('Usage: wafs_blending_0p25 grib2file1 grib2file2 grib2file3') CALL ERREXIT(2) ENDIF CALL GETARG(1,gfile1) CALL GETARG(2,gfile2) CALL GETARG(3,gfile3) call process(trim(gfile1),trim(gfile2),trim(gfile3)) end program main