var/gen_be/gen_be_vertloc.f90

   1 program gen_be_vertloc
   2
   3    use da_gen_be, only : da_eof_decomposition
   4
   5    implicit none
   6
   7    integer, parameter    :: ni = 1000
   8    integer, parameter    :: stdout = 6
   9    character (len=3)     :: cnk                       ! vertical level
  10    real*8, parameter       :: rscale = 0.1
  11    real*8, parameter       :: var_threshold = 0.99
  12
  13    integer               :: i, k, k1, m,nk,nk2
  14    integer               :: nm
  15    integer               :: ktarget
  16    real*8                  :: ni1_inv
  17    real*8                  :: nk_inv, nk3_inv
  18    real*8                  :: kscale, kscale_invsq
  19    real*8                  :: kdist
  20    real*8                  :: r
  21    real*8                  :: totvar, totvar_inv, cumvar
  22    real*8,allocatable      :: rho(:,:)
  23    real*8,allocatable      :: e(:,:)
  24    real,allocatable      :: cov(:,:)
  25    real*8,allocatable      :: eval(:)
  26    real*8,allocatable      :: evec(:,:)
  27    real*8,allocatable      :: v(:)
  28    real*8,allocatable      :: x(:)
  29
  30   cnk=""
  31   call getarg( 1, cnk )
  32   read(cnk,'(i3)')nk2
  33   nk=nk2-1
  34   write(stdout,'(a,i6)') ' vertical level = ', nk2
  35
  36   allocate(rho(1:nk,1:nk))
  37   allocate(e(1:ni,1:nk))
  38   allocate(cov(1:nk,1:nk))
  39   allocate(eval(1:nk))
  40   allocate(evec(1:nk,1:nk))
  41   allocate(v(1:nk))
  42   allocate(x(1:nk))
  43    ni1_inv = 1.0 / real (ni - 1)
  44    nk_inv = 1.0 / real (nk)
  45    nk3_inv = 1.0 / real (nk-3)
  46
  47 !  Specify probability densities:
  48    do k = 1, nk
  49       kscale = 10.0 * real(k) * nk_inv  ! Very simple parametrization of vertical variation of localization scale.
  50 !      kscale = 3.0 + 7.0 * real(k-3) * nk3_inv  ! Very simple parametrization of vertical variation of localization scale.
  51 !      kscale = 10.0
  52       kscale_invsq = 1.0 / ( kscale * kscale )
  53       do k1 = k, nk
  54          kdist = k1 - k
  55          rho(k,k1) = exp ( -real(kdist * kdist) * kscale_invsq )
  56          rho(k1,k) = rho(k,k1)
  57       end do
  58    end do
  59    cov = rho
  60
  61 !  Create random variables:
  62 !   k = 1
  63 !   do i = 1, ni
  64 !      do k1 = 1, nk
  65 !         call da_gauss_noise( r )
  66 !         r = rscale * r
  67 !         e(i,k1) = rho(k,k1) + r
  68 !         write(17,'(i5,f15.5)')k1, e(i,k1)
  69 !      end do
  70 !      k = k + 1
  71 !      if ( k > nk ) k = 1
  72 !   end do
  73
  74 !  Calculate covariance matrix:
  75
  76 !   do k = 1, nk
  77 !      do k1 = k, nk
  78 !         cov(k,k1) = ni1_inv * sum( e(1:ni,k) * e(1:ni,k1) )
  79 !         cov(k,k1) = rho(k,k1)
  80 !         cov(k1,k) = cov(k,k1)
  81 !      end do
  82 !   end do
  83
  84 !------------------------------------------------------------------------------------------------
  85 ! Calculate eigenvectors/values:
  86 !------------------------------------------------------------------------------------------------
  87
  88    call da_eof_decomposition( nk, cov, evec, eval )
  89
  90 !  Eliminate negative eigenvalues:
  91    totvar = sum(eval(1:nk))
  92    do k = 1, nk
  93      if ( eval(k) < 0.0 ) eval(k) = 0.0
  94    end do
  95    totvar_inv = 1.0 / sum(eval(1:nk))
  96    eval(:) = eval(:) * totvar * totvar_inv ! Preserve total variance before -ve values removed.
  97
  98 !  Calculate truncation:
  99    nm = 0
 100    totvar_inv = 1.0 / sum(eval(1:nk))
 101    do k = 1, nk
 102       cumvar = sum(eval(1:k)) * totvar_inv
 103       if ( nm == 0 .and. cumvar >= var_threshold ) nm = k
 104    end do
 105
 106    write(stdout,'(a,f15.5,i6)')' Threshold, truncation = ', var_threshold, nm
 107
 108    open(10, file = 'be.vertloc.dat', status = 'unknown', form='unformatted')
 109    write(10) nk
 110    eval(:) = sqrt(eval(:)) ! Write out L^{1/2} for use in WRF-Var
 111    write(10)eval(1:nk)
 112    write(10)evec(1:nk,1:nk)
 113
 114    close(10)
 115
 116 !------------------------------------------------------------------------------------------------
 117 ! Calculate B = E L E^T = U U^T [..00100..]
 118 !------------------------------------------------------------------------------------------------
 119
 120    ktarget = nk
 121    x(:) = 0.0
 122    x(ktarget) = 1.0
 123
 124 !  v = L^1/2 E^T x:
 125    do m = 1, nm
 126       v(m) = eval(m) * sum(evec(1:nk,m) * x(1:nk))
 127    end do
 128
 129 !  x = E L^1/2 v:
 130    do k = 1, nk
 131       x(k) = sum(evec(k,1:nm) * eval(1:nm) * v(1:nm))
 132 !      write(22,'(i5,2f15.5)')k, rho(ktarget,k), x(k)
 133    end do
 134
 135   deallocate(rho)
 136   deallocate(e)
 137   deallocate(cov)
 138   deallocate(eval)
 139   deallocate(evec)
 140   deallocate(v)
 141   deallocate(x)
 142
 143 end program gen_be_vertloc
 144