#!/bin/ksh 

im=${1}
jm=${2}
lat=0.0
lon=-28.0
Res=${3}
resol=`echo ${Res} / 100 | bc -l | cut -c 1-3`
dlmd=`/gfs/home2/pilotto/Eta/model/scripts/def_eta.scr ${Res}|awk '{print $1}'`
dphd=`/gfs/home2/pilotto/Eta/model/scripts/def_eta.scr ${Res}|awk '{print $2}'`

cat <<EOF> corner.f90
       program CORNERS
!
!
!     *  ROUTINE TO FIND EARTH LATITUDE/LONGITUDE FOR THE CORNER       *
!     *               POINTS OF AN ETA MODEL GRID                      *
!

!-----------------------------------------------------------------------
!                             D I M E N S I O N
!C     & KHL0  (JM),KHH0  (JM), GLAT(IMJM),GLON(IMJM)

	REAL,ALLOCATABLE::GLAT(:),GLON(:)
	INTEGER,ALLOCATABLE::KHL0(:),KHH0(:)

        DATA PI/3.141592654/

        REAL DLMD,DPHD,WBD,SBD,TPH0D,TLM0D,minlat,wlon,maxlat,elon

!C*******************************


       im=${im}
       jm=${jm}
       IMJM=IM*JM-JM/2
       tph0d=${lat}
       tlm0d=${lon}
       dlmd=${dlmd}
       dphd=${dphd}
       res=0${resol}
        WBD=-(float(IM)-1.)*DLMD
        SBD=(-(float(JM)-1.)/2.)*DPHD

      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
	ALLOCATE(KHL0(JM),KHH0(JM))
	ALLOCATE(GLAT(IMJM),GLON(IMJM))

         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
        maxlat=-999.
	minlat=99.
              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)
          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

        if (glat(k) .gt. maxlat) maxlat=glat(k)
	if (glat(k) .lt. minlat) minlat=glat(k)

  200 CONTINUE

	DEALLOCATE(KHL0,KHH0)
	
	
	if (TPH0D .ge. 0) then
        wlon=glon(imjm-im+1)
        elon=glon(imjm)
	else
	wlon=glon(1)
	elon=glon(im)
	endif

!	write(6,*) 'raw lon values (w,e) ', wlon, elon
        if (tlm0d .lt. 0 .and. wlon .gt. 0) wlon=wlon-360.
        if (tlm0d .gt. 0 .and. elon .lt. 0) elon=elon+360.

        rnlat=maxlat+0.1
        slat=minlat-0.1
        west=wlon-0.1
        east=elon+0.1
        write(6,*) 'west= ', west
        write(6,*) 'east= ', east
        write(6,*) 'north= ', rnlat
        write(6,*) 'south= ', slat

        write(6,*) (abs(west-east))/res
	write(6,*) (abs(rnlat-slat))/res
	DEALLOCATE(GLAT,GLON)

      STOP
      END

EOF

f90 -o corners.x corner.f90
#corners.x > corners_${im}X${jm}.out

#rm -f corners.x
rm -f corner.f90

exit