var/da/da_test/da_check_xtoy_adjoint.inc

   1 subroutine da_check_xtoy_adjoint(cv_size, cv, xbx, be, grid, config_flags, iv, y)
   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: Syed RH Rizvi,  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    type (iv_type),             intent(inout) :: iv    ! ob. increment vector.
  20    type (y_type),              intent(inout) :: y     ! y = h (grid%xa)
  21
  22    real                           :: adj_ttl_lhs   ! < y, y >
  23    real                           :: adj_ttl_rhs   ! < x, x_adj >
  24
  25    real                           :: partial_lhs   ! < y, y >
  26    real                           :: partial_rhs   ! < x, x_adj >
  27
  28    real                           :: pertile_lhs   ! < y, y >
  29    real                           :: pertile_rhs   ! < x, x_adj >
  30
  31    real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_u, xa2_v, xa2_t, &
  32                                                  xa2_p, xa2_q, xa2_rh
  33    real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_w
  34    real, dimension(ims:ime, jms:jme)          :: xa2_psfc
  35    real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_qcw, xa2_qci, xa2_qrn, xa2_qsn, xa2_qgr
  36    real, dimension(ims:ime, jms:jme, kms:kme) :: x6a2_u, x6a2_v, x6a2_t, &
  37                                                  x6a2_p, x6a2_q, x6a2_rh
  38    real, dimension(ims:ime, jms:jme, kms:kme) :: x6a2_w
  39    real, dimension(ims:ime, jms:jme)          :: x6a2_psfc
  40    real, dimension(ims:ime, jms:jme, kms:kme) :: x6a2_qcw, x6a2_qci, x6a2_qrn, x6a2_qsn, x6a2_qgr
  41    real, dimension(:,:,:), allocatable :: a_hr_rainc, a_hr_rainnc
  42    integer :: nobwin, i, j, k, fgat_rain
  43    character(len=4) :: filnam
  44    character(len=256) :: timestr
  45    integer :: time_step_seconds
  46    type(x_type) :: shuffle
  47    real             :: subarea, whole_area
  48
  49    if (trace_use) call da_trace_entry("da_check_xtoy_adjoint")
  50
  51    write (unit=stdout, fmt='(/a/)') 'da_check_xtoy_adjoint: Test Results:'
  52
  53    !----------------------------------------------------------------------
  54    ! [1.0] Initialise:
  55    !----------------------------------------------------------------------
  56
  57    partial_lhs = 0.0
  58    pertile_lhs = 0.0
  59
  60 #ifdef A2C
  61   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  62        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
  63      ipe = ipe + 1
  64      ide = ide + 1
  65   end if
  66
  67   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  68        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
  69      jpe = jpe + 1
  70      jde = jde + 1
  71   end if
  72 #endif
  73
  74 #ifdef DM_PARALLEL
  75 #include "HALO_XA.inc"
  76 #endif
  77
  78 #ifdef A2C
  79   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  80        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
  81      ipe = ipe - 1
  82      ide = ide - 1
  83   end if
  84
  85   if ((fg_format==fg_format_wrf_arw_regional  .or. &
  86        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
  87      jpe = jpe - 1
  88      jde = jde - 1
  89   end if
  90 #endif
  91    xa2_u(ims:ime, jms:jme, kms:kme) = grid%xa%u(ims:ime, jms:jme, kms:kme)
  92    xa2_v(ims:ime, jms:jme, kms:kme) = grid%xa%v(ims:ime, jms:jme, kms:kme)
  93    xa2_t(ims:ime, jms:jme, kms:kme) = grid%xa%t(ims:ime, jms:jme, kms:kme)
  94    xa2_p(ims:ime, jms:jme, kms:kme) = grid%xa%p(ims:ime, jms:jme, kms:kme)
  95    xa2_q(ims:ime, jms:jme, kms:kme) = grid%xa%q(ims:ime, jms:jme, kms:kme)
  96    xa2_w(ims:ime, jms:jme, kms:kme) = grid%xa%w(ims:ime, jms:jme, kms:kme)
  97    xa2_rh(ims:ime, jms:jme, kms:kme)= grid%xa%rh(ims:ime, jms:jme, kms:kme)
  98    xa2_psfc(ims:ime, jms:jme)       = grid%xa%psfc(ims:ime, jms:jme)
  99
 100    xa2_qcw = 0.0
 101    xa2_qrn = 0.0
 102    xa2_qci = 0.0
 103    xa2_qsn = 0.0
 104    xa2_qgr = 0.0
 105    if ( cloud_cv_options >= 1 ) then
 106       xa2_qcw(ims:ime, jms:jme, kms:kme) = grid%xa%qcw(ims:ime, jms:jme, kms:kme)
 107       xa2_qrn(ims:ime, jms:jme, kms:kme) = grid%xa%qrn(ims:ime, jms:jme, kms:kme)
 108       if ( cloud_cv_options >= 2 ) then
 109          xa2_qci(ims:ime, jms:jme, kms:kme) = grid%xa%qci(ims:ime, jms:jme, kms:kme)
 110          xa2_qsn(ims:ime, jms:jme, kms:kme) = grid%xa%qsn(ims:ime, jms:jme, kms:kme)
 111          xa2_qgr(ims:ime, jms:jme, kms:kme) = grid%xa%qgr(ims:ime, jms:jme, kms:kme)
 112       end if
 113    end if
 114
 115    x6a2_u = 0.0
 116    x6a2_v = 0.0
 117    x6a2_t = 0.0
 118    x6a2_p = 0.0
 119    x6a2_q = 0.0
 120    x6a2_w = 0.0
 121    x6a2_rh = 0.0
 122    x6a2_psfc = 0.0
 123
 124    x6a2_qcw = 0.0
 125    x6a2_qci = 0.0
 126    x6a2_qrn = 0.0
 127    x6a2_qsn = 0.0
 128    x6a2_qgr = 0.0
 129
 130 #ifdef A2C
 131     if( ite == ide ) &
 132 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))
 133     if( jte == jde ) &
 134 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))
 135 #endif
 136    if (var4d) then
 137 #ifdef VAR4D
 138       call domain_clock_get( grid, current_timestr=timestr )
 139       call da_transfer_xatowrftl(grid, config_flags, 'tl01', timestr)
 140
 141       if ( var4d_lbc ) then
 142          call domain_clock_get (grid, stop_timestr=timestr)
 143          call domain_clock_set( grid, current_timestr=timestr )
 144          grid%u_2 = u6_2 ; grid%v_2 = v6_2; grid%t_2 = t6_2;
 145          grid%w_2 = w6_2 ; grid%mu_2 = mu6_2 ; grid%ph_2 =ph6_2
 146          grid%moist = moist6; grid%p = p6; grid%psfc = psfc6
 147          call da_transfer_wrftoxb(xbx, grid, config_flags)
 148
 149          call da_zero_x(grid%x6a)
 150          shuffle = grid%xa
 151          grid%xa  = grid%x6a
 152          call da_setup_testfield(grid)
 153          grid%xa  = shuffle
 154
 155          x6a2_u(ims:ime, jms:jme, kms:kme) = grid%x6a%u(ims:ime, jms:jme, kms:kme)
 156          x6a2_v(ims:ime, jms:jme, kms:kme) = grid%x6a%v(ims:ime, jms:jme, kms:kme)
 157          x6a2_t(ims:ime, jms:jme, kms:kme) = grid%x6a%t(ims:ime, jms:jme, kms:kme)
 158          x6a2_p(ims:ime, jms:jme, kms:kme) = grid%x6a%p(ims:ime, jms:jme, kms:kme)
 159          x6a2_q(ims:ime, jms:jme, kms:kme) = grid%x6a%q(ims:ime, jms:jme, kms:kme)
 160          x6a2_w(ims:ime, jms:jme, kms:kme) = grid%x6a%w(ims:ime, jms:jme, kms:kme)
 161          x6a2_rh(ims:ime, jms:jme, kms:kme)= grid%x6a%rh(ims:ime, jms:jme, kms:kme)
 162          x6a2_psfc(ims:ime, jms:jme)       = grid%x6a%psfc(ims:ime, jms:jme)
 163
 164          x6a2_qcw(ims:ime, jms:jme, kms:kme) = grid%x6a%qcw(ims:ime, jms:jme, kms:kme)
 165          x6a2_qci(ims:ime, jms:jme, kms:kme) = grid%x6a%qci(ims:ime, jms:jme, kms:kme)
 166          x6a2_qrn(ims:ime, jms:jme, kms:kme) = grid%x6a%qrn(ims:ime, jms:jme, kms:kme)
 167          x6a2_qsn(ims:ime, jms:jme, kms:kme) = grid%x6a%qsn(ims:ime, jms:jme, kms:kme)
 168          x6a2_qgr(ims:ime, jms:jme, kms:kme) = grid%x6a%qgr(ims:ime, jms:jme, kms:kme)
 169
 170          call da_transfer_xatowrftl_lbc(grid, config_flags, 'tl01')
 171
 172          call domain_clock_get( grid, start_timestr=timestr )
 173          call domain_clock_set( grid, current_timestr=timestr )
 174          call da_read_basicstates ( xbx, grid, config_flags, timestr )
 175       else
 176          call da_transfer_xatowrftl_lbc(grid, config_flags, 'tl01')
 177       endif
 178
 179       call da_tl_model
 180
 181       if ( use_rainobs .and. num_fgat_time > 1 ) then
 182          allocate (a_hr_rainc (ims:ime,jms:jme,1:num_fgat_time))
 183          allocate (a_hr_rainnc(ims:ime,jms:jme,1:num_fgat_time))
 184          a_hr_rainc =0.0
 185          a_hr_rainnc=0.0
 186       endif
 187
 188       if (jcdfi_use) then
 189
 190          subarea = SUM ( grid%xb%grid_box_area(its:ite,jts:jte) )
 191          whole_area = wrf_dm_sum_real(subarea)
 192
 193          do j = jms, jme
 194             do k = kms, kme
 195                do i = ims, ime
 196                   pertile_lhs = pertile_lhs - config_flags%jcdfi_penalty * model_grid%jcdfi_u(i,k,j)**2 * &
 197                         grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 198                   pertile_lhs = pertile_lhs - config_flags%jcdfi_penalty * model_grid%jcdfi_v(i,k,j)**2 * &
 199                         grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 200                   pertile_lhs = pertile_lhs - config_flags%jcdfi_penalty * model_grid%jcdfi_t(i,k,j)**2 * &
 201                         (9.81/3.0)**2*grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 202                   pertile_lhs = pertile_lhs - config_flags%jcdfi_penalty * model_grid%jcdfi_p(i,k,j)**2 * &
 203                         (1.0/300.)**2*grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 204                enddo
 205             enddo
 206          enddo
 207
 208          ! We can not calculate the Y just with tile.
 209          partial_lhs = pertile_lhs
 210       endif
 211
 212 #else
 213       write(unit=message(1),fmt='(A)')'Please recompile the code with 4dvar option'
 214       call da_error(__FILE__,__LINE__,message(1:1))
 215 #endif
 216    end if
 217
 218    if ( num_fgat_time > 1 ) then
 219       call domain_clock_get (grid, stop_timestr=timestr)
 220       call domain_clock_set( grid, current_timestr=timestr )
 221       call domain_clock_set (grid, time_step_seconds=-1*var4d_bin)
 222       call domain_clockprint(150, grid, 'get CurrTime from clock,')
 223    endif
 224
 225    fgat_rain = num_fgat_time
 226    do nobwin= num_fgat_time, 1, -1
 227
 228       iv%time = nobwin
 229       iv%info(:)%n1 = iv%info(:)%plocal(iv%time-1) + 1
 230       iv%info(:)%n2 = iv%info(:)%plocal(iv%time)
 231
 232       if (var4d) then
 233 #ifdef VAR4D
 234          call domain_clock_get( grid, current_timestr=timestr )
 235          call da_read_basicstates ( xbx, grid, config_flags, timestr )
 236
 237          write(filnam,'(a2,i2.2)') 'tl',nobwin
 238          call da_transfer_wrftltoxa(grid, config_flags, filnam, timestr)
 239
 240          if ( use_rainobs ) then
 241             if ( num_fgat_time > 1 .and. fgat_rain_flags(nobwin) ) then
 242 !!!            We can not calculate the hourly rainfall in adjoint check, here, we just make sure the adjoint
 243 !!!            algorithm is correct mathmatically. so the amount of forcings doesn't matter
 244                a_hr_rainc (:,:,fgat_rain)=grid%g_rainc(:,:)
 245                a_hr_rainnc(:,:,fgat_rain)=grid%g_rainnc(:,:)
 246             endif
 247          endif
 248 #endif
 249       end if
 250
 251       call da_pt_to_rho_lin(grid)
 252 #ifdef DM_PARALLEL
 253 #include "HALO_XA.inc"
 254 #endif
 255
 256       if (sfc_assi_options == 2) then
 257          call da_transform_xtowtq (grid)
 258 #ifdef DM_PARALLEL
 259 #include "HALO_SFC_XA.inc"
 260 #endif
 261       end if
 262
 263       if (use_ssmt1obs .or. use_ssmt2obs .or. use_gpspwobs .or. &
 264           use_gpsztdobs .or. use_gpsrefobs .or. use_gpsephobs .or. &
 265           use_ssmitbobs .or. use_ssmiretrievalobs) then
 266
 267          ! Now do something for PW
 268          call da_transform_xtotpw(grid)
 269
 270          ! GPS Refractivity:
 271          if (use_gpsrefobs .or. use_gpsztdobs .or. use_gpsephobs) then
 272             call da_transform_xtogpsref_lin(grid)
 273             if (use_gpsztdobs) call da_transform_xtoztd_lin(grid)
 274          end if
 275
 276          if (use_ssmt1obs .or. use_ssmt2obs .or. &
 277              use_ssmitbobs .or. use_ssmiretrievalobs) then
 278             if (global) then
 279               call da_error(__FILE__,__LINE__, &
 280                 (/"grid%xb%speed is not available, see da_transfer_kmatoxb.inc"/))
 281             end if
 282             call da_transform_xtoseasfcwind_lin(grid)
 283          end if
 284
 285          if (use_ssmitbobs) call da_transform_xtotb_lin (grid)
 286
 287 #ifdef DM_PARALLEL
 288 #include "HALO_SSMI_XA.inc"
 289 #endif
 290       end if
 291
 292    ! Compute w increments using Richardson's eqn.
 293
 294    if ( Use_RadarObs ) then
 295       if ( .not. var4d ) call da_uvprho_to_w_lin(grid)
 296
 297       do k=kts,kte
 298          do j=jts,jte
 299             do i=its,ite
 300                grid%xa%wh(i,j,k)=0.5*(grid%xa%w(i,j,k)+grid%xa%w(i,j,k+1))
 301             end do
 302          end do
 303       end do
 304
 305 #ifdef DM_PARALLEL
 306 #include "HALO_RADAR_XA_W.inc"
 307 #endif
 308    end if
 309
 310    if ( cloud_cv_options == 1 )then
 311       ! Partition of hydrometeor increments via warm rain process
 312       call da_moist_phys_lin(grid)
 313    end if
 314
 315       !----------------------------------------------------------------------
 316       ! [2.0] Perform y = Hx transform:
 317       !----------------------------------------------------------------------
 318       call da_transform_xtoy (cv_size, cv, grid, iv, y)
 319
 320 #ifdef VAR4D
 321       if (iv%info(rain)%nlocal > 0 .and. var4d) &
 322          call da_transform_xtoy_rain (grid, iv, y, a_hr_rainc(:,:,nobwin), a_hr_rainnc(:,:,nobwin))
 323 #endif
 324
 325       !----------------------------------------------------------------------
 326       ! [3.0] Calculate LHS of adjoint test equation and
 327       !       Rescale input to adjoint routine :
 328       !----------------------------------------------------------------------
 329
 330       if (iv%info(sound)%nlocal > 0) call da_check_xtoy_adjoint_sound(iv, y, partial_lhs, pertile_lhs)
 331       if (iv%info(sonde_sfc)%nlocal > 0) call da_check_xtoy_adjoint_sonde_sfc (iv, y, partial_lhs, pertile_lhs)
 332       if (iv%info(mtgirs)%nlocal   > 0) call da_check_xtoy_adjoint_mtgirs   (iv, y, partial_lhs, pertile_lhs)
 333       if (iv%info(tamdar)%nlocal   > 0) call da_check_xtoy_adjoint_tamdar   (iv, y, partial_lhs, pertile_lhs)
 334       if (iv%info(tamdar_sfc)%nlocal   > 0) call da_check_xtoy_adjoint_tamdar_sfc(iv, y, partial_lhs, pertile_lhs)
 335       if (iv%info(synop)%nlocal    > 0) call da_check_xtoy_adjoint_synop    (iv, y, partial_lhs, pertile_lhs)
 336       if (iv%info(geoamv)%nlocal   > 0) call da_check_xtoy_adjoint_geoamv   (iv, y, partial_lhs, pertile_lhs)
 337       if (iv%info(polaramv)%nlocal > 0) call da_check_xtoy_adjoint_polaramv (iv, y, partial_lhs, pertile_lhs)
 338       if (iv%info(airep)%nlocal    > 0) call da_check_xtoy_adjoint_airep    (iv, y, partial_lhs, pertile_lhs)
 339       if (iv%info(pilot)%nlocal    > 0) call da_check_xtoy_adjoint_pilot    (iv, y, partial_lhs, pertile_lhs)
 340       if (iv%info(radar)%nlocal    > 0) call da_check_xtoy_adjoint_radar    (iv, y, partial_lhs, pertile_lhs)
 341       if (iv%info(satem)%nlocal    > 0) call da_check_xtoy_adjoint_satem    (iv, y, partial_lhs, pertile_lhs)
 342       if (iv%info(metar)%nlocal    > 0) call da_check_xtoy_adjoint_metar    (iv, y, partial_lhs, pertile_lhs)
 343       if (iv%info(ships)%nlocal    > 0) call da_check_xtoy_adjoint_ships    (iv, y, partial_lhs, pertile_lhs)
 344       if (iv%info(gpspw)%nlocal    > 0) call da_check_xtoy_adjoint_gpspw    (iv, y, partial_lhs, pertile_lhs)
 345       if (iv%info(gpsref)%nlocal   > 0) call da_check_xtoy_adjoint_gpsref   (iv, y, partial_lhs, pertile_lhs)
 346       if (iv%info(gpseph)%nlocal   > 0) call da_check_xtoy_adjoint_gpseph   (iv, y, partial_lhs, pertile_lhs)
 347       if (iv%info(ssmi_tb)%nlocal  > 0) call da_check_xtoy_adjoint_ssmi_tb  (iv, y, partial_lhs, pertile_lhs)
 348       if (iv%info(ssmi_rv)%nlocal  > 0) call da_check_xtoy_adjoint_ssmi_rv  (iv, y, partial_lhs, pertile_lhs)
 349       if (iv%info(ssmt2)%nlocal    > 0) call da_check_xtoy_adjoint_ssmt1    (iv, y, partial_lhs, pertile_lhs)
 350       if (iv%info(ssmt2)%nlocal    > 0) call da_check_xtoy_adjoint_ssmt2    (iv, y, partial_lhs, pertile_lhs)
 351       if (iv%info(qscat)%nlocal    > 0) call da_check_xtoy_adjoint_qscat    (iv, y, partial_lhs, pertile_lhs)
 352       if (iv%info(profiler)%nlocal > 0) call da_check_xtoy_adjoint_profiler (iv, y, partial_lhs, pertile_lhs)
 353       if (iv%info(buoy)%nlocal     > 0) call da_check_xtoy_adjoint_buoy     (iv, y, partial_lhs, pertile_lhs)
 354       if (iv%info(bogus)%nlocal    > 0) call da_check_xtoy_adjoint_bogus    (iv, y, partial_lhs, pertile_lhs)
 355       if (iv%num_inst              > 0) call da_check_xtoy_adjoint_rad      (iv, y, partial_lhs, pertile_lhs)
 356       if (iv%info(rain)%nlocal     > 0) call da_check_xtoy_adjoint_rain     (iv, y, partial_lhs, pertile_lhs)
 357       if (iv%info(pseudo)%nlocal   > 0) call da_check_xtoy_adjoint_pseudo   (iv, y, partial_lhs, pertile_lhs)
 358
 359       !----------------------------------------------------------------------
 360       ! [5.0] Perform adjoint operation:
 361       !----------------------------------------------------------------------
 362       call da_zero_x (grid%xa)
 363
 364       if (use_rainobs .and. num_fgat_time > 1) then
 365          a_hr_rainc(:,:,nobwin) = 0.0
 366          a_hr_rainnc(:,:,nobwin) = 0.0
 367       endif
 368
 369 #ifdef VAR4D
 370       if (iv%info(rain)%nlocal > 0 .and. var4d) then
 371          call da_transform_xtoy_rain_adj (grid, iv, y, a_hr_rainc(:,:,nobwin), a_hr_rainnc(:,:,nobwin))
 372       endif
 373 #endif
 374
 375       call da_transform_xtoy_adj (cv_size, cv, grid, iv, y, grid%xa)
 376
 377 #ifdef A2C
 378     if( ite == ide ) &
 379 print*,__FILE__,jte,' grid%xa%u.grid%xa%u for col= ',ite+1,sum(grid%xa%u(ite+1, jts:jte, kts:kte) * grid%xa%u(ite+1, jts:jte, kts:kte))
 380     if( jte == jde ) &
 381 print*,__FILE__,jte,' grid%xa%v.grid%x%%v for row= ',jte+1,sum(grid%xa%v(its:ite, jte+1, kts:kte) * grid%xa%v(its:ite, jte+1, kts:kte))
 382 #endif
 383
 384    ! Compute w increments using Richardson's eqn.
 385    if ( Use_RadarObs)  then
 386       do k=kts,kte
 387          do j=jts,jte
 388             do i=its,ite
 389                grid%xa%w(i,j,k)=grid%xa%w(i,j,k)+0.5*grid%xa%wh(i,j,k)
 390                grid%xa%w(i,j,k+1)=grid%xa%w(i,j,k+1)+0.5*grid%xa%wh(i,j,k)
 391                grid%xa%wh(i,j,k)=0.0
 392             end do
 393          end do
 394       end do
 395
 396       if ( .not. var4d ) call da_uvprho_to_w_adj(grid)
 397    end if
 398
 399    if ( cloud_cv_options == 1) then
 400       ! Partition of hydrometeor increments via warm rain process
 401       call da_moist_phys_adj(grid)
 402    end if
 403
 404       if (use_ssmt1obs .or. use_ssmt2obs .or. use_gpspwobs .or. &
 405           use_gpsztdobs .or. use_gpsrefobs .or. use_gpsephobs .or. &
 406           use_ssmitbobs .or. use_ssmiretrievalobs) then
 407
 408          if (use_ssmitbobs) call da_transform_xtotb_adj (grid)
 409
 410          ! for PW
 411          call da_transform_xtotpw_adj (grid)
 412
 413          ! GPS Refractivity:
 414          if (use_gpsrefobs .or. use_gpsztdobs .or. use_gpsephobs) then
 415             if (use_gpsztdobs) call da_transform_xtoztd_adj(grid)
 416             call da_transform_xtogpsref_adj (grid)
 417          end if
 418
 419          if (use_ssmt1obs .or. use_ssmt2obs .or. &
 420              use_ssmitbobs .or. use_ssmiretrievalobs) then
 421             if (global) then
 422                call da_error(__FILE__,__LINE__, &
 423                   (/"grid%xb%speed is not available, see da_transfer_kmatoxb.inc"/))
 424             end if
 425             call da_transform_xtoseasfcwind_adj (grid)
 426          end if
 427       end if
 428
 429       ! Now do something for surface variables
 430       if (sfc_assi_options == 2) then
 431          call da_transform_xtowtq_adj (grid)
 432
 433       end if
 434
 435       call da_pt_to_rho_adj (grid)
 436
 437       if (var4d) then
 438 #ifdef VAR4D
 439          grid%g_u_2 = 0.0
 440          grid%g_v_2 = 0.0
 441          grid%g_w_2 = 0.0
 442          grid%g_t_2 = 0.0
 443          grid%g_ph_2 = 0.0
 444          grid%g_p = 0.0
 445          grid%g_mu_2 = 0.0
 446          grid%g_moist = 0.0
 447          grid%g_rainnc = 0.0
 448          grid%g_rainncv = 0.0
 449          grid%g_rainc = 0.0
 450          grid%g_raincv = 0.0
 451
 452          write(unit=filnam,fmt='(a2,i2.2)') 'af',nobwin
 453
 454          if ( use_rainobs ) then
 455             if ( num_fgat_time > 1 .and. fgat_rain_flags(nobwin) ) then
 456 !!!            We can not calculate the hourly rainfall in adjoint check, just ensure the adjoint
 457 !!!            algorithm is right mathmatically. so the amount of forcings doesn't matter
 458                grid%g_rainc(:,:) = grid%g_rainc(:,:) + a_hr_rainc (:,:,fgat_rain)
 459                grid%g_rainnc(:,:) = grid%g_rainnc(:,:) + a_hr_rainnc(:,:,fgat_rain)
 460                a_hr_rainc (:,:,fgat_rain) = 0.0
 461                a_hr_rainnc(:,:,fgat_rain) = 0.0
 462                fgat_rain = fgat_rain - 1
 463             endif
 464          endif
 465
 466          call domain_clock_get( grid, current_timestr=timestr )
 467          call da_transfer_wrftltoxa_adj(grid, config_flags, filnam, timestr)
 468 #endif
 469       end if
 470
 471       if ( nobwin > 1 ) call domain_clockadvance (grid)
 472       call domain_clockprint(150, grid, 'DEBUG Adjoint Check:  get CurrTime from clock,')
 473
 474    end do
 475
 476    if ( num_fgat_time > 1 ) then
 477       call nl_get_time_step ( grid%id, time_step_seconds)
 478       call domain_clock_set (grid, time_step_seconds=time_step_seconds)
 479       call domain_clockprint(150, grid, 'get CurrTime from clock,')
 480    endif
 481
 482    if (var4d) then
 483 #ifdef VAR4D
 484       if (jcdfi_use) then
 485
 486          do j = jms, jme
 487             do k = kms, kme
 488                do i = ims, ime
 489                   model_grid%jcdfi_u(i,k,j) = -1.0 * config_flags%jcdfi_penalty * model_grid%jcdfi_u(i,k,j) * &
 490                         grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 491                   model_grid%jcdfi_v(i,k,j) = -1.0 * config_flags%jcdfi_penalty * model_grid%jcdfi_v(i,k,j) * &
 492                         grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 493                   model_grid%jcdfi_t(i,k,j) = -1.0 * config_flags%jcdfi_penalty * model_grid%jcdfi_t(i,k,j) * &
 494                         (9.81/3.0)**2*grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 495                   model_grid%jcdfi_p(i,k,j) = -1.0 * config_flags%jcdfi_penalty * model_grid%jcdfi_p(i,k,j) * &
 496                         (1.0/300.)**2*grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 497                enddo
 498             enddo
 499          enddo
 500
 501       endif
 502
 503       if ( trajectory_io ) then
 504          ! for memory io, we need to up-side-down the adjoint forcing linked list generated in previous step.
 505          call upsidedown_ad_forcing
 506       endif
 507
 508       call da_ad_model
 509
 510       call da_zero_x(grid%x6a)
 511       if (var4d_lbc) then
 512          call domain_clock_get (grid, stop_timestr=timestr)
 513          call domain_clock_set( grid, current_timestr=timestr )
 514          grid%u_2 = u6_2 ; grid%v_2 = v6_2; grid%t_2 = t6_2;
 515          grid%w_2 = w6_2 ; grid%mu_2 = mu6_2 ; grid%ph_2 =ph6_2
 516          grid%moist = moist6; grid%p = p6; grid%psfc = psfc6
 517          call da_transfer_wrftoxb(xbx, grid, config_flags)
 518
 519          call da_transfer_xatowrftl_adj_lbc(grid, config_flags, 'gr01')
 520
 521          call domain_clock_get( grid, start_timestr=timestr )
 522          call domain_clock_set( grid, current_timestr=timestr )
 523          call da_read_basicstates ( xbx, grid, config_flags, timestr )
 524       end if
 525
 526       call da_zero_x(grid%xa)
 527       call da_transfer_xatowrftl_adj(grid, config_flags, 'gr01')
 528
 529       if ( use_rainobs .and. num_fgat_time > 1 ) then
 530          deallocate (a_hr_rainc)
 531          deallocate (a_hr_rainnc)
 532       endif
 533 #endif
 534    end if
 535
 536
 537    pertile_rhs = sum (grid%xa%u(ims:ime, jms:jme, kms:kme) * xa2_u(ims:ime, jms:jme, kms:kme)) &
 538       + sum (grid%xa%v(ims:ime, jms:jme, kms:kme) * xa2_v(ims:ime, jms:jme, kms:kme))          &
 539       + sum (grid%xa%w(ims:ime, jms:jme, kms:kme) * xa2_w(ims:ime, jms:jme, kms:kme))          &
 540       + sum (grid%xa%t(ims:ime, jms:jme, kms:kme) * xa2_t(ims:ime, jms:jme, kms:kme))          &
 541       + sum (grid%xa%p(ims:ime, jms:jme, kms:kme) * xa2_p(ims:ime, jms:jme, kms:kme))          &
 542       + sum (grid%xa%q(ims:ime, jms:jme, kms:kme) * xa2_q(ims:ime, jms:jme, kms:kme))          &
 543       + sum (grid%xa%rh(ims:ime, jms:jme, kms:kme)* xa2_rh(ims:ime, jms:jme, kms:kme))         &
 544       + sum (grid%xa%psfc(ims:ime, jms:jme) * xa2_psfc(ims:ime, jms:jme))
 545 #ifdef VAR4D
 546    pertile_rhs = pertile_rhs &
 547       + sum (grid%x6a%u(ims:ime, jms:jme, kms:kme) * x6a2_u(ims:ime, jms:jme, kms:kme))        &
 548       + sum (grid%x6a%v(ims:ime, jms:jme, kms:kme) * x6a2_v(ims:ime, jms:jme, kms:kme))        &
 549       + sum (grid%x6a%w(ims:ime, jms:jme, kms:kme) * x6a2_w(ims:ime, jms:jme, kms:kme))        &
 550       + sum (grid%x6a%t(ims:ime, jms:jme, kms:kme) * x6a2_t(ims:ime, jms:jme, kms:kme))        &
 551       + sum (grid%x6a%p(ims:ime, jms:jme, kms:kme) * x6a2_p(ims:ime, jms:jme, kms:kme))        &
 552       + sum (grid%x6a%q(ims:ime, jms:jme, kms:kme) * x6a2_q(ims:ime, jms:jme, kms:kme))        &
 553       + sum (grid%x6a%rh(ims:ime, jms:jme, kms:kme)* x6a2_rh(ims:ime, jms:jme, kms:kme))       &
 554       + sum (grid%x6a%psfc(ims:ime, jms:jme) * x6a2_psfc(ims:ime, jms:jme))
 555 #endif
 556    if ( cloud_cv_options >= 1 ) then
 557       pertile_rhs = pertile_rhs &
 558          + sum (grid%xa%qcw(ims:ime, jms:jme, kms:kme) * xa2_qcw(ims:ime, jms:jme, kms:kme))   &
 559          + sum (grid%xa%qrn(ims:ime, jms:jme, kms:kme) * xa2_qrn(ims:ime, jms:jme, kms:kme))
 560       if ( cloud_cv_options >= 2 ) then
 561          pertile_rhs = pertile_rhs &
 562             + sum (grid%xa%qci(ims:ime, jms:jme, kms:kme) * xa2_qci(ims:ime, jms:jme, kms:kme))&
 563             + sum (grid%xa%qsn(ims:ime, jms:jme, kms:kme) * xa2_qsn(ims:ime, jms:jme, kms:kme))&
 564             + sum (grid%xa%qgr(ims:ime, jms:jme, kms:kme) * xa2_qgr(ims:ime, jms:jme, kms:kme))
 565       end if
 566    end if
 567 #ifdef VAR4D
 568    pertile_rhs = pertile_rhs &
 569       + sum (grid%x6a%qcw(ims:ime, jms:jme, kms:kme) * x6a2_qcw(ims:ime, jms:jme, kms:kme))    &
 570       + sum (grid%x6a%qci(ims:ime, jms:jme, kms:kme) * x6a2_qci(ims:ime, jms:jme, kms:kme))    &
 571       + sum (grid%x6a%qrn(ims:ime, jms:jme, kms:kme) * x6a2_qrn(ims:ime, jms:jme, kms:kme))    &
 572       + sum (grid%x6a%qsn(ims:ime, jms:jme, kms:kme) * x6a2_qsn(ims:ime, jms:jme, kms:kme))    &
 573       + sum (grid%x6a%qgr(ims:ime, jms:jme, kms:kme) * x6a2_qgr(ims:ime, jms:jme, kms:kme))
 574 #endif
 575
 576
 577    !----------------------------------------------------------------------
 578    ! [6.0] Calculate RHS of adjoint test equation:
 579    !----------------------------------------------------------------------
 580
 581    partial_rhs = sum (grid%xa%u(its:ite, jts:jte, kts:kte) * xa2_u(its:ite, jts:jte, kts:kte)) &
 582       + sum (grid%xa%v(its:ite, jts:jte, kts:kte) * xa2_v(its:ite, jts:jte, kts:kte))          &
 583       + sum (grid%xa%w(its:ite, jts:jte, kts:kte+1) * xa2_w(its:ite, jts:jte, kts:kte+1))      &
 584       + sum (grid%xa%t(its:ite, jts:jte, kts:kte) * xa2_t(its:ite, jts:jte, kts:kte))          &
 585       + sum (grid%xa%p(its:ite, jts:jte, kts:kte) * xa2_p(its:ite, jts:jte, kts:kte))          &
 586       + sum (grid%xa%q(its:ite, jts:jte, kts:kte) * xa2_q(its:ite, jts:jte, kts:kte))          &
 587       + sum (grid%xa%rh(its:ite, jts:jte, kts:kte)* xa2_rh(its:ite, jts:jte, kts:kte))         &
 588       + sum (grid%xa%psfc(its:ite, jts:jte) * xa2_psfc(its:ite, jts:jte))
 589 #ifdef VAR4D
 590    partial_rhs = partial_rhs &
 591       + sum (grid%x6a%u(its:ite, jts:jte, kts:kte) * x6a2_u(its:ite, jts:jte, kts:kte))        &
 592       + sum (grid%x6a%v(its:ite, jts:jte, kts:kte) * x6a2_v(its:ite, jts:jte, kts:kte))        &
 593       + sum (grid%x6a%w(its:ite, jts:jte, kts:kte+1) * x6a2_w(its:ite, jts:jte, kts:kte+1))    &
 594       + sum (grid%x6a%t(its:ite, jts:jte, kts:kte) * x6a2_t(its:ite, jts:jte, kts:kte))        &
 595       + sum (grid%x6a%p(its:ite, jts:jte, kts:kte) * x6a2_p(its:ite, jts:jte, kts:kte))        &
 596       + sum (grid%x6a%q(its:ite, jts:jte, kts:kte) * x6a2_q(its:ite, jts:jte, kts:kte))        &
 597       + sum (grid%x6a%rh(its:ite, jts:jte, kts:kte)* x6a2_rh(its:ite, jts:jte, kts:kte))       &
 598       + sum (grid%x6a%psfc(its:ite, jts:jte) * x6a2_psfc(its:ite, jts:jte))
 599 #endif
 600
 601    if ( cloud_cv_options >= 1 ) then
 602       partial_rhs = partial_rhs &
 603          + sum (grid%xa%qcw(its:ite, jts:jte, kts:kte) * xa2_qcw(its:ite, jts:jte, kts:kte))   &
 604          + sum (grid%xa%qrn(its:ite, jts:jte, kts:kte) * xa2_qrn(its:ite, jts:jte, kts:kte))
 605       if ( cloud_cv_options >= 2 ) then
 606          partial_rhs = partial_rhs &
 607             + sum (grid%xa%qci(its:ite, jts:jte, kts:kte) * xa2_qci(its:ite, jts:jte, kts:kte)) &
 608             + sum (grid%xa%qsn(its:ite, jts:jte, kts:kte) * xa2_qsn(its:ite, jts:jte, kts:kte)) &
 609             + sum (grid%xa%qgr(its:ite, jts:jte, kts:kte) * xa2_qgr(its:ite, jts:jte, kts:kte))
 610       end if
 611    end if
 612 #ifdef VAR4D
 613    partial_rhs = partial_rhs &
 614       + sum (grid%x6a%qcw(its:ite, jts:jte, kts:kte) * x6a2_qcw(its:ite, jts:jte, kts:kte)) &
 615       + sum (grid%x6a%qci(its:ite, jts:jte, kts:kte) * x6a2_qci(its:ite, jts:jte, kts:kte)) &
 616       + sum (grid%x6a%qrn(its:ite, jts:jte, kts:kte) * x6a2_qrn(its:ite, jts:jte, kts:kte)) &
 617       + sum (grid%x6a%qsn(its:ite, jts:jte, kts:kte) * x6a2_qsn(its:ite, jts:jte, kts:kte)) &
 618       + sum (grid%x6a%qgr(its:ite, jts:jte, kts:kte) * x6a2_qgr(its:ite, jts:jte, kts:kte))
 619 #endif
 620
 621 #ifdef A2C
 622     if( ite == ide ) then
 623 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))
 624    partial_rhs = partial_rhs &
 625       + sum (grid%xa%u(ite+1, jts:jte, kts:kte) * xa2_u(ite+1, jts:jte, kts:kte))
 626    end if
 627     if( jte == jde ) then
 628 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))
 629    partial_rhs = partial_rhs &
 630       + sum (grid%xa%v(its:ite, jte+1, kts:kte) * xa2_v(its:ite, jte+1, kts:kte))
 631    end if
 632 #endif
 633
 634    !----------------------------------------------------------------------
 635    !  [7.0] Print output:
 636    !----------------------------------------------------------------------
 637    write (unit=stdout, fmt='(A,1pe22.14)') ' Single Domain < y, y     > = ', pertile_lhs
 638    write (unit=stdout, fmt='(A,1pe22.14)') ' Single Domain < x, x_adj > = ', pertile_rhs
 639
 640    adj_ttl_lhs = wrf_dm_sum_real (partial_lhs)
 641    adj_ttl_rhs = wrf_dm_sum_real (partial_rhs)
 642
 643    if (rootproc) then
 644       write(unit=stdout, fmt='(/)')
 645       write (unit=stdout, fmt='(A,1pe22.14)') ' Whole Domain < y, y     > = ', adj_ttl_lhs
 646       write (unit=stdout, fmt='(A,1pe22.14)') ' Whole Domain < x, x_adj > = ', adj_ttl_rhs
 647    end if
 648
 649    write (unit=stdout, fmt='(/a/)') 'da_check_xtoy_adjoint: Test Finished:'
 650    if (trace_use) call da_trace_exit("da_check_xtoy_adjoint")
 651
 652 end subroutine da_check_xtoy_adjoint
 653
 654