var/da/da_wavelet/da_transform_through_wavelet_adj.inc

   1 SUBROUTINE da_transform_through_wavelet_adj(grid,mz,wsd,v,vv)
   2 !
   3 ! Purpose: Horizontal control-var. transform-adjoint vv -> v through dwtai_w.c.
   4 ! Author: Aime' Fournier
   5 !
   6 ! v(:nq,m) = MATMUL(TRANSPOSE(uh(:,:,m)),vv(:ni*nj,m)), m the vert. mode,
   7 ! uh(ij,q,m) = w(q,ij)*wsd(q,m)/p_x(ij) the m'th horiz. CVT matrix,
   8 ! w(q,:ni*nj) the q'th of nq horizontal 2D wavelets,
   9 ! wsd(q,m) the standard deviation of the q'th wavelet coefficient,
  10 ! p_x(ij)**2 the horizontal area in ij.
  11 !
  12 ! Method:
  13 !  3) Diagonal divide by p_x(:ni*nj), cf. gen_be_stage4_regional.
  14 !  2) MatMul by w(:nq,:ni*nj) using dwtai_w();
  15 !  1) Diagonal multiply by wsd(:nq,m);
  16 !
  17  IMPLICIT NONE
  18  TYPE(domain), INTENT( IN)::grid
  19  REAL,         INTENT( IN)::vv(ims:ime,jms:jme,kms:kme)         ! Grid point/EOF equivalent.
  20  INTEGER,      INTENT( IN)::mz                                  ! Vertical truncation.
  21  REAL,         INTENT( IN)::wsd(nij(0,1,2),nij(0,0,2),mz)       ! Wavelet standard deviation.
  22  REAL,         INTENT(OUT)::v(nij(0,0,2)*nij(0,1,2)*mz)         ! Field to be transformed.
  23  INTEGER                  ::dj(0:1),dk(0:1),dv,i,j,m,n
  24  LOGICAL, SAVE            ::call1=.TRUE.
  25
  26 #ifdef WAVELET
  27
  28  v=0.
  29  dj = nij(0,1,0:2:2)                    ! pseudo-eastward lengths.
  30  dk(0) = (nij(0,0,0)-1)*dj(1)+dj(0)     ! horizontal model-var stride.
  31  dk(1) = nij(0,0,2)*dj(1)               ! horizontal control-var size.
  32 !
  33 ! 1D indexing of 3D ("//" & "/" for model- & cv-space):
  34 !             _______________________________________
  35 !            /m+dk1-dj1/ ...  /       / ... /m+dk1-1/
  36 !           /_________/______/_______/_____/_______/
  37 !          /         /      /       /     /       /__
  38 !         /_________/______/_______/_____/_______/  /
  39 !        /m+dk0-dj0/ ...  /m+dk0-1/ ... /       /__/
  40 !       //_______//______//_____//_____/_______/  /__
  41 !      //...    //      //     //     /       /__/  /
  42 !     //_______//______//_____//_____/_______/  /__/
  43 !    //m+dj1  //      //     //     /       /__/  /
  44 !   //_______//______//_____//_____/_______/  /__/
  45 !  //m      // ...  /m+dj0-1/ ... /m+dj1-1/__/  /
  46 ! //_______//______//_____//_____/_______/  /__/
  47 !   //_______//______//_____//_____/_______/  /
  48 !  //m-dk1  //      //     //     /       /__/
  49 ! //_______//______//_____//_____/_______/  /
  50 !   //_______//______//_____//_____/_______/
  51 !  //...    //      //     //     /       /
  52 ! //_______//______//_____//_____/_______/
  53 !
  54  dv = mz*dk(1)                          ! variable 3D length.
  55  n = 1                                  ! n <= 1 + mz:
  56  DO m = 1, 1+dv-dk(1), dk(1)            ! mz vertical-modes loop:
  57     j = m                               ! m <= j <= m+dk0-dj0:
  58     DO i = jts, jts+nij(0,0,0)-1        ! Pseudo-northward loop:
  59        v(j:j+dj(0)-1) = vv(its:ite,i,n)
  60 !      [3.0] Apply box-area factor:
  61        v(j:j+dj(0)-1) = v(j:j+dj(0)-1)/SQRT(grid%xb%grid_box_area(1:dj(0),i))
  62        j = j+dj(1)
  63     ENDDO
  64
  65 !   [2.1] Perform dwtai() in pseudo-northward directions:
  66     DO i = m, m+dj(0)-1                 ! Pseudo-eastward loop:
  67        CALL dwtai_w(namw//CHAR(0),lf,v(i:i+dk(1)-dj(1):dj(1)),ws,nij(:,0,0),nij(:,0,2),nb)
  68     ENDDO
  69
  70 !   [2.0] Perform dwtai() in pseudo-eastward directions:
  71     DO j = m, m+dk(1)-dj(1), dj(1)      ! Pseudo-northward loop:
  72        CALL dwtai_w(namw//CHAR(0),lf,v(j:j+dj(1)-   1       ),ws,nij(:,1,0),nij(:,1,2),nb)
  73     ENDDO
  74 !   [1.0] Multiply by wavelet standard deviations:
  75     v(m:m+dk(1)-1) = RESHAPE(wsd(:,:,n),(/dk(1)/))*v(m:m+dk(1)-1)
  76     n = n+1                             ! Increment mode.
  77  ENDDO                                  ! m (& n) loop.
  78
  79  IF( call1 )THEN
  80     PRINT'(a,": ",ES7.1,"<grid_box_area~",ES7.1,"<",ES7.1)', &
  81        __FILE__,MINVAL(grid%xb%grid_box_area(1:dj(0),jts:jts+nij(0,0,0)-1)), &
  82               SQRT(SUM(grid%xb%grid_box_area(1:dj(0),jts:jts+nij(0,0,0)-1)**2)/(dj(0)*nij(0,0,0))), &
  83                 MAXVAL(grid%xb%grid_box_area(1:dj(0),jts:jts+nij(0,0,0)-1))
  84     call1=.FALSE.
  85  ENDIF
  86
  87 #else
  88     PRINT'(a,": must compile after setenv WAVELET 1")',__FILE__
  89     CALL wrf_abort
  90 #endif
  91
  92 ENDSUBROUTINE da_transform_through_wavelet_adj