orog_utils.F90

Go to the documentation of this file.

 module orog_utils

 implicit none

 private

 real, parameter    :: earth_radius = 6371200.
 real, parameter    :: pi=3.1415926535897931
 real, parameter    :: rad2deg = 180./3.14159265358979
 real, parameter    :: deg2rad = 3.14159265358979/180.

 public :: find_nearest_pole_points
 public :: find_poles
 public :: get_index
 public :: get_lat_angle
 public :: get_lon_angle
 public :: get_xnsum
 public :: get_xnsum2
 public :: get_xnsum3
 public :: inside_a_polygon
 public :: latlon2xyz
 public :: minmax
 public :: remove_isolated_pts
 public :: timef
 public :: transpose_orog
 public :: transpose_mask

 contains

 subroutine minmax(im,jm,a,title,imax,jmax)

 implicit none

 character(len=8), intent(in)   :: title
 
 integer, intent(in)            :: im, jm
 integer, intent(out), optional :: imax, jmax

 real,    intent(in)            :: a(im,jm)

 integer                        :: i, j

 real                           :: rmin,rmax

 rmin=huge(a)
 rmax=-rmin

 if (present(imax) .and. present(jmax)) then
   imax=0
   jmax=0
 endif

 do j=1,jm
   do i=1,im
     if(a(i,j) >= rmax) then
       rmax=a(i,j)
       if (present(imax) .and. present(jmax)) then
         imax = i
         jmax = j
       endif
     endif
     if(a(i,j) <= rmin)rmin=a(i,j)
   enddo
 enddo

 write(6,150) title,rmin,rmax
150 format(' - ',a8,' MIN=',e13.4,2x,'MAX=',e13.4)

 end subroutine minmax

 subroutine latlon2xyz(siz,lon, lat, x, y, z)

 implicit none

 integer, intent(in) :: siz
 real, intent(in)    :: lon(siz), lat(siz)
 real, intent(out)   :: x(siz), y(siz), z(siz)

 integer             :: n

 do n = 1, siz
   x(n) = cos(lat(n))*cos(lon(n))
   y(n) = cos(lat(n))*sin(lon(n))
   z(n) = sin(lat(n))
 enddo

 end subroutine latlon2xyz


 function get_lat_angle(dy)

 implicit none

 real, intent(in)   :: dy

 real               :: get_lat_angle

 get_lat_angle = dy/earth_radius*rad2deg

 end function get_lat_angle


 function get_lon_angle(dx,lat_in)

 implicit none

 real, intent(in)     :: dx, lat_in

 real                 :: get_lon_angle, lat

 lat = lat_in
 if (lat > 89.5) lat = 89.5
 if (lat < -89.5) lat = -89.5

 get_lon_angle = 2*asin( sin(dx/earth_radius*0.5)/cos(lat*deg2rad) )*rad2deg

 end function get_lon_angle


 subroutine transpose_mask(imn, jmn, mask)

 implicit none

 integer, intent(in)       :: imn, jmn
 integer(1), intent(inout) :: mask(imn,jmn)

 integer    :: i, j, it, jt
 integer(1) :: isave

! Transpose from S to N to the NCEP standard N to S.

 do j=1,jmn/2
 do i=1,imn
   jt=jmn - j + 1
   isave = mask(i,j)
   mask(i,j)=mask(i,jt)
   mask(i,jt) = isave
 enddo
 enddo

! Data begins at dateline. NCEP standard is Greenwich.

 do j=1,jmn
 do i=1,imn/2
   it=imn/2 + i
   isave = mask(i,j)
   mask(i,j)=mask(it,j)
   mask(it,j) = isave
 enddo
 enddo

 end subroutine transpose_mask


 subroutine transpose_orog(imn, jmn, glob)

 implicit none

 integer, intent(in)       :: imn, jmn
 integer(2), intent(inout) :: glob(imn,jmn)

 integer    :: i, j, it, jt
 integer(2) :: i2save

! Transpose from S to N to the NCEP standard N to S.

 do j=1,jmn/2
 do i=1,imn
   jt=jmn - j + 1
   i2save = glob(i,j)
   glob(i,j)=glob(i,jt)
   glob(i,jt) = i2save
 enddo
 enddo

! Data begins at dateline. NCEP standard is Greenwich.

 do j=1,jmn
 do i=1,imn/2
   it=imn/2 + i
   i2save = glob(i,j)
   glob(i,j)=glob(it,j)
   glob(it,j) = i2save
 enddo
 enddo

 end subroutine transpose_orog


 function spherical_angle(v1, v2, v3)

 implicit none

 real              :: spherical_angle

 real, parameter   :: epsln30 = 1.e-30

 real, intent(in)  :: v1(3), v2(3), v3(3)
 
 real              :: px, py, pz, qx, qy, qz, ddd

! vector product between v1 and v2

 px = v1(2)*v2(3) - v1(3)*v2(2)
 py = v1(3)*v2(1) - v1(1)*v2(3)
 pz = v1(1)*v2(2) - v1(2)*v2(1)

! vector product between v1 and v3

 qx = v1(2)*v3(3) - v1(3)*v3(2);
 qy = v1(3)*v3(1) - v1(1)*v3(3);
 qz = v1(1)*v3(2) - v1(2)*v3(1);

 ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz);
 if ( ddd <= 0.0 ) then
   spherical_angle = 0.
 else
   ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd);
   if( abs(ddd-1) < epsln30 ) ddd = 1;
   if( abs(ddd+1) < epsln30 ) ddd = -1;
   if ( ddd>1. .or. ddd<-1. ) then
    !FIX to correctly handle co-linear points (angle near pi or 0) */
     if (ddd < 0.) then
       spherical_angle = pi
     else
       spherical_angle = 0.
     endif
   else
     spherical_angle = acos( ddd )
   endif
 endif

 end function spherical_angle


 function inside_a_polygon(lon1, lat1, npts, lon2, lat2)

 implicit none

 logical inside_a_polygon

 real, parameter          :: epsln10 = 1.e-10
 real, parameter          :: epsln8 = 1.e-8
 real, parameter          :: range_check_criteria=0.05

 integer, intent(in)      :: npts

 real, intent(in)         :: lon1, lat1
 real, intent(in)         :: lon2(npts), lat2(npts)

 integer                  :: i, ip1

 real                     :: anglesum, angle
 real                     :: x2(npts), y2(npts), z2(npts)
 real                     :: lon1_1d(1), lat1_1d(1)
 real                     :: x1(1), y1(1), z1(1)
 real                     :: pnt0(3),pnt1(3),pnt2(3)
 real                     :: max_x2,min_x2,max_y2,min_y2,max_z2,min_z2

! first convert to cartesian grid.

 call latlon2xyz(npts,lon2, lat2, x2, y2, z2);
 lon1_1d(1) = lon1
 lat1_1d(1) = lat1
 call latlon2xyz(1,lon1_1d, lat1_1d, x1, y1, z1);
 inside_a_polygon = .false.
 max_x2 = maxval(x2)
 if( x1(1) > max_x2+range_check_criteria ) return
 min_x2 = minval(x2)
 if( x1(1)+range_check_criteria < min_x2 ) return
 max_y2 = maxval(y2)
 if( y1(1) > max_y2+range_check_criteria ) return
 min_y2 = minval(y2)
 if( y1(1)+range_check_criteria < min_y2 ) return
 max_z2 = maxval(z2)
 if( z1(1) > max_z2+range_check_criteria ) return
 min_z2 = minval(z2)
 if( z1(1)+range_check_criteria < min_z2 ) return

 pnt0(1) = x1(1)
 pnt0(2) = y1(1)
 pnt0(3) = z1(1)

 anglesum = 0

 do i = 1, npts
   if(abs(x1(1)-x2(i)) < epsln10 .and.  &
      abs(y1(1)-y2(i)) < epsln10 .and.  &
      abs(z1(1)-z2(i)) < epsln10 ) then ! same as the corner point
     inside_a_polygon = .true.
     return
   endif
   ip1 = i+1
   if(ip1>npts) ip1 = 1
   pnt1(1) = x2(i)
   pnt1(2) = y2(i)
   pnt1(3) = z2(i)
   pnt2(1) = x2(ip1)
   pnt2(2) = y2(ip1)
   pnt2(3) = z2(ip1)
   angle = spherical_angle(pnt0, pnt2, pnt1);
   anglesum = anglesum + angle
 enddo

 if(abs(anglesum-2*pi) < epsln8) then
   inside_a_polygon = .true.
 else
   inside_a_polygon = .false.
 endif

 end function inside_a_polygon

 subroutine find_poles(geolat, nx, ny, i_north_pole, j_north_pole, &
                       i_south_pole, j_south_pole)

 implicit none

 integer, intent(in)  :: nx, ny

 real, intent(in)     :: geolat(nx+1,ny+1)

 integer, intent(out) :: i_north_pole, j_north_pole
 integer, intent(out) :: i_south_pole, j_south_pole

 integer              :: i, j

 real                 :: maxlat, minlat

 print*,'- CHECK IF THE TILE CONTAINS A POLE.'

!--- figure out pole location.

 maxlat = -90
 minlat = 90
 i_north_pole = 0
 j_north_pole = 0
 i_south_pole = 0
 j_south_pole = 0
 do j = 1, ny+1; do i = 1, nx+1
   if( geolat(i,j) > maxlat ) then
      i_north_pole=i
      j_north_pole=j
      maxlat = geolat(i,j)
   endif
   if( geolat(i,j) < minlat ) then
      i_south_pole=i
      j_south_pole=j
      minlat = geolat(i,j)
   endif
 enddo ; enddo

!--- only when maxlat is close to 90. the point is north pole

 if(maxlat < 89.9 ) then
   i_north_pole = 0
   j_north_pole = 0
 endif
 if(minlat > -89.9 ) then
   i_south_pole = 0
   j_south_pole = 0
 endif

 print*, "- MINLAT=", minlat, "MAXLAT=", maxlat
 print*, "- NORTH POLE SUPERGRID INDEX IS ",  &
                i_north_pole, j_north_pole
 print*, "- SOUTH POLE SUPERGRID INDEX IS ",  &
                i_south_pole, j_south_pole

 end subroutine find_poles


 subroutine find_nearest_pole_points(i_north_pole, j_north_pole, &
      i_south_pole, j_south_pole, im, jm, is_north_pole, &
      is_south_pole)

 implicit none

 integer, intent(in)     :: im, jm
 integer, intent(in)     :: i_north_pole, j_north_pole
 integer, intent(in)     :: i_south_pole, j_south_pole

 logical, intent(out)    :: is_north_pole(im,jm)
 logical, intent(out)    :: is_south_pole(im,jm)

 integer                 :: i, j

 print*,'- FIND NEAREST POLE POINTS.'

 is_north_pole=.false.
 is_south_pole=.false.

 if(i_south_pole >0 .and. j_south_pole > 0) then
   if(mod(i_south_pole,2)==0) then ! stretched grid
      do j = 1, jm; do i = 1, im
        if(i==i_south_pole/2 .and. (j==j_south_pole/2  &
                    .or. j==j_south_pole/2+1) ) then
          is_south_pole(i,j) = .true.
          print*, "- SOUTH POLE AT I,J= ", i, j
        endif
      enddo; enddo
   else
      do j = 1, jm; do i = 1, im
        if((i==i_south_pole/2 .or. i==i_south_pole/2+1) &
             .and. (j==j_south_pole/2 .or. &
             j==j_south_pole/2+1) ) then
          is_south_pole(i,j) = .true.
          print*, "- SOUTH POLE AT I,J= ", i, j
        endif
      enddo; enddo
    endif
 endif

 if(i_north_pole >0 .and. j_north_pole > 0) then
   if(mod(i_north_pole,2)==0) then ! stretched grid
      do j = 1, jm; do i = 1, im
        if(i==i_north_pole/2 .and. (j==j_north_pole/2 .or.  &
           j==j_north_pole/2+1) ) then
          is_north_pole(i,j) = .true.
          print*, "- NORTH POLE AT I,J= ", i, j
        endif
      enddo; enddo
    else
      do j = 1, jm; do i = 1, im
        if((i==i_north_pole/2 .or. i==i_north_pole/2+1) &
            .and. (j==j_north_pole/2 .or.  &
                   j==j_north_pole/2+1) ) then
          is_north_pole(i,j) = .true.
          print*, "- NORTH POLE AT I,J= ", i, j
        endif
      enddo; enddo
   endif
 endif

 end subroutine find_nearest_pole_points


 subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol)

 implicit none

 integer, intent(in)       :: im, jm

 real, intent(inout)       :: slm(im,jm)
 real, intent(inout)       :: oro(im,jm)
 real, intent(inout)       :: var(im,jm)
 real, intent(inout)       :: var4(im,jm)
 real, intent(inout)       :: oa(im,jm,4)
 real, intent(inout)       :: ol(im,jm,4)

 integer                   :: i, j, jn, js, k
 integer                   :: iw, ie, wgta, is, ise
 integer                   :: in, ine, inw, isw

 real                      :: slma, oroa, vara, var4a
 real, allocatable         :: oaa(:), ola(:)

!  REMOVE ISOLATED POINTS

 print*,"- REMOVE ISOLATED POINTS."

 allocate (oaa(4),ola(4))

 iso_loop : DO j=2,jm-1
   jn=j-1
   js=j+1
   i_loop : DO i=2,im-1
! Check the points to the 'west' and 'east'.
     iw=i-1
     ie=i+1
     slma=slm(iw,j)+slm(ie,j)
     oroa=oro(iw,j)+oro(ie,j)
     vara=var(iw,j)+var(ie,j)
     var4a=var4(iw,j)+var4(ie,j)
     DO k=1,4
       oaa(k)=oa(iw,j,k)+oa(ie,j,k)
! --- (*j*) fix typo:
       ola(k)=ol(iw,j,k)+ol(ie,j,k)
     ENDDO
     wgta=2
! Check the points to the 'northwest', 'north' and 'northeast'
     in=i
     inw=i-1
     ine=i+1
     slma=slma+slm(inw,jn)+slm(in,jn)+slm(ine,jn)
     oroa=oroa+oro(inw,jn)+oro(in,jn)+oro(ine,jn)
     vara=vara+var(inw,jn)+var(in,jn)+var(ine,jn)
     var4a=var4a+var4(inw,jn)+var4(in,jn)+var4(ine,jn)
     DO k=1,4
       oaa(k)=oaa(k)+oa(inw,jn,k)+oa(in,jn,k)+oa(ine,jn,k)
       ola(k)=ola(k)+ol(inw,jn,k)+ol(in,jn,k)+ol(ine,jn,k)
     ENDDO
     wgta=wgta+3
! Check the points to the 'southwest', 'south' and 'southeast'  
     is=i
     isw=i-1
     ise=i+1
     slma=slma+slm(isw,js)+slm(is,js)+slm(ise,js)
     oroa=oroa+oro(isw,js)+oro(is,js)+oro(ise,js)
     vara=vara+var(isw,js)+var(is,js)+var(ise,js)
     var4a=var4a+var4(isw,js)+var4(is,js)+var4(ise,js)
     DO k=1,4
       oaa(k)=oaa(k)+oa(isw,js,k)+oa(is,js,k)+oa(ise,js,k)
       ola(k)=ola(k)+ol(isw,js,k)+ol(is,js,k)+ol(ise,js,k)
     ENDDO
     wgta=wgta+3
! Take the average of the surrounding the points.
     oroa=oroa/wgta
     vara=vara/wgta
     var4a=var4a/wgta
     DO k=1,4
       oaa(k)=oaa(k)/wgta
       ola(k)=ola(k)/wgta
     ENDDO
! Isolated water point.
     IF(slm(i,j).EQ.0..AND.slma.EQ.wgta) THEN
       print '(" - SEA ",2F8.0," MODIFIED TO LAND",2F8.0,  &
               " AT ",2I8)',oro(i,j),var(i,j),oroa,vara,i,j
       slm(i,j)=1.
       oro(i,j)=oroa
       var(i,j)=vara
       var4(i,j)=var4a
       DO k=1,4
         oa(i,j,k)=oaa(k)
         ol(i,j,k)=ola(k)
       ENDDO
! Isolated land point.
     ELSEIF(slm(i,j).EQ.1..AND.slma.EQ.0.) THEN
       print '(" - LAND",2F8.0," MODIFIED TO SEA ",2F8.0,  &
               " AT ",2I8)',oro(i,j),var(i,j),oroa,vara,i,j
       slm(i,j)=0.
       oro(i,j)=oroa
       var(i,j)=vara
       var4(i,j)=var4a
       DO k=1,4
         oa(i,j,k)=oaa(k)
         ol(i,j,k)=ola(k)
       ENDDO
     ENDIF
   ENDDO i_loop
 ENDDO iso_loop

 deallocate (oaa,ola)

 end subroutine remove_isolated_pts

 subroutine get_index(imn,jmn,npts,lonO,latO,delxn,  &
                jst,jen,ilist,numx)

 implicit none
 integer, intent(in)  :: imn,jmn
 integer, intent(in)  :: npts
 real,    intent(in)  :: lono(npts), lato(npts)
 real,    intent(in)  :: delxn
 integer, intent(out) :: jst,jen
 integer, intent(out) :: ilist(imn)
 integer, intent(out) :: numx

 integer              :: i2, ii, ist, ien
 real                 :: minlat,maxlat,minlon,maxlon
        
 minlat = minval(lato)
 maxlat = maxval(lato)
 minlon = minval(lono)
 maxlon = maxval(lono)
 ist = minlon/delxn+1
 ien = maxlon/delxn+1
 jst = (minlat+90)/delxn+1 
 jen = (maxlat+90)/delxn 
!--- add a few points to both ends of j-direction
 jst = jst - 5
 if(jst<1) jst = 1
 jen = jen + 5
 if(jen>jmn) jen = jmn

!--- when around the pole, just search through all the points.
 if((jst == 1 .OR. jen == jmn) .and. &
    (ien-ist+1 > imn/2) )then
   numx = imn
   do i2 = 1, imn
     ilist(i2) = i2
   enddo
 else if( ien-ist+1 > imn/2 ) then  ! cross longitude = 0
!--- find the minimum that greater than IMN/2 
!--- and maximum that less than IMN/2
   ist = 0
   ien = imn+1
   do i2 = 1, npts
     ii = lono(i2)/delxn+1
     if(ii <0 .or. ii>imn) then
        print*,"- II=",ii,imn,lono(i2),delxn
     endif
     if( ii < imn/2 ) then
       ist = max(ist,ii)
     else if( ii > imn/2 ) then
       ien = min(ien,ii)
     endif
   enddo 
   if(ist<1 .OR. ist>imn) then
     print*, .or."FATAL ERROR: ist<1  ist>IMN"
     call abort()
   endif
   if(ien<1 .OR. ien>imn) then
     print*, .or."FATAL ERROR: iend<1  iend>IMN"
     call abort()
   endif

   numx = imn - ien + 1
   do i2 = 1, numx
     ilist(i2) = ien + (i2-1)           
   enddo
   do i2 = 1, ist
     ilist(numx+i2) = i2
   enddo
   numx = numx+ist
 else
   numx = ien-ist+1
   do i2 = 1, numx
     ilist(i2) = ist + (i2-1)
   enddo
 endif

 end subroutine get_index


 function get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, &
                    glat,zavg,zslm,delxn)

 implicit none

 real                    :: get_xnsum

 integer, intent(in)     :: imn,jmn
 integer, intent(in)     :: zavg(imn,jmn),zslm(imn,jmn)
 real, intent(in)        :: lon1,lat1,lon2,lat2,delxn
 real, intent(in)        :: glat(jmn)

 integer                 :: i, j, ist, ien, jst, jen, i1
 real                    :: oro, height
 real                    :: xland,xwatr,xl1,xs1,slm,xnsum

!-- Figure out ist,ien,jst,jen

 do j = 1, jmn
   if( glat(j) .gt. lat1 ) then
     jst = j
     exit
   endif
 enddo

 do j = 1, jmn
   if( glat(j) .gt. lat2 ) then
     jen = j
     exit
   endif
 enddo

 ist = lon1/delxn + 1
 ien = lon2/delxn
 if(ist .le.0) ist = ist + imn
 if(ien < ist) ien = ien + imn

!--- Compute average oro

 oro = 0.0
 xnsum = 0
 xland = 0
 xwatr = 0
 xl1 = 0
 xs1 = 0
 do j = jst,jen
   do i1 = 1, ien - ist + 1
     i = ist + i1 -1
     if( i .le. 0) i = i + imn
     if( i .gt. imn) i = i - imn
     xland = xland + float(zslm(i,j))
     xwatr = xwatr + float(1-zslm(i,j))
     xnsum = xnsum + 1.
     height = float(zavg(i,j))
     if(height.lt.-990.) height = 0.0
     xl1 = xl1 + height * float(zslm(i,j))
     xs1 = xs1 + height * float(1-zslm(i,j))
   enddo
 enddo

 if( xnsum > 1.) then
   slm = float(nint(xland/xnsum))
   if(slm.ne.0.) then
     oro= xl1 / xland
   else
     oro = xs1 / xwatr
   endif
 endif

 get_xnsum = 0
 do j = jst, jen
 do i1= 1, ien-ist+1
   i = ist + i1 -1
   if( i .le. 0) i = i + imn
   if( i .gt. imn) i = i - imn
   height = float(zavg(i,j))
   if(height.lt.-990.) height = 0.0
   if ( height .gt. oro ) get_xnsum = get_xnsum + 1
 enddo
 enddo

 end function get_xnsum


 subroutine get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn, &
                   glat,zavg,delxn,xnsum1,xnsum2,hc)

 implicit none

 integer, intent(in) :: imn,jmn
 integer, intent(in) :: zavg(imn,jmn)

 real, intent(in)    :: lon1,lat1,lon2,lat2,delxn
 real, intent(in)    :: glat(jmn)
 real, intent(out)   :: xnsum1,xnsum2,hc

 integer             :: i, j, ist, ien, jst, jen, i1

 real                :: height, var
 real                :: xw1,xw2,xnsum

!-- Figure out ist,ien,jst,jen

 do j = 1, jmn
   if( glat(j) .gt. lat1 ) then
     jst = j
     exit
   endif
 enddo

 do j = 1, jmn
   if( glat(j) .gt. lat2 ) then
     jen = j
     exit
   endif
 enddo

 ist = lon1/delxn + 1
 ien = lon2/delxn
 if(ist .le.0) ist = ist + imn
 if(ien < ist) ien = ien + imn

!--- Compute average oro

 xnsum = 0
 xw1 = 0
 xw2 = 0
 do j = jst,jen
   do i1 = 1, ien - ist + 1
     i = ist + i1 -1
     if( i .le. 0) i = i + imn
     if( i .gt. imn) i = i - imn
     xnsum = xnsum + 1.
     height = float(zavg(i,j))
     if(height.lt.-990.) height = 0.0
     xw1 = xw1 + height
     xw2 = xw2 + height ** 2
   enddo
 enddo

 var = sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
 hc = 1116.2 - 0.878 * var
 xnsum1 = 0
 xnsum2 = 0
 do j = jst, jen
   do i1= 1, ien-ist+1
     i = ist + i1 -1
     if( i .le. 0) i = i + imn
     if( i .gt. imn) i = i - imn
     height = float(zavg(i,j))
     if ( height .gt. hc ) xnsum1 = xnsum1 + 1
     xnsum2 = xnsum2 + 1
   enddo
 enddo

 end subroutine get_xnsum2


 subroutine get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn, &
                      glat,zavg,delxn,xnsum1,xnsum2,HC)
 implicit none

 integer, intent(in) :: imn,jmn
 integer, intent(in) :: zavg(imn,jmn)

 real, intent(in)    :: hc, glat(jmn)
 real, intent(in)    :: lon1,lat1,lon2,lat2,delxn
 real, intent(out)   :: xnsum1,xnsum2

 integer             :: i, j, ist, ien, jst, jen, i1

 real                :: height

!-- Figure out ist,ien,jst,jen

! if lat1 or lat 2 is 90 degree. set jst = JMN

 jst = jmn
 jen = jmn
 do j = 1, jmn
   if( glat(j) .gt. lat1 ) then
     jst = j
     exit
   endif
 enddo

 do j = 1, jmn
   if( glat(j) .gt. lat2 ) then
     jen = j
     exit
   endif
 enddo

 ist = lon1/delxn + 1
 ien = lon2/delxn
 if(ist .le.0) ist = ist + imn
 if(ien < ist) ien = ien + imn

 xnsum1 = 0
 xnsum2 = 0
 do j = jst, jen
   do i1= 1, ien-ist+1
     i = ist + i1 -1
     if( i .le. 0) i = i + imn
     if( i .gt. imn) i = i - imn
     height = float(zavg(i,j))
     if ( height .gt. hc ) xnsum1 = xnsum1 + 1
     xnsum2 = xnsum2 + 1
   enddo
 enddo

 end subroutine get_xnsum3

 real function timef()

 implicit none

 character(8)          :: date
 character(10)         :: time
 character(5)          :: zone
 integer,dimension(8)  :: values
 integer               :: total
 real                  :: elapsed

 call date_and_time(date,time,zone,values)
 total=(3600*values(5)) + (60*values(6))+values(7)
 elapsed=float(total) + (1.0e-3*float(values(8)))
 timef=elapsed

 end function timef

 end module orog_utils