program calcmint !****************************************************************************** ! abstract: Compute minimum temperature obs from a given list of stations ! (local fixed file) according to NDFD definition (7P-8A local). ! Local time zones are given in fixed file. ! ! input: ! Previous 24 URMA prepbufr files ! ! output: ! 1 text file with minT observations for given sites. To be read in ! later by GSI. ! ! program history log: ! 2017-04-24 Levine ! !****************************************************************************** implicit none interface subroutine indictb(stnid,buflat,buflon,mtrids,mtrlat,mtrlon,mtrmax,inmtr) implicit none integer(4),parameter::max=99999 integer(4)::i character(len=8),intent(in)::stnid real(4),intent(in)::buflat,buflon real(4),intent(in)::mtrlat(max),mtrlon(max) integer(4),intent(in)::mtrmax character(len=8),intent(in)::mtrids(max) integer(4),intent(out)::inmtr real(4)::latdiff,londiff end subroutine indictb end interface character(len=8)::stnid,ctyp,prov,sprov character(len=23)::namea integer(4),parameter::max=99999 character(len=8)::metars(max),akmetars(max),ships(max),akships(max),mesoa(max),mesob(max),mesoc(max),mesod(max),mesoe(max),mesof(max),mesoak(max),mesohi(max),mesopr(max),mesogu(max) character(len=8)::prva(max),prvb(max),prvc(max),prvd(max),prve(max),prvf(max),prvak(max),prvhi(max),prvpr(max),prvgu(max),prvmtr(max),prvshp(max),akprvmtr(max),akprvshp(max) character(len=8)::sprva(max),sprvb(max),sprvc(max),sprvd(max),sprve(max),sprvf(max),sprvak(max),sprvhi(max),sprvpr(max),sprvgu(max),sprvmtr(max),sprvshp(max),aksprvmtr(max),aksprvshp(max) character(len=20)::metnames(max),shpnames(max),msonames(max),namesak(max) !filled with dummpy arguements real(4)::mtrlat(max),mtrlon(max),shplat(max),shplon(max),akmtrlat(max),akmtrlon(max),akshplat(max),akshplon(max),mtrelv(max),akmtrelv(max),shpelv(max),akshpelv(max) real(4)::mtrtyp(max),mtrdmp(max),shptyp(max),shpdmp(max),akmtrtyp(max),akmtrdmp(max),akshptyp(max),akshpdmp(max) real(4)::msolata(max),msolatb(max),msolatc(max),msolatd(max),msolate(max),msolatf(max),msolatak(max),msolathi(max),msolatpr(max),msolatgu(max) real(4)::msolona(max),msolonb(max),msolonc(max),msolond(max),msolone(max),msolonf(max),msolonak(max),msolonhi(max),msolonpr(max),msolongu(max) real(4)::msodmpa(max),msodmpb(max),msodmpc(max),msodmpd(max),msodmpe(max),msodmpf(max),msodmpak(max),msodmphi(max),msodmppr(max),msodmpgu(max) real(4)::msotypa(max),msotypb(max),msotypc(max),msotypd(max),msotype(max),msotypf(max),msotypak(max),msotyphi(max),msotyppr(max),msotypgu(max) real(4)::msoelva(max),msoelvb(max),msoelvc(max),msoelvd(max),msoelve(max),msoelvf(max),msoelvak(max),msoelvhi(max),msoelvpr(max),msoelvgu(max) real(4),parameter::misreal=9999. !changes to matach with zmiss in MDL sub integer(4),parameter::misint=9999,mxtb=1000000 !urma prepbufr reports contan ~400000 subsets, this should be enough integer(4)::i,j,bufnum,k,m,nstn,idate,iyear,imonth,iday,doy,jdn,nsub,mmax,mtrmax,shpmax,tz,akmtrmax,akshpmax integer(4)::msomaxa,msomaxb,msomaxc,msomaxd,msomaxe,msomaxf,msomaxak,msomaxhi,msomaxpr,msomaxgu,ireadmg,ireadsb,ibfms,nlevs integer(4)::mtrzone(max),msozonea(max),msozoneb(max),msozonec(max),msozoned(max),msozonee(max),msozonef(max),msozoneak(max),msozonehi(max),msozonepr(max),msozonegu(max),shpzone(max),akmtrzone(max),akshpzone(max) character(len=2)::hchar character(len=1)::ns,ew character(len=95)::fname logical::fexist1 real(8)::vars(9),maxs(4,255),tpc(255,20),tobs(2,255,20) real(4)::mtrhrly(max,25),akmtrhrly(max,25),mtrmaxs(max,4),akmtrmaxs(max,4),shpmaxs(max,4),shphrly(max,25),akshphrly(max,25),msomaxs(max,4),rpt(max,25) real(4)::msohrlya(max,25),msohrlyb(max,25),msohrlyc(max,25),msohrlyd(max,25),msohrlye(max,25),msohrlyf(max,25),msohrlyak(max,25),msohrlyhi(max,25),msohrlypr(max,25),msohrlygu(max,25) real(4)::maxt(max) real(4),parameter::pi=3.1415926535897932384626433832795 real(4)::lat,lon,buflat,buflon,buftim,buftmp,bufmax,buftyp,latdif,londif,bufdmp,bufqc,vtcd,dtyp,ddmp,delv,bufelv integer(4)::hours(25) !initialize arrays and dummy variables (variables are defined in comment as initialized) mtrzone(:)=misint !time zone for metar (ADPSFC) sites akmtrzone(:)=misint msozonea(:)=-5 !time zone for mesonets, zone a msozoneb(:)=misint !time zone for mesonets, zone b msozonec(:)=-6 !time zone for mesonets, zone c msozoned(:)=misint !time zone for mesonets, zone d msozonee(:)=misint !time zone for mesonets, zone e msozonef(:)=misint !time zone for mesonets, zone f msozonehi(:)=-9 !time zone for mesonets, zone g msozoneak(:)=misint !time zone for mesonets, zone h msozonepr(:)=-4 msozonegu(:)=5 shpzone(:)=misint !time zone for SFCSHP sites akshpzone(:)=misint mtrlat(:)=misreal !latitude for ADPSFC/METAR sites mtrlon(:)=misreal !longitude for ADPSFC/METAR sites akmtrlat(:)=misreal akmtrlon(:)=misreal msolata(:)=misreal !latitude for mesonets, zone a msolona(:)=misreal !longitude for mesonets, zone a msolatb(:)=misreal !latitude for mesonets, zone b msolonb(:)=misreal !longitude for mesonets, zone b msolatc(:)=misreal !latitude for mesonets, zone c msolonc(:)=misreal !longitude for mesonets, zone c msolatd(:)=misreal !latitude for mesonets, zone d msolond(:)=misreal !longitude for mesonets, zone d msolate(:)=misreal !latitude for mesonets, zone e msolone(:)=misreal !longitude for mesonets, zone e msolatf(:)=misreal !latitude for mesonets, zone f msolonf(:)=misreal msolatpr(:)=misreal msolonpr(:)=misreal msolathi(:)=misreal msolonhi(:)=misreal msolatgu(:)=misreal msolongu(:)=misreal msolatak(:)=misreal msolonak(:)=misreal shplat(:)=misreal !latitude for SHCSHP sites shplon(:)=misreal !longitude for SFCSHP sites akshplat(:)=misreal akshplon(:)=misreal mtrhrly(:,:)=misreal !hourly observations (station number, hour number) for ADPSFC/METAR sites mtrmaxs(:,:)=misreal !6hourly maxT observations from ADPSFC/METAR sites (station number, 1-4: 1=12Z, 2=18Z, 3=00Z, 4=06Z) akmtrmaxs(:,:)=misreal shphrly(:,:)=misreal !hourly observations (station number, hour number) for SFCSHP sites akmtrhrly(:,:)=misreal akshphrly(:,:)=misreal shpmaxs(:,:)=misreal !6hourly maxT observations from shcfhp sites (dummy arguement - will ALWAYS be missing) msomaxs(:,:)=misreal !6hourly maxT observations from mesonet sites (dummy arguement - will ALWAYS be missing) msohrlya(:,:)=misreal !hourly observations (station number, hour number) for MSONET sites, zone a msohrlyb(:,:)=misreal !hourly observations (station number, hour number) for MSONET sites, zone b msohrlyc(:,:)=misreal !hourly observations (station number, hour number) for MSONET sites, zone c msohrlyd(:,:)=misreal !hourly observations (station number, hour number) for MSONET sites, zone d msohrlye(:,:)=misreal !hourly observations (station number, hour number) for MSONET sites, zone e msohrlyf(:,:)=misreal !hourly observations (station number, hour number) for MSONET sites, zone f msohrlypr(:,:)=misreal !hourly observations (station number, hour number) for MSONET sites, zone pr msohrlyhi(:,:)=misreal !hourly observations (station number, hour number) for MSONET sites, zone hi msohrlyak(:,:)=misreal !hourly observations (station number, hour number) for MSONET sites, zone ak msohrlygu(:,:)=misreal maxt(:)=misreal !maximum temperatures (output from calcmaxt subroutine) rpt(:,:)=misreal !reporting time for metar stations (station number, hour number) metars(:)="MISSING " !METAR/ADPSFC stationIDs akmetars(:)="MISSING" mesoa(:)="MISSING " !MSONET stationIDS, zone a mesob(:)="MISSING " !MSONET stationIDS, zone b mesoc(:)="MISSING " !MSONET stationIDS, zone c mesod(:)="MISSING " !MSONET stationIDS, zone d mesoe(:)="MISSING " !MSONET stationIDS, zone e mesof(:)="MISSING " !MSONET stationIDS, zone f mesopr(:)="MISSING " !MSONET stationIDS, zone pr mesohi(:)="MISSING " !MSONET stationIDS, zone hi mesoak(:)="MISSING " mesogu(:)="MISSING " ships(:)="MISSING " !SFCSHP stationIDs akships(:)="MISSING " metnames(:)="MTR NAME GOES HERE " !station description (dummy arguement) shpnames(:)="SHIP NAME GOES HERE " !station description (dummy arguement) msonames(:)="MESO NAME GOES HERE " !station description (dummy arguement) namesak(:)="NAMES GO HERE AK" prva(:)="MISSING " prvb(:)="MISSING " prvc(:)="MISSING " prvd(:)="MISSING " prve(:)="MISSING " prvf(:)="MISSING " prvak(:)="MISSING " prvhi(:)="MISSING " prvpr(:)="MISSING " prvgu(:)="MISSING " prvmtr(:)="ADPSFC " prvshp(:)="SFCSHP " akprvmtr(:)="MISSING " akprvshp(:)="MISSING " sprva(:)="MISSING " sprvb(:)="MISSING " sprvc(:)="MISSING " sprvd(:)="MISSING " sprve(:)="MISSING " sprvf(:)="MISSING " sprvak(:)="MISSING " sprvhi(:)="MISSING " sprvpr(:)="MISSING " sprvgu(:)="MISSING " sprvmtr(:)="MISSING " sprvshp(:)="MISSING " aksprvmtr(:)="MISSING " aksprvshp(:)="MISSING " hours=(/18,19,20,21,22,23,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18/) !first loop through dictionary file to get list of METARS, location and time zone do i=1,max read(11,459,end=10) stnid,delv,dtyp,ddmp,lat,lon,tz metars(i)=stnid mtrlat(i)=lat mtrlon(i)=lon mtrzone(i)=tz mtrtyp(i)=dtyp mtrdmp(i)=ddmp mtrelv(i)=delv enddo 345 format (A8,1X,F6.3,1X,F8.3,1X,I3) 10 continue mtrmax=i-1 print*, "Total number of METARS/ADPSFC stations:",mtrmax close(11) !now the file with "new" Alaskan adpsfc stations do i=1,max read(12,459,end=20) stnid,delv,dtyp,ddmp,lat,lon,tz akmetars(i)=stnid akmtrlat(i)=lat akmtrlon(i)=lon akmtrzone(i)=tz akmtrtyp(i)=dtyp akmtrdmp(i)=ddmp akmtrelv(i)=delv if (tz.gt.-9) then namesak(i)="NAMES GO HERE CN" endif enddo 20 continue akmtrmax=i-1 print*, "Number of AK METAR/ADPSFC stations:",akmtrmax close(12) !now loop through mesonet files, files are already sorted by zone do i=1,max read(13,458,end=40) stnid,prov,delv,dtyp,ddmp,lat,lon,tz mesoa(i)=stnid msolata(i)=lat msolona(i)=lon msozonea(i)=tz prva(i)=prov msotypa(i)=dtyp msoelva(i)=delv msodmpa(i)=ddmp enddo 40 continue msomaxa=i-1 print*, "Total number of mesonet-a stations:",msomaxa close(13) do i=1,max read(14,458,end=41) stnid,prov,delv,dtyp,ddmp,lat,lon,tz !stnid,lat,lon,tz mesob(i)=stnid msolatb(i)=lat msolonb(i)=lon msozoneb(i)=tz prvb(i)=prov msotypb(i)=dtyp msoelvb(i)=delv msodmpb(i)=ddmp enddo 41 continue msomaxb=i-1 print*, "Total number of mesonet-b stations:",msomaxb close(14) do i=1,max read(15,458,end=42) stnid,prov,delv,dtyp,ddmp,lat,lon,tz !stnid,lat,lon,tz mesoc(i)=stnid msolatc(i)=lat msolonc(i)=lon msozonec(i)=tz prvc(i)=prov msotypc(i)=dtyp msoelvc(i)=delv msodmpc(i)=ddmp enddo 42 continue msomaxc=i-1 print*, "Total number of mesonet-c stations:",msomaxc close(15) do i=1,max read(16,458,end=43) stnid,prov,delv,dtyp,ddmp,lat,lon,tz !stnid,lat,lon,tz mesod(i)=stnid msolatd(i)=lat msolond(i)=lon msozoned(i)=tz prvd(i)=prov msotypd(i)=dtyp msoelvd(i)=delv msodmpd(i)=ddmp enddo 43 continue msomaxd=i-1 print*, "Total number of mesonet-d stations:",msomaxd close(16) do i=1,max read(17,458,end=44) stnid,prov,delv,dtyp,ddmp,lat,lon,tz !stnid,lat,lon,tz mesoe(i)=stnid msolate(i)=lat msolone(i)=lon msozonee(i)=tz prve(i)=prov msotype(i)=dtyp msoelve(i)=delv msodmpe(i)=ddmp enddo 44 continue msomaxe=i-1 print*, "Total number of mesonet-e stations:",msomaxe close(17) do i=1,max read(18,458,end=45) stnid,prov,delv,dtyp,ddmp,lat,lon,tz !stnid,lat,lon,tz mesof(i)=stnid msolatf(i)=lat msolonf(i)=lon msozonef(i)=tz prvf(i)=prov msotypf(i)=dtyp msoelvf(i)=delv msodmpf(i)=ddmp enddo 45 continue msomaxf=i-1 print*, "Total number of mesonet-f stations:",msomaxf close(18) do i=1,max read(19,458,end=52) stnid,prov,delv,dtyp,ddmp,lat,lon,tz !stnid,lat,lon,tz mesoak(i)=stnid msolatak(i)=lat msolonak(i)=lon msozoneak(i)=tz prvak(i)=prov msotypak(i)=dtyp msoelvak(i)=delv msodmpak(i)=ddmp enddo 52 continue msomaxak=i-1 print*, "Total number of mesonet-ak stations:",msomaxak close(19) do i=1,max read(20,458,end=53) stnid,prov,delv,dtyp,ddmp,lat,lon,tz !stnid,lat,lon,tz mesopr(i)=stnid msolatpr(i)=lat msolonpr(i)=lon msozonepr(i)=tz prvpr(i)=prov msotyppr(i)=dtyp msoelvpr(i)=delv msodmppr(i)=ddmp enddo 53 continue msomaxpr=i-1 print*, "Total number of mesonet-pr stations:",msomaxpr close(20) do i=1,max read(21,458,end=54) stnid,prov,delv,dtyp,ddmp,lat,lon,tz !stnid,lat,lon,tz mesohi(i)=stnid msolathi(i)=lat msolonhi(i)=lon msozonehi(i)=tz prvhi(i)=prov msotyphi(i)=dtyp msoelvhi(i)=delv msodmphi(i)=ddmp enddo 54 continue msomaxhi=i-1 print*, "Total number of mesonet-hi stations:",msomaxhi close(21) do i=1,max read(22,458,end=55) stnid,prov,delv,dtyp,ddmp,lat,lon,tz !stnid,lat,lon,tz mesogu(i)=stnid msolatgu(i)=lat msolongu(i)=lon msozonegu(i)=tz prvgu(i)=prov msotypgu(i)=dtyp msoelvgu(i)=delv msodmpgu(i)=ddmp enddo 55 continue msomaxgu=i-1 print*, "Total number of mesonet-gu stations:",msomaxgu close(22) do i=1,max read(23,459,end=56) stnid,delv,dtyp,ddmp,lat,lon,tz ships(i)=stnid shplat(i)=lat shplon(i)=lon shpzone(i)=tz shpelv(i)=delv shpdmp(i)=ddmp shptyp(i)=dtyp enddo 56 continue shpmax=i-1 print*, "Total number of ship stations:",shpmax close(23) do i=1,max read(24,459,end=57) stnid,delv,dtyp,ddmp,lat,lon,tz akships(i)=stnid akshplat(i)=lat akshplon(i)=lon akshpzone(i)=tz akshpelv(i)=delv akshpdmp(i)=ddmp akshptyp(i)=dtyp enddo 57 continue akshpmax=i-1 print*, "Total number of ship-ak stations:",akshpmax close(24) !open the first prepbufr file (06Z previous day) and record hourly temperatures vars(:)=misreal maxs(:,:)=misreal call openbf(25,"IN",25) call datelen(10) call ufbqcd(25,'VIRTMP',vtcd) do while (ireadmg(25,ctyp,idate).eq.0) !ADPSFC/METAR procedure: Find station in ADPSFC/METAR dictionary. If in dictionary, see if we have ob !for this station/cycle. If so, record the temp ob as an hourly ob. If there is more than one ob/hour, record the highest ob. !Since this is for the first file (06Z from previous day), any maxT here is irrelevant and therefore ignored if (ctyp.eq. "ADPSFC ") then do while(ireadsb(25).eq.0) vars(:)=misreal !vars is output array from ufbint (R8ARR in BUFRLIB documentation). Reset this array to misreal before reading any new data. call ufbint(25,vars,9,1,nlevs,"SID PRVSTG SPRVSTG YOB XOB RPT T29 TYP ELV")! TOB TQM") if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 418. Nlevs=",nlevs stop endif write(stnid,'(a8)') vars(1) write(prov,'(a8)') vars(2) write(sprov,'(a8)') vars(3) buftyp=vars(8) bufdmp=vars(7) buflat=vars(4) if (vars(5).gt.180) then buflon=vars(5)-360 else buflon=vars(5) endif buftim=vars(6) bufelv=vars(9) if ((buftyp.eq.187.or.buftyp.eq.183.or.buftyp.eq.192.or.buftyp.eq.193).and.(bufdmp.eq.512.or.bufdmp.eq.511)) then tpc(:,:)=misreal call ufbevn(25,tpc,1,255,20,nlevs,'TPC') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 438. Nlevs=",nlevs stop endif tobs(:,:,:)=misreal call ufbevn(25,tobs,2,255,20,nlevs,'TOB TQM') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 444. Nlevs=",nlevs stop endif !initialize buftmp and bufqc to prevent runtime error later buftmp=misreal bufqc=misreal do k=1,nlevs do j=1,20 if (tpc(k,j)==1) then buftmp=tobs(1,k,j) bufqc=tobs(2,k,j) endif if (tpc(k,j)>1000) exit enddo enddo if (buftmp.gt.1000) buftmp = misreal !in case buftmp is missing 10**10 goes to misreal if (buftmp.lt.1000) then !skip if temperature ob is invalid/missing if (bufqc.le.3) then !see if station is in dictionary nstn=0 !nstn=station number in relevant dictionary (0 if not in dictionary), output from indictb subroutine if (buflat.ge.45.and.(buflon.le.-125.or.buflon.ge.155)) then call indictb(stnid,buflat,buflon,akmetars,akmtrlat,akmtrlon,akmtrmax,nstn) if (nstn.gt.0) then if (buftyp.ne.akmtrtyp(nstn)) then akmtrtyp(nstn) = buftyp akmtrdmp(nstn) = bufdmp akmtrelv(nstn) = bufelv endif if (akprvmtr(nstn).eq."MISSING ") then akprvmtr(nstn)=prov aksprvmtr(nstn)=sprov endif if (akmtrhrly(nstn,1).ne.misreal) then if (buftmp.lt.akmtrhrly(nstn,1)) then akmtrhrly(nstn,1)=buftmp endif else !no conflict akmtrhrly(nstn,1)=buftmp endif else call indictb(stnid,buflat,buflon,metars,mtrlat,mtrlon,mtrmax,nstn) if (nstn.gt.0) then if (buftyp.ne.mtrtyp(nstn)) then mtrtyp(nstn) = buftyp mtrdmp(nstn) = bufdmp mtrelv(nstn) = bufelv endif !if there is more than one ob in an hour, we want the one with the highest temp if (prvmtr(nstn).eq."ADPSFC ") then prvmtr(nstn)=prov sprvmtr(nstn)=sprov endif if (mtrhrly(nstn,1).ne.misreal) then if (buftmp.lt.mtrhrly(nstn,1)) then mtrhrly(nstn,1)=buftmp endif else !no conflict mtrhrly(nstn,1)=buftmp endif endif endif endif endif !else ! print*, "WARNING: Invalid ADPSFC temperature report (hour 1) (stnid,type,temp,bufqc):",stnid,buftyp,buftmp,bufqc endif endif !right type enddo !MESONET procedure: First find out which dictionary station is in. Find station in dictionary. If in dictionary, see if we have ob !for this station/cycle. If so, record the temp ob as an hourly ob. If there is more than one ob/hour, record the highest ob. !Recall that raw reported maxTs are not reported for mesonets else if (ctyp.eq."MSONET ") then do while (ireadsb(25).eq.0) vars(:)=misreal call ufbint(25,vars,9,1,nlevs,"SID PRVSTG SPRVSTG YOB XOB RPT T29 TYP ELV")! TOB TQM") if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevs error line 520 hour,nlevs:",i,nlevs stop endif write(stnid,'(a8)') vars(1) write(prov,'(a8)') vars(2) write(sprov,'(a8)') vars(3) !make sure trailing x and a are removed from stationIDS because they just cause confusion if (stnid(8:8).eq."a".or.stnid(8:8).eq."x") then stnid(8:8)=" " endif buftyp=vars(8) bufdmp=vars(7) buflat=vars(4) if (vars(5).gt.180) then buflon=vars(5)-360 else buflon=vars(5) endif buftim=vars(6) bufelv=vars(9) if ((buftyp.eq.188.or.buftyp.eq.195).and.bufdmp.eq.540) then tpc(:,:)=misreal call ufbevn(25,tpc,1,255,20,nlevs,'TPC') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 544. Nlevs=",nlevs stop endif tobs(:,:,:)=misreal call ufbevn(25,tobs,2,255,20,nlevs,'TOB TQM') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 552. Nlevs=",nlevs stop endif !initialize buftmp and bufqc to prevent possible run time error later buftmp=misreal bufqc=misreal do k=1,nlevs do j=1,20 if (tpc(k,j)==1) then buftmp=tobs(1,k,j) bufqc=tobs(2,k,j) endif if (tpc(k,j)>1000) exit enddo enddo if (buftmp.gt.1000) buftmp = misreal if (buftmp.lt.1000) then if (bufqc.le.3) then nstn=0 if (buflon.ge.149.and.buflon.le.162.and.buflat.ge.12.and.buflat.le.17) then ! GUAM call indictb(stnid,buflat,buflon,mesogu,msolatgu,msolongu,msomaxgu,nstn) if (nstn.gt.0) then if (buftyp.ne.msotypgu(nstn)) then msotypgu(nstn) = buftyp msodmpgu(nstn) = bufdmp msoelvgu(nstn) = bufelv endif if (msohrlygu(nstn,1).ne.misreal) then if (buftmp.lt.msohrlygu(nstn,1)) then msohrlygu(nstn,1)=buftmp endif else msohrlygu(nstn,1)=buftmp endif endif elseif (buflon.ge.-149.and.buflon.le.-143.and.buflat.ge.17.and.buflat.le.24) then !HAWAII call indictb(stnid,buflat,buflon,mesohi,msolathi,msolonhi,msomaxhi,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotyphi(nstn)) then msotyphi(nstn) = buftyp msodmphi(nstn) = bufdmp msoelvhi(nstn) = bufelv endif if (msohrlyhi(nstn,1).ne.misreal) then if (buftmp.lt.msohrlya(nstn,1)) then msohrlyhi(nstn,1)=buftmp endif else !otherwise just load the temperature ob msohrlyhi(nstn,1)=buftmp endif endif elseif (buflon.ge.-69.and.buflon.le.-61.and.buflat.ge.16.and.buflat.le.20) then !PUERTO RICO call indictb(stnid,buflat,buflon,mesopr,msolatpr,msolonpr,msomaxpr,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotyppr(nstn)) then msotyppr(nstn) = buftyp msodmppr(nstn) = bufdmp msoelvpr(nstn) = bufelv endif if (msohrlypr(nstn,1).ne.misreal) then if (buftmp.lt.msohrlypr(nstn,1)) then msohrlypr(nstn,1)=buftmp endif else !otherwise just load the temperature ob msohrlypr(nstn,1)=buftmp endif endif elseif (buflat.ge.45.and.(buflon.le.-125.or.buflon.ge.155)) then !ALASKA call indictb(stnid,buflat,buflon,mesoak,msolatak,msolonak,msomaxak,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypak(nstn)) then msotypak(nstn) = buftyp msodmpak(nstn) = bufdmp msoelvak(nstn) = bufelv endif if (msohrlyak(nstn,1).ne.misreal) then if (buftmp.lt.msohrlyak(nstn,1)) then msohrlyak(nstn,1)=buftmp endif else !otherwise just load the temperature ob msohrlyak(nstn,1)=buftmp endif endif elseif (buflon.ge.-82.5.and.buflon.le.-70.and.buflat.ge.25.and.buflat.le.50) then !CONUS A call indictb(stnid,buflat,buflon,mesoa,msolata,msolona,msomaxa,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypa(nstn)) then msotypa(nstn) = buftyp msodmpa(nstn) = bufdmp msoelva(nstn) = bufelv endif if (msohrlya(nstn,1).ne.misreal) then if (buftmp.lt.msohrlya(nstn,1)) then msohrlya(nstn,1)=buftmp endif else !otherwise just load the temperature ob msohrlya(nstn,1)=buftmp endif endif elseif (buflon.ge.-91.and.buflon.le.-82.5.and.buflat.ge.25.and.buflat.le.50) then !CONUS B call indictb(stnid,buflat,buflon,mesob,msolatb,msolonb,msomaxb,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypb(nstn)) then msotypb(nstn) = buftyp msodmpb(nstn) = bufdmp msoelvb(nstn) = bufelv endif if (msohrlyb(nstn,1).ne.misreal) then if (buftmp.lt.msohrlyb(nstn,1)) then msohrlyb(nstn,1)=buftmp endif else !otherwise just load the temperature ob msohrlyb(nstn,1)=buftmp endif endif elseif (buflon.ge.-99.5.and.buflon.le.-91.and.buflat.ge.25.and.buflat.le.50) then !CONUS C call indictb(stnid,buflat,buflon,mesoc,msolatc,msolonc,msomaxc,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypc(nstn)) then msotypc(nstn) = buftyp msodmpc(nstn) = bufdmp msoelvc(nstn) = bufelv endif if (msohrlyc(nstn,1).ne.misreal) then if (buftmp.lt.msohrlya(nstn,1)) then msohrlyc(nstn,1)=buftmp endif else !otherwise just load the temperature ob msohrlyc(nstn,1)=buftmp endif endif elseif (buflon.ge.-118.5.and.buflon.le.-99.5.and.buflat.ge.25.and.buflat.le.50) then !CONUS D call indictb(stnid,buflat,buflon,mesod,msolatd,msolond,msomaxd,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypd(nstn)) then msotypd(nstn) = buftyp msodmpd(nstn) = bufdmp msoelvd(nstn) = bufelv endif if (msohrlyd(nstn,1).ne.misreal) then if (buftmp.lt.msohrlyd(nstn,1)) then msohrlyd(nstn,1)=buftmp endif else !otherwise just load the temperature ob msohrlyd(nstn,1)=buftmp endif endif elseif (((buflon.ge.-70.and.buflon.le.0).or.buflon.le.-118.5).and.buflat.ge.25.and.buflat.le.50) then !CONUS E call indictb(stnid,buflat,buflon,mesoe,msolate,msolone,msomaxe,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotype(nstn)) then msotype(nstn) = buftyp msodmpe(nstn) = bufdmp msoelve(nstn) = bufelv endif if (msohrlye(nstn,1).ne.misreal) then if (buftmp.lt.msohrlye(nstn,1)) then msohrlye(nstn,1)=buftmp endif else !otherwise just load the temperature ob msohrlye(nstn,1)=buftmp endif endif else !CONUS F call indictb(stnid,buflat,buflon,mesof,msolatf,msolonf,msomaxf,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypf(nstn)) then msotypf(nstn) = buftyp msodmpf(nstn) = bufdmp msoelvf(nstn) = bufelv endif if (msohrlyf(nstn,1).ne.misreal) then if (buftmp.lt.msohrlyf(nstn,1)) then msohrlyf(nstn,1)=buftmp endif else !otherwise just load the temperature ob msohrlyf(nstn,1)=buftmp endif endif !if in dictionary endif endif !else !no warnings for mesonets, too many don't report temperature ! print*, "WARNING: Invalid mesonet temperature report (hour 1) (stnid,prov,type,temp,bufqc):",stnid,prov,buftyp,buftmp,bufqc endif endif enddo !SHIP PROCEDURE: Like mesonet, but with only one dictionary! else if (ctyp.eq."SFCSHP ") then do while (ireadsb(25).eq.0) vars(:)=misreal call ufbint(25,vars,9,1,nlevs,"SID PRVSTG SPRVSTG YOB XOB RPT T29 TYP ELV")! TOB TQM") if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevs error line 751 hour,nlevs:",i,nlevs stop endif write(stnid,'(a8)') vars(1) write(prov,'(a8)') vars(2) write(sprov,'(a8)') vars(3) buftyp=vars(8) bufdmp=vars(7) buflat=vars(4) if (vars(5).gt.180) then buflon=vars(5)-360 else buflon=vars(5) endif buftim=vars(6) bufelv=vars(9) if (buftyp.eq.180.or.buftyp.eq.183.or.buftyp.eq.194) then tpc(:,:)=misreal call ufbevn(25,tpc,1,255,20,nlevs,'TPC') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 770. Nlevs=",nlevs stop endif tobs(:,:,:)=misreal call ufbevn(25,tobs,2,255,20,nlevs,'TOB TQM') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 776. Nlevs=",nlevs stop endif !initialize buftmp and bufqc to prevent runtime error later buftmp=misreal bufqc=misreal do k=1,nlevs do j=1,20 if (tpc(k,j)==1) then buftmp=tobs(1,k,j) bufqc=tobs(2,k,j) endif if (tpc(k,j)>1000) exit enddo enddo if (buftmp.gt.1000) buftmp = misreal if (buftmp.lt.1000) then !if ((buftyp.ge.180.and.buftyp.le.199).and.(bufdmp.eq.531.or.bufdmp.eq.532.or.bufdmp.eq.561.or.bufdmp.eq.562).and.bufqc.le.3) then !if valid ship/buoy ob by type, dump, qc mark if (bufqc.le.3) then nstn=0 if (buflat.ge.45.and.(buflon.le.-125.or.buflon.ge.155)) then call indictb(stnid,buflat,buflon,akships,akshplat,akshplon,akshpmax,nstn) if (nstn.gt.0) then if (buftyp.ne.akshptyp(nstn)) then akshptyp(nstn) = buftyp akshpdmp(nstn) = bufdmp akshpelv(nstn) = bufelv endif !if more than one ob/hour, take max ob from that hour, just overwrite the old ob. if (akprvshp(nstn).eq."MISSING ") then akprvshp(nstn)=prov aksprvshp(nstn)=sprov endif if (akshphrly(nstn,1).ne.misreal) then if (buftmp.lt.akshphrly(nstn,1)) then akshphrly(nstn,1)=buftmp endif else !otherwise just load the temperature ob akshphrly(nstn,1)=buftmp endif endif else call indictb(stnid,buflat,buflon,ships,shplat,shplon,shpmax,nstn) if (nstn.gt.0) then if (buftyp.ne.shptyp(nstn)) then shptyp(nstn) = buftyp shpdmp(nstn) = bufdmp shpelv(nstn) = bufelv endif !if more than one ob/hour, take max ob from that hour, just overwrite the old ob. if (prvshp(nstn).eq."SFCSHP ") then prvshp(nstn)=prov sprvshp(nstn)=sprov endif if (shphrly(nstn,1).ne.misreal) then if (buftmp.lt.shphrly(nstn,1)) then shphrly(nstn,1)=buftmp endif else !otherwise just load the temperature ob akshphrly(nstn,1)=buftmp endif endif endif endif !no warning for sfcshp obs since many don't report temperature !else ! print*, "WARNING: Invalid SFCSHP temperature report (hour 1) (stnid,type,temp,qc):",stnid,buftyp,buftmp,bufqc endif endif enddo endif enddo call closbf(25) print*, "Read in first prepbufr file!" allhours: do i=2,25 !get hour number and open, read through the prepbufr file for that hour bufnum=i+24 write(hchar,"(I2.2)") hours(i) !pull the obs from the prepbufr file call openbf(bufnum,"IN",bufnum) call datelen(10) call ufbqcd(bufnum,'VIRTMP',vtcd) print*, "Hour,vtcd=",hours(i),vtcd do while (ireadmg(bufnum,ctyp,idate).eq.0) if (ctyp.eq."ADPSFC ") then do while(ireadsb(bufnum).eq.0) vars(:)=misreal call ufbint (bufnum,vars,9,1,nlevs,"SID PRVSTG SPRVSTG YOB XOB RPT T29 TYP ELV")! TOB TQM") if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevs error line 862 hour,nlevs:",i,nlevs stop endif write(stnid,'(a8)') vars(1) write(prov,'(a8)') vars(2) write(sprov,'(a8)') vars(3) buftyp=vars(8) bufdmp=vars(7) buflat=vars(4) if (vars(5).gt.180) then buflon=vars(5)-360 else buflon=vars(5) endif buftim=vars(6) bufelv=vars(9) if ((buftyp.eq.187.or.buftyp.eq.183.or.buftyp.eq.193.or.buftyp.eq.192).and.(bufdmp.eq.512.or.bufdmp.eq.511)) then bufmax=misreal !initialize bufmax !check for 6-hour maxes only if we are in an hour which we know a priori will have a 6-hour max if (hours(i).eq.12.or.hours(i).eq.18.or.hours(i).eq.0.or.hours(i).eq.6) then maxs(:,:)=misreal call ufbint (bufnum,maxs,4,255,nlevs,".DTHMXTM MXTM .DTHMITM MITM") if (nlevs.ne.1.and.nlevs.ne.0) then print*, "WARNING: maxT warning line 887 hour,nlevs:",hours(i),nlevs !stop endif !get 6-hour minimum temperature if available maxdo: do k=1,255 !if dthmxtm is missing value, that means we have no max, so set to missing if (ibfms(maxs(3,k).eq.1)) then bufmax=misreal exit maxdo !if dthmxtm is 6, that means we have our proper value. Collect and store elseif (maxs(3,k).eq.6) then if (ibfms(maxs(4,k)).eq.0) then bufmax=maxs(4,k)-273.15 else bufmax=misreal endif exit maxdo endif enddo maxdo endif tpc(:,:)=misreal call ufbevn(bufnum,tpc,1,255,20,nlevs,'TPC') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 910. Nlevs=",nlevs stop endif tobs(:,:,:)=misreal call ufbevn(bufnum,tobs,2,255,20,nlevs,'TOB TQM') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 916. Nlevs=",nlevs stop endif !initialize buftmp and bufqc to prevent runtime error later buftmp=misreal bufqc=misreal getv: do k=1,nlevs do j=1,20 if (tpc(k,j)==1) then buftmp=tobs(1,k,j) bufqc=tobs(2,k,j) endif if (tpc(k,j)>1000) exit getv enddo enddo getv if (buftmp.gt.1000) buftmp=misreal !if buftmp is missing 10**10 goes to misreal if (buftmp.lt.1000.or.bufmax.lt.1000) then !ensure temp obs are valid if (bufqc.le.3) then nstn=0 if (buflat.ge.45.and.(buflon.le.-125.or.buflon.ge.155)) then call indictb(stnid,buflat,buflon,akmetars,akmtrlat,akmtrlon,akmtrmax,nstn) if (nstn.gt.0) then if (buftyp.ne.akmtrtyp(nstn)) then akmtrtyp(nstn) = buftyp akmtrdmp(nstn) = bufdmp akmtrelv(nstn) = bufelv endif if (akprvmtr(nstn).eq."MISSING ") then akprvmtr(nstn)=prov aksprvmtr(nstn)=sprov endif !record 6hour max if available if (bufmax.ne.misreal) then if (hours(i).eq.0) then if (akmtrmaxs(nstn,1).eq.misreal) then akmtrmaxs(nstn,1)=bufmax!-273.15 elseif (bufmax.ne.akmtrmaxs(nstn,1)) then print*, "WARNING: Multiple mint (12Z) (skipping mint1):",stnid,bufmax,akmtrmaxs(nstn,1) !stop akmtrmaxs(nstn,1)=misreal else print*, "WARNING: Multiple same min (12Z):",stnid,bufmax endif elseif (hours(i).eq.6) then if (akmtrmaxs(nstn,2).eq.misreal) then akmtrmaxs(nstn,2)=bufmax!-273.15 elseif (bufmax.ne.akmtrmaxs(nstn,2)) then print*, "WARNING: Multiple mint (18Z) (skipping min2):",stnid,bufmax,akmtrmaxs(nstn,2) !stop akmtrmaxs(nstn,2)=misreal else print*, "WARNING: Multiple same min (18Z):",stnid,bufmax endif elseif (hours(i).eq.12) then if (akmtrmaxs(nstn,3).eq.misreal) then akmtrmaxs(nstn,3)=bufmax!-273.15 elseif (bufmax.ne.akmtrmaxs(nstn,3)) then print*, "WARNING: Multiple mint (00Z) (skipping min3):",stnid,bufmax,akmtrmaxs(nstn,3) !stop akmtrmaxs(nstn,3)=misreal else print*, "WARNING: Multiple same min (00Z):",stnid,bufmax endif elseif (hours(i).eq.18) then if (akmtrmaxs(nstn,4).eq.misreal) then akmtrmaxs(nstn,4)=bufmax!-273.15 elseif (bufmax.ne.akmtrmaxs(nstn,4)) then print*, "WARNING: Multiple mint (06Z) (skipping min4):",stnid,bufmax,akmtrmaxs(nstn,4) !stop akmtrmaxs(nstn,4)=misreal else print*, "WARNING: Multiple same min (06Z):",stnid,bufmax endif endif endif if (akmtrhrly(nstn,i).ne.misreal) then if (buftmp.gt.akmtrhrly(nstn,i)) then akmtrhrly(nstn,i)=buftmp endif else akmtrhrly(nstn,i)=buftmp endif endif else call indictb(stnid,buflat,buflon,metars,mtrlat,mtrlon,mtrmax,nstn) if (nstn.gt.0) then if (buftyp.ne.mtrtyp(nstn)) then mtrtyp(nstn) = buftyp mtrdmp(nstn) = bufdmp mtrelv(nstn) = bufelv endif !record 6hour max if available if (prvmtr(nstn).eq."ADPSFC ") then prvmtr(nstn)=prov sprvmtr(nstn)=sprov endif if (bufmax.ne.misreal) then if (hours(i).eq.0) then if (mtrmaxs(nstn,1).eq.misreal) then mtrmaxs(nstn,1)=bufmax!-273.15 elseif (bufmax.ne.mtrmaxs(nstn,1)) then print*, "WARNING: Multiple mint (12Z) (skipping mint1):",stnid,bufmax,mtrmaxs(nstn,1) !stop mtrmaxs(nstn,1)=misreal else print*, "WARNING: Multiple same min (12Z):",stnid,bufmax endif elseif (hours(i).eq.6) then if (mtrmaxs(nstn,2).eq.misreal) then mtrmaxs(nstn,2)=bufmax!-273.15 elseif (bufmax.ne.mtrmaxs(nstn,2)) then print*, "WARNING: Multiple mint (18Z) (skipping min2):",stnid,bufmax,mtrmaxs(nstn,2) !stop mtrmaxs(nstn,2)=misreal else print*, "WARNING: Multiple same min (18Z):",stnid,bufmax endif elseif (hours(i).eq.12) then if (mtrmaxs(nstn,3).eq.misreal) then mtrmaxs(nstn,3)=bufmax!-273.15 elseif (bufmax.ne.mtrmaxs(nstn,3)) then print*, "WARNING: Multiple mint (00Z) (skipping min3):",stnid,bufmax,mtrmaxs(nstn,3) !stop mtrmaxs(nstn,3)=misreal else print*, "WARNING: Multiple same min (00Z):",stnid,bufmax endif elseif (hours(i).eq.18) then if (mtrmaxs(nstn,4).eq.misreal) then mtrmaxs(nstn,4)=bufmax!-273.15 elseif (bufmax.ne.mtrmaxs(nstn,4)) then print*, "WARNING: Multiple mint (06Z) (skipping min4):",stnid,bufmax,mtrmaxs(nstn,4) !stop mtrmaxs(nstn,4)=misreal else print*, "WARNING: Multiple same min (06Z):",stnid,bufmax endif endif endif !conflict situation: look at tob and compare to hourly ob currently in. Stick with the larger one if (mtrhrly(nstn,i).ne.misreal) then if (buftmp.lt.mtrhrly(nstn,i)) then mtrhrly(nstn,i)=buftmp endif else mtrhrly(nstn,i)=buftmp endif endif endif endif !else ! print*, "WARNING: Invalid ADPSFC temperature (and maxT) report hour (stnid,type,hour,temp,maxT,qc):",stnid,buftyp,hours(i),buftmp,bufmax,bufqc endif endif enddo elseif (ctyp.eq."MSONET ") then do while(ireadsb(bufnum).eq.0) vars(:)=misreal call ufbint (bufnum,vars,9,1,nlevs,"SID PRVSTG SPRVSTG YOB XOB RPT T29 TYP ELV") if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevs error line 1072 hour,nlevs:",i,nlevs stop endif write(stnid,'(a8)') vars(1) write(prov,'(a8)') vars(2) write(sprov,'(a8)') vars(3) !make sure trailing x and a are removed from stationIDS because they just cause confusion if (stnid(8:8).eq."a".or.stnid(8:8).eq."x") then stnid(8:8)=" " endif buftyp=vars(8) bufdmp=vars(7) buflat=vars(4) if (vars(5).gt.180) then buflon=vars(5)-360 else buflon=vars(5) endif buftim=vars(6) bufelv=vars(9) if ((buftyp.eq.188.or.buftyp.eq.195).and.bufdmp.eq.540) then tpc(:,:)=misreal call ufbevn(bufnum,tpc,1,255,20,nlevs,'TPC') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 1098. Nlevs=",nlevs stop endif tobs(:,:,:)=misreal call ufbevn(bufnum,tobs,2,255,20,nlevs,'TOB TQM') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 1104. Nlevs=",nlevs stop endif !initialize buftmp and bufqc to prevent runtime error later buftmp=misreal bufqc=misreal do k=1,nlevs do j=1,20 if (tpc(k,j)==1) then buftmp=tobs(1,k,j) bufqc=tobs(2,k,j) endif if (tpc(k,j)>1000) exit enddo enddo if (buftmp.lt.1000) then if (bufqc.le.3) then nstn=0 if (buflon.ge.149.and.buflon.le.162.and.buflat.ge.12.and.buflat.le.17) then ! GUAM call indictb(stnid,buflat,buflon,mesogu,msolatgu,msolongu,msomaxgu,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypgu(nstn)) then msotypgu(nstn) = buftyp msodmpgu(nstn) = bufdmp msoelvgu(nstn) = bufelv endif if (msohrlygu(nstn,i).ne.misreal) then if (buftmp.lt.msohrlygu(nstn,i)) then !.and.buftmp.ne.misreal) then msohrlygu(nstn,i)=buftmp endif else !otherwise just load the temperature ob msohrlygu(nstn,i)=buftmp endif endif !if in dictionary elseif (buflon.ge.-149.and.buflon.le.-143.and.buflat.ge.17.and.buflat.le.24) then !HAWAII call indictb(stnid,buflat,buflon,mesohi,msolathi,msolonhi,msomaxhi,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotyphi(nstn)) then msotyphi(nstn) = buftyp msodmphi(nstn) = bufdmp msoelvhi(nstn) = bufelv endif if (msohrlyhi(nstn,i).ne.misreal) then if (buftmp.lt.msohrlyhi(nstn,i)) then !.and.buftmp.ne.misreal) then msohrlyhi(nstn,i)=buftmp endif else !otherwise just load the temperature ob msohrlyhi(nstn,i)=buftmp endif endif !if in dictionary elseif (buflon.ge.-69.and.buflon.le.-61.and.buflat.ge.16.and.buflat.le.20) then !PUERTO RICO call indictb(stnid,buflat,buflon,mesopr,msolatpr,msolonpr,msomaxpr,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotyppr(nstn)) then msotyppr(nstn) = buftyp msodmppr(nstn) = bufdmp msoelvpr(nstn) = bufelv endif if (msohrlypr(nstn,i).ne.misreal) then if (buftmp.lt.msohrlypr(nstn,i)) then !.and.buftmp.ne.misreal) then msohrlypr(nstn,i)=buftmp endif else !otherwise just load the temperature ob msohrlypr(nstn,i)=buftmp endif endif !if in dictionary elseif (buflat.ge.45.and.(buflon.le.-125.or.buflon.ge.155)) then !ALASKA call indictb(stnid,buflat,buflon,mesoak,msolatak,msolonak,msomaxak,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypak(nstn)) then msotypak(nstn) = buftyp msodmpak(nstn) = bufdmp msoelvak(nstn) = bufelv endif if (msohrlyak(nstn,i).ne.misreal) then if (buftmp.lt.msohrlyak(nstn,i)) then !.and.buftmp.ne.misreal) then msohrlyak(nstn,i)=buftmp endif else !otherwise just load the temperature ob msohrlyak(nstn,i)=buftmp endif endif !if in dictionary elseif (buflon.ge.-82.5.and.buflon.le.-70.and.buflat.ge.25.and.buflat.le.50) then !CONUS A call indictb(stnid,buflat,buflon,mesoa,msolata,msolona,msomaxa,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypa(nstn)) then msotypa(nstn) = buftyp msodmpa(nstn) = bufdmp msoelva(nstn) = bufelv endif if (msohrlya(nstn,i).ne.misreal) then if (buftmp.lt.msohrlya(nstn,i)) then !.and.buftmp.ne.misreal) then msohrlya(nstn,i)=buftmp endif else !otherwise just load the temperature ob msohrlya(nstn,i)=buftmp endif endif !if in dictionary elseif (buflon.ge.-91.and.buflon.le.-82.5.and.buflat.ge.25.and.buflat.le.50) then !CONUS B call indictb(stnid,buflat,buflon,mesob,msolatb,msolonb,msomaxb,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypb(nstn)) then msotypb(nstn) = buftyp msodmpb(nstn) = bufdmp msoelvb(nstn) = bufelv endif if (msohrlyb(nstn,i).ne.misreal) then if (buftmp.lt.msohrlyb(nstn,i)) then msohrlyb(nstn,i)=buftmp endif else !otherwise just load the temperature ob msohrlyb(nstn,i)=buftmp endif endif !if in dictionary elseif (buflon.ge.-99.5.and.buflon.le.-91.and.buflat.ge.25.and.buflat.le.50) then !CONUS C call indictb(stnid,buflat,buflon,mesoc,msolatc,msolonc,msomaxc,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypc(nstn)) then msotypc(nstn) = buftyp msodmpc(nstn) = bufdmp msoelvc(nstn) = bufelv endif if (msohrlyc(nstn,i).ne.misreal) then if (buftmp.lt.msohrlyc(nstn,i)) then msohrlyc(nstn,i)=buftmp endif else !otherwise just load the temperature ob msohrlyc(nstn,i)=buftmp endif endif !if in dictionary elseif (buflon.ge.-118.5.and.buflon.le.-99.5.and.buflat.ge.25.and.buflat.le.50) then !CONUS D call indictb(stnid,buflat,buflon,mesod,msolatd,msolond,msomaxd,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypd(nstn)) then msotypd(nstn) = buftyp msodmpd(nstn) = bufdmp msoelvd(nstn) = bufelv endif if (msohrlyd(nstn,i).ne.misreal) then if (buftmp.lt.msohrlyd(nstn,i)) then msohrlyd(nstn,i)=buftmp endif else !otherwise just load the temperature ob msohrlyd(nstn,i)=buftmp endif endif !if in dictionary elseif (((buflon.ge.-70.and.buflon.le.0).or.buflon.le.-118.5).and.buflat.ge.25.and.buflat.le.50) then !CONUS E call indictb(stnid,buflat,buflon,mesoe,msolate,msolone,msomaxe,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotype(nstn)) then msotype(nstn) = buftyp msodmpe(nstn) = bufdmp msoelve(nstn) = bufelv endif if (msohrlye(nstn,i).ne.misreal) then if (buftmp.lt.msohrlye(nstn,i)) then msohrlye(nstn,i)=buftmp endif else !otherwise just load the temperature ob msohrlye(nstn,i)=buftmp endif endif !if in dictionary else !CONUS F call indictb(stnid,buflat,buflon,mesof,msolatf,msolonf,msomaxf,nstn) if (nstn.gt.0) then !if in dictionary if (buftyp.ne.msotypf(nstn)) then msotypf(nstn) = buftyp msodmpf(nstn) = bufdmp msoelvf(nstn) = bufelv endif if (msohrlyf(nstn,i).ne.misreal) then if (buftmp.lt.msohrlyf(nstn,i)) then !.and.buftmp.ne.misreal) then msohrlyf(nstn,i)=buftmp endif else !otherwise just load the temperature ob msohrlyf(nstn,i)=buftmp endif endif !if in dictionary endif endif !no warnings for mesonets - too many don't report temperature !else ! print*, "WARNING: Invalid MSONET temperature report hour (stnid,prov,type,hour,temp,qc):",stnid,prov,buftyp,hours(i),buftmp,bufqc endif endif enddo elseif (ctyp.eq."SFCSHP ") then do while(ireadsb(bufnum).eq.0) vars(:)=misreal call ufbint (bufnum,vars,9,1,nlevs,"SID PRVSTG SPRVSTG YOB XOB RPT T29 TYP ELV")! TOB TQM") if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevs error line 1298 hour,nlevs:",i,nlevs stop endif write(stnid,'(a8)') vars(1) write(prov,'(a8)') vars(2) write(sprov,'(a8)') vars(3) buftyp=vars(8) bufdmp=vars(7) buflat=vars(4) if (vars(5).gt.180) then buflon=vars(5)-360 else buflon=vars(5) endif buftim=vars(6) bufelv=vars(9) if ((buftyp.ge.180.and.buftyp.le.199).and.(bufdmp.eq.531.or.bufdmp.eq.532.or.bufdmp.eq.561.or.bufdmp.eq.562)) then tpc(:,:)=misreal call ufbevn(bufnum,tpc,1,255,20,nlevs,'TPC') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 1318. Nlevs=",nlevs stop endif tobs(:,:,:)=misreal call ufbevn(bufnum,tobs,2,255,20,nlevs,'TOB TQM') if (nlevs.ne.1) then print*, "FATAL ERROR: Nlevels error line 1324. Nlevs=",nlevs stop endif !initialize buftmp and bufqc to prevent runtime error later buftmp=misreal bufqc=misreal do k=1,nlevs do j=1,20 if (tpc(k,j)==1) then buftmp=tobs(1,k,j) bufqc=tobs(2,k,j) endif if (tpc(k,j)>1000) exit enddo enddo if (buftmp.lt.1000) then if (bufqc.le.3) then nstn=0 if (buflat.ge.45.and.(buflon.le.-125.or.buflon.ge.155)) then call indictb(stnid,buflat,buflon,akships,akshplat,akshplon,akshpmax,nstn) if (nstn.gt.0) then if (buftyp.ne.akshptyp(nstn)) then akshptyp(nstn) = buftyp akshpdmp(nstn) = bufdmp akshpelv(nstn) = bufelv endif if (akprvshp(nstn).eq."MISSING ") then akprvshp(nstn)=prov aksprvshp(nstn)=sprov endif if (akshphrly(nstn,i).ne.misreal) then if (buftmp.lt.akshphrly(nstn,i)) then akshphrly(nstn,i)=buftmp endif else !otherwise just load the temperature ob akshphrly(nstn,i)=buftmp endif endif else call indictb(stnid,buflat,buflon,ships,shplat,shplon,shpmax,nstn) if (nstn.gt.0) then if (buftyp.ne.shptyp(nstn)) then shptyp(nstn) = buftyp shpdmp(nstn) = bufdmp shpelv(nstn) = bufelv endif if (prvshp(nstn).eq."SFCSHP ") then prvshp(nstn)=prov sprvshp(nstn)=sprov endif if (shphrly(nstn,i).ne.misreal) then if (buftmp.lt.shphrly(nstn,i)) then shphrly(nstn,i)=buftmp endif else !otherwise just load the temperature ob shphrly(nstn,i)=buftmp endif endif endif endif !no warning message for sfcshp since many do not report temperature !else ! print*, "WARNING: Invalid SFCSHP temperature report hour (stnid,type,hour,temp,qc):",stnid,buftyp,hours(i),buftmp,bufqc endif endif enddo endif enddo close(bufnum) enddo allhours 125 format(A8,2X,6(F6.2,2X)) print*, "Finished reding in obs - now calculating mimT for METARS" print*, "While calcing mins: mtrmax=",mtrmax call maxmin(71,mtrhrly,mtrmaxs,mtrzone,metars,metnames,maxt,max,4,25,1,2,mtrmax) do i=1,mtrmax if (maxt(i).ne.misreal) then write(61,356) metars(i),prvmtr(i),sprvmtr(i),mtrtyp(i),mtrdmp(i),mtrlat(i),mtrlon(i),mtrelv(i),mtrzone(i),maxt(i) endif enddo call maxmin(72,akmtrhrly,akmtrmaxs,akmtrzone,akmetars,namesak,maxt,max,4,25,1,2,akmtrmax) do i=1,akmtrmax if (maxt(i).ne.misreal) then write(61,356) akmetars(i),akprvmtr(i),aksprvmtr(i),akmtrtyp(i),akmtrdmp(i),akmtrlat(i),akmtrlon(i),akmtrelv(i),akmtrzone(i),maxt(i) endif enddo call maxmin(73,shphrly,shpmaxs,shpzone,ships,shpnames,maxt,max,4,25,1,2,shpmax) print*, "While calcing mins: shpmax=",shpmax do i=1,shpmax if (maxt(i).ne.misreal) then write(61,356) ships(i),prvshp(i),sprvshp(i),shptyp(i),shpdmp(i),shplat(i),shplon(i),shpelv(i),shpzone(i),maxt(i) endif enddo call maxmin(74,akshphrly,shpmaxs,akshpzone,akships,namesak,maxt,max,4,25,1,2,akshpmax) do i=1,akshpmax if (maxt(i).ne.misreal) then write(61,356) akships(i),akprvshp(i),aksprvshp(i),akshptyp(i),akshpdmp(i),akshplat(i),akshplon(i),akshpelv(i),akshpzone(i),maxt(i) endif enddo call maxmin(75,msohrlya,msomaxs,msozonea,mesoa,msonames,maxt,max,4,25,1,2,msomaxa) print*, "While calcing mins: msomaxa=",msomaxa do i=1,msomaxa if (maxt(i).ne.misreal) then write(61,356) mesoa(i),prva(i),sprva(i),msotypa(i),msodmpa(i),msolata(i),msolona(i),msoelva(i),msozonea(i),maxt(i) endif enddo print*, "While calcing mins: msomaxb=",msomaxb call maxmin(76,msohrlyb,msomaxs,msozoneb,mesob,msonames,maxt,max,4,25,1,2,msomaxb) do i=1,msomaxb if (maxt(i).ne.misreal) then write(61,356) mesob(i),prvb(i),sprvb(i),msotypb(i),msodmpb(i),msolatb(i),msolonb(i),msoelvb(i),msozoneb(i),maxt(i) endif enddo print*, "While calcing mins: msomaxc=",msomaxc call maxmin(77,msohrlyc,msomaxs,msozonec,mesoc,msonames,maxt,max,4,25,1,2,msomaxc) do i=1,msomaxc if (maxt(i).ne.misreal) then write(61,356) mesoc(i),prvc(i),sprvc(i),msotypc(i),msodmpc(i),msolatc(i),msolonc(i),msoelvc(i),msozonec(i),maxt(i) endif enddo print*, "While calcing mins: msomaxd=",msomaxd call maxmin(78,msohrlyd,msomaxs,msozoned,mesod,msonames,maxt,max,4,25,1,2,msomaxd) do i=1,msomaxd if (maxt(i).ne.misreal) then write(61,356) mesod(i),prvd(i),sprvd(i),msotypd(i),msodmpd(i),msolatd(i),msolond(i),msoelvd(i),msozoned(i),maxt(i) endif enddo print*, "While calcing mins: msomaxe=",msomaxe call maxmin(79,msohrlye,msomaxs,msozonee,mesoe,msonames,maxt,max,4,25,1,2,msomaxe) do i=1,msomaxe if (maxt(i).ne.misreal) then write(61,356) mesoe(i),prve(i),sprve(i),msotype(i),msodmpe(i),msolate(i),msolone(i),msoelve(i),msozonee(i),maxt(i) endif enddo print*, "While calcing mins: msomaxf=",msomaxf call maxmin(80,msohrlyf,msomaxs,msozonef,mesof,msonames,maxt,max,4,25,1,2,msomaxf) do i=1,msomaxf if (maxt(i).ne.misreal) then write(61,356) mesof(i),prvf(i),sprvf(i),msotypf(i),msodmpf(i),msolatf(i),msolonf(i),msoelvf(i),msozonef(i),maxt(i) endif enddo print*, "While calcing mins: msomaxpr=",msomaxpr call maxmin(81,msohrlypr,msomaxs,msozonepr,mesopr,msonames,maxt,max,4,25,1,2,msomaxpr) do i=1,msomaxpr if (maxt(i).ne.misreal) then write(61,356) mesopr(i),prvpr(i),sprvpr(i),msotyppr(i),msodmppr(i),msolatpr(i),msolonpr(i),msoelvpr(i),msozonepr(i),maxt(i) endif enddo print*, "While calcing mins: msomaxhi=",msomaxhi call maxmin(82,msohrlyhi,msomaxs,msozonehi,mesohi,msonames,maxt,max,4,25,1,2,msomaxhi) do i=1,msomaxhi if (maxt(i).ne.misreal) then write(61,356) mesohi(i),prvhi(i),sprvhi(i),msotyphi(i),msodmphi(i),msolathi(i),msolonhi(i),msoelvhi(i),msozonehi(i),maxt(i) endif enddo print*, "While calcing mins: msomaxgu=",msomaxgu call maxmin(83,msohrlygu,msomaxs,msozonegu,mesogu,msonames,maxt,max,4,25,1,2,msomaxgu) do i=1,msomaxgu if (maxt(i).ne.misreal) then write(61,356) mesogu(i),prvgu(i),sprvgu(i),msotypgu(i),msodmpgu(i),msolatgu(i),msolongu(i),msoelvgu(i),msozonegu(i),maxt(i) endif enddo print*, "While calcing mins: msomaxak=",msomaxak call maxmin(84,msohrlyak,msomaxs,msozoneak,mesoak,namesak,maxt,max,4,25,1,2,msomaxak) do i=1,msomaxak if (maxt(i).ne.misreal) then write(61,356) mesoak(i),prvak(i),sprvak(i),msotypak(i),msodmpak(i),msolatak(i),msolonak(i),msoelvak(i),msozoneak(i),maxt(i) endif enddo close(61) 315 format (A8,3X,3(F8.3,3X),F8.3,3X,F8.3,3X,I3,3X,5(F8.3,3X)) 355 format (2(A8,3X),3(F8.3,3X),F8.3,3X,F8.3,3X,I3,3X,F8.3) 356 format (3(A8,3X),3(F8.3,3X),F8.3,3X,F8.3,3X,I3,3X,F8.3) 456 format (A8,3X,3(F7.2,3X)) 457 format (A8,3X,2(F7.2),3X,A8) 458 format (A8,3X,A8,3X,5(F7.2,3X),I3) 459 format (A8,3X,5(F7.2,3X),I3) end program calcmint