var/da/da_wavelet/da_transform_through_wavelet.inc

   1  SUBROUTINE da_transform_through_wavelet(grid,mz,wsd,v,vv)
   2 !
   3 ! Purpose: Horizontal control-variable transform v -> vv through idwtai_w.c.
   4 ! Author: Aime' Fournier
   5 !
   6 ! vv(:ni*nj,m) = MATMUL(uh(:,:,m),v(:nq,m)), m the vertical 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 !  1) Diagonal multiply by wsd(:nq,m);
  14 !  2) MatMul by Transpose(w(:nq,:ni*nj)) using idwtai_w();
  15 !  3) Diagonal divide by p_x(:ni*nj), cf. gen_be_stage4_regional.
  16 !
  17  IMPLICIT NONE
  18  TYPE(domain),  INTENT( IN)::grid
  19  REAL,          INTENT(OUT)::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( IN)::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  REAL                      ::u(nij(0,1,2),nij(0,0,2))           ! Since INTENT(IN)::v in caller.
  26
  27 #ifdef WAVELET
  28
  29  vv=0.
  30  dj = nij(0,1,0:2:2)                    ! pseudo-eastward lengths.
  31  dk(0) = (nij(0,0,0)-1)*dj(1)+dj(0)     ! horizontal model-var stride.
  32  dk(1) = nij(0,0,2)*dj(1)               ! horizontal control-var size.
  33 !
  34 ! 1D indexing of 3D ("//" & "/" for model- & cv-space):
  35 !             _______________________________________
  36 !            /m+dk1-dj1/ ...  /       / ... /m+dk1-1/
  37 !           /_________/______/_______/_____/_______/
  38 !          /         /      /       /     /       /__
  39 !         /_________/______/_______/_____/_______/  /
  40 !        /m+dk0-dj0/ ...  /m+dk0-1/ ... /       /__/
  41 !       //_______//______//_____//_____/_______/  /__
  42 !      //...    //      //     //     /       /__/  /
  43 !     //_______//______//_____//_____/_______/  /__/
  44 !    //m+dj1  //      //     //     /       /__/  /
  45 !   //_______//______//_____//_____/_______/  /__/
  46 !  //m      // ...  /m+dj0-1/ ... /m+dj1-1/__/  /
  47 ! //_______//______//_____//_____/_______/  /__/
  48 !   //_______//______//_____//_____/_______/  /
  49 !  //m-dk1  //      //     //     /       /__/
  50 ! //_______//______//_____//_____/_______/  /
  51 !   //_______//______//_____//_____/_______/
  52 !  //...    //      //     //     /       /
  53 ! //_______//______//_____//_____/_______/
  54 !
  55  dv = mz*dk(1)                          ! variable 3D length.
  56  n = 1                                  ! n <= 1 + mz:
  57  DO m = 1, 1+dv-dk(1), dk(1)            ! mz vertical-modes loop:
  58     u=RESHAPE(v(m:m+dk(1)-1),nij(0,1:0:-1,2))
  59 !   [1.0] Multiply by wavelet standard deviations:
  60     u=wsd(:,:,n)*u
  61 !   [2.0] Perform idwtai() in pseudo-eastward directions:
  62     DO j = 1, nij(0,0,2)                ! Pseudo-northward loop:
  63        CALL idwtai_w(namw//CHAR(0),lf,u(:,j),ws,nij(:,1,0),nij(:,1,1),nij(:,1,2),nb)
  64     ENDDO
  65
  66 !   [2.1] Perform idwtai() in pseudo-northward directions:
  67     DO i = 1, dj(0)                     ! Pseudo-eastward loop:
  68        CALL idwtai_w(namw//CHAR(0),lf,u(i,:),ws,nij(:,0,0),nij(:,0,1),nij(:,0,2),nb)
  69     ENDDO
  70
  71     DO i = jts, jts+nij(0,0,0)-1        ! Pseudo-northward loop:
  72 !      [3.0] Apply box-area factor:
  73        u(:dj(0),i) = u(:dj(0),i)/SQRT(grid%xb%grid_box_area(1:dj(0),i))
  74        vv(its:ite,i,n) = u(:dj(0),i)
  75     ENDDO
  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