subroutine ogw_topo


	include 'com_topo.fc'

      parameter (itlon=1801,itlat=901) 
      dimension D(itlon,itlat),
     *	mu1(igit,nb),gamma(igit,nb),
     *	alon(itlon), itopo(itlon),alat(itlat),temp(5) 
	integer A
	real dEx1,dx1,dEy,dy1,MM,LL,dis
	real lon,lat,sum,num,sum2,D,itopo,alon,alat
	real mu1,temp,temps
	CHARACTER*20 fileint
	 
		alonl = 0.  
		alonr = 354.325   
		alatt = 87.5 
		alatb = -87.5

		dis= 0.2
	A=6400000.

 	igs1 = 180. / nb
	igs2 = 360. / igit 
	fileint = 'et2.dat'	 
      open ( 201,file=fileint)

!     compute the latitude and longitude points
      do i = 1,itlon
        alon(i) = (float(i)-1.0)*2.0/10.0 
      enddo
      do j = 1,itlat
        alat(j) = -90.0 + (float(j) -1.0)*2.0/10.0  
      enddo
      		               
!     compute boundary points      
	ilonl = nint((alonl + 0.)/igs2) + 1
      ilonr = nint((alonr + 0.)/igs2) + 1
      ilatb = nint((alatb + 90.)/igs1) !  + 1
      ilatt = nint((alatt + 90.)/igs1) ! + 1
!	reading topographic data
      do 45 j=1,itlat
       read(201,'(F15.7)')(itopo(i),i=1,itlon)
        
		do 45 i=1,itlon
      if(i.le.901) then  
       D(i+900,j)=itopo(i)
	else
	 D(i-901,j)=itopo(i)
!	D(i,j)=itopo(i)
	endif	
 45   continue 
 !	reading surface wind from ncep


!----------------------read latitude points in a row---------------------------------------     
	do 46 j=1,nb
	 lat=alatb+(j-ilatb)*igs1
!---------------------read longitudes to each latitude-------------------------------------
	do 47 i=1,igit
	  lon=alonl+(i-ilonl)*igs2
	
	  num=0
	  sum=0
	  sum2=0
	  sum4=0
	  sum5=0
		  dy1=A*dis*3.1415926535/180
		  dx1=dy1*cos(lat*3.1415926535/180)

	do 48 k=2,itlat-1
	if (alat(k).gt.(lat+igs1/2.).or.alat(k).lt.(lat-igs1/2.)) goto 48
      do 49 l=2,itlon-1
	if (alon(l).gt.(lon+igs2/2.).or.alon(l).lt.(lon-igs2/2.)) goto 49
  
!---------computing average elevations and dephts ----------------
	   sum=sum+D(l,k)
	   num=num+1
	   sum5=sum5+D(l,k)*D(l,k)

!-------------computing dE/dx & dE/dy1------------
			dEy=(D(l,k+1)-D(l,k-1))/(2*dy1)
              dEx1=(D(l+1,k)-D(l-1,k))/(2.*dx1)
		   !-------------end computing----------------------
		   LL=(dEx1**2.-dEy**2.)/2.
		   MM=dEx1*dEy
			   sum2=sum2+MM
			   sum4=sum4+LL
	
 49	continue
 48	continue
!-----------everage koeffitients--------------------------------------
           E(i,j)=sum/num
!	print*,E(i,j),sum,num
	       mu1(i,j)=sqrt(sum5/num-(sum/num)**2)
	       mu1(i,j)=mu1(i,j) + E(i,j)/2.
           MM=sum2/num
		   LL=sum4/num
       if (LL.eq.0.) teta(i,j)=0. 
       if (LL.ne.0.) teta(i,j)=atan(MM/LL)/2.
		dxx=dx1*cos(teta(i,j))+dy1*sin(teta(i,j))
		dyy=dy1*cos(teta(i,j))-dx1*sin(teta(i,j))

	
!	write(211,*) mu(i,j)
!---------------comuting sigma & gamma---------------------------------	 
	 num1=0; sum3=0; sum4=0

	do 52 k=2,itlat-1
	if (alat(k).gt.(lat+igs1/2.).or.alat(k).lt.(lat-igs1/2.)) go to 52
      do 53 l=2,itlon-1
	if (alon(l).gt.(lon+igs2/2.).or.alon(l).lt.(lon-igs2/2.)) go to 53
	   
	   num1=num1+1
	   dExx=(D(l+1,k)-D(l-1,k))/(2*dxx)
	   dEyy=(D(l,k+1)-D(l,k-1))/(2*dyy)
       sum3=sum3+dExx*dExx
	   sum4=sum4+dEyy*dEyy

 53	continue
 52	continue
           sigma(i,j)=sqrt(sum3/num1)
		   if (sum3.eq.0.) then
		   gamma(i,j)=0
		   else
		   gamma(i,j)=sqrt(sum4/sum3)
		   end if

!---------------end comuting sigma & gamma---------------------------------	 

!---------------------------------------------------------------------------
      B(i,j)=1-0.18*gamma(i,j)-0.04*(gamma(i,j))**2 
      C(i,j)=0.48*gamma(i,j)+0.3*(gamma(i,j))**2
!	IF(B(i,j).gt.0) then 
!	write(211,*) i,j,'B=',B(i,j),'C=',C(i,j),'top'
!	ENDIF 
!	write(211,*)E(i,j)

 47	continue
 46	continue
       do i=1,5 
	temp(i) = 0.
       end do 
	do 557 j=1,nb
	do 558 i=1,igit
	temp = 0 
	count1 = 3 
	count2 = 2

		if(i.eq.1) then 
 			temp(1)=mu1(i,j)
			temp(2)=mu1(igit,j)
			temp(3)=mu1(i+1,j)
		elseif(i.eq.igit) then
 			temp(1)=mu1(i,j)
			temp(2)=mu1(i-1,j)
			temp(3)=mu1(1,j)
		else
 			temp(1)=mu1(i,j)
			temp(2)=mu1(i-1,j)
			temp(3)=mu1(i+1,j)
		endif
	if(j.ne.1.and.j.ne.nb)	then
 		if(i.eq.1) then 
			temp(4)=mu1(i,j-1)
			temp(5)=mu1(i,j+1)
		elseif(i.eq.igit) then
			temp(4)=mu1(i,j-1)
			temp(5)=mu1(i,j+1)
		else
			temp(4)=mu1(i,j-1)
			temp(5)=mu1(i,j+1)
		endif
	count1 = 5
	count2 = 3
	endif
	if(j.eq.1)	then
 		if(i.eq.1) then 
			temp(4)=mu1(i,2)
		elseif(i.eq.igit) then
			temp(4)=mu1(i,2)
		else
			temp(4)=mu1(i,2)
		endif
	count1 = 4
	count2 = 2
	endif
	if(j.eq.nb)	then
 		if(i.eq.1) then 
			temp(4)=mu1(i,j-1)
		elseif(i.eq.igit) then
			temp(4)=mu1(i,j-1)
		else
			temp(4)=mu1(i,j-1)
		endif
	count1 = 4
	count2 = 2
	endif

      pass = 1
 1001 continue
      sorted = 1
      do 1005 l = 1,count1-pass
        if(temp(l) .gt. temp(l+1)) then
          temps = temp(l)
          temp(l) = temp(l+1)
          temp(l+1) = temps
          sorted = 0
        endif
 1005 continue
      pass = pass +1
      if(sorted .eq. 0) goto 1001 
	mu(i,j) =  temp(count2)
!	write(211,*) 'mu=',mu(i,j),i,j
 558	continue
 557	continue

c	print*,'topo!!!!!!!!!!!!!!!!!!!!!'
	end subroutine ogw_topo