       PROGRAM CORNERS
C
C     *  ROUTINE TO PRINT OUT EARTH LATITUDE/LONGITUDE FOR THE CORNER  *
C     *               POINTS OF AN ETA MODEL GRID                      *
C
                             P A R A M E T E R
     & (IM=61,JM=99
     &, IMJM=IM*JM-JM/2)
C-----------------------------------------------------------------------
	real*8 dlmd,dphd,wbd,sbd
                             D I M E N S I O N
     & KHL0  (JM),KHH0  (JM), GLAT(IMJM),GLON(IMJM)
                             D A T A
     & TLM0D/-100./,
     & TPH0D/40./,
     & PI/3.141592654/

	DLMD = .4
	DPHD = .4
C
C  set dlat for the output lat-lon grid (spacing in degrees lat and lon)
C
	dlat=0.50

C**************************************************************************
C**************************************************************************
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC 
CCC SHOULDN'T NEED TO CHANGE BELOW THIS POINT CCCCCCCCCC
CCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

	WBD=-(float(IM)-1.)*DLMD
	SBD=(-(float(JM)-1.)/2.)*DPHD

	write(6,*) ' '
	write(6,*) '**********************************************'
	write(6,*) ' '
	write(6,*) 'IM= ', IM, 'JM= ',JM
	write(6,*) 'DLMD= ', DLMD
	write(6,*) 'DPHD= ', DPHD
	write(6,*) 'WBD= ',WBD
	write(6,*) 'SBD= ',SBD
	write(6,*) 'grid centered at ',TPH0D, TLM0D

      DTR=PI/180.
      TPH0=TPH0D*DTR
      WB=WBD*DTR
      SB=SBD*DTR
      DLM=DLMD*DTR
      DPH=DPHD*DTR
      TDLM=DLM+DLM
      TDPH=DPH+DPH
C
      STPH0=SIN(TPH0)
      CTPH0=COS(TPH0)
C
         DO 100 J=1,JM
      KHL0(J)=IM*(J-1)-(J-1)/2+1
      KHH0(J)=IM*J-J/2
C     WRITE(6,9999) J, KHL0(J), KHH0(J)
C9999 FORMAT(2X,3(I10,1X))
  100 CONTINUE
C--------------GEOGRAPHIC LAT AND LONG OF TLL GRID POINTS---------------
              TPH=SB-DPH
              DO 200 J=1,JM
              KHL=KHL0(J)
              KHH=KHH0(J)
C
              TLM=WB-TDLM+MOD(J+1,2)*DLM
              TPH=TPH+DPH
              STPH=SIN(TPH)
              CTPH=COS(TPH)
C
          DO 200 K=KHL,KHH
      TLM=TLM+TDLM
      SPH=CTPH0*STPH+STPH0*CTPH*COS(TLM)
      GLAT(K)=ASIN(SPH)
      CLM=CTPH*COS(TLM)/(COS(GLAT(K))*CTPH0)-TAN(GLAT(K))*TAN(TPH0)
          IF(CLM.GT.1.)      CLM=1.
      FACT=1.
          IF(TLM.GT.0.)      FACT=-1.
      GLON(K)=(-TLM0D*DTR+FACT*ACOS(CLM))/DTR

Cmp	at this point GLON is in DEGREES WEST
	if (GLON(K) .lt. 0) GLON(K)=GLON(K)+360.
	if (GLON(K) .gt. 360.) GLON(K)=GLON(K)-360.
	if (GLON(K) .lt. 180) GLON(K)=-GLON(K)         ! make WH negative
	if (GLON(K) .gt. 180) GLON(K)=360.-GLON(K)     ! make EH 

	GLAT(K)=GLAT(K)/DTR

  200 CONTINUE

C
	write(6,1020) GLAT(1),GLON(1)
	write(6,1021) GLAT(IMJM),GLON(IMJM)
	write(6,*) ' '
	write(6,*) '**********************************************'
C	write(6,*) ' '

	wlon=int(glon(1))+1
	slat=int(TPH0D+SBD)+1
	delta=abs(TLM0D-wlon)
	elon=int(TLM0D+delta)-1
	nlat=int(GLAT(IMJM))-1
Cmp	nlat=int(TPH0D-SBD)-1

Cmp	search the lat values corresponding to wlon

	do K=1,imjm
	dif=wlon-glon(k)
	if (dif .lt. 1 .and. dif .gt. 0.) then
C	write(6,*) 'possible nlat value... ', glat(k)
	tnlat=glat(k)
	endif
	enddo

	nlat=int(tnlat)-1

	jdim=((nlat-slat)/dlat)+1
	idim=((elon-wlon)/dlat)+1
	islat=slat*1000
	iwlon=wlon*1000
	nlat=nlat*1000
	ielon=elon*1000
	idlat=(dlat*1000)
	open (unit=23,file='outjob_input_lat',form='formatted',
     +	access='sequential',status='unknown')
	write(23,299) IDIM,JDIM,islat,iwlon,nlat,ielon,idlat,idlat
  299   format('255 0 ',2(I3,x),I6,x,I7,x,'128 ',I6,x,I7,x,2(I3,x),'64')
        close(23)


Cmp	now compute the values needed for a lambert conic conformal 
Cmp	projection

	IDIM=IM*2-1
	JDIM=JM

	idx=INT(DPHD*111*1000)
	idy=idx

	LATONE=INT(GLAT(1)*1000)
	LONONE=INT(GLON(1)*1000)
	
		
	open (unit=24,file='outjob_input_lmbc',form='formatted',
     +	access='sequential',status='unknown')
	write(24,301) IDIM,JDIM,latone,lonone,int(TLM0D*1000),idx,idy,
     +	int(TPH0D*1000),INT(TPH0D*1000)

  301	format('255 3 ',2(I3,x),I6,x,I7,x,'8 ',I7,x,2(I6,x),'0 64',
     +	2(I6,x))
	
 1010 FORMAT(2X,I8,3(F12.5,1X))
 1015 FORMAT(2X,I8,5(F12.5,1X))
 1020 format(' LOWER LEFT  POINT= ',f10.5,x,f10.5)
 1021 FORMAT(' UPPER RIGHT POINT= ',F10.5,X,F10.5)
  
C
      STOP
      END
