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(lightning)%nlocal> 0) call da_check_xtoy_adjoint_lightning(iv, y, partial_lhs, pertile_lhs)
 342       if (iv%info(satem)%nlocal    > 0) call da_check_xtoy_adjoint_satem    (iv, y, partial_lhs, pertile_lhs)
 343       if (iv%info(metar)%nlocal    > 0) call da_check_xtoy_adjoint_metar    (iv, y, partial_lhs, pertile_lhs)
 344       if (iv%info(ships)%nlocal    > 0) call da_check_xtoy_adjoint_ships    (iv, y, partial_lhs, pertile_lhs)
 345       if (iv%info(gpspw)%nlocal    > 0) call da_check_xtoy_adjoint_gpspw    (iv, y, partial_lhs, pertile_lhs)
 346       if (iv%info(gpsref)%nlocal   > 0) call da_check_xtoy_adjoint_gpsref   (iv, y, partial_lhs, pertile_lhs)
 347       if (iv%info(gpseph)%nlocal   > 0) call da_check_xtoy_adjoint_gpseph   (iv, y, partial_lhs, pertile_lhs)
 348       if (iv%info(ssmi_tb)%nlocal  > 0) call da_check_xtoy_adjoint_ssmi_tb  (iv, y, partial_lhs, pertile_lhs)
 349       if (iv%info(ssmi_rv)%nlocal  > 0) call da_check_xtoy_adjoint_ssmi_rv  (iv, y, partial_lhs, pertile_lhs)
 350       if (iv%info(ssmt2)%nlocal    > 0) call da_check_xtoy_adjoint_ssmt1    (iv, y, partial_lhs, pertile_lhs)
 351       if (iv%info(ssmt2)%nlocal    > 0) call da_check_xtoy_adjoint_ssmt2    (iv, y, partial_lhs, pertile_lhs)
 352       if (iv%info(qscat)%nlocal    > 0) call da_check_xtoy_adjoint_qscat    (iv, y, partial_lhs, pertile_lhs)
 353       if (iv%info(profiler)%nlocal > 0) call da_check_xtoy_adjoint_profiler (iv, y, partial_lhs, pertile_lhs)
 354       if (iv%info(buoy)%nlocal     > 0) call da_check_xtoy_adjoint_buoy     (iv, y, partial_lhs, pertile_lhs)
 355       if (iv%info(bogus)%nlocal    > 0) call da_check_xtoy_adjoint_bogus    (iv, y, partial_lhs, pertile_lhs)
 356       if (iv%num_inst              > 0) call da_check_xtoy_adjoint_rad      (iv, y, partial_lhs, pertile_lhs)
 357       if (iv%info(rain)%nlocal     > 0) call da_check_xtoy_adjoint_rain     (iv, y, partial_lhs, pertile_lhs)
 358       if (iv%info(pseudo)%nlocal   > 0) call da_check_xtoy_adjoint_pseudo   (iv, y, partial_lhs, pertile_lhs)
 359
 360       !----------------------------------------------------------------------
 361       ! [5.0] Perform adjoint operation:
 362       !----------------------------------------------------------------------
 363       call da_zero_x (grid%xa)
 364
 365       if (use_rainobs .and. num_fgat_time > 1) then
 366          a_hr_rainc(:,:,nobwin) = 0.0
 367          a_hr_rainnc(:,:,nobwin) = 0.0
 368       endif
 369
 370 #ifdef VAR4D
 371       if (iv%info(rain)%nlocal > 0 .and. var4d) then
 372          call da_transform_xtoy_rain_adj (grid, iv, y, a_hr_rainc(:,:,nobwin), a_hr_rainnc(:,:,nobwin))
 373       endif
 374 #endif
 375
 376       call da_transform_xtoy_adj (cv_size, cv, grid, iv, y, grid%xa)
 377
 378 #ifdef A2C
 379     if( ite == ide ) &
 380 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))
 381     if( jte == jde ) &
 382 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))
 383 #endif
 384
 385    ! Compute w increments using Richardson's eqn.
 386    if ( Use_RadarObs)  then
 387       do k=kts,kte
 388          do j=jts,jte
 389             do i=its,ite
 390                grid%xa%w(i,j,k)=grid%xa%w(i,j,k)+0.5*grid%xa%wh(i,j,k)
 391                grid%xa%w(i,j,k+1)=grid%xa%w(i,j,k+1)+0.5*grid%xa%wh(i,j,k)
 392                grid%xa%wh(i,j,k)=0.0
 393             end do
 394          end do
 395       end do
 396
 397       if ( .not. var4d ) call da_uvprho_to_w_adj(grid)
 398    end if
 399
 400    if ( cloud_cv_options == 1) then
 401       ! Partition of hydrometeor increments via warm rain process
 402       call da_moist_phys_adj(grid)
 403    end if
 404
 405       if (use_ssmt1obs .or. use_ssmt2obs .or. use_gpspwobs .or. &
 406           use_gpsztdobs .or. use_gpsrefobs .or. use_gpsephobs .or. &
 407           use_ssmitbobs .or. use_ssmiretrievalobs) then
 408
 409          if (use_ssmitbobs) call da_transform_xtotb_adj (grid)
 410
 411          ! for PW
 412          call da_transform_xtotpw_adj (grid)
 413
 414          ! GPS Refractivity:
 415          if (use_gpsrefobs .or. use_gpsztdobs .or. use_gpsephobs) then
 416             if (use_gpsztdobs) call da_transform_xtoztd_adj(grid)
 417             call da_transform_xtogpsref_adj (grid)
 418          end if
 419
 420          if (use_ssmt1obs .or. use_ssmt2obs .or. &
 421              use_ssmitbobs .or. use_ssmiretrievalobs) then
 422             if (global) then
 423                call da_error(__FILE__,__LINE__, &
 424                   (/"grid%xb%speed is not available, see da_transfer_kmatoxb.inc"/))
 425             end if
 426             call da_transform_xtoseasfcwind_adj (grid)
 427          end if
 428       end if
 429
 430       ! Now do something for surface variables
 431       if (sfc_assi_options == 2) then
 432          call da_transform_xtowtq_adj (grid)
 433
 434       end if
 435
 436       call da_pt_to_rho_adj (grid)
 437
 438       if (var4d) then
 439 #ifdef VAR4D
 440          grid%g_u_2 = 0.0
 441          grid%g_v_2 = 0.0
 442          grid%g_w_2 = 0.0
 443          grid%g_t_2 = 0.0
 444          grid%g_ph_2 = 0.0
 445          grid%g_p = 0.0
 446          grid%g_mu_2 = 0.0
 447          grid%g_moist = 0.0
 448          grid%g_rainnc = 0.0
 449          grid%g_rainncv = 0.0
 450          grid%g_rainc = 0.0
 451          grid%g_raincv = 0.0
 452
 453          write(unit=filnam,fmt='(a2,i2.2)') 'af',nobwin
 454
 455          if ( use_rainobs ) then
 456             if ( num_fgat_time > 1 .and. fgat_rain_flags(nobwin) ) then
 457 !!!            We can not calculate the hourly rainfall in adjoint check, just ensure the adjoint
 458 !!!            algorithm is right mathmatically. so the amount of forcings doesn't matter
 459                grid%g_rainc(:,:) = grid%g_rainc(:,:) + a_hr_rainc (:,:,fgat_rain)
 460                grid%g_rainnc(:,:) = grid%g_rainnc(:,:) + a_hr_rainnc(:,:,fgat_rain)
 461                a_hr_rainc (:,:,fgat_rain) = 0.0
 462                a_hr_rainnc(:,:,fgat_rain) = 0.0
 463                fgat_rain = fgat_rain - 1
 464             endif
 465          endif
 466
 467          call domain_clock_get( grid, current_timestr=timestr )
 468          call da_transfer_wrftltoxa_adj(grid, config_flags, filnam, timestr)
 469 #endif
 470       end if
 471
 472       if ( nobwin > 1 ) call domain_clockadvance (grid)
 473       call domain_clockprint(150, grid, 'DEBUG Adjoint Check:  get CurrTime from clock,')
 474
 475    end do
 476
 477    if ( num_fgat_time > 1 ) then
 478       call nl_get_time_step ( grid%id, time_step_seconds)
 479       call domain_clock_set (grid, time_step_seconds=time_step_seconds)
 480       call domain_clockprint(150, grid, 'get CurrTime from clock,')
 481    endif
 482
 483    if (var4d) then
 484 #ifdef VAR4D
 485       if (jcdfi_use) then
 486
 487          do j = jms, jme
 488             do k = kms, kme
 489                do i = ims, ime
 490                   model_grid%jcdfi_u(i,k,j) = -1.0 * config_flags%jcdfi_penalty * model_grid%jcdfi_u(i,k,j) * &
 491                         grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 492                   model_grid%jcdfi_v(i,k,j) = -1.0 * config_flags%jcdfi_penalty * model_grid%jcdfi_v(i,k,j) * &
 493                         grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 494                   model_grid%jcdfi_t(i,k,j) = -1.0 * config_flags%jcdfi_penalty * model_grid%jcdfi_t(i,k,j) * &
 495                         (9.81/3.0)**2*grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 496                   model_grid%jcdfi_p(i,k,j) = -1.0 * config_flags%jcdfi_penalty * model_grid%jcdfi_p(i,k,j) * &
 497                         (1.0/300.)**2*grid%xb%grid_box_area(i,j)/whole_area*grid%xb%dnw(k)
 498                enddo
 499             enddo
 500          enddo
 501
 502       endif
 503
 504       if ( trajectory_io ) then
 505          ! for memory io, we need to up-side-down the adjoint forcing linked list generated in previous step.
 506          call upsidedown_ad_forcing
 507       endif
 508
 509       call da_ad_model
 510
 511       call da_zero_x(grid%x6a)
 512       if (var4d_lbc) then
 513          call domain_clock_get (grid, stop_timestr=timestr)
 514          call domain_clock_set( grid, current_timestr=timestr )
 515          grid%u_2 = u6_2 ; grid%v_2 = v6_2; grid%t_2 = t6_2;
 516          grid%w_2 = w6_2 ; grid%mu_2 = mu6_2 ; grid%ph_2 =ph6_2
 517          grid%moist = moist6; grid%p = p6; grid%psfc = psfc6
 518          call da_transfer_wrftoxb(xbx, grid, config_flags)
 519
 520          call da_transfer_xatowrftl_adj_lbc(grid, config_flags, 'gr01')
 521
 522          call domain_clock_get( grid, start_timestr=timestr )
 523          call domain_clock_set( grid, current_timestr=timestr )
 524          call da_read_basicstates ( xbx, grid, config_flags, timestr )
 525       end if
 526
 527       call da_zero_x(grid%xa)
 528       call da_transfer_xatowrftl_adj(grid, config_flags, 'gr01')
 529
 530       if ( use_rainobs .and. num_fgat_time > 1 ) then
 531          deallocate (a_hr_rainc)
 532          deallocate (a_hr_rainnc)
 533       endif
 534 #endif
 535    end if
 536
 537
 538    pertile_rhs = sum (grid%xa%u(ims:ime, jms:jme, kms:kme) * xa2_u(ims:ime, jms:jme, kms:kme)) &
 539       + sum (grid%xa%v(ims:ime, jms:jme, kms:kme) * xa2_v(ims:ime, jms:jme, kms:kme))          &
 540       + sum (grid%xa%w(ims:ime, jms:jme, kms:kme) * xa2_w(ims:ime, jms:jme, kms:kme))          &
 541       + sum (grid%xa%t(ims:ime, jms:jme, kms:kme) * xa2_t(ims:ime, jms:jme, kms:kme))          &
 542       + sum (grid%xa%p(ims:ime, jms:jme, kms:kme) * xa2_p(ims:ime, jms:jme, kms:kme))          &
 543       + sum (grid%xa%q(ims:ime, jms:jme, kms:kme) * xa2_q(ims:ime, jms:jme, kms:kme))          &
 544       + sum (grid%xa%rh(ims:ime, jms:jme, kms:kme)* xa2_rh(ims:ime, jms:jme, kms:kme))         &
 545       + sum (grid%xa%psfc(ims:ime, jms:jme) * xa2_psfc(ims:ime, jms:jme))
 546 #ifdef VAR4D
 547    pertile_rhs = pertile_rhs &
 548       + sum (grid%x6a%u(ims:ime, jms:jme, kms:kme) * x6a2_u(ims:ime, jms:jme, kms:kme))        &
 549       + sum (grid%x6a%v(ims:ime, jms:jme, kms:kme) * x6a2_v(ims:ime, jms:jme, kms:kme))        &
 550       + sum (grid%x6a%w(ims:ime, jms:jme, kms:kme) * x6a2_w(ims:ime, jms:jme, kms:kme))        &
 551       + sum (grid%x6a%t(ims:ime, jms:jme, kms:kme) * x6a2_t(ims:ime, jms:jme, kms:kme))        &
 552       + sum (grid%x6a%p(ims:ime, jms:jme, kms:kme) * x6a2_p(ims:ime, jms:jme, kms:kme))        &
 553       + sum (grid%x6a%q(ims:ime, jms:jme, kms:kme) * x6a2_q(ims:ime, jms:jme, kms:kme))        &
 554       + sum (grid%x6a%rh(ims:ime, jms:jme, kms:kme)* x6a2_rh(ims:ime, jms:jme, kms:kme))       &
 555       + sum (grid%x6a%psfc(ims:ime, jms:jme) * x6a2_psfc(ims:ime, jms:jme))
 556 #endif
 557    if ( cloud_cv_options >= 1 ) then
 558       pertile_rhs = pertile_rhs &
 559          + sum (grid%xa%qcw(ims:ime, jms:jme, kms:kme) * xa2_qcw(ims:ime, jms:jme, kms:kme))   &
 560          + sum (grid%xa%qrn(ims:ime, jms:jme, kms:kme) * xa2_qrn(ims:ime, jms:jme, kms:kme))
 561       if ( cloud_cv_options >= 2 ) then
 562          pertile_rhs = pertile_rhs &
 563             + sum (grid%xa%qci(ims:ime, jms:jme, kms:kme) * xa2_qci(ims:ime, jms:jme, kms:kme))&
 564             + sum (grid%xa%qsn(ims:ime, jms:jme, kms:kme) * xa2_qsn(ims:ime, jms:jme, kms:kme))&
 565             + sum (grid%xa%qgr(ims:ime, jms:jme, kms:kme) * xa2_qgr(ims:ime, jms:jme, kms:kme))
 566       end if
 567    end if
 568 #ifdef VAR4D
 569    pertile_rhs = pertile_rhs &
 570       + sum (grid%x6a%qcw(ims:ime, jms:jme, kms:kme) * x6a2_qcw(ims:ime, jms:jme, kms:kme))    &
 571       + sum (grid%x6a%qci(ims:ime, jms:jme, kms:kme) * x6a2_qci(ims:ime, jms:jme, kms:kme))    &
 572       + sum (grid%x6a%qrn(ims:ime, jms:jme, kms:kme) * x6a2_qrn(ims:ime, jms:jme, kms:kme))    &
 573       + sum (grid%x6a%qsn(ims:ime, jms:jme, kms:kme) * x6a2_qsn(ims:ime, jms:jme, kms:kme))    &
 574       + sum (grid%x6a%qgr(ims:ime, jms:jme, kms:kme) * x6a2_qgr(ims:ime, jms:jme, kms:kme))
 575 #endif
 576
 577
 578    !----------------------------------------------------------------------
 579    ! [6.0] Calculate RHS of adjoint test equation:
 580    !----------------------------------------------------------------------
 581
 582    partial_rhs = sum (grid%xa%u(its:ite, jts:jte, kts:kte) * xa2_u(its:ite, jts:jte, kts:kte)) &
 583       + sum (grid%xa%v(its:ite, jts:jte, kts:kte) * xa2_v(its:ite, jts:jte, kts:kte))          &
 584       + sum (grid%xa%w(its:ite, jts:jte, kts:kte+1) * xa2_w(its:ite, jts:jte, kts:kte+1))      &
 585       + sum (grid%xa%t(its:ite, jts:jte, kts:kte) * xa2_t(its:ite, jts:jte, kts:kte))          &
 586       + sum (grid%xa%p(its:ite, jts:jte, kts:kte) * xa2_p(its:ite, jts:jte, kts:kte))          &
 587       + sum (grid%xa%q(its:ite, jts:jte, kts:kte) * xa2_q(its:ite, jts:jte, kts:kte))          &
 588       + sum (grid%xa%rh(its:ite, jts:jte, kts:kte)* xa2_rh(its:ite, jts:jte, kts:kte))         &
 589       + sum (grid%xa%psfc(its:ite, jts:jte) * xa2_psfc(its:ite, jts:jte))
 590 #ifdef VAR4D
 591    partial_rhs = partial_rhs &
 592       + sum (grid%x6a%u(its:ite, jts:jte, kts:kte) * x6a2_u(its:ite, jts:jte, kts:kte))        &
 593       + sum (grid%x6a%v(its:ite, jts:jte, kts:kte) * x6a2_v(its:ite, jts:jte, kts:kte))        &
 594       + sum (grid%x6a%w(its:ite, jts:jte, kts:kte+1) * x6a2_w(its:ite, jts:jte, kts:kte+1))    &
 595       + sum (grid%x6a%t(its:ite, jts:jte, kts:kte) * x6a2_t(its:ite, jts:jte, kts:kte))        &
 596       + sum (grid%x6a%p(its:ite, jts:jte, kts:kte) * x6a2_p(its:ite, jts:jte, kts:kte))        &
 597       + sum (grid%x6a%q(its:ite, jts:jte, kts:kte) * x6a2_q(its:ite, jts:jte, kts:kte))        &
 598       + sum (grid%x6a%rh(its:ite, jts:jte, kts:kte)* x6a2_rh(its:ite, jts:jte, kts:kte))       &
 599       + sum (grid%x6a%psfc(its:ite, jts:jte) * x6a2_psfc(its:ite, jts:jte))
 600 #endif
 601
 602    if ( cloud_cv_options >= 1 ) then
 603       partial_rhs = partial_rhs &
 604          + sum (grid%xa%qcw(its:ite, jts:jte, kts:kte) * xa2_qcw(its:ite, jts:jte, kts:kte))   &
 605          + sum (grid%xa%qrn(its:ite, jts:jte, kts:kte) * xa2_qrn(its:ite, jts:jte, kts:kte))
 606       if ( cloud_cv_options >= 2 ) then
 607          partial_rhs = partial_rhs &
 608             + sum (grid%xa%qci(its:ite, jts:jte, kts:kte) * xa2_qci(its:ite, jts:jte, kts:kte)) &
 609             + sum (grid%xa%qsn(its:ite, jts:jte, kts:kte) * xa2_qsn(its:ite, jts:jte, kts:kte)) &
 610             + sum (grid%xa%qgr(its:ite, jts:jte, kts:kte) * xa2_qgr(its:ite, jts:jte, kts:kte))
 611       end if
 612    end if
 613 #ifdef VAR4D
 614    partial_rhs = partial_rhs &
 615       + sum (grid%x6a%qcw(its:ite, jts:jte, kts:kte) * x6a2_qcw(its:ite, jts:jte, kts:kte)) &
 616       + sum (grid%x6a%qci(its:ite, jts:jte, kts:kte) * x6a2_qci(its:ite, jts:jte, kts:kte)) &
 617       + sum (grid%x6a%qrn(its:ite, jts:jte, kts:kte) * x6a2_qrn(its:ite, jts:jte, kts:kte)) &
 618       + sum (grid%x6a%qsn(its:ite, jts:jte, kts:kte) * x6a2_qsn(its:ite, jts:jte, kts:kte)) &
 619       + sum (grid%x6a%qgr(its:ite, jts:jte, kts:kte) * x6a2_qgr(its:ite, jts:jte, kts:kte))
 620 #endif
 621
 622 #ifdef A2C
 623     if( ite == ide ) then
 624 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))
 625    partial_rhs = partial_rhs &
 626       + sum (grid%xa%u(ite+1, jts:jte, kts:kte) * xa2_u(ite+1, jts:jte, kts:kte))
 627    end if
 628     if( jte == jde ) then
 629 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))
 630    partial_rhs = partial_rhs &
 631       + sum (grid%xa%v(its:ite, jte+1, kts:kte) * xa2_v(its:ite, jte+1, kts:kte))
 632    end if
 633 #endif
 634
 635    !----------------------------------------------------------------------
 636    !  [7.0] Print output:
 637    !----------------------------------------------------------------------
 638    write (unit=stdout, fmt='(A,1pe22.14)') ' Single Domain < y, y     > = ', pertile_lhs
 639    write (unit=stdout, fmt='(A,1pe22.14)') ' Single Domain < x, x_adj > = ', pertile_rhs
 640
 641    adj_ttl_lhs = wrf_dm_sum_real (partial_lhs)
 642    adj_ttl_rhs = wrf_dm_sum_real (partial_rhs)
 643
 644    if (rootproc) then
 645       write(unit=stdout, fmt='(/)')
 646       write (unit=stdout, fmt='(A,1pe22.14)') ' Whole Domain < y, y     > = ', adj_ttl_lhs
 647       write (unit=stdout, fmt='(A,1pe22.14)') ' Whole Domain < x, x_adj > = ', adj_ttl_rhs
 648    end if
 649
 650    write (unit=stdout, fmt='(/a/)') 'da_check_xtoy_adjoint: Test Finished:'
 651    if (trace_use) call da_trace_exit("da_check_xtoy_adjoint")
 652
 653 end subroutine da_check_xtoy_adjoint
 654
 655