program extrae
!	
!	reads the new CRU global gridded dataset 1901-2000 0.5x0.5 degree
!	cells prepared by Tim Mitchell.
!	written by pjk-trl 4/2003
!	
	integer    month,maxlon,maxlat,iu,ou,mxyr,mnyr
	parameter (mxyr=2000,mnyr=1901,month=12,iu=13,ou=14)
	parameter (maxlon=720,maxlat=360,gridskip=100)
	integer   yskip1,yskip2,yskip3
	integer   yrmin,yrmax,vartype,gridx,gridy,numgrid
	integer   ydata(month)
	integer   p,j,k
	real      indxlon(720),indxlat(360),lat,lon
	real      mxlon,mnlon,mxlat,mnlat
	character*27 infl,outfl,title
	character*11 dattype,fmt
	character*17 rfmt
	character*80 titles(5)
!	
	logical ok1
!	
	data fmt /'(12i5)'/
!	
	write(*,'(1x,''Program CRUts2b PJK-TRL-LDEO 5/2003'',/)')
	write(*,'(/1x,''  This program serves gridded 0.5x0.5 degree LAND AREA ONLY'')')
	write(*,'(1x,''gridded data from the most recent CRU TS 2.0 data set'')')
	write(*,'(1x,''This dataset is an update of the previous CRU TS 1.0 data'')')
	write(*,'(1x,''set which ends in 1995-96 and includes WET-DAY information.'')')
	write(*,'(1x,''  This data set was put together by Dr. Tim Mitchel CRU-UEA'')')
	write(*,'(1x,''The proper citation, at this time, for this data set is:'')')
	write(*,'(/1x,'' Mitchell, T.D., et al., 2003: A comprehensive set of climate'')')
	write(*,'(1x,''senarios for Europe and the globe. (in preparation)'')')	
	write(*,'(/,4x,''Select the Variable you want to extract:'',/)')
	write(*,'(/,1x,''    Variable Type                Name    Units     Length'')')  
	write(*,'(1x,''[1] Precipitation                pre     mm        1901-2000'')')
	write(*,'(1x,''[2] Mean temperature             tmp     degC*10   1901-2000'')')
	write(*,'(1x,''[3] Diurnal temperature range    dtr     degC*10   1901-2000'')')
	write(*,'(1x,''[4] Vapour pressure              vap     hPa*10    1901-2000'')')
	write(*,'(1x,''[5] Cloud cover                  cld     oktas*10  1901-2000'')')
	write(*,'(1x,''[0] Exit Program'')')
	write(*,'(/,1x,''Enter a value from 0 to 5 ==> '',$)')
	read(*,'(i1)')vartype
	select case (vartype)
		case(1)
		    dattype='CRUts2.pre.'
		    title='glo.pre. mm'
		    rfmt='RFMT=(1x,i4,12i5)'
		    infl='obs.globe.lan.1901-2000.pre'
		case(2)
		    dattype='CRUts2.tmp.'
		    title='glo.tmp. degC*10'
		    rfmt='RFMT=(1x,i4,12i5)'
		    infl='obs.globe.lan.1901-2000.tmp'
		case(3)
		    dattype='CRUts2.dtr.'
		    title='glo.dtr. degC*10'
		    rfmt='RFMT=(1x,i4,12i5)'
		    infl='obs.globe.lan.1901-2000.dtr'
		case(4)
		    dattype='CRUts2.vap.'
		    title='glo.vap. hPa*10'
		    rfmt='RFMT=(1x,i4,12i5)'
		    infl='obs.globe.lan.1901-2000.vap'
		case(5)
		    dattype='CRUts2.cld.'
		    title='glo.cld. oktas*10'
		    rfmt='RFMT=(1x,i4,12i5)'
		    infl='obs.globe.lan.1901-2000.cld'
		case(0)
		    write(*,'(1x,/''Program Terminating'')')
		    stop
	end select
	outfl=dattype//'TSF'

        open(unit=iu,file=infl,status='old')
	open(unit=ou,file=outfl,status='unknown')
	ok1=.true.
!
5	write(*,'(/,4x,''Select the Geographical Range of data to extract'')')
	write(*,'(/,2x,''  Enter the range of Latitude and Longitude'')')
	write(*,'(2x,''from which you would like to extract your data.'')')
	write(*,'(2x,''Enter values as REAL numbers to the nearest 0.5'')')
	write(*,'(2x,''degree. EXAMPLE a Latitude entered as -89.5 defines '')')
	write(*,'(2x,''cells bounded on the bottom by 90.0 degrees SOUTH.'')')
	write(*,'(2x,''Similarly a longitude entered as 180.0 are the cells '')')
	write(*,'(2x,''bounded to the right by -179.5 degrees WEST).'')')
	write(*,'(/,3x,''  All LONGITUDES WEST or [LEFT] of Greenwich  '')')
	write(*,'(2x,''  are NEGATIVE and EAST or [RIGHT] of Greenwich  '')')
	write(*,'(2x,''  are POSITIVE'')')
	write(*,'(2x,''    All LATITUDES SOUTH or [BELOW] the equator are'')')
	write(*,'(2x,''  NEGATIVE and NORTH latitudes [ABOVE] the equator are'')')
	write(*,'(2x,''  POSITIVE.'')')
	write(*,'(15x,'' THIS IS VERY IMPORTANT!'')')
	write(*,'(/,2x,''  Each grid cell value is the observation for that  '')')
	write(*,'(2x,''cells mid-point. Example: for the cell 90.0 Lat. -179.5'')')
	write(*,'(2x,''Lon. The actual grid point value given is centered '')')
	write(*,'(2x,''over (89.75,-179.75).''/)')
	Pause' HIT RETURN to continue'
!
10	write(*,'(8x,''Enter Max Latitude  [TOP]   ==> '',$)')
	read(*,'(f6.1)')mxlat
	if((mxlat.gt.90.0).or.(mxlat.lt.-90.0))then
		  write(*,'(1x,''Invalid latitude: '',f6.1)')mxlat
		  write(*,'(/,8x,''Valid range is 90.0 <-> -90.0'')')
		  go to 10
	endif
20	write(*,'(8x,''Enter Min Latitude  [BOTTOM]==> '',$)')
	read(*,'(f6.1)')mnlat
	if(mnlat.gt.mxlat)then
		  write(*,'(/,8x,''Your Min. Lat is greater than your Max.Lat'')')
		  write(*,'(8x,'' MAX. LAT =>  '',f6.1)')mxlat
		  write(*,'(8x,''Valid range is'',f6.1,'' -> -90.0'')')mxlat
		  go to 20
	elseif(mnlat.lt.-90.0)then
		  write(*,'(1x,''Invalid latitude: '',f6.1)')mnlat
		  write(*,'(/,8x,''Valid range is'',f6.1,'' -> -90.0'')')mxlat
		  go to 20
	endif
30	write(*,'(8x,''Enter Max Longitude [RIGHT] ==> '',$)')
	read(*,'(f6.1)')mxlon
	if((mxlon.gt.180.0).or.(mxlon.lt.-180.0))then
		write(*,'(1x,''Invalid longitude: '',f6.1)')mxlon
		write(*,'(/,8x,''Valid range is 180.0 <-> -180.0'')')
		go to 30
	endif
40	write(*,'(8x,''Enter Min Longitude [LEFT]  ==> '',$)')
	read(*,'(f6.1)')mnlon
	if((mnlon.gt.180.0).or.(mnlon.lt.-180.0))then
		  write(*,'(1x,''Invalid longitude: '',f6.1)')mnlon
		  write(*,'(/,8x,''Valid range is 180.0 <-> -180.0'')')
		  go to 40
	endif
	write(*,'(1x,''Values Entered are:'')')
	write(*,'(1x,''mxlon= '',f6.1,'' mnlon= '',f6.1,'' mxlat= '',f6.1,'' mnlat= '',f6.1)')mxlon,mnlon,mxlat,mnlat
	write(*,'(1x,''Is this ok? [T]rue or [F]alse ==> '',$)')
	read(*,*)ok1
	if(.not.ok1)goto 5
!	
	write(*,'(/1x,''Enter the temporal range to select data from'')')
	write(*,'(2x,''Valid years are '',i4, ''-'',i4)')mnyr,mxyr
50	write(*,'(/,2x,''Enter the First Year Of Data you want to extract.'')')
	write(*,'(2x,'' FYOD ==> '',$)')
	read(*,'(i4)')yrmin
	if((yrmin.lt.mnyr).or.(yrmin.gt.mxyr))then
		  write(*,'(/1x,''*WARNING* FYOD must be > '',i4,'' and < '',i4)')mnyr,mxyr
		  goto 50
	endif
60	write(*,'(/2x,''Enter the Last Year Of Data you want to extract.'')')
	write(*,'(2x,'' LYOD ==> '',$)')
	read(*,'(i4)')yrmax
	if(yrmax.lt.yrmin)then
		  write(*,'(1x,''*WARNING* LYOD must be greater than'',i4)')yrmin
		  write(*,'(2x,''Valid range = '',i4,''-'',i4 )')yrmin,mxyr
		  goto 60	
	elseif(yrmax.gt.mxyr)then
		  write(*,'(1x,''FYOD cannot exceed '',i4)')mxyr
		  write(*,'(/2x,''Valid range = '',i4,'' to '',i4)')yrmin,mxyr
		  goto 60
	endif
	write(*,'(1x,''Values Entered are:'')')
	write(*,'(1x,''yrmin= '',i5,'' yrmax= '',i5)')yrmin,yrmax
	write(*,'(1x,''Is this ok? [T]rue or [F]alse ==> '',$)')
	read(*,*)ok1
	if(.not.ok1)goto 60
!	
	yskip1=yrmin-mnyr       !number line to skip to first year of data
	yskip2=yrmax-yrmin+1     !number of lines to read data from
	yskip3=mxyr-yrmax       !number of lines to skip to end of cell
!	
	indxlon(1)=-179.5 !180.0	
	indxlat(1)=-90.0
	do j=2,maxlon
	   indxlon(j)=indxlon(j-1)+0.5
	end do
	do j=2,maxlat
	   indxlat(j)=indxlat(j-1)+0.5
	end do
	numgrid=0
!	
	do j=1,5
           read(iu,'(a)')titles(j)
	   write(ou,'(1x,a)')titles(j)
	end do
70      read(iu,'(9x,i4,1x,i4)',err=999,end=100)gridx,gridy
        if((indxlon(gridx).ge.mnlon).and.(indxlon(gridx).le.mxlon))then
		     if((indxlat(gridy).ge.mnlat).and.(indxlat(gridy).le.mxlat))then
		     	        numgrid=numgrid+1
!	        		write(*,'(1x,''cell found'')')
!              			write(*,'(1x,''gridx= '',i5,''gridy= '',i4)')gridx,gridy
!              			write(*,'(1x,''indxlon(gridx)= '',f8.3,'' indxlat(gridy) = '',f8.3)')
!     				*             indxlon(gridx),indxlat(gridy)
!	        		write(*,'(1x,''yskip1 = '',i4)')yskip1
!	        		write(*,'(1x,''yskip2 = '',i4)')yskip2
!	        		write(*,'(1x,''yskip3 = '',i4)')yskip3
		 		if(indxlon(gridx).lt.0.0)then
		 		   	lon=indxlon(gridx)+0.25
		 		else
		 			lon=indxlon(gridx)-0.25
				end if
				if(indxlat(gridy).lt.0.0)then
		 			lat=indxlat(gridy)+0.25
				else
		 			lat=indxlat(gridy)-0.25
				end if
			        write(ou,'(1x,2i4,2f8.2,2x,a,1x,a)')gridx,gridy,lon,lat,rfmt,title
		 		do p=1,yskip1
			    		read(iu,*,end=110)
	         		end do
		 		k=0
		 		do p=1,yskip2
		            		read(iu,fmt,err=998,end=120)(ydata(j),j=1,month)
			    		write(ou,'(1x,i4,1x,12i5)')yrmin+k,(ydata(j),j=1,month)
			    		k=k+1
		 		end do
		 		do p=1,yskip3
		        	    	read(iu,*,end=130)
	  	 		end do	
			else
	         		do p=1,gridskip
		            		read(iu,*,end=140)
                 		end do
		 		goto 70
			endif
	 else
	    do p=1,gridskip
	          read(iu,*,end=140)
	    end do
	    goto 70
	endif
	    goto 70
!
100	write(*,'(/1x,''PROGRAM TERMINATING NORMALLY'')')
	write(*,'(1x,''End Of File '',a,'' Reached'')')infl
	write(*,'(1x,i10,'' Grids Found'')')numgrid
	write(*,'(1x,''Output is in file '',a)')outfl
	close(iu)
	close(ou)
	stop
110	write(*,'(/1x,''End Of File '',a,'' Reached'')')infl
	write(*,'(1x,''Program Terminating '')')
	write(*,'(1x,''yskip1= '',i4)')yskip1
	close(iu)
	close(ou)
	stop
120	write(*,'(/1x,''End Of File '',a,'' Reached'')')infl
	write(*,'(1x,''Program Terminating '')')
	write(*,'(1x,''yskip2= '',i4)')yskip2
	close(iu)
	close(ou)
	stop
130	write(*,'(/1x,''End Of File '',a,'' Reached'')')infl
	write(*,'(1x,''Program Terminating '')')
	write(*,'(1x,''yskip3= '',i4)')yskip3
	close(iu)
	close(ou)
	stop
140	write(*,'(/1x,''End Of File '',a,'' Reached'')')infl
	write(*,'(1x,''Program Terminating '')')
	write(*,'(1x,''gridskip= '',i4)')gridskip
	close(iu)
	close(ou)
	stop
998	write(*,'(/1x,''Error Reading Grid Values'')')
	write(*,'(1x,''Program Terminating Without Completions'')')
	close(iu)
	close(ou)
	stop	           
999	write(*,'(/1x,''Error Reading Grid reference numbers'')')
	write(*,'(1x,''Program Terminating Without Completions'')')
	close(iu)
	close(ou)
	stop	           
	end