var/da/da_test/da_check_vvtovp_adjoint.inc

   1 subroutine da_check_vvtovp_adjoint(grid, ne, xb, be, vv, vp)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Test Vv to Vp routine and adjoint for compatibility.
   5    !
   6    ! Method:  Standard adjoint test: < Vp, Vp > = < Vv_adj, Vv >.
   7    !---------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    type (domain), intent(in)         :: grid
  12    integer, intent(in)               :: ne    ! Ensemble size.
  13    type (xb_type), intent(in)        :: xb    ! first guess (local).
  14    type (be_type), intent(in)        :: be    ! background error structure.
  15    type (vp_type), intent(inout)     :: vv    ! CV(i,j,m).
  16    type (vp_type), intent(inout)     :: vp    ! CV(i,j,k)
  17
  18    real                              :: adj_par_lhs ! < x, x >
  19    real                              :: adj_par_rhs ! < v_adj, v >
  20    real                              :: adj_sum_lhs ! < x, x >
  21    real                              :: adj_sum_rhs ! < v_adj, v >
  22
  23    real                              :: vv2_v1(ims:ime,jms:jme,kms:kme)
  24    real                              :: vv2_v2(ims:ime,jms:jme,kms:kme)
  25    real                              :: vv2_v3(ims:ime,jms:jme,kms:kme)
  26    real                              :: vv2_v4(ims:ime,jms:jme,kms:kme)
  27    real                              :: vv2_v5(ims:ime,jms:jme,kms:kme)
  28 !  real                              :: vv2_alpha(ims:ime,jms:jme,kts:kte,1:ne)
  29    real                              :: vv2_alpha(ims_int:ime_int,jms_int:jme_int,kts_int:kte_int,1:ne)
  30
  31    if (trace_use) call da_trace_entry("da_check_vvtovp_adjoint")
  32
  33    if (cv_options == 3 ) then
  34       write(unit=stdout, fmt='(/a,i2,a/)') 'cv_options =',cv_options, &
  35                      '   no da_check_vvtovp_adjoint check...'
  36       goto 1235
  37    end if
  38
  39    !----------------------------------------------------------------------
  40    ! [1.0] Initialise:
  41    !----------------------------------------------------------------------
  42
  43    write(unit=stdout, fmt='(/a/)') 'da_check_vvtovp_adjoint: Test Results:'
  44
  45    !----------------------------------------------------------------------
  46    ! [2.0] Perform Vp = U_v Vv transform:
  47    !----------------------------------------------------------------------
  48
  49    call da_vertical_transform(grid, 'u', be, &
  50                                xb % vertical_inner_product, &
  51                                vv, vp)
  52
  53    !----------------------------------------------------------------------
  54    ! [3.0] Calculate LHS of adjoint test equation:
  55    !----------------------------------------------------------------------
  56
  57    adj_par_lhs = sum(vp % v1(its:ite,jts:jte,kts:kte)**2) * inv_typ_vp1_sumsq &
  58                + sum(vp % v2(its:ite,jts:jte,kts:kte)**2) * inv_typ_vp2_sumsq &
  59                + sum(vp % v3(its:ite,jts:jte,kts:kte)**2) * inv_typ_vp3_sumsq &
  60                + sum(vp % v4(its:ite,jts:jte,kts:kte)**2) * inv_typ_vp4_sumsq &
  61                + sum(vp % v5(its:ite,jts:jte,kts:kte)**2) * inv_typ_vp5_sumsq
  62
  63    if (be % ne > 0) then
  64       adj_par_lhs = adj_par_lhs + &
  65 !        sum(vp % alpha(its:ite,jts:jte,kts:kte,1:be%ne)**2) * inv_typ_vpalpha_sumsq
  66          sum(vp % alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:be%ne)**2) * inv_typ_vpalpha_sumsq
  67    end if
  68
  69    !----------------------------------------------------------------------
  70    ! [4.0] Rescale input to adjoint routine:
  71    !----------------------------------------------------------------------
  72
  73    vp % v1(its:ite,jts:jte,kts:kte) = vp % v1(its:ite,jts:jte,kts:kte) * &
  74       inv_typ_vp1_sumsq
  75    vp % v2(its:ite,jts:jte,kts:kte) = vp % v2(its:ite,jts:jte,kts:kte) * &
  76       inv_typ_vp2_sumsq
  77    vp % v3(its:ite,jts:jte,kts:kte) = vp % v3(its:ite,jts:jte,kts:kte) * &
  78       inv_typ_vp3_sumsq
  79    vp % v4(its:ite,jts:jte,kts:kte) = vp % v4(its:ite,jts:jte,kts:kte) * &
  80       inv_typ_vp4_sumsq
  81    vp % v5(its:ite,jts:jte,kts:kte) = vp % v5(its:ite,jts:jte,kts:kte) * &
  82       inv_typ_vp5_sumsq
  83
  84    if (be % ne > 0) then
  85 !     vp % alpha(its:ite,jts:jte,kts:kte,1:be%ne) = &
  86 !        vp % alpha(its:ite,jts:jte,kts:kte,1:be%ne) * inv_typ_vpalpha_sumsq
  87       vp % alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:be%ne) = &
  88          vp % alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:be%ne) * inv_typ_vpalpha_sumsq
  89    end if
  90
  91    !----------------------------------------------------------------------
  92    ! [5.0] Perform adjoint operation:
  93    !----------------------------------------------------------------------
  94
  95    vv2_v1(its:ite,jts:jte,1:be%v1%mz) = vv % v1(its:ite,jts:jte,1:be%v1%mz)
  96    vv2_v2(its:ite,jts:jte,1:be%v2%mz) = vv % v2(its:ite,jts:jte,1:be%v2%mz)
  97    vv2_v3(its:ite,jts:jte,1:be%v3%mz) = vv % v3(its:ite,jts:jte,1:be%v3%mz)
  98    vv2_v4(its:ite,jts:jte,1:be%v4%mz) = vv % v4(its:ite,jts:jte,1:be%v4%mz)
  99    vv2_v5(its:ite,jts:jte,1:be%v5%mz) = vv % v5(its:ite,jts:jte,1:be%v5%mz)
 100
 101    if (be % ne > 0) then
 102 !     vv2_alpha(its:ite,jts:jte,kts:kte,1:be%ne) = vv % alpha(its:ite,jts:jte,kts:kte,1:be%ne)
 103       vv2_alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:be%ne) = &
 104            vv % alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:be%ne)
 105    end if
 106
 107    call da_vertical_transform(grid, 'u_adj', be, &
 108                                xb % vertical_inner_product, &
 109                                vv, vp)
 110
 111    !----------------------------------------------------------------------
 112    ! [6.0] Calculate RHS of adjoint test equation:
 113    !----------------------------------------------------------------------
 114
 115    adj_par_rhs = 0.0
 116    if (be % v1 % mz > 0) &
 117       adj_par_rhs = sum(vv % v1(its:ite,jts:jte,1:be%v1%mz) * &
 118          vv2_v1(its:ite,jts:jte,1:be%v1%mz)) + adj_par_rhs
 119    if (be % v2 % mz > 0) &
 120       adj_par_rhs = sum(vv % v2(its:ite,jts:jte,1:be%v2%mz) * &
 121          vv2_v2(its:ite,jts:jte,1:be%v2%mz)) + adj_par_rhs
 122    if (be % v3 % mz > 0) &
 123       adj_par_rhs = sum(vv % v3(its:ite,jts:jte,1:be%v3%mz) * &
 124          vv2_v3(its:ite,jts:jte,1:be%v3%mz)) + adj_par_rhs
 125    if (be % v4 % mz > 0) &
 126       adj_par_rhs = sum(vv % v4(its:ite,jts:jte,1:be%v4%mz) * &
 127          vv2_v4(its:ite,jts:jte,1:be%v4%mz)) + adj_par_rhs
 128    if (be % v5 % mz == 1) &
 129       adj_par_rhs = sum(vv % v5(its:ite,jts:jte,1:be%v5%mz) * &
 130          vv2_v5(its:ite,jts:jte,1:be%v5%mz)) + adj_par_rhs
 131 !  if (be % ne > 0) &
 132 !     adj_par_rhs = sum(vv % alpha(its:ite,jts:jte,kts:kte,1:be%ne) * &
 133 !        vv2_alpha(its:ite,jts:jte,kts:kte,1:be%ne)) + adj_par_rhs
 134    if (be % ne > 0) &
 135       adj_par_rhs = sum(vv % alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:be%ne) * &
 136          vv2_alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:be%ne)) + adj_par_rhs
 137
 138    !----------------------------------------------------------------------
 139    ! [7.0] Print output:
 140    !----------------------------------------------------------------------
 141
 142    write(unit=stdout, fmt='(a,1pe22.14)') &
 143         'Single domain < vp,     vp > = ', adj_par_lhs, &
 144         'Single domain < Vv_adj, Vv > = ', adj_par_rhs
 145
 146    adj_sum_lhs = wrf_dm_sum_real(adj_par_lhs)
 147    adj_sum_rhs = wrf_dm_sum_real(adj_par_rhs)
 148
 149
 150    if (rootproc) then
 151       write(unit=stdout, fmt='(/)')
 152       write(unit=stdout, fmt='(a,1pe22.14)') &
 153          'Whole  Domain: < Vp, Vp >     = ', adj_sum_lhs, &
 154          'Whole  Domain: < Vv_adj, Vv > = ', adj_sum_rhs
 155    end if
 156
 157    vv % v1(its:ite,jts:jte,1:be%v1%mz) = vv2_v1(its:ite,jts:jte,1:be%v1%mz)
 158    vv % v2(its:ite,jts:jte,1:be%v2%mz) = vv2_v2(its:ite,jts:jte,1:be%v2%mz)
 159    vv % v3(its:ite,jts:jte,1:be%v3%mz) = vv2_v3(its:ite,jts:jte,1:be%v3%mz)
 160    vv % v4(its:ite,jts:jte,1:be%v4%mz) = vv2_v4(its:ite,jts:jte,1:be%v4%mz)
 161    vv % v5(its:ite,jts:jte,1:be%v5%mz) = vv2_v5(its:ite,jts:jte,1:be%v5%mz)
 162
 163    if (be % ne > 0) then
 164 !     vv % alpha(its:ite,jts:jte,kts:kte,1:be%ne) = vv2_alpha(its:ite,jts:jte,kts:kte,1:be%ne)
 165       vv % alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:be%ne) = &
 166                  vv2_alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:be%ne)
 167    end if
 168
 169    write(unit=stdout, fmt='(/a/)') 'da_check_vvtovp_adjoint: Test Finished.'
 170
 171 1235 continue
 172
 173    if (trace_use) call da_trace_exit("da_check_vvtovp_adjoint")
 174
 175 end subroutine da_check_vvtovp_adjoint
 176
 177