var/da/da_test/da_check_vtox_adjoint_chem.inc

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