!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine calc_shvort_llgrid(uwnd, vwnd, lat, lon, nlon, nlat, shvort)

implicit none

real, parameter     :: pid    = atan(1.0) / 45.0
real, parameter     :: Rearth = 6378388.

integer, intent(in) :: nlon, nlat
real, intent(in)    :: uwnd(nlon,nlat), vwnd(nlon,nlat), lat(nlat), lon(nlon)
real, intent(out)   :: shvort(nlon,nlat)

integer  :: i, j, i1, i2, j1, j2

real     :: vlen, ivec, jvec, v1, v2, u1, u2, dlon

do i = 1, nlon

  i1 = i-1
  if ( i1 < 1    )  i1 = i1 + nlon
  i2 = i+1
  if ( i2 > nlon )  i2 = i2 - nlon

  dlon = lon(i2)-lon(i1)
  if ( dlon >  180. )  dlon = dlon - 360.0
  if ( dlon < -180. )  dlon = dlon + 360.0

  do j = 1, nlat

    j1 = max(j-1,   1)
    j2 = min(j+1,nlat)

    vlen = sqrt(uwnd(i,j)*uwnd(i,j) + vwnd(i,j)*vwnd(i,j))
    ivec = uwnd(i,j) / vlen  ;  jvec = vwnd(i,j) / vlen    

    v1 = (ivec * uwnd(i1,j) + jvec * vwnd(i1,j)) * jvec
    v2 = (ivec * uwnd(i2,j) + jvec * vwnd(i2,j)) * jvec

    u1 = (ivec * uwnd(i,j1) + jvec * vwnd(i,j1)) * ivec
    u2 = (ivec * uwnd(i,j2) + jvec * vwnd(i,j2)) * ivec

    shvort(i,j) = (v2 - v1) / max(dlon*pid*Rearth*cos(pid*lat(j)), 1.0) - &
                  (u2 - u1) / ((lat(j2)-lat(j1))*pid*Rearth)

  end do

end do 

return
end subroutine calc_shvort_llgrid

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine calc_vort_llgrid(uwnd, vwnd, lat, lon, nlon, nlat, vort)

implicit none

real, parameter     :: pid    = atan(1.0) / 45.0
real, parameter     :: Rearth = 6378388.

integer, intent(in) :: nlon, nlat
real, intent(in)    :: uwnd(nlon,nlat), vwnd(nlon,nlat), lat(nlat), lon(nlon)
real, intent(out)   :: vort(nlon,nlat)

integer :: i, j, i1, i2, j1, j2

real    :: dlon

do i = 1, nlon
  i1 = i-1
  if ( i1 < 1 )  i1 = i1 + nlon
  i2 = i+1
  if ( i2 > nlon )  i2 = i2 - nlon

  dlon = lon(i2)-lon(i1)
  if ( dlon >  180. )  dlon = dlon - 360.0
  if ( dlon < -180. )  dlon = dlon + 360.0

  do j = 1, nlat
    vort(i,j) = (vwnd(i2,j)-vwnd(i1,j)) / max(dlon*pid*Rearth*cos(pid*lat(j)),1.0)
  end do
end do

do j = 1, nlat
  j1 = max(j-1, 1)
  j2 = min(j+1,nlat)
  do i = 1, nlon
    vort(i,j) = vort(i,j) - (uwnd(i,j2)-uwnd(i,j1)) / &
                             ((lat(j2)-lat(j1))*Rearth*pid)
  end do
end do

return
end subroutine calc_vort_llgrid


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine circulation_llgrid(vort, lat, lon, circrad, nlon, nlat, nens, circ)

implicit none

real, parameter     :: pid    = atan(1.0) / 45.0
real, parameter     :: Rearth = 6378388.

integer, intent(in) :: nlon, nlat, nens
real, intent(in)    :: vort(nens,nlon,nlat), lon(nlon), lat(nlat), circrad
real, intent(out)   :: circ(nens,nlon,nlat)

integer :: deli, delj, i, j, n, ii, jj, la1, la2, nla, ii1s, ii1e, ii2s, ii2e

real    :: dX, garea, gc_dist, dlat, dlon, asum(nlon,nlat)


dlat = abs(lat(2)-lat(1))
dlon = abs(lon(2)-lon(1))

delj = ceiling(circrad * 1000.0 / (Rearth*pid*dlat))

circ(:,:,:) = 0.0  ;  asum(:,:) = 0.0

do j = 1, nlat

  la1 = max(j-delj,    1)
  la2 = min(j+delj, nlat)
  nla = la2 - la1 + 1

  dX = Rearth
  do jj = la1, la2
    dX = min(abs(dlon*Rearth*pid*cos(pid*lat(jj))),dX)
  end do
  deli = min(ceiling(circrad / dX * 1000.0), nlon/2)

  do i = 1, nlon

    if ( deli == nlon/2 ) then
      ii1s = 1
      ii1e = nlon
      ii2s = 0
      ii2e = -1
    else
      ii1s = i - deli
      ii1e = i + deli
      ii2s = 0
      ii2e = -1
      if ( ii1s < 1 ) then
        ii2s = ii1s + nlon
        ii2e = nlon
        ii1s = 1
      end if
      if ( ii1e > nlon ) then
        ii2s = 1
        ii2e = ii1e - nlon
        ii1e = nlon
      end if
    end if

    do jj = la1, la2

      do ii = ii1s, ii1e

        if ( gc_dist(lat(j),lon(i),lat(jj),lon(ii))*0.001 <= circrad ) then
          garea       = max(dlon*pid*Rearth*cos(pid*lat(jj))*dlat*pid*Rearth,0.001)
          circ(:,i,j) = circ(:,i,j) + vort(:,ii,jj) * garea
          asum(i,j)   = asum(i,j) + garea
        end if

      end do

      do ii = ii2s, ii2e

        if ( gc_dist(lat(j),lon(i),lat(jj),lon(ii))*0.001 <= circrad ) then
          garea       = max(dlon*pid*Rearth*cos(pid*lat(jj))*dlat*pid*Rearth,0.001)
          circ(:,i,j) = circ(:,i,j) + vort(:,ii,jj) * garea
          asum(i,j)   = asum(i,j) + garea
        end if

      end do

    end do

  end do

end do

do n = 1, nens
  do i = 1, nlon  ;  do j = 1, nlat
    circ(n,i,j) = circ(n,i,j) / asum(i,j)
  end do  ;  end do
end do

return
end subroutine circulation_llgrid


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine circulation_llgrid_5d(vort, lat, lon, circrad, nlon, nlat, nens, ntimes, circ)

implicit none

real, parameter     :: pid    = atan(1.0) / 45.0
real, parameter     :: Rearth = 6378388.

integer, intent(in) :: nlon, nlat, nens, ntimes
real, intent(in)    :: vort(nens,ntimes,nlon,nlat), lon(nlon), lat(nlat), circrad
real, intent(out)   :: circ(nens,ntimes,nlon,nlat)

integer :: deli, delj, i, j, n, t, ii, jj, la1, la2, nla, ii1s, ii1e, ii2s, ii2e

real    :: dX, garea, gc_dist, dlat, dlon, asum(nlon,nlat)


dlat = abs(lat(2)-lat(1))
dlon = abs(lon(2)-lon(1))

delj = ceiling(circrad * 1000.0 / (Rearth*pid*dlat))

circ(:,:,:,:) = 0.0  ;  asum(:,:) = 0.0

do j = 1, nlat

  la1 = max(j-delj,    1)
  la2 = min(j+delj, nlat)
  nla = la2 - la1 + 1

  dX = Rearth
  do jj = la1, la2
    dX = min(abs(dlon*Rearth*pid*cos(pid*lat(jj))),dX)
  end do
  deli = min(ceiling(circrad / dX * 1000.0), nlon/2)

  do i = 1, nlon

    if ( deli == nlon/2 ) then
      ii1s = 1
      ii1e = nlon
      ii2s = 0
      ii2e = -1
    else
      ii1s = i - deli
      ii1e = i + deli
      ii2s = 0
      ii2e = -1
      if ( ii1s < 1 ) then
        ii2s = ii1s + nlon
        ii2e = nlon
        ii1s = 1
      end if
      if ( ii1e > nlon ) then
        ii2s = 1
        ii2e = ii1e - nlon
        ii1e = nlon
      end if
    end if

    do jj = la1, la2

      do ii = ii1s, ii1e

        if ( gc_dist(lat(j),lon(i),lat(jj),lon(ii))*0.001 <= circrad ) then
          garea       = max(dlon*pid*Rearth*cos(pid*lat(jj))*dlat*pid*Rearth,0.001)
          circ(:,:,i,j) = circ(:,:,i,j) + vort(:,:,ii,jj) * garea
          asum(i,j)   = asum(i,j) + garea
        end if

      end do

      do ii = ii2s, ii2e

        if ( gc_dist(lat(j),lon(i),lat(jj),lon(ii))*0.001 <= circrad ) then
          garea       = max(dlon*pid*Rearth*cos(pid*lat(jj))*dlat*pid*Rearth,0.001)
          circ(:,:,i,j) = circ(:,:,i,j) + vort(:,:,ii,jj) * garea
          asum(i,j)   = asum(i,j) + garea
        end if

      end do

    end do

  end do

end do

do n = 1, nens
  do t = 1, ntimes
    do i = 1, nlon  ;  do j = 1, nlat
      circ(n,t,i,j) = circ(n,t,i,j) / asum(i,j)
    end do  ;  end do
  end do
end do

return
end subroutine circulation_llgrid_5d


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


function gc_dist(lat1,lon1,lat2,lon2)

implicit none

real, parameter     :: pid    = atan(1.0) / 45.0
real, parameter     :: Rearth = 6378388.

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

real :: gc_dist

gc_dist = Rearth * acos( sin(lat1*pid) * sin(lat2*pid) + &
                         cos(lat1*pid) * cos(lat2*pid) * &
                         cos((lon2-lon1)*pid) )

return
end function gc_dist