var/da/da_test/da_check_vptox_adjoint.inc

   1 subroutine da_check_vptox_adjoint(grid, ne, be, ep, vp, cv_size)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Test Vp to X routine and adjoint for compatibility.
   5    !
   6    ! Method:  Standard adjoint test: < x, x > = < v_adj, v >.
   7    !---------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    type (domain), intent(inout)     :: grid
  12
  13    integer, intent(in)              :: ne   ! Ensemble size.
  14    type (be_type), intent(in)       :: be   ! background errors.
  15    type (ep_type), intent(in)       :: ep   ! Ensemble perturbation type.
  16    type (vp_type),intent(inout)     :: vp   ! grdipt/level cv (local)
  17    integer, intent(in)              :: cv_size
  18
  19 ! control variable
  20    real                             :: cv(1:cv_size), cv_2(1:cv_size)
  21
  22    real                             :: adj_par_lhs ! < x, x >
  23    real                             :: adj_par_rhs ! < v_adj, v >
  24    real                             :: adj_sum_lhs ! < x, x >
  25    real                             :: adj_sum_rhs ! < v_adj, v >
  26    real                             :: vp2_v1(ims:ime,jms:jme,kms:kme)
  27    real                             :: vp2_v2(ims:ime,jms:jme,kms:kme)
  28    real                             :: vp2_v3(ims:ime,jms:jme,kms:kme)
  29    real                             :: vp2_v4(ims:ime,jms:jme,kms:kme)
  30    real                             :: vp2_v5(ims:ime,jms:jme,kms:kme)
  31
  32    real                             :: vp2_v6(ims:ime,jms:jme,kms:kme)
  33    real                             :: vp2_v7(ims:ime,jms:jme,kms:kme)
  34    real                             :: vp2_v8(ims:ime,jms:jme,kms:kme)
  35    real                             :: vp2_v9(ims:ime,jms:jme,kms:kme)
  36    real                             :: vp2_v10(ims:ime,jms:jme,kms:kme)
  37    real                             :: vp2_v11(ims:ime,jms:jme,kms:kme)
  38
  39 !  real                             :: vp2_alpha(ims:ime,jms:jme,kms:kme,1:ne)
  40    real                             :: vp2_alpha(ims_int:ime_int,jms_int:jme_int,kms_int:kme_int,1:ne)
  41
  42    if (trace_use) call da_trace_entry("da_check_vptox_adjoint")
  43
  44    !-------------------------------------------------------------------------
  45    ! [1.0] Initialise:
  46    !-------------------------------------------------------------------------
  47
  48
  49    call da_zero_x(grid%xa)
  50
  51    vp2_v1(:,:,:) = vp % v1(:,:,:)
  52    vp2_v2(:,:,:) = vp % v2(:,:,:)
  53
  54    call da_psichi_to_uv(vp % v1, vp % v2, grid%xb % coefx, &
  55                         grid%xb % coefy , grid%xa % u, grid%xa % v)
  56
  57    adj_par_lhs = sum(grid%xa % u(its:ite,jts:jte,:)**2) / typical_u_rms**2
  58    adj_par_lhs = sum(grid%xa % v(its:ite,jts:jte,:)**2) / typical_v_rms**2 + &
  59       adj_par_lhs
  60
  61    grid%xa % u(:,:,:) = grid%xa % u(:,:,:) / typical_u_rms**2
  62    grid%xa % v(:,:,:) = grid%xa % v(:,:,:) / typical_v_rms**2
  63
  64    vp%v1(:,:,:)=0.0
  65    vp%v2(:,:,:)=0.0
  66
  67    call da_psichi_to_uv_adj(grid%xa % u, grid%xa % v, grid%xb % coefx,   &
  68                              grid%xb % coefy, vp % v1, vp % v2)
  69
  70    adj_par_rhs = sum(vp % v1(its:ite,jts:jte,:) * vp2_v1(its:ite,jts:jte,:))
  71    adj_par_rhs = sum(vp % v2(its:ite,jts:jte,:) * vp2_v2(its:ite,jts:jte,:)) + &
  72       adj_par_rhs
  73
  74       write(unit=stdout, fmt='(/a/)') &
  75           'da_check_da_psichi_to_uv_adjoint: Test Results:'
  76       write(unit=stdout, fmt='(/)')
  77       write(unit=stdout, fmt='(a,1pe22.14)') &
  78           'Single Domain: < u_v,     u_v         > = ', adj_par_lhs, &
  79           'Single Domain: < psi_chi, psi_chi_adj > = ', adj_par_rhs
  80
  81    adj_sum_lhs = wrf_dm_sum_real(adj_par_lhs)
  82    adj_sum_rhs = wrf_dm_sum_real(adj_par_rhs)
  83
  84    if (rootproc) then
  85
  86       write(unit=stdout, fmt='(/)')
  87       write(unit=stdout, fmt='(a,1pe22.14)') &
  88           'Whole  Domain: < u_v,     u_v         > = ', adj_sum_lhs, &
  89           'Whole  Domain: < psi_chi, psi_chi_adj > = ', adj_sum_rhs
  90    end if
  91       write(unit=stdout, fmt='(/a/)') &
  92           'da_check_da_psichi_to_uv_adjoint: Test Finished:'
  93
  94    vp%v1(:,:,:) = vp2_v1(:,:,:)
  95    vp%v2(:,:,:) = vp2_v2(:,:,:)
  96
  97    call da_zero_x(grid%xa)
  98
  99    vp2_v1(:,:,:) = vp % v1(:,:,:)
 100    vp2_v2(:,:,:) = vp % v2(:,:,:)
 101    vp2_v3(:,:,:) = vp % v3(:,:,:)
 102    vp2_v4(:,:,:) = vp % v4(:,:,:)
 103    vp2_v5(:,:,:) = vp % v5(:,:,:)
 104
 105    if ( cloud_cv_options >= 2 ) then
 106       vp2_v6(:,:,:) = vp % v6(:,:,:)
 107       vp2_v7(:,:,:) = vp % v7(:,:,:)
 108       vp2_v8(:,:,:) = vp % v8(:,:,:)
 109       vp2_v9(:,:,:) = vp % v9(:,:,:)
 110       vp2_v10(:,:,:)= vp % v10(:,:,:)
 111    end if
 112
 113    if ( use_cv_w ) then
 114       vp2_v11(:,:,:)= vp % v11(:,:,:)
 115    end if
 116
 117    !hcl if (be % ne > 0) vp2_alpha(:,:,:,:) = vp % alpha(:,:,:,:)
 118
 119    !-------------------------------------------------------------------------
 120    ! [2.0] Perform x = U vp transform:
 121    !-------------------------------------------------------------------------
 122
 123    if ( cv_options == 3 ) then
 124
 125       call random_number(cv(:))
 126       cv(:) = cv(:) - 0.5
 127       cv_2 = cv
 128
 129       call da_apply_be( be, cv, vp, grid)
 130       call da_transform_bal( vp, be, grid)
 131
 132    else
 133
 134       call da_transform_vptox(grid, vp, be, ep)
 135
 136    end if
 137
 138    !-------------------------------------------------------------------------
 139    ! [3.0] Calculate LHS of adjoint test equation:
 140    !-------------------------------------------------------------------------
 141
 142    !  grid%xa % u(:,:,:) = 0.0
 143    !  grid%xa % v(:,:,:) = 0.0
 144    !  grid%xa % t(:,:,:) = 0.0
 145    !  grid%xa % q(:,:,:) = 0.0
 146    !  grid%xa%psfc(:,:) = 0.0
 147
 148    !  grid%xa % p(:,:,:) = 0.0
 149    !  grid%xa % rho(:,:,:) = 0.0
 150    !  grid%xa % w(:,:,:) = 0.0
 151    !  grid%xa % wh(:,:,:) = 0.0
 152    !  grid%xa % rh(:,:,:) = 0.0
 153    !  grid%xa % qt(:,:,:) = 0.0
 154    !  grid%xa % qcw(:,:,:) = 0.0
 155    !  grid%xa % qrn(:,:,:) = 0.0
 156
 157    adj_par_lhs = sum(grid%xa%u(its:ite,jts:jte,:)**2)/typical_u_rms**2
 158    adj_par_lhs = sum(grid%xa%v(its:ite,jts:jte,:)**2)/typical_v_rms**2 + adj_par_lhs
 159    adj_par_lhs = sum(grid%xa%t(its:ite,jts:jte,:)**2)/typical_t_rms**2 + adj_par_lhs
 160    if ( (use_radar_rf .or. crtm_cloud) .and. (cloud_cv_options == 1) ) then
 161       adj_par_lhs = sum(grid%xa%qt(its:ite,jts:jte,:)**2)/typical_q_rms**2 + adj_par_lhs
 162    else
 163       adj_par_lhs = sum(grid%xa%q(its:ite,jts:jte,:)**2)/typical_q_rms**2 + adj_par_lhs
 164    end if
 165    adj_par_lhs = sum(grid%xa%psfc(its:ite,jts:jte)**2)/typical_p_rms**2 + adj_par_lhs
 166
 167    adj_par_lhs = sum(grid%xa%p(its:ite,jts:jte,:)**2)/typical_p_rms**2 + adj_par_lhs
 168    adj_par_lhs = sum(grid%xa%rho(its:ite,jts:jte,:)**2)/typical_rho_rms**2 + &
 169       adj_par_lhs
 170
 171    if (use_radarobs) then
 172       adj_par_lhs = adj_par_lhs &
 173                 + sum(grid%xa % wh (its:ite, jts:jte, kts:kte)**2)/typical_w_rms**2
 174    else
 175       adj_par_lhs = adj_par_lhs &
 176                 + sum(grid%xa % w  (its:ite, jts:jte, kts:kte)**2)/typical_w_rms**2
 177    end if
 178
 179    if (cv_options_hum == cv_options_hum_relative_humidity) then
 180       adj_par_lhs = sum(grid%xa % rh(its:ite,jts:jte,:)**2) / &
 181          typical_rh_rms**2 + adj_par_lhs
 182    end if
 183
 184    if ( cloud_cv_options >= 2 ) then
 185       adj_par_lhs = sum(grid%xa % qcw(its:ite,jts:jte,kts:kte)**2) / &
 186            typical_qcw_rms**2 + adj_par_lhs
 187       adj_par_lhs = sum(grid%xa % qrn(its:ite,jts:jte,kts:kte)**2) / &
 188            typical_qrn_rms**2 + adj_par_lhs
 189       adj_par_lhs = sum(grid%xa % qci (its:ite,jts:jte,kts:kte)**2) / &
 190            typical_qci_rms**2 + adj_par_lhs
 191       adj_par_lhs = sum(grid%xa % qsn(its:ite,jts:jte,kts:kte)**2) / &
 192            typical_qsn_rms**2 + adj_par_lhs
 193       adj_par_lhs = sum(grid%xa % qgr (its:ite,jts:jte,kts:kte)**2) / &
 194            typical_qgr_rms**2 + adj_par_lhs
 195    end if
 196
 197    if (use_radar_rf .or. crtm_cloud) then
 198       if ( cloud_cv_options == 1 ) then
 199          adj_par_lhs = sum(grid%xa % qcw(its:ite,jts:jte,kts:kte)**2) / &
 200             typical_qcw_rms**2 + adj_par_lhs
 201          adj_par_lhs = sum(grid%xa % qrn(its:ite,jts:jte,kts:kte)**2) / &
 202             typical_qrn_rms**2 + adj_par_lhs
 203       end if
 204    end if
 205
 206    !-------------------------------------------------------------------------
 207    ! [4.0] Rescale input to adjoint routine:
 208    !-------------------------------------------------------------------------
 209
 210    grid%xa % u(:,:,:) = grid%xa % u(:,:,:) / typical_u_rms**2
 211    grid%xa % v(:,:,:) = grid%xa % v(:,:,:) / typical_v_rms**2
 212    grid%xa % t(:,:,:) = grid%xa % t(:,:,:) / typical_t_rms**2
 213    if ( (use_radar_rf .or. crtm_cloud) .and. (cloud_cv_options == 1) ) then
 214       grid%xa % qt(:,:,:) = grid%xa % qt(:,:,:) / typical_q_rms**2
 215    else
 216       grid%xa % q(:,:,:) = grid%xa % q(:,:,:) / typical_q_rms**2
 217    end if
 218    grid%xa%psfc(:,:) = grid%xa%psfc(:,:) / typical_p_rms**2
 219
 220    grid%xa % p(:,:,:) = grid%xa % p(:,:,:) / typical_p_rms**2
 221    grid%xa % rho(:,:,:) = grid%xa % rho(:,:,:) / typical_rho_rms**2
 222
 223    if (use_radarobs) then
 224       grid%xa %wh(:,:,:) = grid%xa %wh(:,:,:) / typical_w_rms**2
 225       grid%xa % w(:,:,:) = 0.0
 226    else
 227       grid%xa %w (:,:,:) = grid%xa %w (:,:,:) / typical_w_rms**2
 228    end if
 229
 230    if (cv_options_hum == cv_options_hum_relative_humidity) then
 231       grid%xa % rh(:,:,:) = grid%xa % rh(:,:,:) / typical_rh_rms**2
 232    end if
 233
 234    if ( cloud_cv_options >= 2 ) then
 235       grid%xa % qcw(:,:,:) = grid%xa % qcw(:,:,:) / typical_qcw_rms**2
 236       grid%xa % qrn(:,:,:) = grid%xa % qrn(:,:,:) / typical_qrn_rms**2
 237       grid%xa % qci(:,:,:) = grid%xa % qci(:,:,:) / typical_qci_rms**2
 238       grid%xa % qsn(:,:,:) = grid%xa % qsn(:,:,:) / typical_qsn_rms**2
 239       grid%xa % qgr(:,:,:) = grid%xa % qgr(:,:,:) / typical_qgr_rms**2
 240    end if
 241
 242    if (use_radar_rf .or. crtm_cloud) then
 243       if ( cloud_cv_options == 1 ) then
 244          grid%xa % qcw(:,:,:) = grid%xa % qcw(:,:,:) / typical_qcw_rms**2
 245          grid%xa % qrn(:,:,:) = grid%xa % qrn(:,:,:) / typical_qrn_rms**2
 246       end if
 247    end if
 248
 249    !-------------------------------------------------------------------------
 250    ! [5.0] Perform adjoint operation:
 251    !-------------------------------------------------------------------------
 252
 253    call da_zero_vp_type (vp)
 254
 255    if (cv_options == 3 ) then
 256
 257       cv = 0.0
 258
 259       call da_transform_bal_adj( vp, be, grid)
 260       call da_apply_be_adj( be, cv, vp, grid)
 261
 262    else
 263
 264        call da_transform_vptox_adj(grid, vp, be, ep)
 265
 266    end if
 267
 268    !-------------------------------------------------------------------------
 269    ! [6.0] Calculate RHS of adjoint test equation:
 270    !-------------------------------------------------------------------------
 271
 272    adj_par_rhs = sum(vp % v1(its:ite,jts:jte,:) * vp2_v1(its:ite,jts:jte,:))
 273    adj_par_rhs = sum(vp % v2(its:ite,jts:jte,:) * vp2_v2(its:ite,jts:jte,:)) + &
 274       adj_par_rhs
 275    adj_par_rhs = sum(vp % v3(its:ite,jts:jte,:) * vp2_v3(its:ite,jts:jte,:)) + &
 276       adj_par_rhs
 277    adj_par_rhs = sum(vp % v4(its:ite,jts:jte,:) * vp2_v4(its:ite,jts:jte,:)) + &
 278       adj_par_rhs
 279    adj_par_rhs = sum(vp % v5(its:ite,jts:jte,:) * vp2_v5(its:ite,jts:jte,:)) + &
 280       adj_par_rhs
 281
 282    if ( cloud_cv_options >= 2 ) then
 283       adj_par_rhs = sum(vp % v6(its:ite,jts:jte,:) * vp2_v6(its:ite,jts:jte,:)) + &
 284          adj_par_rhs
 285       adj_par_rhs = sum(vp % v7(its:ite,jts:jte,:) * vp2_v7(its:ite,jts:jte,:)) + &
 286          adj_par_rhs
 287       adj_par_rhs = sum(vp % v8(its:ite,jts:jte,:) * vp2_v8(its:ite,jts:jte,:)) + &
 288          adj_par_rhs
 289       adj_par_rhs = sum(vp % v9(its:ite,jts:jte,:) * vp2_v9(its:ite,jts:jte,:)) + &
 290          adj_par_rhs
 291       adj_par_rhs = sum(vp % v10(its:ite,jts:jte,:)* vp2_v10(its:ite,jts:jte,:)) + &
 292          adj_par_rhs
 293    end if
 294
 295    if ( use_cv_w ) then
 296       adj_par_rhs = sum(vp % v11(its:ite,jts:jte,:)* vp2_v11(its:ite,jts:jte,:)) + &
 297          adj_par_rhs
 298    end if
 299
 300    !hcl if (be % ne > 0) then
 301    !hcl    adj_par_rhs = sum(vp % alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,:) * &
 302    !hcl       vp2_alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,:)) + adj_par_rhs
 303    !hcl end if
 304
 305    if ( cv_options == 3 ) adj_par_rhs = sum (cv_2*cv)
 306
 307    !-------------------------------------------------------------------------
 308    ! [7.0] Print output:
 309    !-------------------------------------------------------------------------
 310
 311    adj_sum_lhs = wrf_dm_sum_real(adj_par_lhs)
 312    adj_sum_rhs = wrf_dm_sum_real(adj_par_rhs)
 313
 314    write(unit=stdout, fmt='(/a/)') 'da_check_vptox_adjoint: Test Results:'
 315       write(unit=stdout, fmt='(/)')
 316       write(unit=stdout, fmt='(a,1pe22.14)') &
 317          'Single Domain: < x, x >       = ', adj_par_lhs, &
 318          'Single Domain: < vp_adj, vp > = ', adj_par_rhs
 319
 320    if (rootproc) then
 321       write(unit=stdout, fmt='(/)')
 322       write(unit=stdout, fmt='(a,1pe22.14)') &
 323          'Whole  Domain: < x, x >       = ', adj_sum_lhs, &
 324          'Whole  Domain: < vp_adj, vp > = ', adj_sum_rhs
 325    end if
 326
 327    vp % v1(:,:,:) = vp2_v1(:,:,:)
 328    vp % v2(:,:,:) = vp2_v2(:,:,:)
 329    vp % v3(:,:,:) = vp2_v3(:,:,:)
 330    vp % v4(:,:,:) = vp2_v4(:,:,:)
 331    vp % v5(:,:,:) = vp2_v5(:,:,:)
 332
 333    if ( cloud_cv_options >= 2 ) then
 334       vp % v6(:,:,:) = vp2_v6(:,:,:)
 335       vp % v7(:,:,:) = vp2_v7(:,:,:)
 336       vp % v8(:,:,:) = vp2_v8(:,:,:)
 337       vp % v9(:,:,:) = vp2_v9(:,:,:)
 338       vp % v10(:,:,:)= vp2_v10(:,:,:)
 339    end if
 340
 341    if ( use_cv_w ) then
 342       vp % v11(:,:,:)= vp2_v11(:,:,:)
 343    end if
 344
 345    !hcl if (be % ne > 0) vp % alpha(:,:,:,:) = vp2_alpha(:,:,:,:)
 346
 347    write(unit=stdout, fmt='(/a/)') 'da_check_vptox_adjoint: Test Finished:'
 348
 349    if (trace_use) call da_trace_exit("da_check_vptox_adjoint")
 350
 351 end subroutine da_check_vptox_adjoint
 352
 353