| 1 |
c*** ensemble square root filter that uses an optional cut-off |
| 2 |
c*** radius and boost factor. |
| 3 |
|
| 4 |
subroutine gregfilt_loc(xens,yo,iobsloc,ngp,mobs,Rs,nens,nx,ny) |
| 5 |
|
| 6 |
implicit none |
| 7 |
|
| 8 |
! Arguments |
| 9 |
integer, intent(in) :: nens, mobs, ngp, nx, ny |
| 10 |
real*8, intent(inout) :: xens(ngp,nens) |
| 11 |
real*8, intent(in) :: yo(mobs), Rs(mobs) |
| 12 |
real*8, intent(in) :: iobsloc(mobs) |
| 13 |
|
| 14 |
! Local Variables |
| 15 |
integer :: xob(mobs), yob(mobs) |
| 16 |
integer :: ind, k, j, i, r2, kk, jj, ko, jo, kj, g |
| 17 |
real*8 :: PHT(ngp), Ks(ngp), Khat(ngp) |
| 18 |
real*8 :: xp(ngp), xa(ngp), zp(ngp,nens), R(mobs) |
| 19 |
real*8 :: HPHT, alpha, boost |
| 20 |
|
| 21 |
! Filter Stuff |
| 22 |
integer :: d, rad, rad2 |
| 23 |
real*8 :: dr, rrad, filt |
| 24 |
real*8, external :: cov_factor |
| 25 |
|
| 26 |
c*** observational standard deviation |
| 27 |
R = sqrt(Rs) |
| 28 |
|
| 29 |
c*** set cutoff radius |
| 30 |
rad = 100 |
| 31 |
rad2 = 2*rad |
| 32 |
r2 = rad*rad |
| 33 |
rrad = float(rad) |
| 34 |
|
| 35 |
c*** set inflation factor |
| 36 |
boost = 1.00 |
| 37 |
! boost = 1.05 |
| 38 |
! boost = 1.10 |
| 39 |
|
| 40 |
c*** rename ensemble matrix |
| 41 |
zp = xens |
| 42 |
|
| 43 |
!*** Find the initial ensemble mean |
| 44 |
do j = 1, ngp |
| 45 |
xp(j) = sum(zp(j,:))/float(nens) |
| 46 |
end do |
| 47 |
|
| 48 |
!*** Apply inflation factor to initial ensemble |
| 49 |
do k = 1, nens |
| 50 |
zp(:,k) = boost*(zp(:,k) - xp) + xp |
| 51 |
end do |
| 52 |
|
| 53 |
!*** Find the xob and yob arrays from H |
| 54 |
c do j = 1, ngp |
| 55 |
do k = 1, mobs |
| 56 |
c if ( H(k,j) == 1. ) then |
| 57 |
c xob(k) = mod(j-1,nx) + 1 |
| 58 |
c yob(k) = (j-1)/nx + 1 |
| 59 |
xob(k) = mod(iobsloc(k)-1,nx) + 1 |
| 60 |
yob(k) = (iobsloc(k)-1)/nx + 1 |
| 61 |
c end if |
| 62 |
end do |
| 63 |
c end do |
| 64 |
|
| 65 |
!*** Now process each observation sequentially abiding by cut-off radius |
| 66 |
do j = 1, mobs |
| 67 |
|
| 68 |
ind = nx*( yob(j) - 1 ) + xob(j) |
| 69 |
|
| 70 |
!*** Find PH' first |
| 71 |
PHT = 0. |
| 72 |
do jj = yob(j)-rad2, yob(j)+rad2 |
| 73 |
do kk = xob(j)-rad2, xob(j)+rad2 |
| 74 |
jo = jj |
| 75 |
ko = kk |
| 76 |
|
| 77 |
!*** Point is within block of radius, but it may not be within the |
| 78 |
! basin boundaries |
| 79 |
if ( ko>0 .and. ko<=nx .and. jo>0 .and. jo<=ny ) then |
| 80 |
|
| 81 |
!*** Since we've sequestered a square of side 2*rad and the |
| 82 |
! cut-off radius assumes a circle, we need to check to |
| 83 |
! make sure the point we're considering is actually |
| 84 |
! within the circle. |
| 85 |
|
| 86 |
d = (ko - xob(j))**2 + (jo - yob(j))**2 |
| 87 |
|
| 88 |
dr = sqrt( float( d ) ) |
| 89 |
|
| 90 |
!*** The element of interest in the 1D vector according to addresses |
| 91 |
! kk and jj is: |
| 92 |
kj = nx*(jo-1) + ko |
| 93 |
|
| 94 |
!*** Evaluate the filter coefficient based on distance from center d |
| 95 |
! filt = cov_factor(dr,rrad) |
| 96 |
filt = 1. |
| 97 |
|
| 98 |
!*** Now contribute to PHT sum |
| 99 |
do g = 1, nens |
| 100 |
PHT(kj) = PHT(kj) + filt*(zp(kj,g) - xp(kj))* |
| 101 |
& (zp(ind,g) - xp(ind)) |
| 102 |
end do |
| 103 |
|
| 104 |
end if |
| 105 |
end do |
| 106 |
end do |
| 107 |
PHT = PHT/float(nens - 1) |
| 108 |
|
| 109 |
!*** Now find HPH' from PH'. Because of cut-off radius, there is a |
| 110 |
! (good) chance that HPH' will be zero. |
| 111 |
HPHT = PHT(ind) |
| 112 |
|
| 113 |
!*** Evaluate Ks |
| 114 |
Ks = PHT/( HPHT + Rs(j) ) |
| 115 |
|
| 116 |
!*** Update all effected elements in the mean |
| 117 |
xa = xp + Ks*( yo(j) - xp(ind) ) |
| 118 |
|
| 119 |
!*** Now update all ensemble members as perturbations about mean |
| 120 |
alpha = 1./( 1. + sqrt( Rs(j)/( HPHT + Rs(j) ) ) ) |
| 121 |
Khat = alpha*Ks |
| 122 |
|
| 123 |
do g = 1, nens |
| 124 |
zp(:,g) = ((zp(:,g) - xp) - Khat*( zp(ind,g) - |
| 125 |
& xp(ind) )) + xa |
| 126 |
end do |
| 127 |
|
| 128 |
!*** Use analysis ensemble as the background for the next observation |
| 129 |
xp = xa |
| 130 |
|
| 131 |
end do |
| 132 |
|
| 133 |
c print*, 'EnSRF:: Done Mobs Loop' |
| 134 |
|
| 135 |
xens = zp |
| 136 |
|
| 137 |
return |
| 138 |
|
| 139 |
end subroutine gregfilt_loc |
| 140 |
|
| 141 |
!---------------------------------------------------------------------------- |
| 142 |
c*** distance filter |
| 143 |
function cov_factor(z_in, c) |
| 144 |
|
| 145 |
implicit none |
| 146 |
|
| 147 |
double precision :: cov_factor |
| 148 |
double precision, intent(in) :: z_in, c |
| 149 |
double precision :: z, r |
| 150 |
|
| 151 |
z = abs(z_in) |
| 152 |
r = z / c |
| 153 |
|
| 154 |
if ( z >= 2*c ) then |
| 155 |
cov_factor = 0. |
| 156 |
else if ( z >= c .and. z < 2*c ) then |
| 157 |
cov_factor = r**5/12. - r**4/2. + r**3 * 5./8. + |
| 158 |
& r**2 * 5./3. - 5*r + 4. - (2 * c) / (3 * z) |
| 159 |
else |
| 160 |
cov_factor = r**5 * (-1./4.) + r**4/2. + r**3 * |
| 161 |
& 5./8. - r**2 * 5./3. + 1. |
| 162 |
end if |
| 163 |
|
| 164 |
end function cov_factor |
| 165 |
|
| 166 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph