wrfv2_fire/var/da/da_vtox_transforms/da_vertical_transform.inc

   1 subroutine da_vertical_transform(grid, string, be, vertical_wgt, vv, vp)
   2
   3    !---------------------------------------------------------------------
   4    ! Purpose: TBD
   5    !---------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (domain),    intent(in)    :: grid
  10    character(len=*), intent(in)    :: string      ! Character operation
  11    type (be_type),   intent(in)    :: be          ! Background error structure.
  12    real,             intent(in)    :: vertical_wgt(ims:ime,jms:jme,kms:kme) ! Weighting.
  13    type (vp_type),   intent(inout) :: vv          ! CV in gridpt/EOF space.
  14    type (vp_type),   intent(inout) :: vp          ! CV in gridpt/level space.
  15
  16    integer                         :: j, m, n     ! Loop counters.
  17    real                            :: alpha_stddev_inv ! 1/ sigma_alpha
  18    real                            :: size_inv    ! 1 / size.
  19    real                            :: alpha_me, alpha_ms, alpha_sd ! Alpha statistics.
  20
  21    if (trace_use) call da_trace_entry("da_vertical_transform")
  22
  23    select case(string)
  24
  25    case ('u');
  26
  27       !-------------------------------------------------------------------
  28       ! [1.0] Perform vp(i,j,k) = E L^{1/2} vv(i,j,m) transform:
  29       !-------------------------------------------------------------------
  30
  31       if (be % v1 % mz > 0) then
  32          call da_transform_vvtovp (grid, be % v1 % evec, be % v1 % val, vertical_wgt, &
  33             vv % v1, vp % v1, be % v1 % mz, kte)
  34       else
  35          vp % v1(its:ite,jts:jte,kts:kte) = 0.0
  36       end if
  37
  38       if (be % v2 % mz > 0) then
  39          call da_transform_vvtovp (grid, be % v2 % evec, be % v2 % val, vertical_wgt, &
  40             vv % v2, vp % v2, be % v2 % mz, kte)
  41       else
  42          vp % v2(its:ite,jts:jte,kts:kte) = 0.0
  43       end if
  44
  45       if (be % v3 % mz > 0) then
  46          call da_transform_vvtovp (grid, be % v3 % evec, be % v3 % val, vertical_wgt, &
  47             vv % v3, vp % v3, be % v3 % mz, kte)
  48       else
  49          vp % v3(its:ite,jts:jte,kts:kte) = 0.0
  50       end if
  51
  52       if (be % v4 % mz > 0) then
  53          call da_transform_vvtovp (grid, be % v4 % evec, be % v4 % val, vertical_wgt, &
  54             vv % v4, vp % v4, be % v4 % mz, kte)
  55       else
  56          vp % v4(its:ite,jts:jte,kts:kte) = 0.0
  57       end if
  58
  59       if (be % v5 % mz > 0) then
  60          if (global) then
  61             vp % v5(its:ite,jts:jte,1) = vv % v5(its:ite,jts:jte,1)
  62          else
  63          call da_transform_vvtovp (grid, be % v5 % evec, be % v5 % val, vertical_wgt, &
  64             vv % v5, vp % v5, be % v5 % mz, kts)
  65          end if
  66       else
  67          vp % v5(its:ite,jts:jte,kts:kts) = 0.0
  68       end if
  69
  70       if ( be % ne > 0 .and. be % alpha % mz > 0 ) then
  71          do n = 1, be % ne
  72             call da_transform_vvtovp (grid, be % alpha % evec, be % alpha % val, vertical_wgt, &
  73                                        vv % alpha(:,:,:,n), vp % alpha(:,:,:,n), be % alpha % mz, kte)
  74          end do
  75
  76 !        Calculate alpha standard deviation diagnostic:
  77 !         size_inv = 1.0 / ( (ite-its+1) * ( jte-jts+1) * be % ne * be % alpha % mz )
  78 !         alpha_me = sum(vp % alpha(its:ite,jts:jte,:)) * size_inv
  79 !         alpha_ms = sum(vp % alpha(its:ite,jts:jte,:) * vp % alpha(its:ite,jts:jte,:)) * &
  80 !                    size_inv
  81 !         alpha_sd = sqrt( alpha_ms - alpha_me * alpha_me )
  82 !         write(6,'(a,f15.5)')' Alpha std. dev = ', alpha_sd
  83       end if
  84
  85    case ('u_inv');
  86
  87       !-------------------------------------------------------------------
  88       ! [2.0] Perform vv(i,j,m) = L^{-1/2} E^T vp(i,j,k) transform:
  89       !-------------------------------------------------------------------
  90
  91       if (be % v1 % mz > 0) then
  92          call da_transform_vptovv (be % v1 % evec, be % v1 % val, vertical_wgt, &
  93             vp % v1, vv % v1, be % v1 % mz, kds,kde, ims,ime, jms,jme, kms,kme, &
  94             its,ite, jts,jte, kts,kte)
  95       end if
  96
  97       if (be % v2 % mz > 0) then
  98          call da_transform_vptovv (be % v2 % evec, be % v2 % val, vertical_wgt, &
  99             vp % v2, vv % v2, be % v2 % mz, kds,kde, ims,ime, jms,jme, kms,kme, &
 100             its,ite, jts,jte, kts,kte)
 101       end if
 102
 103       if (be % v3 % mz > 0) then
 104          call da_transform_vptovv (be % v3 % evec, be % v3 % val, vertical_wgt, &
 105             vp % v3, vv % v3, be % v3 % mz, kds,kde, ims,ime, jms,jme, kms,kme, &
 106             its,ite, jts,jte, kts,kte)
 107       end if
 108
 109       if (be % v4 % mz > 0) then
 110          call da_transform_vptovv (be % v4 % evec, be % v4 % val, vertical_wgt, &
 111             vp % v4, vv % v4, be % v4 % mz, kds,kde, ims,ime, jms,jme, kms,kme, &
 112             its,ite, jts,jte, kts,kte)
 113       end if
 114
 115       if (be % v5 % mz > 0) then
 116          if (global) then
 117             vv % v5(its:ite,jts:jte,1) = vp % v5(its:ite,jts:jte,1)
 118          else
 119             call da_transform_vptovv (be % v5 % evec, be % v5 % val, vertical_wgt, &
 120                vp % v5, vv % v5, be % v5 % mz, kds,kde, ims,ime, jms,jme, kms,kme, &
 121                its,ite, jts,jte, kts,kte)
 122          end if
 123       end if
 124
 125       if ( be % ne > 0 .and. be % alpha % mz > 0 ) then
 126          do n = 1, be % ne
 127             call da_transform_vptovv (be % alpha % evec, be % alpha % val, vertical_wgt, &
 128                                       vp % alpha(:,:,:,n), vv % alpha(:,:,:,n), be % alpha % mz, kds,kde, &
 129                                       ims,ime, jms,jme, kms,kme, its,ite, jts,jte, kts,kte)
 130          end do
 131       end if
 132
 133    case ('u_adj');
 134
 135       !-------------------------------------------------------------------
 136       ! [3.0] Perform vv_adj = U_{v}^{T} vp_adj transform:
 137       !-------------------------------------------------------------------
 138
 139       if (be % v1 % mz > 0) then
 140          call da_transform_vvtovp_adj (grid, be % v1 % evec, be % v1 % val, vertical_wgt, &
 141             vp % v1, vv % v1, be % v1 % mz, kte)
 142       end if
 143
 144       if (be % v2 % mz > 0) then
 145          call da_transform_vvtovp_adj (grid, be % v2 % evec, be % v2 % val, vertical_wgt, &
 146             vp % v2, vv % v2, be % v2 % mz, kte)
 147       end if
 148
 149       if (be % v3 % mz > 0) then
 150          call da_transform_vvtovp_adj (grid, be % v3 % evec, be % v3 % val, vertical_wgt, &
 151             vp % v3, vv % v3, be % v3 % mz, kte)
 152       end if
 153
 154       if (be % v4 % mz > 0) then
 155          call da_transform_vvtovp_adj (grid, be % v4 % evec, be % v4 % val, vertical_wgt, &
 156             vp % v4, vv % v4, be % v4 % mz, kte)
 157       end if
 158
 159       if (be % v5 % mz > 0) then
 160          if (global) then
 161             vv % v5(its:ite,jts:jte,1) = vp % v5(its:ite,jts:jte,1)
 162          else
 163             call da_transform_vvtovp_adj (grid, be % v5 % evec, be % v5 % val, vertical_wgt, &
 164                vp % v5, vv % v5, be % v5 % mz, kts)
 165          end if
 166       end if
 167
 168       if ( be % ne > 0 .and. be % alpha % mz > 0 ) then
 169          do n = 1, be % ne
 170             call da_transform_vvtovp_adj (grid, be % alpha % evec, be % alpha % val, vertical_wgt, &
 171                                            vp % alpha(:,:,:,n), vv % alpha(:,:,:,n), be % alpha % mz, kte)
 172          end do
 173       end if
 174
 175    case default;
 176
 177       call da_error(__FILE__,__LINE__, &
 178          (/"Invalid da_vertical_transform option "//string/))
 179
 180    end select
 181
 182    if (trace_use) call da_trace_exit("da_vertical_transform")
 183
 184 end subroutine da_vertical_transform
 185
 186