var/da/da_test/da_check_vtox_adjoint.inc

   1 subroutine da_check_vtox_adjoint(grid, cv_size, xbx, be, ep, cv1, vv, vp)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Test V 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)             :: cv_size
  14    type (xbx_type),intent(in)      :: xbx    ! For header & non-grid arrays.
  15    type (be_type), intent(in)      :: be     ! background error structure.
  16    type (ep_type), intent(in)      :: ep     ! Ensemble perturbation structure.
  17    real, intent(in)                :: cv1(1:cv_size) ! control variable (local).
  18    type (vp_type), intent(inout)   :: vv     ! Grdipt/EOF CV.
  19    type (vp_type), intent(inout)   :: vp     ! Grdipt/level CV.
  20
  21    real                   :: cv2(1:cv_size)    ! control variable (local).
  22    real                   :: adj_par_lhs ! < x, x >
  23    real                   :: adj_par_rhs ! < x, x >
  24    real                   :: adj_sum_lhs ! < x, x >
  25    real                   :: adj_sum_rhs ! < v_adj, v >
  26
  27    if (trace_use) call da_trace_entry("da_check_vtox_adjoint")
  28
  29    write(unit=stdout, fmt='(/a/)') &
  30       'da_check_vtox_adjoint: Adjoint Test Results:'
  31
  32    !----------------------------------------------------------------------
  33    ! [1.0] Initialise:
  34    !----------------------------------------------------------------------
  35
  36    cv2(:) = 0.0
  37
  38    !----------------------------------------------------------------------
  39    ! [2.0] Perform x = U v transform:
  40    !----------------------------------------------------------------------
  41
  42    call da_zero_x (grid%xa)
  43
  44    call da_transform_vtox(grid,cv_size, xbx, be, ep, cv1, vv, vp)
  45    call da_transform_xtoxa(grid)
  46
  47    !----------------------------------------------------------------------
  48    ! [3.0] Calculate LHS of adjoint test equation:
  49    !----------------------------------------------------------------------
  50
  51    adj_par_lhs = sum(grid%xa % u(its:ite, jts:jte, kts:kte)**2) / typical_u_rms**2 &
  52                + sum(grid%xa % v(its:ite, jts:jte, kts:kte)**2) / typical_v_rms**2 &
  53                + sum(grid%xa % p(its:ite, jts:jte, kts:kte)**2) / typical_p_rms**2 &
  54                + sum(grid%xa % t(its:ite, jts:jte, kts:kte)**2) / typical_t_rms**2 &
  55                + sum(grid%xa % q(its:ite, jts:jte, kts:kte)**2) / typical_q_rms**2 &
  56                + sum(grid%xa % rho(its:ite,jts:jte,kts:kte)**2)/ typical_rho_rms**2 &
  57                + sum(grid%xa % psfc(its:ite, jts:jte)**2) / typical_p_rms**2
  58
  59    if (cv_options_hum == cv_options_hum_relative_humidity) then
  60       adj_par_lhs = adj_par_lhs &
  61               + sum(grid%xa % rh(its:ite, jts:jte, kts:kte)**2) / typical_rh_rms**2
  62    end if
  63 !
  64    if (use_radar_rf .or. use_radar_rhv .or.  crtm_cloud ) then
  65       adj_par_lhs = adj_par_lhs &
  66              + sum(grid%xa % qcw(its:ite, jts:jte, kts:kte)**2)/typical_qcw_rms**2 &
  67              + sum(grid%xa % qrn(its:ite, jts:jte, kts:kte)**2)/typical_qrn_rms**2
  68       if ( cloud_cv_options /= 1 ) then
  69          adj_par_lhs = adj_par_lhs &
  70              + sum(grid%xa % qci(its:ite, jts:jte, kts:kte)**2)/typical_qci_rms**2 &
  71              + sum(grid%xa % qsn(its:ite, jts:jte, kts:kte)**2)/typical_qsn_rms**2 &
  72              + sum(grid%xa % qgr(its:ite, jts:jte, kts:kte)**2)/typical_qgr_rms**2
  73       end if
  74    end if
  75
  76    if (use_radarobs) then
  77       adj_par_lhs = adj_par_lhs &
  78          + sum(grid%xa % wh (its:ite, jts:jte, kts:kte)**2)/typical_w_rms**2
  79    else
  80       adj_par_lhs = adj_par_lhs &
  81          + sum(grid%xa % w  (its:ite, jts:jte, kts:kte)**2)/typical_w_rms**2
  82    end if
  83
  84    !-------------------------------------------------------------------------
  85    ! [4.0] Rescale input to adjoint routine:
  86    !-------------------------------------------------------------------------
  87
  88    grid%xa % u(:,:,:) = grid%xa % u(:,:,:) / typical_u_rms**2
  89    grid%xa % v(:,:,:) = grid%xa % v(:,:,:) / typical_v_rms**2
  90    grid%xa % p(:,:,:) = grid%xa % p(:,:,:) / typical_p_rms**2
  91    grid%xa % t(:,:,:) = grid%xa % t(:,:,:) / typical_t_rms**2
  92    grid%xa % q(:,:,:) = grid%xa % q(:,:,:) / typical_q_rms**2
  93    grid%xa % rho(:,:,:) = grid%xa % rho(:,:,:) / typical_rho_rms**2
  94
  95    grid%xa % psfc(:,:) = grid%xa % psfc(:,:) / typical_p_rms**2
  96
  97    if (cv_options_hum == cv_options_hum_relative_humidity) then
  98       grid%xa % rh(:,:,:) = grid%xa % rh(:,:,:) / typical_rh_rms**2
  99    end if
 100
 101 !
 102    if (use_radar_rf .or. use_radar_rhv .or. crtm_cloud ) then
 103       grid%xa % qcw(:,:,:) = grid%xa % qcw(:,:,:) / typical_qcw_rms**2
 104       grid%xa % qrn(:,:,:) = grid%xa % qrn(:,:,:) / typical_qrn_rms**2
 105       if ( cloud_cv_options /= 1 ) then
 106          grid%xa % qci(:,:,:) = grid%xa % qci(:,:,:) / typical_qci_rms**2
 107          grid%xa % qsn(:,:,:) = grid%xa % qsn(:,:,:) / typical_qsn_rms**2
 108          grid%xa % qgr(:,:,:) = grid%xa % qgr(:,:,:) / typical_qgr_rms**2
 109       end if
 110    end if
 111
 112    if (use_radarobs) then
 113       grid%xa %wh(:,:,:) = grid%xa %wh(:,:,:) / typical_w_rms**2
 114       grid%xa % w(:,:,:) = 0.0
 115    else
 116       grid%xa %w (:,:,:) = grid%xa %w (:,:,:) / typical_w_rms**2
 117    end if
 118
 119    !-------------------------------------------------------------------------
 120    ! [5.0] Perform adjoint operation:
 121    !-------------------------------------------------------------------------
 122
 123    call da_transform_xtoxa_adj(grid)
 124    !as of V3.9, da_zero_vp_type(vp) is commented out in the da_transform_vtox_adj
 125    call da_zero_vp_type(vp) !as of V3.9, zero out vp here
 126    call da_transform_vtox_adj(grid, cv_size, xbx, be, ep, vp, vv, cv2)
 127
 128    !-------------------------------------------------------------------------
 129    ! [6.0] Calculate RHS of adjoint test equation:
 130    !-------------------------------------------------------------------------
 131
 132    adj_par_rhs = sum(cv1(1:cv_size) * cv2(1:cv_size))
 133
 134    !-------------------------------------------------------------------------
 135    ! [7.0] Print output:
 136    !-------------------------------------------------------------------------
 137
 138    if (.not. global ) then
 139     if( num_procs == 1) then
 140       write(unit=stdout, fmt='(/)')
 141       write(unit=stdout, fmt='(a,1pe22.14)') &
 142          'Single Domain: < x, x >     = ', adj_par_lhs, &
 143          'Single Domain: < v_adj, v > = ', adj_par_rhs
 144     else
 145       write(unit=stdout, fmt='(/a/,a/)')&
 146         'It is Multi Processor Run: ',&
 147         'For Single Domain: da_check_vtox_adjoint Test: Not Performed'
 148     endif
 149    end if
 150
 151    adj_sum_lhs = wrf_dm_sum_real(adj_par_lhs)
 152
 153    if (global) then
 154       adj_sum_rhs = adj_par_rhs
 155    else
 156       adj_sum_rhs = wrf_dm_sum_real(adj_par_rhs)
 157    end if
 158
 159    if (rootproc) then
 160       write(unit=stdout, fmt='(/)')
 161       write(unit=stdout, fmt='(a,1pe22.14)') &
 162          'Whole  Domain: < x, x >     = ', adj_sum_lhs, &
 163          'Whole  Domain: < v_adj, v > = ', adj_sum_rhs
 164    end if
 165
 166    write(unit=stdout, fmt='(/a/)') 'da_check_vtox_adjoint: Finished'
 167
 168    if (trace_use) call da_trace_exit("da_check_vtox_adjoint")
 169
 170 end subroutine da_check_vtox_adjoint
 171
 172