program degrib2model_driver c implicit none c integer nx,ny,nz,kpds(200),datsav(9),n !Degrib grid dimensions 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 read(5,'(a)') gribfile read(5,'(a)') outdir c n=index(gribfile,' ')-1 write(6,*) 'first call of GET_GDS' CALL GET_GDS(gribfile(1:n),datsav,kpds) nx=datsav(1) ny=datsav(2) if (nx .eq. 360 .and. ny .eq. 181) then nz=26 elseif (nx .eq. 144 .and. ny .eq. 73) then nz=17 else write(6,*) 'what grid is this??? ', nx,ny write(6,*) ' ' write(6,*) 'I do not know the vertical levels for this ' write(6,*) 'grid...please add some code to dgeta2model.f' write(6,*) ' ' STOP 999 endif c Cold call degrib2model(degrib_dir,degrib_file,outdir,nx,ny,nz) call degrib2model(gribfile,outdir,nx,ny,nz) c end c c=============================================================================== c Cold subroutine degrib2model(degrib_dir,degrib_file,outdir,nx,ny,nz) subroutine degrib2model(gribfile,outdir,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,ifcsthr . ,count,kgds(200) c real ht(nx,ny,nz) !Isobaric heights (m) . ,tp(nx,ny,nz) !Isobaric temps (K) . ,th(nx,ny,nz) !Isobaric theta (K) . ,uw(nx,ny,nz) !Isobaric u-wind (m/s) . ,vw(nx,ny,nz) !Isobaric v-wind (m/s) . ,rh(nx,ny,nz) !Isobaric rh,mr (%,kg/kg) . ,pr(nz),pri(nz) !Isobaric pressures (mb) . ,lat1,lat2,lon0,sw(2),ne(2) . ,xe,mrsat,esat,es . ,avnlevs(26),reanlevs(17) . ,slp(nx,ny,12) !slp(i,j,1)=EMSL;slp(i,j,2)=PMSL . ,temp(nx,ny) !slp(i,j,3)=PSFC;slp(i,j,4)=ZSFC !slp(i,j,5 & 6) are 0-10 cm STC and SMC !slp(i,j,7 & 8) are 10-200 cm STC and SMC REAL, ALLOCATABLE:: htout(:,:,:), + uwout(:,:,:),vwout(:,:,:),rhout(:,:,:), + slpout(:,:,:) c character*255 gribfile,outdir,outfile,gdsfile character*2 gproj character*11 atime character*7 model data AVNLEVS/1000,975,950,925,900,850,800,750,700,650, + 600,550,500,450,400,350,300,250,200,150, + 100, 70, 50, 30, 20, 10/ data REANLEVS/1000,925,850,700,600,500, + 400,300,250,200,150,100, + 70, 50, 30, 20, 10/ c common /estab/esat(15000:45000),es(15000:45000) c_______________________________________________________________________________ c c *** Fill pressure levels. c Cmp nsfcfld will be 12 although 4 of these are set to missing for AVN... nsfcfld=12 Cmp write(6,*) 'using hardwired stuff! ', nx,ny,nz do k=1,nz if (nz .eq. 26) then pr(k)=AVNLEVS(k) elseif (nz .eq. 17) then pr(k)=reanlevs(k) endif C write(6,*) 'pr(K): ', pr(K) enddo c c *** Read in degrib data. c n=index(gribfile,' ')-1 write(6,*) 'call read_degrib' call read_degrib(gribfile(1:n) . ,nx*ny,nz,pr,ht,tp,rh,uw,vw,slp,atime) write(6,*) 'return read_degrib' c c *** Check for any missing data. c 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 Cmp Cmp Handle missing surface data.... Cmp C C For reanalysis runs these data will not be used. Just bogus C to hold a place. C write(6,*) 'going to use nsfcfld= ', nsfcfld C If rean put bogus data into slp if (nz .eq. 17) then ! Chou if (REAN) then % deveria controlar via variavel logica write(6,*) ' reading soil parm from REAN file, skipping bogus data' Chou do K=5,8 Chou do j=1,ny Chou do i=1,nx Chou if (mod(k,2) .eq. 0) then Chou slp(i,j,k)=0.14 Chou else Chou slp(i,j,k)=275. Chou endif Chou enddo Chou enddo Chou enddo Chou endif Chou do k=5,nsfcfld if (slp(1,1,k) .eq. -99999.) then write(6,*) 'SURFACE DATA MISSING... ', K if (k.eq.9.or.k.eq.11) then write(6,*) 'filling soil temp data...' do j=1,ny do i=1,nx slp(i,j,k)=slp(i,j,7) enddo enddo elseif(k.eq.10.or.k.eq.12.) then write(6,*) 'filling soil moisture data...' do j=1,ny do i=1,nx slp(i,j,k)=slp(i,j,8) enddo enddo endif endif Cmp if (mod(k,2) .eq. 0) then do j=1,ny do i=1,nx if (slp(i,j,k) .eq. 0) then slp(i,j,k)=0.14 endif enddo enddo else do j=1,ny do i=1,nx if (slp(i,j,k) .eq. 0) then slp(i,j,k)=273.15 endif enddo enddo endif Cmp 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 C 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. If reanalysis data, skip this part Cmp if (nz .eq. 17) goto 969 do k=1,nz do j=1,ny do i=1,nx temp(i,j)=th(i,j,k) enddo enddo do J=1,ny do i=1,nx th(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=uw(i,j,k) enddo enddo do J=1,ny do i=1,nx uw(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=vw(i,j,k) enddo enddo do J=1,ny do i=1,nx vw(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=ht(i,j,k) enddo enddo if (k.eq.1) write(6,*) 'flipping isobaric data - ' do J=1,ny do i=1,nx ht(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=rh(i,j,k) enddo enddo do J=1,ny do i=1,nx rh(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx C temp(i,j)=ex(i,j,k) enddo enddo do J=1,ny do i=1,nx C ex(i,j,k)=temp(i,ny-J+1) enddo enddo C----------------------------------------- enddo Cmp generalize eventually! do k=1,nsfcfld C write(6,*) 'flipping surface fields... ' do j=1,ny do i=1,nx temp(i,j)=slp(i,j,k) enddo enddo do J=1,ny do i=1,nx slp(i,j,k)=temp(i,ny-J+1) enddo enddo enddo 969 continue CC CC CC Add bogus column if reanalysis data CC CC if (nx .eq. 144) then ALLOCATE(HTOUT(145,ny,nz),UWOUT(145,ny,nz),VWOUT(145,ny,nz)) ALLOCATE(RHOUT(145,ny,nz),SLPOUT(145,ny,nsfcfld)) DO J=1,NY DO I=1,NX DO K=1,NSFCFLD SLPOUT(I,J,K)=SLP(I,J,K) ENDDO DO K=1,NZ HTOUT(I,J,K)=HT(I,J,K) UWOUT(I,J,K)=UW(I,J,K) VWOUT(I,J,K)=VW(I,J,K) RHOUT(I,J,K)=RH(I,J,K) ENDDO ENDDO DO K=1,NSFCFLD SLPOUT(145,J,K)=SLP(1,J,K) ENDDO DO K=1,NZ HTOUT(145,J,K)=HT(1,J,K) UWOUT(145,J,K)=UW(1,J,K) VWOUT(145,J,K)=VW(1,J,K) RHOUT(145,J,K)=RH(1,J,K) ENDDO ENDDO endif c c c *** Create output file name. c n=index(outdir,' ')-1 Cmp make ETA_avn to differentiate from files based on eta grib input if (nz .eq. 26) then model='ETA_avn' elseif (nz .eq. 17) then model='ETA_rnl' nx=145 endif 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 gproj='LL' 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' if (nx .eq. 145) then KGDS(8)=0 KGDS(2)=145 endif 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') C C MODIFY GDS TO REFLECT FLIPPING OF DATA IN N-S SENSE C if (nz .eq. 26) then kgds(4)= -kgds(4) kgds(7)= -kgds(7) endif write(6,*) 'writing kgds ', (kgds(I),i=1,14) write(14) kgds close(14) write(6,*) 'past 14 close' C C end new stuff c c *** Write isobaric upper air data. c write(6,*) 'starting writes to unit1' if (nx .ne. 145) then write(1) uw write(1) vw write(1) ht write(1) rh write(1) pr write(1) slp else WRITE(1) UWOUT WRITE(1) VWOUT WRITE(1) HTOUT WRITE(1) RHOUT WRITE(1) PR WRITE(1) SLPOUT endif if (allocated(uwout)) deallocate(uwout) if (allocated(vwout)) deallocate(vwout) if (allocated(htout)) deallocate(htout) if (allocated(rhout)) deallocate(rhout) if (allocated(slpout)) deallocate(slpout) Cmp write(6,*) 'writing these sfc values to the ETA_avn file...' do k=1,nsfcfld write(6,*) 'field, value ', k, slp(nx/2,ny/2,k) enddo do k=1,nsfcfld write(6,*) 'field, value ', k, slp(3*nx/4,3*ny/4,k) enddo c close(1) c return end c c=============================================================================== c subroutine read_degrib(etafile,nxny,nz,pr . ,ht,tp,rh,uw,vw,slp,atime) c implicit none c integer nxny,nz,i c real pr(nz),ht(nxny,nz),tp(nxny,nz) . ,rh(nxny,nz),uw(nxny,nz),vw(nxny,nz) . ,dummy,slp(nxny,12) c integer kgds(200),kpds(50),len,kerr . ,lenpds,lenkgds,nwords,kpdsl . ,j,k,datsav(9) . ,IRETO,IRET1,JPDS(200),JGDS(200),KNUM logical BITMAP(nxny) c character*(*) etafile character(LEN=1):: pds(50) character(LEN=11):: atime c_______________________________________________________________________________ c len=index(etafile//' ',' ')-1 c c c *** Fill time stamp (assume all records are for same time). c write(6,*) 'calling get_gds in read_degrib ',etafile call get_gds(etafile(1:len),datsav,kpds) write(6,*) 'return calling get_gds in read_degrib' if (kpds(8) .eq. 100) kpds(8)=0 write(atime,'(i2.2,i2.2,i2.2,i2.2,i3.3)') . kpds(8),kpds(9),kpds(10),kpds(11),kpds(14) c c *** Fill a missing value flag into first space of each variable. c do k=1,nz ht(1,k)=-99999. tp(1,k)=-99999. rh(1,k)=-99999. uw(1,k)=-99999. vw(1,k)=-99999. enddo Cmp initialize surface fields to -99999. so the interp code can handle Cmp appropriately write(6,*) 'nxny,datsav(1)*datsav(2): ', nxny,datsav(1)*datsav(2) do k=1,12 do j=1,datsav(1)*datsav(2) slp(j,k)=-99999. enddo enddo Cmp c c *** Now put the data into the corresponding arrays. c C write(6,*) 'here2 : calling baopen ', etafile C call baopen(11,etafile,IRETO) write(6,*) 'here3' if (IRETO .ne. 0) write(6,*) 'BAOPEN TROUBLE!!!! ', IRETO jpds=-1 jgds=-1 jpds(5)=81 jpds(6)=1 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,1),IRET1) write(6,*) 'first GETGB, IRET1= ', IRET1 if (IRET1 .eq. 0) then write(6,*) 'LAND/SEA READ!!!!! ' write(6,*) (slp(j,1),j=nwords/2,nwords/2+5) endif jpds(5)=2 jpds(6)=102 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,2),IRET1) write(6,*) 'PMSL READ!!!!! ' write(6,*) (slp(j,2),j=nwords/2,nwords/2+8) jpds(5)=1 jpds(6)=1 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,3),IRET1) write(6,*) 'PSFC READ!!!!! ' write(6,*) (slp(j,3),j=nwords/2,nwords/2+8) jpds(5)=7 jpds(6)=1 jpds(7)=0 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,4),IRET1) write(6,*) 'ZSFC READ!!!!! ' write(6,*) (slp(j,4),j=nwords/2,nwords/2+8) Cmp SOIL FIELDS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C AVN stores soil T in 11, ETA in 85 jpds(5)=11 jpds(6)=112 jpds(7)=10 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,5),IRET1) if (IRET1.eq.0) then write(6,*) 'found soil temp at ',jpds(7) write(6,*) (slp(j,5),j=nwords/4,nwords/4+8) write(6,*) (slp(j,5),j=nwords/2,nwords/2+8) endif jpds(7)=2600 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,7),IRET1) if (IRET1.eq.0) then write(6,*) 'found soil temp at ',jpds(7) write(6,*) (slp(j,7),j=nwords/2,nwords/2+8) write(6,*) (slp(j,7),j=nwords/2,nwords/2+8) endif jpds(7)=10340 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,9),IRET1) if (IRET1.eq.0) then write(6,*) 'found soil temp at ',jpds(7) write(6,*) (slp(j,9),j=nwords/2,nwords/2+5) endif jpds(7)=25800 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,11),IRET1) if (IRET1.eq.0) then write(6,*) 'found soil temp at ',jpds(7) write(6,*) (slp(j,11),j=nwords/4,nwords/4+5) endif C 10-200 layer jpds(7)=2760 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,7),IRET1) if (IRET1.eq.0) then write(6,*) 'found soil temp at ',jpds(7) write(6,*) (slp(j,7),j=nwords/4,nwords/4+8) endif C MOISTURE jpds(5)=144 jpds(7)=10 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,6),IRET1) if (IRET1.eq.0) then write(6,*) 'found soil wet at ',jpds(7) write(6,*) (slp(j,6),j=nwords/4,nwords/4+8) endif jpds(7)=2600 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,8),IRET1) if (IRET1.eq.0) then write(6,*) 'found soil wet at ',jpds(7) write(6,*) (slp(j,8),j=nwords/2,nwords/2+8) endif jpds(7)=10340 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,10),IRET1) if (IRET1.eq.0) then write(6,*) 'found soil wet at ',jpds(7) write(6,*) (slp(j,10),j=nwords/2,nwords/2+5) endif jpds(7)=25800 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,12),IRET1) if (IRET1.eq.0) then write(6,*) 'found soil wet at ',jpds(7) write(6,*) (slp(j,12),j=nwords/2,nwords/2+5) endif C 10-200 layer jpds(7)=2760 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,slp(1,8),IRET1) if (IRET1.eq.0) then write(6,*) 'found soil wet at ',jpds(7) write(6,*) (slp(j,8),j=nwords/4,nwords/4+5) endif C ISOBARIC DATA jpds(6)=100 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,ht(1,k),IRET1) 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,tp(1,k),IRET1) 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,rh(1,k),IRET1) 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,uw(1,k),IRET1) 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,vw(1,k),IRET1) if (IRET1 .ne. 0) write(6,*) ' AT LEVEL ', jpds(7) , jpds(5) write(6,*) 'Z,T,Q,U,V ', ht(nxny/2,k),tp(nxny/2,k), + rh(nxny/2,k),uw(nxny/2,k),vw(nxny/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*4 nx,ny,nz,i,j,k,kk c real*4 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=============================================================================== subroutine get_gds(etafile,gdsinfo,kpds) character*(*) etafile character*1 pds(50) Ctst integer*4 kgds(200),kpds(50),len,kerr integer*4 kgds(200),kpds(200),len,kerr . ,lenpds,lenkgds,nwords,kpdsl . ,j,k,gdsinfo(9) . ,gdsav,IRETO,JGDS(200),JPDS(200) real tmp(70000) LOGICAL*1 BITMAP(70000) nxny=360*181 JPDS=-1 JGDS=-1 len=index(etafile//' ',' ')-1 write(6,*) 'want to open unit 11 for ', etafile(1:len) call baopen(11,etafile(1:len),IRETO) write(6,*)'AQUI!!' 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 call getgb(11,0,nxny,0,JPDS,JGDS,nwords,KNUM,KPDS,KGDS, & BITMAP,tmp,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) write(6,*) gdsinfo return print *,'GETGDS PROBLEM' stop end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine get_fullgds(etafile,nx,ny,kgds) character*(*) etafile character*1 pds(50) integer*4 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 C write(6,*) 'calling getgb ' C write(6,*) 'jpds = ', jpds C 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