subroutine nvap_byteswap(bytebuff, isize) byte bytebuff(1440), btemp integer*2 isize c do i = 1,isize ifb = (i-1)*2 + 1 isb = ifb + 1 btemp = bytebuff(ifb) bytebuff(ifb) = bytebuff(isb) bytebuff(isb) = btemp enddo return end c c------------------------------------------------------------------------------------ c subroutine nvap_globalave(glb) c c******************************************************************** c This routine averages the whole globe, the southern hemisphere * c and the northern hemisphere zones into averages weighted by * c earth area. * c * c written by Dave Randel cira/csu/metsat * c on march 20 1992/ modified 8/2003 * c * c * c output parameters glb(3) nh,sh,global averages * c * c******************************************************************** c real fact(360) integer nsx(3),nex(3), ysize_2 real glb(3) c c--------common block which holds the data and header information--------- c real nvdata(720,360), headr(4), zindef integer*2 headi(11) character label*40 common /nvapdata/ nvdata, headr, zindef, headi, label c c**** define NH and SH zones to average c cmy fix nsx(1)=1 nsx(2)=2 nsx(3)=3 nex(1)=1 nex(2)=2 nex(3)=3 cend my fix ysize_2 = nint( float(headi(11)) / 2.) if(nsx(1).eq. 0 .and. nex(1).eq.0) then nsx(1) = 1 nsx(2) = ysize_2 + 1 nsx(3) = 1 nex(1) = ysize_2 nex(2) = headi(11) nex(3) = headi(11) endif c c*** define area scaling factors ******** c theta = headr(3) + (0.5 * -headr(2)) !starting NH at 90 degrees qcons = 3.14159/180. do j = 1,ysize_2 fact(j)=(sin(qcons*theta)-sin(qcons*(theta+headr(2))) ) theta = theta + headr(2) enddo do j = ysize_2+1,headi(11) fact(j) = fact(headi(11) + 1 - j) enddo c c*** sum up zones applying weighting ******** c do m = 1,3 sum = 0. glb(m) = 0. do j = nsx(m),nex(m) do i = 1,headi(10) if(nvdata(i,j) .ne. zindef) then sum = sum + fact(j) endif enddo enddo c c*** calc averages ************************ c do j = nsx(m),nex(m) rmult = fact(j) / sum do i = 1,headi(10) if(nvdata(i,j).ne. zindef) then glb(m) = (nvdata(i,j) * rmult) + glb(m) endif enddo enddo enddo c return end c c-------------------------------------------------------------------------- c subroutine nvap_read(flname,igrd,istat) C C******************************************************************** C C This routine reads the nvap data files. It will work with either c 1 degree data (1988 - 1999) or 1/2 degree data (2000 - ??). The c input data file is opened with direct access in order to jump c down the file to other grids without having to read all the bytes. c c This can cause some problems between operating systems since in c fortran a record length must be specified in the OPEN statement c for direct access files. This recl may be in bytes OR bytes/4 c depending on operating systems or compiler settings. c c The switch to byte swap data needs to be turned on if working in UNIX c or other big_endian system (remove comment character in byte swap c line near the bottom of this subroutine) c C written by: Dave Randel CIRA/STC-METSAT C August, 2003 C C INPUT PARAMETERS flname* : input filename c igrd : grid number in file C C OUTPUT PARAMETERS istat : status code (0 = normal) C C******************************************************************** C character flname*(*) integer*2 datatc,ixsize,iysize, igrd, iy2k1, iy2k2 integer*2 intbuff(720) integer*4 idt,istat, ishead C real*4 offset,scale, datbuff(720) cmy fix character blank1*1, blank6*6,buff*504, bufft*38 character header_pad*1296 cend my fix C C--------COMMON BLOCK WHICH HOLDS THE DATA AND HEADER INFORMATION--------- c --------INCLUDE THIS IN YOUR CALLING PROGRAM --------------- C real nvdata(720,360), headr(4), zindef integer*2 headi(11) character label*40 common /nvapdata/ nvdata, headr, zindef, headi, label c c------------------------------------------------------------------------- c istat = -1 C C**OPEN FILE TO GET SIZE INFORMATION FROM HEADER C cmy fix open(unit=101,file=flname,form='formatted', open(unit=101,file=flname,form='unformatted', > access='direct',recl = 38, status='old',err=90) read(101,rec=1,err=91) bufft read(bufft,1) ixsize,iysize 1 format( 30x, i4, i4) close(unit=101) C C**SET RECORD LENGTH AND OPEN FILE------------------------------ C if(ixsize .eq. 360) then ! 1 degree nvap data irecl = 720 / 4 ! bytes / 4 ***this may need to be just bytes ipadl = 720 - 144 ! pad bytes in header idatal = 720 ! length in bytes of data (360*2) endif if(ixsize .eq. 720) then ! 1/2 degree nvap data irecl = 1440 / 4 ! bytes / 4 ***this may need to be just bytes ipadl = 1440 - 144 ! pad bytes in header idatal = 1440 ! length in bytes of data (720*2) endif iheadl = 144 ! ascii header length c c --- open statement for reading on a WINDOWS PC: depending on the platform c used, integer*2 data may need to have the bytes swapped. You might be c able to do this with a switch in the open statement. c open(unit=101,file=flname,form='unformatted', > access='direct',status='old',recl = irecl, err=90) C**COMPUTE STARTING RECORD FROM INPUT GRID NUMBER AND FIELD SIZE cmy fix set this variable c igrd=1 ishead = (igrd-1) * (iysize+1) + 1 C C**READ IN REQUESTED GRID HEADER (HEADER FORMAT SPECIFIED AT BOTTOM OF THIS FILE) C c2---6--------------------------------40------------------------------71 read(101,rec=ishead,err=92) buff(1:iheadl), header_pad(1:ipadl) read(buff,11)tapec,(headi(J),J=1,3), blank1, iy2k1, iy2k2, > (headi(j),j=4,11), > (headr(j),j=1,4), > offset, scale, zindef, > blank6, label C 11 FORMAT( a4, i3, i3, i1, a1, i2, i2, > i2, i3, i2, i2, i3, i2, > i4, i4, f6.2, f6.2, f7.3, f8.3, > e11.5, e11.5, e11.5, > a6, a40) C C**CHANGE 2 CHARACTER YEAR TO 4 CHARACTER YEAR IF NECESSARY C IF(iy2k1.eq.0) then headi(4) = 19*100 + headi(4) headi(7) = 19*100 + headi(7) else headi(4) = iy2k1*100 + headi(4) headi(7) = iy2k2*100 + headi(7) endif c WRITE(6,*)'READ_NVAP: START, END YEAR = ',HEADI(4),HEADI(7) C c**INITIALIZE DATA ARRAY WITH ZINDEF VALUES c do i = 1,360 do j = 1,720 nvdata(j,i) = zindef enddo enddo C C**READ IN ALL DATA RECORDS in GRID and apply scale and offset C do i = 1,iysize read(101,rec=ishead+i,err=93) (intbuff(j),j=1,ixsize) c c--------this swap necessary for 88-95 nvap ssmi (lwp,pwc,clw) grids cmyfix if(headi(3) .eq. 7) call nvap_byteswap(intbuff, ixsize) cmyfix c c------------------------------------------------------------------ c-------- if necessary (UNIX, MacOS...) call byte swapping subroutine c call nvap_byteswap(intbuff, ixsize) c c------------------------------------------------------------------ c------------------------------------------------------------------ do j = 1,ixsize if(float(intbuff(j)) .eq. zindef) then nvdata(j,i) = zindef else nvdata(j,i) = (float(intbuff(j)) * scale) + offset endif enddo enddo C C**CLOSE FILE AND RETURN C istat = 0 !sucessful close(unit=101) return C C**ERROR FOUND C 90 write(6,*)'ERROR opening file' istat = -1 return 91 write(6,*)'ERROR reading first header ' istat = -2 close(unit=101) return 92 write(6,*)'ERROR reading requested header rec', ishead istat = -3 close(unit=101) return 93 write(6,*)'ERROR reading requested data' istat = -4 close(unit=101) return end c c------------------------------------------------------------------------ c subroutine nvap_read_dsb(flname,ifld,dsb,istat) c-------------------------------------------------------------------- c c This routine reads in the NVAP (1/2 degree data) Data Source c Bit file passing it back to the calling program in dsb(720,360). c c ** if on UNIX system or other big_endian OS, remove comment near bottom c ** of this routine to swap bytes c c written by Dave Randel STC/METSAT c February 2003 c-------------------------------------------------------------------- c integer*2 dsb(720,360), intbuff(720), istat character tapec*4,flname*50 integer*2 datatc integer*2 xsize,ysize,iy2k1,iy2k2 real offset,scale character blank1*1, blank6*6, buff*144 C--------COMMON BLOCK WHICH HOLDS THE DATA AND HEADER INFORMATION--------- c real nvdata(720,360), headr(4), zindef integer*2 headi(11) character label*40 common /nvapdata/ nvdata, headr, zindef, headi, label c c------------------------------------------------------------------------- c istat = -1 C C**SET RECORD LENGTH AND OPEN FILE xsize = 720 ysize = 360 idatal = 1440 irecl = idatal / 4 !bytes divided by 4...may have to be changed to bytes open(unit=104,file=flname,form='unformatted', > access='direct',status='old',recl = irecl,err=90) C**COMPUTE STARTING RECORD FROM INPUT FIELD NUMBER AND FIELD SIZE ishead = (ifld-1) * (ysize+1) + 1 C C**READ IN REQUESTED FIELD HEADER C read(104,rec=ishead,err=91) buff read(buff,11)tapec,(headi(j),j=1,3),blank1,iy2k1,iy2k2, > (headi(j),j=4,11), > (headr(j),j=1,4), > offset,scale,zindef, > blank6,label C 11 format( a4, i3, i3, i1, a1, i2, i2, > i2, i3, i2, i2, i3, i2, > i4, i4, f6.2, f6.2, f7.3, f8.3, > e11.5, e11.5, e11.5, > a6, a40) c c**CHANGE 2 CHARACTER YEAR TO 4 CHARACTER YEAR c if(iy2k1.eq.0) then headi(4) = 19*100 + headi(4) headi(7) = 19*100 + headi(7) else headi(4) = iy2k1*100 + headi(4) headi(7) = iy2k2*100 + headi(7) endif C c**Initialize data array c do j = 1,360 do i = 1,720 nvdata(i,j) = ZINDEF enddo enddo C C**READ IN DATA RECORDS C do irec = 1,ysize read(104,rec=ishead+irec,err=92) (intbuff(j),j=1,xsize) c-------- if necessary (UNIX, MacOS...) call byte swapping subroutine c c call nvap_byteswap(intbuff, xsize) c c------------------------------------------------------------------ c do j = 1,xsize dsb(j,irec) = intbuff(j) enddo enddo c c**close file and return c istat = 0 close(unit=104) return C C**error found C 90 write(6,*)'error opening file' istat = -1 return 91 write(6,*)'error reading requested header rec', ishead istat = -2 close(unit=104) return 92 write(6,*)'error reading requested data' istat = -3 close(unit=104) return end c c------------------------------------------------------------------------------------ c subroutine nvap_write(flname, igrd, istat) C C******************************************************************** C C This routine writes out a NVAP formatted grid file to: flname c This version of write will calc a scale and offset to apply to c the data and write out data as a int*2 field. c It will write out at grid location number igrd. c c written by dave randel metsat/cira/csu c in march, 1992 (modified 8/2003) c c input parameters flname* : data filename c igrd : field number c output parameters istat : status code c c******************************************************************** c character tapec*4 real offset,scale character blank6*6, blank1*1 character buff*504 integer*2 dataint(720,360), int2, iy2k1, iy2k2, isy, iey character flname*(*) logical ex character sta*3 c c--------COMMON BLOCK WHICH HOLDS THE DATA AND HEADER INFORMATION--------- c real nvdata(720,360), headr(4), zindef integer*2 headi(11) character label*40 common /nvapdata/ nvdata, headr, zindef, headi, label c c------------------------------------------------------------------------- data blank1,blank6,tapec/' ',' ','ccda'/ c istat = -1 rnewzind = -32768 c c**set size and record length parameters from size specified in header c ixsize = headi(10) if(ixsize .eq. 360) then ! 1 degree nvap data irecl = 720 / 4 ! bytes / 4 ***this may need to be just bytes ipadl = 720 - 144 ! pad bytes in header idatal = 720 ! length in bytes of data (360*2) endif if(ixsize .eq. 720) then ! 1/2 degree nvap data irecl = 1440 / 4 ! bytes / 4 ***this may need to be just bytes ipadl = 1440 - 144 ! pad bytes in header idatal = 1440 ! length in bytes of data (720*2) endif iheadl = 144 ! ascii header length c c**open file for writing c open(unit=102,file=flname,form='unformatted',recl=irecl, > access='direct',status='unknown', err=90) c*** get scale and offset factors, apply to data creating dataint c rmin = 1.0e10 rmax = -1.0e10 igood = 0 do j = 1,headi(11) do i = 1,headi(10) if(nvdata(i,j) .ne. zindef) then if(nvdata(i,j) .gt. rmax) rmax = nvdata(i,j) if(nvdata(i,j) .lt. rmin) rmin = nvdata(i,j) igood = 1 endif enddo enddo c if(igood .eq. 0) then !all points in field zindef do j = 1,headi(11) do i = 1,headi(10) dataint(i,j) = rnewzind scale = 1.0 offset = 0.0 enddo enddo go to 20 endif c if(rmin .eq. rmax) then !all good points in field identical scale=(rmax)/65534. offset = rmax + 32767.* scale else scale=(rmax-rmin)/65534. !normal field with different values offset = rmin + 32767.* scale endif do j = 1,headi(11) do i = 1,headi(10) if(nvdata(i,j) .ne. zindef) then if(scale .ne. 0.0) then dataint(i,j) = nint((nvdata(i,j)-offset) / scale) else dataint(i,j) = rnewzind endif else dataint(i,j) = rnewzind endif enddo enddo 20 continue c c** fill up ascii header buffer c if(headi(4) .gt. 1999) then iy2k1 = 20 else iy2k1 = 19 endif if(headi(7) .gt. 1999) then iy2k2 = 20 else iy2k2 = 19 endif isy = headi(4) - (iy2k1*100) iey = headi(7) - (iy2k2*100) c write(buff,11,err=91) tapec,(headi(j),j=1,3),blank1,iy2k1,iy2k2, > isy,(headi(j),j=5,6),iey,(headi(j),j=8,11), > (headr(j),j=1,4), > offset,scale,rnewzind, > blank6,label 11 format( a4, i3, i3, i1, a1, i2, i2, > i2, i3, i2, i2, i3, i2, > i4, i4, f6.2, f6.2, f7.3, f8.3, > e11.5, e11.5, e11.5, > a6, a40) c** write header and field data to file c ihrec = (igrd-1) * (headi(11)+1) + 1 write(102,rec=ihrec,err=92) buff(1:iheadl) c do i = 1,headi(11) write(102,rec=ihrec+i,err=93)(dataint(j,i),j=1,headi(10)) enddo c c istat = 0 close(unit=102) return c 90 write(6,*)'error opening output file ',flname(1:40) close(unit=102) pause stop90 91 write(6,*)'error writing data into ascii buffer' write(6,*)' chk 4 character year ' close(unit=102) pause stop91 92 write(6,*)'error writing output header, field number: ',igrd close(unit=102) pause stop92 93 write(6,*)'error writing output data, field number: ',igrd close(unit=102) pause stop93 end c c------------------------------------------------------------------------------------ c c NVAP DATA FILE HEADER and GRID FORMAT c c Each grid has it's own 144 byte ASCII header with either 1296 (1440-144) padded c bytes for the 1/2 degree data, or 576 (720-144) for 1 degree data. Following the c header the grid data there are either (720*360) or (360*180) 2 byte integers. c Multiply these integers by the header scale factor, and add the offset. (SEE c CODE ABOVE). Values that equal the ZINDEF value should not be scaled. Grid c value (1,1) is centered at either (0.25lon, 89.75lat) for 1/2 degree data or c (0.5lon, 89.5lat) for the 1 degree data. c c BYTE NUMBER c POSITION ofBytes TYPE c------------------------------------------------------------------------------------ c 1-4 (4) (A) Format code (CCDA = CIRA Climate Data Archive) c 5-7 (3) (I) HEADI(1) data source code c 8-10 (3) (I) HEADI(2) parameter number or code c 11-11 (1) (I) HEADI(3) data type (1=real*4, 2=int*4, 3=int*2, 4=byte) c non-VMS (5=real*4, 6=int*4, 7=int*2) c 12-12 (1) blank c 13-14 (2) (I) y2k start time century (either 19 or 20) if blank then defaults to 19 c 15-16 (2) (I) y2k end time century (either 19 or 20) if blank then defaults to 19 c c 17-18 (2) (I) HEADI(4) start year c 19-21 (3) (I) HEADI(5) start jul day c 22-23 (2) (I) HEADI(6) start hour c 24-25 (2) (I) HEADI(7) end year c 26-28 (3) (I) HEADI(8) end jul day c 29-30 (2) (I) HEADI(9) end hour c c 31-34 (4) (I) HEADI(10) xsize number horizontal data points c 35-38 (4) (I) HEADI(11) ysize number vertical data points c 39-44 (6) (F) HEADR(1) dx horizontal degree spacing (+ for W to E) c 45-50 (6) (F) HEADR(2) dy vertical degree spacing (- for N to S) c 51-57 (7) (F) HEADR(3) center gridbox latitude of first data location ( + for NH ) c 58-65 (8) (F) HEADR(4) center gridbox longitude of western most data location c 66-76 (11) (E) offset for data (add to data after applying scale) c 77-87 (11) (E) scale factor for data (multiply data by this) c 88-98 (11) (E) ZINDEF indefinite value c 99-104 (6) blanks c 105-144 (40) (A) LABEL character label (40 bytes) c c c--------------------------------------------------------------------------------------