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