helmholtz.f90

集合卡尔曼滤波(EnKF) 数据同化方法可以避免了EKF 中协方差演变方程预报过程中出现的计算不准确和关于协方差矩阵的大量数据的存储问题,最主要的是可以有效的控制估计误差方差的增长,改善预报的效果。

!! Copyright (C) 2008 Pavel Sakov
!! 
!! This file is part of EnKF-Matlab. EnKF-Matlab is a free software. See 
!! LICENSE for details.

! File:           helmholtz.f90
!
! Created:        31/08/2007
!
! Last modified:  08/02/2008
!
! Author:         Pavel Sakov
!                 CSIRO Marine and Atmospheric Research
!                 NERSC
!
! Purpose:        Fortran code for QG model. Solver for Helmholtz equation.
!
! Description:    
!
! Revisions:
module helmholtz_mod

  use utils_mod

  implicit none

  save

  public helmholtz

  integer, private :: iq
  real(8), allocatable, dimension(:) :: q
  integer, parameter, private :: LEN20 = 20
  integer, dimension(LEN20), private :: nst, imx, jmx
  real(8), dimension(LEN20), private :: h

contains

  subroutine helmholtz(G, F, zk2, NX1, NY1, M, H1, NU1, NU2, NCYC, NX, MY, P)
    !
    ! Solves helmholtz equation:
    ! (\nabla^2 - zk2) P  = F, with Dirichlet B.C. The B.C. and the initial
    ! guess for the solution are given in the array G.
    !
    ! The dimension of the arrays should be:
    ! (NX, MY) = (p * 2^(M - 1) + 1, q * 2^(M - 1) + 1).
    ! p and q are small integers, equal to the number of grid boxes in the
    ! coarsest grid used in the multi-grid solution.
    !
    ! Note that number of boxes on a side of the grid is equal to the number
    ! of grid points on this side - 1.
    !
    ! In the code below p = NX1; q = NX2.                                       **
    ! M is the finest level of the multigrid. M=1 means use (p,q)     **
    ! boxes; M=2 means use twice as many boxes on each lide of grid.  **
    !
    ! H1 is the delx for the coarsest grid. (NU1,NU2) should be (1,2). **
    ! NCYCL should be 3-5 for a bad initial guess, and 1 for a good   **
    ! initial guess.                                                  **
    ! The dimension of the work array q should be NX * MY * 3.          **
    !
    real(8), intent(in), dimension(NX, MY) :: G, F
    real(8), intent(in) :: zk2
    integer, intent(in) :: NX1, NY1, M
    real(8), intent(in) :: H1
    integer, intent(in) :: NU1, NU2, NCYC, NX, MY
    real(8), intent(out), dimension(NX, MY) :: P

    integer :: ir
    integer :: k, ksq
    real(8) :: wu
    integer :: ic, km

    allocate(q(NX * MY * 3))

    iq = 1
    do k = 1, M
       ksq = 2 ** (k - 1)
       call grdfn(k, NX1 * ksq + 1, NY1 * ksq + 1, H1 / ksq)
       call grdfn(k + M, NX1 * ksq + 1, NY1 * ksq + 1, H1 / ksq)
    end do
    wu = 0.0d0
    call putf1(M, G, 0)
    call putf1(2 * M, F, 2)
    do ic = 1, NCYC
      do km = 1, M
         k = 1 + M - km
         if (k /= M) then
            call putz(k)
         end if
         do ir = 1, NU1
            call relax(k, k + M, wu, M, zk2)
         end do
         if (k > 1) then
            call rescal(k, k + M, k + M - 1, zk2)
         end if
      end do
      do k = 1, M
         do ir = 1, NU2
            call relax(k, k + M, wu, M, zk2)
         end do
         if (k < M) then
            call intadd(k, k + 1)
         end if
      end do
   end do
   call getsol(M, P)
   deallocate(q)
 end subroutine helmholtz

 subroutine grdfn(k, M, N, hh)
   integer, intent(in) :: k, M, N
   real(8), intent(in) :: hh

   nst(k) = iq
   imx(k) = M
   jmx(k) = N
   h(k) = hh
   iq = iq + M * N
 end subroutine grdfn

 subroutine key(k, ist, M, N, hh)
   integer, intent(in) :: k
   integer, dimension(:), intent(out) :: ist
   integer, intent(out) :: M, N
   real(8), intent(out) :: hh

   integer :: is, i

   M = imx(k)
   N = jmx(k)
   is = nst(k) - N - 1
   do i = 1, M
      is = is + N
      ist(i) = is
   end do
   hh = h(k)
 end subroutine key

 subroutine putf1(k, F, nh)
   integer, intent(in) :: k
   real(8), dimension(:, :), intent(in) :: F
   integer, intent(in) :: nh

   integer, dimension(200) :: ist

   integer :: ii, jj
   real(8) :: hh, hh2
   integer :: i, j

   call key (k, ist, ii, jj, hh)
   hh2 = hh ** nh
   do i = 1, ii
      do j = 1, jj
         q(ist(i) + j) = F(i, j) * hh2
      end do
   end do
 end subroutine putf1

 subroutine getsol(k, P)
   integer, intent(in) :: k
   real(8), dimension(:, :), intent(out) :: P

   integer, dimension(200) :: ist
   integer :: ii, jj
   real(8) :: hh
   integer :: i, j

   call key(k, ist, ii, jj, hh)
   do i = 1, ii
      do j = 1, jj
         P(i, j) = q(ist(i) + j)
      end do
   end do
 end subroutine getsol

 subroutine putz(k)
   integer, intent(in) :: k

   integer, dimension(200) :: ist
   integer :: ii, jj
   real(8) :: hh
   integer :: i, j

   call key( k, ist, ii, jj, hh)
   do  i = 1, ii
      do j = 1, jj
         q(ist(i) + j) = 0.0d0
      end do
   end do
 end subroutine putz

 subroutine relax(k, krhs, wu, M, zk2)
   integer, intent(in) :: k
   integer, intent(in) :: krhs
   real(8), intent(inout) :: wu
   integer, intent(in) :: M
   real(8), intent(in) :: zk2

   integer, dimension(200) :: ist, irhs
   integer :: ii, jj
   real(8) :: hh
   integer:: i1, j1, ir, iq, im, ip
   real(8) diag
   integer :: i, j
   real(8) :: a

   call key(k, ist, ii, jj, hh)
   call key(krhs, irhs, ii, jj, hh)
   diag = 4.0d0 + zk2 * hh ** 2
   i1 = ii - 1
   j1 = jj - 1
   do i = 2, i1
      ir = irhs(i)
      iq = ist(i)
      im = ist(i - 1)
      ip = ist(i + 1)
      do j = 2, j1
         a = q(ir + j) - q(iq + j + 1) - q(iq + j - 1) - q(im + j) - q(ip + j)
         q(iq + j) = - a / diag
      end do
   end do
   wu = wu + 4.0d0 ** (k - M)
 end subroutine relax

 subroutine intadd(kc, kf)
   integer, intent(in) :: kc
   integer, intent(in) :: kf

   integer, dimension(200) :: istc, istf
   real(8) :: hc, hf
   integer :: iic, jjc, iif, jjf
   integer :: ic, jc, if, jf
   integer :: ifo, ifm, ico, icm
   real(8) :: a, am

   call key(kc, istc, iic, jjc, hc)
   call key(kf, istf, iif, jjf, hf)
   do ic = 2, iic
      if = 2 * ic - 1
      jf = 1
      ifo = istf(if)
      ifm = istf(if - 1)
      ico = istc(ic)
      icm = istc(ic - 1)
      do jc = 2, jjc
         jf = jf + 2
         a = 0.5d0 * (q(ico + jc) + q(ico + jc - 1))
         am = 0.5d0 * (q(icm + jc) + q(icm + jc - 1))
         q(ifo+jf) = q(ifo + jf) + q(ico + jc)
         q(ifm + jf) = q(ifm + jf) + 0.5d0 * (q(ico + jc) + q(icm + jc))
         q(ifo + jf - 1) = q(ifo + jf - 1) + a
         q(ifm + jf - 1)  =  q(ifm + jf - 1) + 0.5d0 * (a + am)
      end do
   end do
 end subroutine intadd

 subroutine rescal(kf, krf, krc, zk2)
   integer, intent(in) :: kf, krf, krc
   real(8), intent(in) :: zk2

   integer, dimension(200) :: iuf, irf, irc
   integer :: iif, jjf, iic, jjc, iic1, jjc1
   integer :: if, ic, jf, jc
   integer :: icr, ifo, ifm, ifr, ifp
   real(8) :: hf, hc
   real(8) :: s
   real(8) :: diag

   call key(kf, iuf, iif, jjf, hf)
   call key(krf, irf, iif, jjf, hf)
   call key(krc, irc, iic, jjc, hc)

   diag = 4.0d0 + zk2 * hf ** 2
   iic1 = iic - 1
   jjc1 = jjc - 1
   do ic = 2, iic1
      icr = irc(ic)
      if = 2 * ic - 1
      jf = 1
      ifr = irf(if)
      ifo = iuf(if)
      ifm = iuf(if - 1)
      ifp = iuf(if + 1)
      do jc = 2, jjc1
         jf = jf + 2
         s = q(ifo + jf + 1) + q(ifo + jf - 1) + q(ifm + jf) + q(ifp + jf)
         q(icr + jc) = 4.0d0 * (q(ifr + jf) - s + diag * q(ifo + jf))
      end do
   end do
 end subroutine rescal

end module helmholtz_mod