program CRUts2
c
c	reads the new CRU global gridded dataset 1901-2000 0.5x0.5 degree
c	cells prepared by Tim Mitchell.
c	written by pjk-trl 4/2003
c
	integer    month,maxlon,maxlat,iu,ou1,mxyr,mnyr
	parameter (mxyr=2000,mnyr=1901,month=12,iu=13,ou1=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)
c
	logical ok1
c
	data fmt /'(12i5)'/
c
	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,''If you want land AND sea gridded data you must use the'')')
	write(*,'(1x,''CRU TS 1.0 data set which ends in 1995-96'')')
	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',readonly)
	open(unit=ou,file=outfl,status='unknown')
	ok1=.true.
c
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'
c
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
c
	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(string,'(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(string,'(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
c
	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
c
	indxlon(1)=-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
c
	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
c	        write(*,'(1x,''cell found'')')
c               write(*,'(1x,''gridx= '',i5,''gridy= '',i4)')gridx,gridy
c               write(*,'(1x,''indxlon(gridx)= '',f8.3,'' indxlat(gridy) = '',f8.3)')
c     *             indxlon(gridx),indxlat(gridy)
c	        write(*,'(1x,''yskip1 = '',i4)')yskip1
c	        write(*,'(1x,''yskip2 = '',i4)')yskip2
c	        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.1,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
c
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