program degrib2model_driver c implicit none c integer nx,ny,nz,kpds(200),datsav(9),n !Degrib grid dimensions integer filemax,noct,I parameter (filemax=8) integer octnum(filemax),jm1,jp1,ip1,ip2,ip3,ip4,latcov,loncov, + readin c character(LEN=255):: gribfile,outdir c real esat,es c common /estab/esat(15000:45000),es(15000:45000) c_______________________________________________________________________________ c c *** Initialize tables. c call es_ini c c *** Read file names. c Cmp Cmp all in same file....read in desired octant numbers???? Cmp read(5,'(a)') gribfile read(5,'(a)') outdir do I=1,filemax+1 C write(6,*) 'reading for I= ', I read(5,23) readin if (readin .eq. -9) then noct=I-1 goto 1103 endif octnum(I)=readin enddo 23 format(I2) 1103 continue c write(6,*) 'want to work with the following octants: ' do I=1,noct write(6,*) octnum(I) enddo n=index(gribfile,' ')-1 IP1=0 IP2=0 IP3=0 IP4=0 JP1=0 JM1=0 do I=1,noct C CALL GET_GDS(gribfile(1:n),octnum(I),datsav,kpds) if (octnum(I) .le. 40) JP1=1 if (octnum(I) .gt. 40) JM1=1 if (octnum(I) .eq. 37 .or. octnum(I) .eq. 41) IP1=1 if (octnum(I) .eq. 38 .or. octnum(I) .eq. 42) IP2=1 if (octnum(I) .eq. 39 .or. octnum(I) .eq. 43) IP3=1 if (octnum(I) .eq. 40 .or. octnum(I) .eq. 44) IP4=1 C datsav(3) SW lat C datsav(4) SW lon enddo latcov=JM1+JP1 loncov=(IP1+IP2+IP3+IP4) write(6,*) 'octants selected are ', loncov, ' in longitude ' write(6,*) 'and ', latcov, ' in latitude' C write(6,*) 'calculated NX x NY is : ', loncov*73, latcov*73 nz=11 C nx=loncov*73 C ny=latcov*73 nx=loncov*73-(loncov-1) ny=latcov*73-(latcov-1) write(6,*) 'nx= ny= ', nx,ny C stop call degrib2model(gribfile,outdir,octnum,noct,nx,ny,nz) c end c c=============================================================================== c subroutine degrib2model(gribfile,outdir,octnum,noct,nx,ny,nz) c implicit none c real cp,kappa,dlat,dlon,sum,avg parameter(cp=1004.,kappa=287./1004.) c integer nx,ny,nz,i,j,k,l,n,it,ip,jp,nsfcfld . ,iyear,imonth,iday,ijulday,julday_laps,ifcsthr . ,count,jinput,noct,octnum(noct),swlatsv,swlonsv integer swlat(noct),swlon(noct),III,datsav(9),II integer londif,latdif,istart,jstart,ival,jval,kpds(200) integer KGDS(200) c real ht(nx,ny,nz),htin(73,73,nz) !Isobaric heights (m) . ,tp(nx,ny,nz),tpin(73,73,nz) !Isobaric temps (K) . ,th(nx,ny,nz),thin(73,73,nz) !Isobaric theta (K) . ,uw(nx,ny,nz),uwin(73,73,nz) !Isobaric u-wind (m/s) . ,vw(nx,ny,nz),vwin(73,73,nz) !Isobaric v-wind (m/s) . ,rh(nx,ny,nz),rhin(73,73,nz) !Isobaric rh,mr (%,kg/kg) . ,ex(nx,ny,nz) !Isobaric Exner . ,pr(nz),pri(nz) !Isobaric pressures (mb) . ,lat1,lat2,lon0,sw(2),ne(2) . ,xe,mrsat,esat,es . ,avnlevs(11) . ,temp(nx,ny) . ,slp(nx,ny,12) c character(LEN=255):: gribfile,outdir,outfile,gdsfile character(LEN=2):: gproj character(LEN=11):: atime character(LEN=8):: model data AVNLEVS/1000,850,700,500,400,300,250,200,150,100,70/ c common /estab/esat(15000:45000),es(15000:45000) c___________________________________________________________________________ Cmp Loop over whole subroutine n=index(gribfile,' ')-1 do III=1,noct CALL GET_GDS(gribfile(1:n),octnum(III),datsav,kpds) swlat(III)=datsav(3) swlon(III)=datsav(4) enddo do I=1,noct swlatsv=99999 if (swlat(I) .eq. -90000) then swlatsv=-90000 goto 97 endif enddo C lat not 90S, so must be 0 swlatsv=0 97 continue Cmp logic here is BAD! C C Need to be able to wrap from Grid 40 (120W at west) to grid C 37 (30W at west). C do I=1,noct if (swlon(I) .eq. -30000) then swlonsv=-30000 goto 98 endif enddo do I=1,noct if (swlon(I) .eq. 60000) then swlonsv=60000 goto 98 endif enddo do I=1,noct if (swlon(I) .eq. 150000) then swlonsv=150000 goto 98 endif enddo swlonsv=-120000 98 continue C C C at this point have initial guess of "westernmost" point, C but not sure if it is the "left" most point on the final grid if (swlonsv .eq. -30000) then do III=1,noct if (swlon(III) .eq. -120000) then write(6,*) 'special case...to 120W' swlonsv=-120000 do II=1,noct if (swlon(II) .eq. 150000) then write(6,*) 'special case...to 150E' swlonsv=150000 endif enddo endif enddo endif do III=1,noct jinput=octnum(III) c c *** Fill pressure levels. c C write(6,*) 'using hardwired stuff! ', nx,ny,nz do k=1,nz pr(k)=AVNLEVS(k) C write(6,*) 'pressure value of ', pr(k) enddo c c *** Read in degrib data. c Cmp Cmp "thinned" WAFS grids have 73 rows, and when filled are 73 X 73 Cmp However, in thinned status each octant has 3447 gridpoints Cmp n=index(gribfile,' ')-1 call read_degrib(gribfile(1:n) . ,3447,nz,pr,htin,tpin,rhin,uwin,vwin,jinput,atime) C C Find destination londif=swlon(III)-swlonsv latdif=swlat(III)-swlatsv if (londif .lt. 0) londif=londif+360000 if ( londif .eq. 0) ISTART=1 C if ( londif.eq.90000) ISTART=74 if ( londif.eq.90000) ISTART=73 C if ( londif.eq.180000) ISTART=147 if ( londif.eq.180000) ISTART=145 C if ( londif.eq.270000) ISTART=146+73 if ( londif.eq.270000) ISTART=145+72 if ( latdif .eq. 0) JSTART=1 C if ( latdif .eq. 90000) JSTART=74 if ( latdif .eq. 90000) JSTART=73 C C C Put the "in" arrays into their final places for the isobaric fields C do K=1,NZ do J=JSTART,JSTART+73-1 do I=ISTART,ISTART+73-1 ival=(I-istart)+1 jval=(J-jstart)+1 ht(I,J,K)=htin(ival,jval,K) tp(I,J,K)=tpin(ival,jval,K) rh(I,J,K)=rhin(ival,jval,K) uw(I,J,K)=uwin(ival,jval,K) vw(I,J,K)=vwin(ival,jval,K) enddo enddo enddo c c *** Check for any missing data. c if (III .eq. NOCT) THEN do k=1,nz if (ht(1,1,k) .eq. -99999.) then print *,'Height data missing at level: ',pr(k) stop elseif (tp(1,1,k) .eq. -99999.) then print *,'Temperature data missing at level: ',pr(k) stop elseif (rh(1,1,k) .eq. -99999.) then print *,'RH data missing at level: ',pr(k) print *,'Calling RH patch.' call rh_fix(nx,ny,nz,rh,pr) elseif (uw(1,1,k) .eq. -99999.) then print *,'U-wind data missing at level: ',pr(k) stop elseif (vw(1,1,k) .eq. -99999.) then print *,'V-wind data missing at level: ',pr(k) stop endif enddo c c *** Convert 3d temp to theta. c *** Compute Exner function. c *** Convert 3d rh to mr. c do k=1,nz pri(k)=1./pr(k) enddo c do k=1,nz do j=1,ny do i=1,nx th(i,j,k)=tp(i,j,k)*(1000.*pri(k))**kappa ex(i,j,k)=cp*tp(i,j,k)/th(i,j,k) it=tp(i,j,k)*100 it=min(45000,max(15000,it)) xe=esat(it) mrsat=0.00622*xe/(pr(k)-xe) rh(i,j,k)=rh(i,j,k)*mrsat enddo enddo enddo Cmp Cmp for some reason AVN data appears to be "upside down" in the N-S Cmp sense. Flip the arrays of data. Cmp C do k=1,nz C do j=1,ny C do i=1,nx C temp(i,j)=th(i,j,k) C enddo C enddo C do J=1,ny C do i=1,nx C th(i,j,k)=temp(i,ny-J+1) C enddo C enddo C====================================== C do j=1,ny C do i=1,nx C temp(i,j)=uw(i,j,k) C enddo C enddo C do J=1,ny C do i=1,nx C uw(i,j,k)=temp(i,ny-J+1) C enddo C enddo C====================================== C do j=1,ny C do i=1,nx C temp(i,j)=vw(i,j,k) C enddo C enddo C do J=1,ny C do i=1,nx C vw(i,j,k)=temp(i,ny-J+1) C enddo C enddo C====================================== C do j=1,ny C do i=1,nx C temp(i,j)=ht(i,j,k) C enddo C enddo C if (k.eq.1) write(6,*) 'flipping isobaric data - ' C do J=1,ny C do i=1,nx C ht(i,j,k)=temp(i,ny-J+1) C enddo C enddo C====================================== C do j=1,ny C do i=1,nx C temp(i,j)=rh(i,j,k) C enddo C enddo C do J=1,ny C do i=1,nx C rh(i,j,k)=temp(i,ny-J+1) C enddo C enddo C====================================== C do j=1,ny C do i=1,nx C temp(i,j)=ex(i,j,k) C enddo C enddo C do J=1,ny C do i=1,nx C ex(i,j,k)=temp(i,ny-J+1) C enddo C enddo C----------------------------------------- C enddo c C print *,'ua :',ht(nx/2,ny/2,nz),th(nx/2,ny/2,nz),rh(nx/2,ny/2,nz) C . ,uw(nx/2,ny/2,nz),vw(nx/2,ny/2,nz) C print *,'ua :',ht(nx/2,ny/2,1),th(nx/2,ny/2,1),rh(nx/2,ny/2,1) C . ,uw(nx/2,ny/2,1),vw(nx/2,ny/2,1) c c *** Create output file name. c n=index(outdir,' ')-1 model='ETA_wafs' outfile=outdir(1:n)//'/'//atime//'.'//model n=index(outfile,' ')-1 print *,model,' data ---> ',outfile(1:n) open(1,file=outfile(1:n),status='unknown',form='unformatted') c c *** Write header stuff. c ip=1 jp=1 Cmp nsfcfld will be 12 although 4 of these are set to missing for AVN... nsfcfld=12 Cmp is the AVN 'LL' or 'CE'?????? gproj='LL' ! write(1) nx,ny,nz,nx,ny,ip,jp,nsfcfld,gproj write(1) nz C new stuff for GDS file C n=index(gribfile,' ')-1 write(6,*) 'calling get_fullgds with ', + gribfile(1:n) call get_fullgds(gribfile(1:n),nx,ny,kgds) write(6,*) 'back from get_fullgds' Cmoved n=index(outdir,' ')-1 C gdsfile=outdir(1:n)//'/'//'gdsinfo.'//model C C n=index(gdsfile,' ')-1 C write(6,*) 'GDS written to ', gdsfile(1:n) C open(unit=14,file=gdsfile(1:n),form='unformatted', C + access='sequential') C C write(6,*) 'writing kgds ', (kgds(I),i=1,14) C write(14) kgds C Cmoved close(14) C C end new stuff C C C this header info is octant-dependent...also needs work for patched C together octants C C lat1=-90.0 C lat2=0.0 C lon0=0.0 C sw(1)=-90 C sw(2)=0. C ne(1)=90. C ne(2)=359. lat1=swlatsv/1000. lat2=0. lon0=swlonsv/1000. sw(1)=swlatsv/1000. sw(2)=swlonsv/1000. Cmp We know the SW corner coordinates, we know that each octant is 90 deg Cmp by 90 deg, and we know how many octants we are dealing with. Calc Cmp the NE corner write(6,*) 'calculating the corner points assuming + this is SW corner: ' write(6,*) swlatsv,swlonsv write(6,*) 'and nx and ny is this: ', nx,ny ne(1)= ( swlatsv + int(ny/73.+0.5)*90000 )/1000. ne(2)= ( swlonsv + int(nx/73.+0.5)*90000 )/1000. if (ne(2) .gt. 180.) ne(2)=ne(2)-360. write(6,*) 'using these values for the corners!!!!' write(6,*) 'SW ', sw(1),sw(2) write(6,*) 'NE ', ne(1),ne(2) dlat=1.25 dlon=1.25 Cmp write(1) nx,ny,nz,lat1,lat2,lon0,sw,ne ! write(1) nx,ny,nz,lat1,lon0,dlat,dlon c c *** Write isobaric upper air data. c Cmp write a dummy value for WAFS slp field to prevent Cmp read errors. slp=-9999. write(1) uw write(1) vw C write(1) th write(1) ht write(1) rh C write(1) ex write(1) pr write(1) slp Cmp n=index(outdir,' ')-1 gdsfile=outdir(1:n)//'/'//'gdsinfo.'//model n=index(gdsfile,' ')-1 write(6,*) 'GDS written to ', gdsfile(1:n) open(unit=14,file=gdsfile(1:n),form='unformatted', + access='sequential') kgds(2)=nx kgds(3)=ny kgds(4)=sw(1)*1000 kgds(5)=sw(2)*1000 kgds(7)=ne(1)*1000 kgds(8)=ne(2)*1000 write(6,*) 'writing kgds ', (kgds(I),i=1,14) write(14) kgds close(14) C write(6,*) 'writing these values to the ETA_wafs file...' do k=1,nz write(6,*) 'level, value ', k, ht(nx/2,ny/2,k) C write(6,*) 'level, value ', k, th(nx/2,ny/2,k) C write(6,*) 'level, value ', k, uw(nx/2,ny/2,k) C write(6,*) 'level, value ', k, vw(nx/2,ny/2,k) C write(6,*) 'level, value ', k, ex(nx/2,ny/2,k) enddo close(1) c endif enddo C end loop III return end c c=============================================================================== c subroutine read_degrib(etafile,nxny,nz,pr . ,ht,tp,rh,uw,vw,jinput,atime) c implicit none c integer nxny,nz,i,nx,ny parameter(nx=73,ny=73) c real pr(nz),ht(nx,ny,nz),tp(nx,ny,nz) . ,rh(nx,ny,nz),uw(nx,ny,nz),vw(nx,ny,nz) . ,dummy,waftmp(nxny) c integer kgds(200),kpds(200),len,kerr . ,lenpds,lenkgds,nwords,kpdsl . ,ijulday,julday_laps,j,k,datsav(9) . ,IRETO,IRET1,JPDS(200),JGDS(200),KNUM,jinput . ,HEMFLAG logical BITMAP(nxny) c character*(*) etafile character(LEN=1):: pds(200) character(LEN=11):: atime c_______________________________________________________________________________ c len=index(etafile//' ',' ')-1 c c c *** Fill time stamp (assume all records are for same time). c call get_gds(etafile(1:len),jinput,datsav,kpds) C ijulday=julday_laps(kpds(10),kpds(9),kpds(8)) if (kpds(8) .eq. 100) kpds(8)=0 write(atime,'(4(i2.2),i3.3)') . kpds(8),kpds(9),kpds(10),kpds(11),kpds(14) write(6,*) 'atime= ', atime c c *** Fill a missing value flag into first space of each variable. c do k=1,nz ht(1,1,k)=-99999. tp(1,1,k)=-99999. rh(1,1,k)=-99999. uw(1,1,k)=-99999. vw(1,1,k)=-99999. enddo c c *** Now put the data into the corresponding arrays. c call baopen(11,etafile,IRETO) if (IRETO .ne. 0) write(6,*) 'BAOPEN TROUBLE!!!! ', IRETO jpds=-1 jgds=-1 jpds(3)=jinput C ISOBARIC DATA C C C ************************************************** C C NOTE::::NEED TO DISTINGUISH BETWEEN NH AND SH OCTANTS C AS THEY GIVE A DIFFERENT ARGUMENT TO W3FT33!!!!!!!!!!!!!! C C ************************************************** C C jpds(6)=100 if (jpds(3).le.40) HEMFLAG=1 if (jpds(3).ge.41) HEMFLAG=-1 do k=1,nz jpds(7)=nint(pr(k)) jpds(5)=7 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,waftmp,IRET1) if (IRET1 .eq. 0) call w3ft33(waftmp,ht(1,1,k),HEMFLAG) if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5) jpds(5)=11 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,waftmp,IRET1) if (IRET1 .eq. 0) call w3ft33(waftmp,tp(1,1,k),HEMFLAG) if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5) jpds(5)=52 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,waftmp,IRET1) if (IRET1 .eq. 0) call w3ft33(waftmp,rh(1,1,k),HEMFLAG) if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5) jpds(5)=33 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,waftmp,IRET1) if (IRET1 .eq. 0) call w3ft33(waftmp,uw(1,1,k),HEMFLAG) if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5) jpds(5)=34 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,waftmp,IRET1) if (IRET1 .eq. 0) call w3ft33(waftmp,vw(1,1,k),HEMFLAG) if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5) C write(6,*) 'Z,T,Q,U,V ', ht(nx/2,ny/2,k),tp(nx/2,ny/2,k), C + rh(nx/2,ny/2,k),uw(nx/2,ny/2,k),vw(nx/2,ny/2,k) enddo c c *** Normal finish. c 1000 continue return c c *** Premature end of file. c 1100 continue print *,'Premature end of file.' print *,'Abort...' stop c end c c=============================================================================== c subroutine rh_fix(nx,ny,nz,rh,pr) c implicit none c integer nx,ny,nz,i,j,k,kk c real rh(nx,ny,nz),pr(nz) c_______________________________________________________________________________ c c *** Fix bottom levels if necessary. c if (rh(1,1,1) .eq. -99999.) then do k=2,nz if (rh(1,1,k) .ne. -99999.) then do kk=k-1,1,-1 print *,'Copying',nint(pr(kk+1)),' mb to' . , nint(pr(kk)),' mb.' do j=1,ny do i=1,nx rh(i,j,kk)=rh(i,j,kk+1) enddo enddo enddo goto 10 endif enddo print *,'RH patch did not work.' stop endif c c *** Fix upper levels if necessary. c 10 continue if (rh(1,1,nz) .eq. -99999.) then do k=nz-1,1,-1 if (rh(1,1,k) .ne. -99999.) then do kk=k+1,nz print *,'Copying',nint(pr(kk-1)),' mb to' . , nint(pr(kk)),' mb.' do j=1,ny do i=1,nx rh(i,j,kk)=rh(i,j,kk-1) enddo enddo enddo goto 20 endif enddo endif c 20 continue do k=1,nz if (rh(1,1,k) .eq. -99999.) then print *,'RH patch did not work.' stop endif enddo c return end c c=============================================================================== c subroutine es_ini c common /estab/esat(15000:45000),es(15000:45000) c c *** Create tables of the saturation vapour pressure with up to c two decimal figures of accuraccy: c do it=15000,45000 t=it*0.01 p1 = 11.344-0.0303998*t p2 = 3.49149-1302.8844/t c1 = 23.832241-5.02808*alog10(t) esat(it) = 10.**(c1-1.3816E-7*10.**p1+ . 8.1328E-3*10.**p2-2949.076/t) es(it) = 610.78*exp(17.269*(t-273.16)/(t-35.86)) enddo c return end c c=============================================================================== c function julday_laps(day,month,year) c implicit none c integer julday_laps . ,day,month,year . ,ndays(12),i c data ndays/31,28,31,30,31,30,31,31,30,31,30,31/ c_______________________________________________________________________________ c julday_laps=0 do i=1,month-1 julday_laps=julday_laps+ndays(i) enddo julday_laps=julday_laps+day if (mod(year,4) .eq. 0 .and. month .gt. 2) . julday_laps=julday_laps+1 return c end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine get_gds(etafile,jinput,gdsinfo,kpds) character*(*) etafile character(LEN=1):: pds(50) Ctst integer*4 kgds(200),kpds(50),len,kerr integer kgds(200),kpds(200),len,kerr . ,lenpds,lenkgds,nwords,kpdsl . ,ijulday,julday_laps,j,k,gdsinfo(9),jinput . ,gdsav,IRETO,JGDS(200),JPDS(200) real tmp(6000) LOGICAL BITMAP(6000) nxny=73*73 JPDS=-1 JGDS=-1 len=index(etafile//' ',' ')-1 call baopen(11,etafile(1:len),IRETO) C write(6,*) 'BAOPEN in get_gds: ', IRETO if (IRETO .ne. 0) then print *,'Error opening unit=11, file name = ',etafile(1:len) . ,' iostat = ',kerr stop endif jpds(5)=7 jpds(6)=100 jpds(7)=500 Cmp set octant jpds(3)=jinput Cmp C write(6,*) 'calling getgb' call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,tmp,IRET1) C write(6,*) 'return from getgb ', IRET1 gdsinfo(1)=KGDS(2) gdsinfo(2)=KGDS(3) gdsinfo(3)=KGDS(4) gdsinfo(4)=KGDS(5) gdsinfo(5)=KGDS(7) gdsinfo(6)=KGDS(8) gdsinfo(7)=KGDS(9) gdsinfo(8)=KGDS(12) gdsinfo(9)=KGDS(13) Cmp write(6,*) gdsinfo return print *,'GETGDS PROBLEM' stop end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE W3FT33(AIN,OUT,NSFLAG) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: W3FT33 THICKEN THINNED WAFS GRIB GRID 37-44 C PRGMMR: RALPH PETTERSON ORG: W/NMCXX DATE: 94-11-13 C C ABSTRACT: SUBROUTINE THICKENS ONE THINNED WAFS GRIB GRID TO A C REAL ARRAY OF 5329 NUMBERS (73,73) 1.25 DEGREE GRID. C C PROGRAM HISTORY LOG: C 94-??-?? RALPH PETERSON C 94-11-07 R.E.JONES ADD DOC BLOCK, CHANGE CALL TO 3 C PARAMETERS. REPLACE COS WITH TABLE C LOOKUP. C 95-06-02 RALPH PETERSON CHANGES TO CORRECT MISS-POSITION C BETWEEN + OR - 8.75 N/S. C 95-06-03 R.E.JONES CHANGES SO 8 ROWS WITH 73 VALUES C ARE NOT THICKENED, 10% FASTER. C C USAGE: CALL W3FT33(AIN, OUT, NSFLAG) C INPUT ARGUMENT LIST: C AIN - REAL 3447 WORD ARRAY WITH UNPACKED THINNED WAFS C GRIB TYPE 37-44. C NSFLAG - INTEGER = 1 AIN IS WAFS GRIB GRID 37-40 N. HEMI. C = -1 AIN IS WAFS GRIB GRID 41-44 S. HEMI. C C OUTPUT ARGUMENT LIST: C OUT - REAL (73,73) WORD ARRAY WITH THICKENED WAFS GRIB C GRID 37-44. C C REMARKS: THE POLE POINT FOR U AND V WIND COMPONENTS WILL HAVE ONLY C ONE POINT. IF YOU NEED THE POLE ROW CORRECTED SEE PAGE 9 SECTION C 1 IN OFFICE NOTE 388. YOU NEED BOTH U AND V TO MAKE THE C CORRECTION. C C ATTRIBUTES: C LANGUAGE: SiliconGraphics 5.2 FORTRAN 77 C MACHINE: SiliconGraphics IRIS-4D/25, 35, INDIGO, Indy C C$$$ C PARAMETER (NX=73,NY=73) PARAMETER (NIN=3447) C REAL AIN(*) REAL OUT(NX,NY) C INTEGER IPOINT(NX) C SAVE C DATA IPOINT/ & 73, 73, 73, 73, 73, 73, 73, 73, 72, 72, 72, 71, 71, 71, 70, & 70, 69, 69, 68, 67, 67, 66, 65, 65, 64, 63, 62, 61, 60, 60, & 59, 58, 57, 56, 55, 54, 52, 51, 50, 49, 48, 47, 45, 44, 43, & 42, 40, 39, 38, 36, 35, 33, 32, 30, 29, 28, 26, 25, 23, 22, & 20, 19, 17, 16, 14, 12, 11, 9, 8, 6, 5, 3, 2/ C NXM = NX - 1 FNXM = FLOAT(NXM) C C TEST FOR GRIDS (37-40) C IF (NSFLAG.GT.0) THEN C C DO NOT THICKEN 8 ROWS WITH 73 VALUES, MOVE DATA C TO OUT ARRAY. GRIDS (37-40) N. C IS = 0 DO J = 1,8 DO I = 1,NX IS = IS + 1 OUT(I,J) = AIN(IS) END DO END DO C IE = NX * 8 DO J = 9,NY NPOINT = IPOINT(J) IS = IE + 1 IE = IS + NPOINT - 1 DPTS = (FLOAT(NPOINT)-1.) / FNXM PW = 1.0 PE = PW + DPTS OUT(1,J) = AIN(IS) VALW = AIN(IS) VALE = AIN(IS+1) DVAL = (VALE-VALW) DO I = 2,NXM WGHT = PE -FLOAT(IFIX(PE)) OUT(I,J) = VALW + WGHT * DVAL PW = PE PE = PE + DPTS IF (IFIX(PW).NE.IFIX(PE)) THEN IS = IS + 1 VALW = VALE VALE = AIN(IS+1) DVAL = (VALE - VALW) END IF END DO OUT(NX,J) = AIN(IE) END DO C ELSE C C DO NOT THICKEN 8 ROWS WITH 73 VALUES, MOVE DATA C TO OUT ARRAY. GRIDS (41-44) S. C IS = NIN - (8 * NX) DO J = 66,NY DO I = 1,NX IS = IS + 1 OUT(I,J) = AIN(IS) END DO END DO C IE = 0 DO J = 1,65 NPOINT = IPOINT(74-J) IS = IE + 1 IE = IS + NPOINT - 1 DPTS = (FLOAT(NPOINT)-1.) / FNXM PW = 1.0 PE = PW + DPTS OUT(1,J) = AIN(IS) VALW = AIN(IS) VALE = AIN(IS+1) DVAL = (VALE-VALW) DO I = 2,NXM WGHT = PE -FLOAT(IFIX(PE)) OUT(I,J) = VALW + WGHT * DVAL PW = PE PE = PE + DPTS IF (IFIX(PW).NE.IFIX(PE)) THEN IS = IS + 1 VALW = VALE VALE = AIN(IS+1) DVAL = (VALE - VALW) END IF END DO OUT(NX,J) = AIN(IE) END DO END IF C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine get_fullgds(etafile,nx,ny,kgds) character*(*) etafile character(LEN=1) pds(50) integer kgds(200),kpds(200),len,kerr,jpds(200) . ,lenpds,lenkgds,nwords,kpdsl,jgds(200) . ,j,k,KNUM,nx,ny logical bitmap(nx*ny) real tmp(nx*ny) write(6,*) 'inside get_fullgds...', etafile nxny=nx*ny len=index(etafile//' ',' ')-1 jpds=-1 jgds=-1 jpds(5)=11 jpds(6)=100 jpds(7)=500 write(6,*) 'calling getgb ' write(6,*) 'jpds = ', jpds write(6,*) 'jgds = ', jgds call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,tmp,IRET1) write(6,*) 'return from getgb ', IRET1 if (IRET1 .ne. 0) then print *,'Error getting GDS in get_fullgds ', IRET1 stop endif write(6,*) 'leaving get_fullgds' return end