	program ERsstv3_8
c	Extended Reconstructed Sea Surface Temperatures, version 2 (ERSST.v3), beginning January
c	1854 to December 2008.
c	pjk, trl-ldeo 1/2008
c
	integer iyrbeg,iyrend,ou,ou2,iu,iu2,mi,mj,mxyr,month,mxfl,land,fyrmon
	real    di,dj,aln1,alt1
	parameter(iu=23,iu2=25,ou=24,ou2=26,mxyr=300)
	parameter(mi=180,mj=89,di=2.0,dj=2.0,aln1=0.0,alt1=-88.0,land=-9999)
	parameter(month=12)
	integer p,j,k,i,l,mm,m,wcnt
	integer mon,yr,styr,stpyr,yrnum,yrskip,lyod,lmod,tlyod,tmod,imonth
	integer iyw(mi,mj,mxyr*month)
	integer latskp,strlat,stplat,latnum,tmonths,datype
	integer strlon1,strlon2,stplon1,stplon2,lonskp1,lonskp2
	real    tlat(mj),tlon(mi),mnlat,mxlat,mnlon,mxlon
	character*40 outfl1,outfl2,rfmt,infl1,infl2(mxyr)
	character*60 head1,head2
	character*60 listfl
	character*3 cmonth(12),tmonth1,tmonth2
	
	logical flok(mxyr),flfail,ok,split
	data cmonth/'JAN','FEB','MAR','APR','MAY','JUN',
     * 'JUL','AUG','SEP','OCT','NOV','DEC'/
	data outfl1 /'ERsst.v5b.tsf'/
	data infl1 /'ERsst.dat'/
	data head1 /'ERsst Anomalies degrees C*100 RFMT=(i5,12f6.2)'/
	data head2 /'ERsst ESD degrees C RFMT=(i5,12f6.2)'/
c      

ccc
	write(*,'(//1x,''Program ERsst SERVER *************************************************'')')
	write(*,'(1x,''ERSST data are an improved extended reconstruction of SST.'')')
	write(*,'(/1x,''All fields are computed on a 2-degree grid. The grid has 180 zonal points and 89 '')')
	write(*,'(1x,''meridional points. The 2-degree cells are centered on 0E, 2E, ...,and 88S,86S .... '')')
	write(*,'(1x,''Over land a value of -99.99 is assigned.  In addition any SST mean value less than '')')
	write(*,'(1x,''-1.8C was set to -1.8C.'')')
	write(*,'(//1x,''Program ERsst will prompt you for a range of years and a geographic region'')')
	write(*,'(1x,''from which to extract monthly, grid-cell data.  The resulting grid cell infor- '')')
	write(*,'(1x,''mation is saved in a new file called ERsst.TSF. This output file is REPLACED  '')') 
	write(*,'(1x,''with each run of the program.'')') 
	write(*,'(/1x,''The proper citation for these data is:'')')
	write(*,'(/1x,''  Smith, T. M., and R. W. Reynolds, 2004: Improved extended reconstruction of '')')
	write(*,'(1x,''         SST (1854-1997). J. Climate, 17, 2466-2477'')')
	write(*,'(1x,''pjk trl-ldeo Jan.2008; update May, 2019'')')
	write(*,'(1x,''_____________________________________________________________________________''/)')
	pause' Hit RETURN to continue'
c
c
c check to see if ERsst.v3b.1854-2009.dat exists and if not try to build it.
	datype=1
	
	flok(1)=.false.
	inquire(file=infl1,exist=flok(1))
c
	if(flok(1))then
	   write(*,'(/2x,''> Found file: '',a,'' <'')')infl1
	   open(unit=iu,file=infl1,status='old',readonly)
	   read(iu,*)
	   read(iu,'(29x,a3,1x,i4,4x,a3,1x,i4)')tmonth1,iyrbeg,tmonth2,lyod
	   do i=1,month
	      if(tmonth1.eq.cmonth(i))imonth=i
	      if(tmonth2.eq.cmonth(i))lmod=i
	   enddo
	   tmod=lmod;tlyod=lyod
	   rewind(iu)
	   goto 10
	else
	  write(*,'(/1x,''Program ERsst needs to build a composite data file from the annual source files.  '')')
	  write(*,'(/1x,''The first time the program is opened it will search for these files then open '')')
	  write(*,'(1x,''each one and combine them into one large file of monthly sst data called '')')
	  write(*,'(1x,''ERsst.dat .   '')')
	  write(*,'(1x,''After creating this new sst data file the original decadal files may be deleted.  '')')
	  write(*,'(1x,''For all future runs of the program only the file ERsst.dat is required.'')')
	  write(*,'(1x,''To concatonate all the annual ERSST grid files into one global data file for all years'')')
	  write(*,'(/1x,''1) Create a file containing the names of all the annual files each on a seperate line.'')')
	  write(*,'(3x,''*** NOTE :> The column of file names must be written in accending order by year. ***'')')
	  write(*,'(/1x,''2) Launch this program a second time and when you get to the following prompt enter the '')')
	  write(*,'(1x,''list file name. The program will concatonate all the annual files into one.'')')
	  write(*,'(1x,''This process may take some time depending on your machine. You may want to go out and '')')
	  write(*,'(1x,''get a coffee.'')')
	  write(*,'(1x,''Once done you will not have to do this until adding an annual update.'',/)')
	  write(*,'(1x,''Enter [1] to CONTINUE'')')
	  write(*,'(1x,''Enter [2] to QUIT   > '',$)')
	  read(*,*)i
	  if(i.eq.2)then
	     write(*,'(1x,''Program Terminating'')')
	     stop
	  else
c	    read list file names store in infl3
	    write(*,'(/1x,''Enter ERsst list file name > '',$)')
	    read(*,*)listfl
	    open(unit=iu,file=listfl,status='old')
	    mxfl=1
1	    read(iu,*,end=2,err=199)infl2(mxfl)
              mxfl=mxfl+1
           goto 1
           
	  endif
c
2	  close(iu); flfail=.false.;mxfl=mxfl-1

	  write(*,'(/1x,''Building ERsst.dat file ________________________________________'')')
	  do p=1,mxfl
	     flok(p)=.false.
	     write(*,'(1x,''Searching for file : '',a)')infl2(p)
	     inquire(file=infl2(p),exist=flok(p))
	     if(.not.flok(p))then
	       write(*,'(6x,''WARNING: File > '',a,1x,''<*DOES NOT EXIST*>'')')infl2(p)
              flfail=.true.
	     else
	       write(*,'(15x,''File : '',a,'' *FOUND*'',)')infl2(p)
	       if(p.eq.1)read(infl2(p),'(9x,i4)')iyrbeg
	     endif
	  end do
c	  
	  if(flfail)then
	      write(*,'(/1x,''You are missing one or more of the listed ERsst data files. Goto'')')
	      write(*,'(1x,''http://www.jisao.washington.edu/data_sets/ersst/index3.html#data download '')')
	      write(*,'(1x,''the all the neccessary ersst files, decompress them, and start over. '',/)')
	      write(*,'(/1x,''Program Terminating'')')
	      stop
	  else
	      write(*,'(/1x,''The program will now try to make a single composite data file from the '',i2)')mxfl
	      write(*,'(1x,''decadal files for further use.  This may take some time so relax.'',/)')
c concatonate all annual files and save data in  ERsst.dat
	      write(*,*)' '
	      lyod=iyrbeg
	      mon=1
	      do p=1,mxfl  !decade files
	        open(unit=iu2,file=infl2(p),status='old',readonly)
	        yrnum=10
	        if(p.eq.1)yrnum=6  !first file has only 6 years of data
	        do j=1,yrnum   !years
	          do k=1,month   !month
	            lmod=k
	            do l=1,mi
	              read(iu2,'(90i6)',end=5,err=6)(iyw(l,mm,mon),mm=1,mj)
	            end do
	            mon=mon+1
	          end do
	          lyod=lyod+1
	        end do
5	        close(iu2)
	        write(*,'(1x,''Processed file : '',a)')infl2(p)
	      end do
	      mon=mon-1
	      lmod=lmod-1
	      if(lmod.eq.0)then
	        tmod=12
	        tlyod=lyod-1
	      else
	        tmod=lmod
	        tlyod=lyod
	      endif
	      open(unit=iu,file=infl1,status='unknown')
	      write(iu,'(1x,''Extended Reconstructed Sea Surface Temperatures'')')
	      write(iu,'(1x,''Version 5 (ERSST.v5b), from '',a3,1x,i4,'' to '',a3,1x,i4)')cmonth(1),iyrbeg,cmonth(tmod),tlyod
	      write(iu,'(1x,''High-frequency SST anomalies *100'')')
	      l=iyrbeg
	      mm=1
	      do k=1,mon
	        write(iu,'(1x,2i4)')l,mm
	        do p=1,mi
	            write(iu,'(90i6)')(iyw(p,j,k),j=1,mj)
	        end do
	        if(mod(k,12).eq.0)then
	          l=l+1
	          mm=1
	        else
	          mm=mm+1
	        endif
	      end do
	      close(iu)
c
	       open(unit=iu,file=infl1,status='old',readonly)
	       write(*,'(/1x,''The program has sucessfully created file '',a30)')infl1
	      write(*,'(/1x,''The time span of these data is '',i4,'' to '',i2,''/'',i4)')iyrbeg,tmod,tlyod
	      write(*,'(/1x,''After this run you may delete the original the ERsst annual files.'')') 
	      write(*,'(1x,''**Note**: With every annual update you will need to repeat this entire '')')
	      write(*,'(1x,''process requiring these '',i4,'' files in addition to the update.'')')mxfl
	      goto 10
c
6	      write(*,'(/1x,''Read Error in file '',a,2i4,'' line '',i4)')infl2(p),lyod,k,l
	      pause' Hit Return to STOP'
	      stop
            endif
	endif
c	
10	write(*,'(/1x,''Select years to extract __________________________________________________'')')
	write(*,'(/3x,''Valid years are '',i4,'' and '',i4)')iyrbeg,tlyod
	write(*,'(/3x,''Enter first year of data to extract  : '',$)')
	read(*,'(i4)',err=10)styr
	if(styr.lt.iyrbeg.or.styr.gt.tlyod)then
	   write(*,'(3x,''Wooah..Snowy! between '',i4,'' and '',i4)')iyrbeg,tlyod
	   pause' Hit RETURN to continue'
	   goto 10
	end if
20	write(*,'(/3x,''Enter last year of data to extract   : '',$)')
	read(*,'(i4)',err=20)stpyr
	if(stpyr.lt.styr.or.stpyr.gt.tlyod)then
	   write(*,'(3x,''Wooah..Snowy! Between '',i4,'' and '',i4)')styr,tlyod
	   pause'   Hit RETURN to continue'
	   goto 20
	end if
c
30	write(*,'(/1x,''Select geographic region to extract_______________________________________'',/)')
	write(*,'(1x,''  This data set is spatially shaped like a rectangle where the'')')
	write(*,'(3x,''top left corner grid cell, mid-point is defined by the coordinates'')')
	write(*,'(3x,''0.0E Lon./90.0N Lat. and the bottom right coordinates are'')')
	write(*,'(3x,''-1.0W Lon./-89.0 Lat.  (see diagram below)'')')
	write(*,'(//8x,''  0.0E/88.0N         178.0E|-180.0W         -2.0W/88.0N'')')
	write(*,'(8x,''      :--------------------|---------------------:'')')
	write(*,'(8x,''      :                    |                     :'')')
	write(*,'(8x,''      :                    |                     :'')')
	write(*,'(8x,''      :                    |                     :'')')
	write(*,'(8x,''      :            Our Lovely Planet             :'')')
	write(*,'(8x,''      :-----           -equator-            -----:'')')
	write(*,'(8x,''      :                    |                     :'')')
	write(*,'(8x,''      :                    |                     :'')')
	write(*,'(8x,''      :                    |                     :'')')
	write(*,'(8x,''      :                    |                     :'')')
	write(*,'(8x,''      :---------------(dateline)-----------------:'')')
	write(*,'(8x,''  0.0E/-88.0S                               -2.0W/-88.0S'')')
	write(*,'(/3x,''All latitudes SOUTH of the equator (0.0) are negative and NORTH'')')
	write(*,'(3x,''of the equator are positive. All longitudes WEST of the dateline'')')
	write(*,'(3x,''are positive and EAST of the dateline are negative.'')')
	write(*,'(/3x,''Minumum latitudes (bottom) are values most negative and Maximum '')')
	write(*,'(3x,''Latitudes (top) are values most positve.'')')
	write(*,'(3x,''Maximum longitudes are values most EAST (right) and MAY be negative.'')')
	write(*,'(3x,''Minimum longitudes are values most WEST (left).'')')
	write(*,'(3x,''NOTE: All valus must be entered as real numbers with'')')
	write(*,'(3x,''their decimal fraction e.g., -100.25 , 151.0 or 30.5, -40. , 21.'')')
c
35	write(*,'(/13x,''Enter MAX. Latutude  [most north]    TOP: '',$)')
	read(*,'(f8.1)',err=35)mxlat
40	write(*,'(13x,''Enter MIN. Latitude  [most south] BOTTOM: '',$)')
	read(*,'(f8.1)',err=40)mnlat
	if(mxlat.lt.mnlat)then
	   write(*,'(/13x,''MAX. LAT. must be greater than MIN. LAT.'')')
	   write(*,'(13x,''MIN. LAT. entered = '',f8.2,'' MAX. LAT. = '',f8.2)')mnlat,mxlat
	   pause' Hit RETURN to continue'
	   goto 35
	 endif
45	write(*,'(/13x,''Enter MIN. Longitude  [most west]   LEFT: '',$)')
	read(*,'(f8.1)',err=45)mnlon
50	write(*,'(13x,''Enter Max. Longitude  [most east]  RIGHT: '',$)')
	read(*,'(f8.1)',err=50)mxlon
	write(*,'(/13x,''Area selected is :'')')
	 write(*,'(/4x,''    Lon.    Lat                                      Lon.     Lat.'')')
	 write(*,'(4x,f8.2,''/'',f7.2,''  :------------------------------:'',f7.2,''/'',f6.2)')
     *           mnlon,mxlat,mxlon,mxlat
	 write(*,'(22x,'':                              :'')')
	 write(*,'(22x,'':                              :'')')
	 write(*,'(22x,'':        Your Selection        :'')')
	 write(*,'(22x,'':                              :'')')
	 write(*,'(22x,'':                              :'')')
	 write(*,'(4x,f8.2,''/'',f7.2,''  :------------------------------:'',f7.2,''/'',f6.2)')
     *           mnlon,mnlat,mxlon,mnlat
	 write(*,'(/13x,''Is this ok? Enter [T]rue or [F]alse : '',$)')
	 read(*,*)ok
	 if(.not.ok)then
	   goto 30
	 endif
c make arrays of lats and lons
	tlon(1)=aln1
	do i=2,90
	  tlon(i)=tlon(i-1)+di
	end do
	tlon(91)=-180.0
	i=91
	do i=92,mi
	  tlon(i)=tlon(i-1)+di
	end do
c
	tlat(1)=alt1
	do j=2,mj
	  tlat(j)=tlat(j-1)+dj
	end do
c prepare write array
	do p=1,mi
	  do j=1,mj
	    do k=1,(iyrend-iyrbeg+1)*12
	      iyw(p,j,k)=-9999
	    end do
	  end do
	end do
c compute year skp and year num
	yrnum=stpyr-styr+1
	yrskip=styr-iyrbeg
	if((stpyr.eq.lyod).and.(lmod.ne.12))then
	  tmonths=((yrnum-1)*12)+lmod  !total month with partial last year
	else
	  tmonths=yrnum*12
	endif
c	
c compute latitude bounds indices
	split=.false.
	strlon1=0;stplon1=0;strlon2=0;stplon2=0;lonskp1=0;lonskp2=0
	latskp=0;strlat=0;stplat=0;latnum=0
	do p=1,mj
	  if((mnlat.ge.tlat(p)).and.(mnlat.le.tlat(p+1)))then
	    latskp=p-1
	    strlat=p
	  endif
c	  if((mxlat.gt.tlat(p)).and.(mxlat.le.tlat(p+1)))then
	  if((mxlat.ge.tlat(p)).and.(mxlat.le.tlat(p+1)))then
	    stplat=p   !+1
	  endif
	end do
	latnum=(stplat-strlat)+1
c	write(*,'(/1x,''mnlat mxlat latskp strlat stplat latnum'')')
c	write(*,*)mnlat,mxlat,latskp,strlat,stplat,latnum
c	write(*,'(1x,''tlat(strlat) tlat(stplat)'')')
c	write(*,*)tlat(strlat),tlat(stplat)
c
c compute longitude bounds indices
	if(mnlon.lt.0.0.and.mxlon.ge.0.0)then  !case -mnlon and +mxlon split=true span prime meridian
	   split=.true.
	   strlon1=1
	   stplon2=mi
	   do p=1,90
	     if((mxlon.gt.tlon(p)).and.(mxlon.le.tlon(p+1)))stplon1=p+1
	   end do
	   do p=91,mi
	     if((mnlon.ge.tlon(p)).and.(mnlon.lt.tlon(p+1)))strlon2=p
	   end do
	   lonskp1=(strlon2-stplon1)-1
c	   write(*,'(/1x,''split strlon1 stplon1 strlon2 stplon2 lonskp1'')')
c	   write(*,*)split ,strlon1, stplon1 ,strlon2 ,stplon2 ,lonskp1
c	   write(*,'(1x,''tlon(strlon1) tlon(stplon1) tlon(strlon2) tlon(stplon2) '')')
c	   write(*,*)tlon(strlon1), tlon(stplon1), tlon(strlon2), tlon(stplon2)
	elseif((mnlon.ge.0.0).and.(mxlon.lt.0.0))then  !case +mnlon and -mxlon
	   do p=1,90
	     if((mnlon.ge.tlon(p)).and.(mnlon.le.tlon(p+1)))strlon1=p
	   end do
	   do p=91,mi
	     if((mxlon.gt.tlon(p)).and.(mxlon.le.tlon(p+1)))stplon1=p+1
	   end do
	   lonskp1=strlon1-1
	   lonskp2=mi-stplon1
c	   write(*,'(/1x,''split strlon1 stplon1 strlon2 stplon2 lonskp1 lonskp2'')')
c	   write(*,*)split ,strlon1, stplon1 ,strlon2 ,stplon2 ,lonskp1,lonskp2
c	   write(*,'(1x,''tlon(strlon1) tlon(stplon1) tlon(strlon2) tlon(stplon2) '')')
c	   write(*,*)tlon(strlon1), tlon(stplon1), tlon(strlon2), tlon(stplon2)
	elseif((mnlon.ge.0.0).and.(mxlon.ge.0.0))then   !case both positive
	   do p=1,90
	     if((mnlon.ge.tlon(p)).and.(mnlon.le.tlon(p+1)))strlon1=p
	     if((mxlon.ge.tlon(p)).and.(mxlon.le.tlon(p+1)))stplon1=p
	   end do
	   lonskp1=strlon1-1
	   lonskp2=mi-stplon1
c	   write(*,'(/1x,''split strlon1 stplon1 strlon2 stplon2 lonskp1 lonskp2'')')
c	   write(*,*)split ,strlon1, stplon1 ,strlon2 ,stplon2 ,lonskp1,lonskp2
c	   write(*,'(1x,''tlon(strlon1) tlon(stplon1) tlon(strlon2) tlon(stplon2) '')')
c	   write(*,*)tlon(strlon1), tlon(stplon1), tlon(strlon2), tlon(stplon2)
	elseif((mnlon.lt.0.0).and.(mxlon.lt.0.0))then   !case both negative
	   do p=91,mi
	     if((mnlon.ge.tlon(p)).and.(mnlon.lt.tlon(p+1)))strlon1=p
	     if((mxlon.gt.tlon(p)).and.(mxlon.le.tlon(p+1)))stplon1=p+1
	   end do
	   lonskp1=strlon1-1
	   lonskp2=mi-stplon1
c	   write(*,'(/1x,''split strlon1 stplon1 strlon2 stplon2 lonskp1 lonskp2'')')
c	   write(*,*)split ,strlon1, stplon1 ,strlon2 ,stplon2 ,lonskp1,lonskp2
c	   write(*,'(1x,''tlon(strlon1) tlon(stplon1) tlon(strlon2) tlon(stplon2) '')')
c	   write(*,*)tlon(strlon1), tlon(stplon1), tlon(strlon2), tlon(stplon2)
	endif
c
	rfmt='(   x,   (i6))'
	write(rfmt(2:4),'(i3)')latskp*6
	write(rfmt(7:9),'(i3)')latnum
c	write(*,*)rfmt
c	pause'before reading'
c begin reading from file
c  skip down to proper year
	do p=1,3
	  read(iu,*)
	end do
	do p=1,yrskip*month
	  read(iu,'(1x,2i4)')yr,mon
c	  write(*,'(1x,''skipping ''2i4)')yr,mon
	  do j=1,mi
	    read(iu,*)
	  end do
	end do
c start big loop to read monthly ersst data
	do k=1,tmonths   !yrnum*month
	    read(iu,'(1x,2i4)',err=195)yr,mon
c	    write(*,'(1x,''reading '',2i4)')yr,mon
	    if(split)then
	      do p=strlon1,stplon1
	        read(iu,rfmt,end=150,err=196)(iyw(p,j,k),j=strlat,stplat)
c	        write(*,rfmt)(iyw(p,j,k),j=strlat,stplat)
	      end do
	      do p=1,lonskp1   ! skip to next longitude block to read negative longitudes west of 0.0
	        read(iu,*)
	      end do
	      do p=strlon2,stplon2
	        read(iu,rfmt,end=150,err=197)(iyw(p,j,k),j=strlat,stplat)
c	        write(*,rfmt)(iyw(p,j,k),j=strlat,stplat)
	      end do
	    else
	      do p=1,lonskp1   ! if lonskp1 .gt.0 then skip down to proper longitude if required
	        read(iu,*)
	      end do
	      do p=strlon1,stplon1
	        read(iu,rfmt,end=150,err=198)(iyw(p,j,k),j=strlat,stplat)
c	        write(*,rfmt)(iyw(p,j,k),j=strlat,stplat)
	      end do
	      do p=1,lonskp2   !if lonskp2 .gt.0 skip down to next monthly block if required
	        read(iu,*)
	      end do
	    endif
	end do
c
c write retrieved data to file.
c   open output files
	open(unit=ou2,file='ERsst_v3b.xy',status='unknown')
	if(datype.eq.1)then
	  open(unit=ou,file=outfl1,status='unknown')
	  write(ou,'(1x,''Extended Reconstructed Sea Surface Temperatures'')')
	  write(ou,'(1x,''Version 3 (ERSST.v3b). '',i4,'' to '',i4)')styr,stpyr
	  write(ou,'(1x,''High-frequency SST anomalies *100'')')
	elseif(datype.eq.2)then
	  open(unit=ou,file=outfl2,status='unknown')
	  write(ou,'(1x,''Extended Reconstructed Sea Surface Temperatures'')')
	  write(ou,'(1x,''Version 3 (ERSST.v3b). '',i4,'' to '',i4)')styr,stpyr
	  write(ou,'(1x,''Error Standard Deviation degrees C'')')
	endif
c
	write(*,'(/1x,''Writing Selected Grid Cell data _______________________________________'')')
	wcnt=1
	do j=strlat,stplat
	  if(.not.split)then
	    do p=strlon1,stplon1
	       if(iyw(p,j,1).eq.land)goto 100
	          if(datype.eq.1)then
	             write(ou,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head1
	             write(ou2,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head1
	             write(*,'(1x,i5,1x,2f7.2,1x,a)')wcnt,tlon(p),tlat(j),head1
	          elseif(datype.eq.2)then
	             write(ou,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head2
	             write(ou2,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head2
	             write(*,'(1x,i5,1x,2f7.2,1x,a)')wcnt,tlon(p),tlat(j),head2
	          endif
	          m=1
	          mm=12
	          do i=styr,stpyr
	             write(ou,'(i4,1x,12i6)')i,(iyw(p,j,k), k=m,mm)
	             m=mm+1
	             mm=m+11
	          end do
	          wcnt=wcnt+1
100	       continue	       
	    end do
	  else
	    do p=strlon1,stplon1
	       if(iyw(p,j,1).eq.land)goto 110
	          if(datype.eq.1)then
	             write(ou,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head1
	             write(ou2,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head1
	             write(*,'(1x,i5,1x,2f7.2,1x,a)')wcnt,tlon(p),tlat(j),head1
	          elseif(datype.eq.2)then
	             write(ou,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head2
	             write(ou2,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head2
	             write(*,'(1x,i5,1x,2f7.2,1x,a)')wcnt,tlon(p),tlat(j),head2
	          endif
	          m=1
	          mm=12
	          do i=styr,stpyr
	             write(ou,'(i4,1x,12i6)')i,(iyw(p,j,k), k=m,mm)
	             m=mm+1
	             mm=m+11
	          end do
	          wcnt=wcnt+1
110		continue
	    end do
	    do p=strlon2,stplon2
	       if(iyw(p,j,1).eq.land)goto 120
	          if(datype.eq.1)then
	             write(ou,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head1
	             write(ou2,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head1
	             write(*,'(1x,i5,1x,2f7.2,1x,a)')wcnt,tlon(p),tlat(j),head1
	          elseif(datype.eq.2)then
	             write(ou,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head2
	             write(ou2,'(1x,2f7.2,1x,a)')tlon(p),tlat(j),head2
	             write(*,'(1x,i5,1x,2f7.2,1x,a)')wcnt,tlon(p),tlat(j),head2
	          endif
	          m=1
	          mm=12
	          do i=styr,stpyr
	             write(ou,'(i4,1x,12i6)')i,(iyw(p,j,k), k=m,mm)
	             m=mm+1
	             mm=m+11
	          end do
	          wcnt=wcnt+1
120		continue	       
	    end do
	  endif
	end do
	write(*,'(1x,''_______________________________________________________________________'',/)')
	wcnt=wcnt-1
c
c
150	close(iu)
	close(ou)
	close(ou2)
	write(*,'(/1x,''Total number of grid cells found = '',i5)')wcnt
	write(*,'(1x,''Normal Termination of program'',/)')
	pause' Hit RETURN to quit'
	stop
195    write(*,'(/1x,''Error reading year and month'')')
	pause' Hit RETURN to stop'
	stop
196    write(*,'(/1x,''Error reading data 1'')')
	pause' Hit RETURN to stop'
	stop
197    write(*,'(/1x,''Error reading data 2'')')
	pause' Hit RETURN to stop'
	stop
198    write(*,'(/1x,''Error reading data 3'')')
	pause' Hit RETURN to stop'
	stop
199    write(*,'(/1x,''Error reading listfile'')')
	pause' Hit RETURN to stop'
	stop
	end
	
