var/da/da_setup_structures/da_get_alpha_vertloc.inc

   1 subroutine da_get_alpha_vertloc (xb, alpha_val, alpha_evec)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: To get vertical localization information for alpha control variable
   5    ! History:
   6    !    07/2019 - initial version
   7    !-----------------------------------------------------------------------
   8
   9    implicit none
  10
  11    type (xb_type), intent(in)  :: xb               ! first guess structure
  12    real*8, intent(inout)       :: alpha_val(:)     ! alpha_vertloc eigenvalues
  13    real*8, intent(inout)       :: alpha_evec(:,:)  ! alpha_vertloc eigenvectors
  14
  15    character(len=filename_len) :: filename
  16    logical                     :: fexist
  17    integer                     :: iunit, ier
  18    integer                     :: kz, nk, k, k1, i
  19    real*8                      :: kscale, kscale_invsq, kdist
  20    real*8                      :: totvar, totvar_inv
  21    real*8                      :: cutoff
  22    real,   allocatable         :: cov(:,:)
  23    real*8, allocatable         :: rho(:,:)
  24    real*8, allocatable         :: eval(:)
  25    real*8, allocatable         :: evec(:,:)
  26    real*8, allocatable         :: p_w(:)  ! pressure at full levels
  27    real*8, allocatable         :: ln_p_w(:)
  28    real*8, allocatable         :: dlnp(:,:)
  29    real*8, allocatable         :: kcutoff(:), kcutoff_smth(:)
  30
  31    kz = xb%mkz
  32
  33    if ( size(alpha_val(:)) /= kz ) then
  34       ! this should not happen because dimension check is already handled in
  35       ! the parent subroutine da_setup_be_regional
  36       write(unit=message(1),fmt='(a,a)')  &
  37             'Vertical dimension of the assimilation domain do not match that in ', trim(filename)
  38       call da_error(__FILE__,__LINE__,message(1:1))
  39    end if
  40
  41    filename = 'be.vertloc.dat'
  42
  43    select case ( alpha_vertloc_opt )
  44
  45       case ( 1 )  ! reading in from be.vertloc.dat
  46
  47          inquire (file=trim(filename), exist=fexist)
  48          if ( .not. fexist ) then
  49             write(unit=message(1),fmt='(a,a)') trim(filename), ' not found for reading vertloc info.'
  50             write(unit=message(2),fmt='(a)')   'Either make sure it exists in the working directory, or'
  51             write(unit=message(3),fmt='(a)')   'set alpha_vertloc_opt=2 (recommended) in namelist.input (&wrfvar16) to have WRFDA calculate log-P based vertloc.'
  52             call da_error(__FILE__,__LINE__,message(1:3))
  53          end if
  54          call da_get_unit(iunit)
  55          open(unit=iunit,file=trim(filename), status='old', form='unformatted')
  56
  57          write(unit=message(1),fmt='(a)') ''
  58          write(unit=message(2),fmt='(a,a)') 'alpha_vertloc_opt=1: reading alpha_vertloc info from ', trim(filename)
  59          call da_message(message(1:2))
  60
  61          read (iunit, iostat=ier) nk
  62          if (ier /= 0) then
  63             write (unit=message(1),fmt='(a,a)') 'Error in reading ', trim(filename)
  64             call da_error(__FILE__,__LINE__,message(1:1))
  65          end if
  66          if ( nk /= kz ) then
  67             write(unit=message(1),fmt='(a,a)')  &
  68                'Vertical dimension of the assimilation domain do not match that in ', trim(filename)
  69             write(unit=message(2),fmt='(3x,a,a,i4)') 'in ', trim(filename), ': nk = ', nk
  70             write(unit=message(3),fmt='(3x,a,i4)')   'in fg:        , kz = ', kz
  71             call da_error(__FILE__,__LINE__,message(1:4))
  72          end if
  73          read (iunit) alpha_val(1:nk)
  74          read (iunit) alpha_evec(1:nk,1:nk)
  75          close(iunit)
  76          call da_free_unit(iunit)
  77
  78       case ( 2 ) ! log-P method
  79
  80          allocate(p_w(1:kz+1))
  81          allocate(ln_p_w(1:kz+1))
  82          allocate(dlnp(1:kz,1:kz))
  83          allocate(kcutoff_smth(1:kz))
  84          allocate(kcutoff(1:kz))
  85
  86          ! empirical settings
  87          cutoff   = 1.0/2.71828 !1/e
  88
  89          do k = 1, kz+1
  90             p_w(k) = xb%znw(k)*(base_pres - xb%ptop) + xb%ptop
  91             ln_p_w(k) = log(max(p_w(k),0.0001))
  92          end do
  93          kcutoff(:) = 1.0  ! initialize
  94          do k = 1, kz
  95             do k1 = k+1, kz
  96                dlnp(k,k1) = abs(ln_p_w(k)-ln_p_w(k1))
  97                if ( dlnp(k,k1) <= cutoff ) then
  98                   kcutoff(k) = k1-k+1
  99                end if
 100             end do
 101          end do
 102          ! smoothing
 103          kcutoff_smth(:) = kcutoff(:)
 104          do i = 1, 2 ! two passes
 105             do k = 2, kz-1
 106                kcutoff_smth(k)  = kcutoff(k)+ 0.25 * ( kcutoff(k-1) + kcutoff(k+1) -2.0*kcutoff(k) )
 107             end do
 108             kcutoff(:) = kcutoff_smth(:)
 109          end do
 110          deallocate(kcutoff_smth)
 111          deallocate(dlnp)
 112          deallocate(ln_p_w)
 113          deallocate(p_w)
 114
 115          allocate(rho(1:kz,1:kz))
 116          allocate(cov(1:kz,1:kz))
 117          allocate(eval(1:kz))
 118          allocate(evec(1:kz,1:kz))
 119
 120          ! specify probability densities:
 121          do k = 1, kz
 122             kscale = kcutoff(k)
 123             kscale_invsq = 1.0 / ( kscale * kscale )
 124             do k1 = k, kz
 125                kdist = k1 - k
 126                rho(k,k1) = exp ( -1.0 * real(kdist * kdist) * kscale_invsq )
 127                rho(k1,k) = rho(k,k1)
 128             end do
 129          end do
 130          cov = rho
 131
 132          deallocate(kcutoff)
 133
 134          ! calculate eigenvectors/values:
 135
 136          call da_eof_decomposition( kz, cov, evec, eval )
 137
 138          ! eliminate negative eigenvalues
 139          totvar = sum(eval(1:kz))
 140          do k = 1, kz
 141             if ( eval(k) < 0.0 ) eval(k) = 0.0
 142          end do
 143          totvar_inv = 1.0 / sum(eval(1:kz))
 144          eval(:) = eval(:) * totvar * totvar_inv ! preserve total variance before -ve values removed.
 145
 146          eval(:) = sqrt(eval(:))
 147
 148          alpha_val(:)    = eval(:)
 149          alpha_evec(:,:) = evec(:,:)
 150
 151          write(unit=message(1),fmt='(a)') ''
 152          write(unit=message(2),fmt='(a,a)') 'alpha_vertloc_opt=2: calculate and write out alpha_vertloc info to ', trim(filename)
 153          call da_message(message(1:2))
 154
 155          ! writing it out for diagnostics
 156          call da_get_unit(iunit)
 157          open(unit=iunit,file=trim(filename), status='replace', form='unformatted')
 158          write(iunit) kz
 159          write(iunit) alpha_val(1:kz)
 160          write(iunit) alpha_evec(1:kz,1:kz)
 161          close(iunit)
 162          call da_free_unit(iunit)
 163
 164          deallocate(evec)
 165          deallocate(eval)
 166          deallocate(cov)
 167          deallocate(rho)
 168
 169       case default
 170
 171          write (unit=message(1),fmt='(a,i3)') ' Unknown option alpha_vertloc_opt = ', alpha_vertloc_opt
 172          write (unit=message(2),fmt='(a)')    ' Please set alpha_vertloc_opt = 1 or alpha_vertloc_opt = 2.'
 173          call da_error(__FILE__,__LINE__,message(1:2))
 174
 175    end select
 176
 177 end subroutine da_get_alpha_vertloc