var/da/da_test/da_check_dynamics_adjoint.inc

   1 subroutine da_check_dynamics_adjoint(cv_size, cv, xbx, be, grid, config_flags)
   2
   3    !--------------------------------------------------------------------------
   4    ! Purpose: Test observation operator transform and adjoint for compatibility.
   5    !
   6    ! Method:  Standard adjoint test: < y, y > = < x, x_adj >.
   7    !    Updated for Analysis on Arakawa-C grid
   8    !    Author: Xiaowen Tang,  MMM/ESSL/NCAR,  Date: 10/22/2008
   9    !---------------------------------------------------------------------------
  10
  11    implicit none
  12
  13    integer, intent(in)                       :: cv_size ! Size of cv array.
  14    type (be_type),             intent(in)    :: be    ! background error structure.
  15    real, intent(inout)                       :: cv(1:cv_size)   ! control variables.
  16    type (xbx_type),            intent(inout) :: xbx   ! Header & non-gridded vars.
  17    type (domain),              intent(inout) :: grid
  18    type(grid_config_rec_type), intent(inout) :: config_flags
  19
  20    real                           :: adj_ttl_lhs   ! < y, y >
  21    real                           :: adj_ttl_rhs   ! < x, x_adj >
  22
  23    real                           :: partial_lhs   ! < y, y >
  24    real                           :: partial_rhs   ! < x, x_adj >
  25
  26    real                           :: pertile_lhs   ! < y, y >
  27    real                           :: pertile_rhs   ! < x, x_adj >
  28
  29    real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_u, xa2_v, xa2_t, &
  30                                                  xa2_p, xa2_q, xa2_rh
  31    real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_w
  32    real, dimension(ims:ime, jms:jme)          :: xa2_psfc
  33    real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_qcw, xa2_qci, xa2_qrn, xa2_qsn, xa2_qgr
  34    real, dimension(ims:ime, jms:jme, kms:kme) :: x6a2_u, x6a2_v, x6a2_t, &
  35                                                  x6a2_p, x6a2_q, x6a2_rh
  36    real, dimension(ims:ime, jms:jme, kms:kme) :: x6a2_w
  37    real, dimension(ims:ime, jms:jme)          :: x6a2_psfc
  38    real, dimension(ims:ime, jms:jme, kms:kme) :: x6a2_qcw, x6a2_qci, x6a2_qrn, x6a2_qsn, x6a2_qgr
  39    real, dimension(ims:ime, jms:jme, kms:kme) :: tempv
  40
  41    integer :: nobwin, i, j, k, fgat_rain
  42    character(len=4) :: filnam
  43    character(len=256) :: timestr
  44    integer :: time_step_seconds
  45    type(x_type) :: shuffle
  46    real             :: subarea, whole_area
  47
  48    if (trace_use) call da_trace_entry("da_check_dynamics_adjoint")
  49
  50    write (unit=stdout, fmt='(/a/)') 'da_check_dynamics_adjoint: Test Results:'
  51
  52    !----------------------------------------------------------------------
  53    ! [1.0] Initialise:
  54    !----------------------------------------------------------------------
  55
  56    partial_lhs = 0.0
  57    pertile_lhs = 0.0
  58
  59 #ifdef A2C
  60   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  61        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
  62      ipe = ipe + 1
  63      ide = ide + 1
  64   end if
  65
  66   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  67        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
  68      jpe = jpe + 1
  69      jde = jde + 1
  70   end if
  71 #endif
  72
  73 #ifdef DM_PARALLEL
  74 #include "HALO_XA.inc"
  75 #endif
  76
  77 #ifdef A2C
  78   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  79        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
  80      ipe = ipe - 1
  81      ide = ide - 1
  82   end if
  83
  84   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  85        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
  86      jpe = jpe - 1
  87      jde = jde - 1
  88   end if
  89 #endif
  90
  91    call da_transform_vtox(grid, cv_size, xbx, be, grid%ep, cv, grid%vv, grid%vp)
  92    call da_transform_xtoxa(grid)
  93
  94    xa2_u(ims:ime, jms:jme, kms:kme) = grid%xa%u(ims:ime, jms:jme, kms:kme)
  95    xa2_v(ims:ime, jms:jme, kms:kme) = grid%xa%v(ims:ime, jms:jme, kms:kme)
  96    xa2_t(ims:ime, jms:jme, kms:kme) = grid%xa%t(ims:ime, jms:jme, kms:kme)
  97    xa2_p(ims:ime, jms:jme, kms:kme) = grid%xa%p(ims:ime, jms:jme, kms:kme)
  98    xa2_q(ims:ime, jms:jme, kms:kme) = grid%xa%q(ims:ime, jms:jme, kms:kme)
  99    xa2_w(ims:ime, jms:jme, kms:kme) = grid%xa%w(ims:ime, jms:jme, kms:kme)
 100    xa2_rh(ims:ime, jms:jme, kms:kme)= grid%xa%rh(ims:ime, jms:jme, kms:kme)
 101    xa2_psfc(ims:ime, jms:jme)       = grid%xa%psfc(ims:ime, jms:jme)
 102
 103    xa2_qcw = 0.0
 104    xa2_qrn = 0.0
 105    xa2_qci = 0.0
 106    xa2_qsn = 0.0
 107    xa2_qgr = 0.0
 108    if ( cloud_cv_options >= 1 ) then
 109       xa2_qcw(ims:ime, jms:jme, kms:kme) = grid%xa%qcw(ims:ime, jms:jme, kms:kme)
 110       xa2_qrn(ims:ime, jms:jme, kms:kme) = grid%xa%qrn(ims:ime, jms:jme, kms:kme)
 111       if ( cloud_cv_options >= 2 ) then
 112          xa2_qci(ims:ime, jms:jme, kms:kme) = grid%xa%qci(ims:ime, jms:jme, kms:kme)
 113          xa2_qsn(ims:ime, jms:jme, kms:kme) = grid%xa%qsn(ims:ime, jms:jme, kms:kme)
 114          xa2_qgr(ims:ime, jms:jme, kms:kme) = grid%xa%qgr(ims:ime, jms:jme, kms:kme)
 115       end if
 116    end if
 117
 118    x6a2_u = 0.0
 119    x6a2_v = 0.0
 120    x6a2_t = 0.0
 121    x6a2_p = 0.0
 122    x6a2_q = 0.0
 123    x6a2_w = 0.0
 124    x6a2_rh = 0.0
 125    x6a2_psfc = 0.0
 126
 127    x6a2_qcw = 0.0
 128    x6a2_qci = 0.0
 129    x6a2_qrn = 0.0
 130    x6a2_qsn = 0.0
 131    x6a2_qgr = 0.0
 132
 133 #ifdef A2C
 134     if( ite == ide ) &
 135 print*,__FILE__,jte,' xa2_u.xa2_u for col= ',ite+1,sum(xa2_u(ite+1, jts:jte, kts:kte) * xa2_u(ite+1, jts:jte, kts:kte))
 136     if( jte == jde ) &
 137 print*,__FILE__,jte,' xa2_v.xa2_v for row= ',jte+1,sum(xa2_v(its:ite, jte+1, kts:kte) * xa2_v(its:ite, jte+1, kts:kte))
 138 #endif
 139
 140    if ( num_fgat_time > 1 ) then
 141       call domain_clock_get (grid, stop_timestr=timestr)
 142       call domain_clock_set( grid, current_timestr=timestr )
 143       call domain_clock_set (grid, time_step_seconds=-1*var4d_bin)
 144       call domain_clockprint(150, grid, 'get CurrTime from clock,')
 145    endif
 146
 147    fgat_rain = num_fgat_time
 148    do nobwin= num_fgat_time, 1, -1
 149       !----------------------------------------------------------------------
 150       ! [1.0] Perform y = Hx transform:
 151       !----------------------------------------------------------------------
 152       ! call da_uv_to_divergence(grid%xb, grid%xa%u, grid%xa%v, tempv)
 153       call da_divergence_constraint(grid%xb, grid%xa%u, grid%xa%v, tempv)
 154       partial_lhs = partial_lhs + SUM( tempv(ims:ime, jms:jme, kms:kme) * tempv(ims:ime, jms:jme, kms:kme) )
 155       pertile_lhs = pertile_lhs + SUM( tempv(its:ite, jts:jte, kts:kte) * tempv(its:ite, jts:jte, kts:kte) )
 156
 157       !----------------------------------------------------------------------
 158       ! [5.0] Perform adjoint operation:
 159       !----------------------------------------------------------------------
 160       call da_zero_x (grid%xa)
 161       ! call da_uv_to_divergence_adj(grid, grid%xa%u, grid%xa%v, tempv)
 162       call da_divergence_constraint_adj(grid, grid%xa%u, grid%xa%v, tempv)
 163
 164       if ( nobwin > 1 ) call domain_clockadvance (grid)
 165       call domain_clockprint(150, grid, 'DEBUG Adjoint Check:  get CurrTime from clock,')
 166
 167    end do
 168
 169    if ( num_fgat_time > 1 ) then
 170       call nl_get_time_step ( grid%id, time_step_seconds)
 171       call domain_clock_set (grid, time_step_seconds=time_step_seconds)
 172       call domain_clockprint(150, grid, 'get CurrTime from clock,')
 173    endif
 174
 175
 176    pertile_rhs = sum (grid%xa%u(ims:ime, jms:jme, kms:kme) * xa2_u(ims:ime, jms:jme, kms:kme)) &
 177       + sum (grid%xa%v(ims:ime, jms:jme, kms:kme) * xa2_v(ims:ime, jms:jme, kms:kme))          &
 178       + sum (grid%xa%w(ims:ime, jms:jme, kms:kme) * xa2_w(ims:ime, jms:jme, kms:kme))          &
 179       + sum (grid%xa%t(ims:ime, jms:jme, kms:kme) * xa2_t(ims:ime, jms:jme, kms:kme))          &
 180       + sum (grid%xa%p(ims:ime, jms:jme, kms:kme) * xa2_p(ims:ime, jms:jme, kms:kme))          &
 181       + sum (grid%xa%q(ims:ime, jms:jme, kms:kme) * xa2_q(ims:ime, jms:jme, kms:kme))          &
 182       + sum (grid%xa%rh(ims:ime, jms:jme, kms:kme)* xa2_rh(ims:ime, jms:jme, kms:kme))         &
 183       + sum (grid%xa%psfc(ims:ime, jms:jme) * xa2_psfc(ims:ime, jms:jme))
 184 #ifdef VAR4D
 185    pertile_rhs = pertile_rhs &
 186       + sum (grid%x6a%u(ims:ime, jms:jme, kms:kme) * x6a2_u(ims:ime, jms:jme, kms:kme))        &
 187       + sum (grid%x6a%v(ims:ime, jms:jme, kms:kme) * x6a2_v(ims:ime, jms:jme, kms:kme))        &
 188       + sum (grid%x6a%w(ims:ime, jms:jme, kms:kme) * x6a2_w(ims:ime, jms:jme, kms:kme))        &
 189       + sum (grid%x6a%t(ims:ime, jms:jme, kms:kme) * x6a2_t(ims:ime, jms:jme, kms:kme))        &
 190       + sum (grid%x6a%p(ims:ime, jms:jme, kms:kme) * x6a2_p(ims:ime, jms:jme, kms:kme))        &
 191       + sum (grid%x6a%q(ims:ime, jms:jme, kms:kme) * x6a2_q(ims:ime, jms:jme, kms:kme))        &
 192       + sum (grid%x6a%rh(ims:ime, jms:jme, kms:kme)* x6a2_rh(ims:ime, jms:jme, kms:kme))       &
 193       + sum (grid%x6a%psfc(ims:ime, jms:jme) * x6a2_psfc(ims:ime, jms:jme))
 194 #endif
 195    if ( cloud_cv_options >= 1 ) then
 196       pertile_rhs = pertile_rhs &
 197          + sum (grid%xa%qcw(ims:ime, jms:jme, kms:kme) * xa2_qcw(ims:ime, jms:jme, kms:kme))   &
 198          + sum (grid%xa%qrn(ims:ime, jms:jme, kms:kme) * xa2_qrn(ims:ime, jms:jme, kms:kme))
 199       if ( cloud_cv_options >= 2 ) then
 200          pertile_rhs = pertile_rhs &
 201             + sum (grid%xa%qci(ims:ime, jms:jme, kms:kme) * xa2_qci(ims:ime, jms:jme, kms:kme))&
 202             + sum (grid%xa%qsn(ims:ime, jms:jme, kms:kme) * xa2_qsn(ims:ime, jms:jme, kms:kme))&
 203             + sum (grid%xa%qgr(ims:ime, jms:jme, kms:kme) * xa2_qgr(ims:ime, jms:jme, kms:kme))
 204       end if
 205    end if
 206 #ifdef VAR4D
 207    pertile_rhs = pertile_rhs &
 208       + sum (grid%x6a%qcw(ims:ime, jms:jme, kms:kme) * x6a2_qcw(ims:ime, jms:jme, kms:kme))    &
 209       + sum (grid%x6a%qci(ims:ime, jms:jme, kms:kme) * x6a2_qci(ims:ime, jms:jme, kms:kme))    &
 210       + sum (grid%x6a%qrn(ims:ime, jms:jme, kms:kme) * x6a2_qrn(ims:ime, jms:jme, kms:kme))    &
 211       + sum (grid%x6a%qsn(ims:ime, jms:jme, kms:kme) * x6a2_qsn(ims:ime, jms:jme, kms:kme))    &
 212       + sum (grid%x6a%qgr(ims:ime, jms:jme, kms:kme) * x6a2_qgr(ims:ime, jms:jme, kms:kme))
 213 #endif
 214
 215
 216    !----------------------------------------------------------------------
 217    ! [6.0] Calculate RHS of adjoint test equation:
 218    !----------------------------------------------------------------------
 219
 220    partial_rhs = sum (grid%xa%u(its:ite, jts:jte, kts:kte) * xa2_u(its:ite, jts:jte, kts:kte)) &
 221       + sum (grid%xa%v(its:ite, jts:jte, kts:kte) * xa2_v(its:ite, jts:jte, kts:kte))          &
 222       + sum (grid%xa%w(its:ite, jts:jte, kts:kte+1) * xa2_w(its:ite, jts:jte, kts:kte+1))      &
 223       + sum (grid%xa%t(its:ite, jts:jte, kts:kte) * xa2_t(its:ite, jts:jte, kts:kte))          &
 224       + sum (grid%xa%p(its:ite, jts:jte, kts:kte) * xa2_p(its:ite, jts:jte, kts:kte))          &
 225       + sum (grid%xa%q(its:ite, jts:jte, kts:kte) * xa2_q(its:ite, jts:jte, kts:kte))          &
 226       + sum (grid%xa%rh(its:ite, jts:jte, kts:kte)* xa2_rh(its:ite, jts:jte, kts:kte))         &
 227       + sum (grid%xa%psfc(its:ite, jts:jte) * xa2_psfc(its:ite, jts:jte))
 228 #ifdef VAR4D
 229    partial_rhs = partial_rhs &
 230       + sum (grid%x6a%u(its:ite, jts:jte, kts:kte) * x6a2_u(its:ite, jts:jte, kts:kte))        &
 231       + sum (grid%x6a%v(its:ite, jts:jte, kts:kte) * x6a2_v(its:ite, jts:jte, kts:kte))        &
 232       + sum (grid%x6a%w(its:ite, jts:jte, kts:kte+1) * x6a2_w(its:ite, jts:jte, kts:kte+1))    &
 233       + sum (grid%x6a%t(its:ite, jts:jte, kts:kte) * x6a2_t(its:ite, jts:jte, kts:kte))        &
 234       + sum (grid%x6a%p(its:ite, jts:jte, kts:kte) * x6a2_p(its:ite, jts:jte, kts:kte))        &
 235       + sum (grid%x6a%q(its:ite, jts:jte, kts:kte) * x6a2_q(its:ite, jts:jte, kts:kte))        &
 236       + sum (grid%x6a%rh(its:ite, jts:jte, kts:kte)* x6a2_rh(its:ite, jts:jte, kts:kte))       &
 237       + sum (grid%x6a%psfc(its:ite, jts:jte) * x6a2_psfc(its:ite, jts:jte))
 238 #endif
 239    if ( cloud_cv_options >= 1 ) then
 240       partial_rhs = partial_rhs &
 241          + sum (grid%xa%qcw(its:ite, jts:jte, kts:kte) * xa2_qcw(its:ite, jts:jte, kts:kte))   &
 242          + sum (grid%xa%qrn(its:ite, jts:jte, kts:kte) * xa2_qrn(its:ite, jts:jte, kts:kte))
 243       if ( cloud_cv_options >= 2 ) then
 244          partial_rhs = partial_rhs &
 245             + sum (grid%xa%qci(its:ite, jts:jte, kts:kte) * xa2_qci(its:ite, jts:jte, kts:kte)) &
 246             + sum (grid%xa%qsn(its:ite, jts:jte, kts:kte) * xa2_qsn(its:ite, jts:jte, kts:kte)) &
 247             + sum (grid%xa%qgr(its:ite, jts:jte, kts:kte) * xa2_qgr(its:ite, jts:jte, kts:kte))
 248       end if
 249    end if
 250 #ifdef VAR4D
 251    partial_rhs = partial_rhs &
 252       + sum (grid%x6a%qcw(its:ite, jts:jte, kts:kte) * x6a2_qcw(its:ite, jts:jte, kts:kte)) &
 253       + sum (grid%x6a%qci(its:ite, jts:jte, kts:kte) * x6a2_qci(its:ite, jts:jte, kts:kte)) &
 254       + sum (grid%x6a%qrn(its:ite, jts:jte, kts:kte) * x6a2_qrn(its:ite, jts:jte, kts:kte)) &
 255       + sum (grid%x6a%qsn(its:ite, jts:jte, kts:kte) * x6a2_qsn(its:ite, jts:jte, kts:kte)) &
 256       + sum (grid%x6a%qgr(its:ite, jts:jte, kts:kte) * x6a2_qgr(its:ite, jts:jte, kts:kte))
 257 #endif
 258
 259 #ifdef A2C
 260     if( ite == ide ) then
 261 print*,__FILE__,' contribution from ',ite+1,' col for U : ',sum (grid%xa%u(ite+1, jts:jte, kts:kte) * xa2_u(ite+1, jts:jte, kts:kte))
 262    partial_rhs = partial_rhs &
 263       + sum (grid%xa%u(ite+1, jts:jte, kts:kte) * xa2_u(ite+1, jts:jte, kts:kte))
 264    end if
 265     if( jte == jde ) then
 266 print*,__FILE__,' contribution from ',jte+1,' row for V : ',sum(grid%xa%v(its:ite, jte+1, kts:kte) * xa2_v(its:ite, jte+1, kts:kte))
 267    partial_rhs = partial_rhs &
 268       + sum (grid%xa%v(its:ite, jte+1, kts:kte) * xa2_v(its:ite, jte+1, kts:kte))
 269    end if
 270 #endif
 271
 272    !----------------------------------------------------------------------
 273    !  [7.0] Print output:
 274    !----------------------------------------------------------------------
 275    write (unit=stdout, fmt='(A,1pe22.14)') ' Single Domain < y, y     > = ', pertile_lhs
 276    write (unit=stdout, fmt='(A,1pe22.14)') ' Single Domain < x, x_adj > = ', pertile_rhs
 277
 278    adj_ttl_lhs = wrf_dm_sum_real (partial_lhs)
 279    adj_ttl_rhs = wrf_dm_sum_real (partial_rhs)
 280
 281    if (rootproc) then
 282       write(unit=stdout, fmt='(/)')
 283       write (unit=stdout, fmt='(A,1pe22.14)') ' Whole Domain < y, y     > = ', adj_ttl_lhs
 284       write (unit=stdout, fmt='(A,1pe22.14)') ' Whole Domain < x, x_adj > = ', adj_ttl_rhs
 285    end if
 286
 287    write (unit=stdout, fmt='(/a/)') 'da_check_dynamics_adjoint: Test Finished:'
 288    if (trace_use) call da_trace_exit("da_check_dynamics_adjoint")
 289
 290 end subroutine da_check_dynamics_adjoint
 291
 292