wrfv2_fire/var/da/da_radiance/da_transform_xtoy_rttov_adj.inc

   1 subroutine da_transform_xtoy_rttov_adj ( iv, jo_grad_y, jo_grad_x )
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: transform gradient from obs space to model grid space.
   5    !
   6    ! METHOD:  jo_grad_x = H^T jo_grad_y =  - H^T R^-1 ( d - H delta_x )
   7    !           1. input gradient in obs space and reference state of RTTOV
   8    !           2. call adjoint of RTM
   9    !           3. adjoint of interpolation from model grid to obs loc
  10    !---------------------------------------------------------------------------
  11
  12    implicit none
  13
  14 !#include "rttov_v2q.interface"
  15
  16    type (x_type), intent(inout)   :: jo_grad_x !
  17    type (y_type),  intent(in)     :: jo_grad_y ! H' delta_x
  18    type (iv_type), intent(in)     :: iv        ! O-B structure.
  19
  20 #ifdef RTTOV
  21
  22    integer                        :: i, j, k,c  ! Index dimension.
  23    integer                        :: nlevels ! Number of obs levels.
  24    integer                        :: num_rad  ! Number of radiance obs
  25 #ifdef DM_PARALLEL
  26    integer :: status(mpi_status_size) ! MPI status.
  27 #endif
  28
  29    real, allocatable              :: model_mr(:)
  30
  31
  32    integer            :: nchan, n
  33
  34    type(con_vars_type),  allocatable :: con_vars(:), con_vars_ad(:)
  35    type(aux_vars_type),  allocatable :: aux_vars(:)
  36
  37    ! evenly distributed versions
  38    type(con_vars_type),  allocatable  :: d_con_vars(:), d_con_vars_ad(:)
  39    type(aux_vars_type),  allocatable  :: d_aux_vars(:)
  40    real,    allocatable :: d_tb(:,:)
  41    integer :: d_num_rad, l_num_rad,e,s,len,nl
  42    real, allocatable :: temp_t(:,:),temp_q(:,:), temp_ps(:), rtemp(:,:)
  43    real, allocatable :: temp_tb(:,:)
  44    integer, allocatable :: itemp(:,:)
  45
  46    if (iv%num_inst < 1) return
  47
  48    if (trace_use) call da_trace_entry("da_transform_xtoy_rttov_adj")
  49
  50    do i = 1, iv%num_inst                 ! loop for sensor
  51
  52       if (iv%instid(i)%num_rad_glo < 1) cycle
  53
  54       num_rad  = iv%instid(i)%num_rad
  55       nlevels  = iv%instid(i)%nlevels
  56       nchan    = iv%instid(i)%nchan
  57
  58       if (num_rad > 0) then
  59
  60          allocate (model_mr(1:nlevels))
  61
  62          allocate (con_vars(num_rad))
  63          allocate (con_vars_ad(num_rad))
  64          allocate (aux_vars(num_rad))
  65
  66          !---------------------------------------------------------------
  67          ! [1.0] assign tb = R^-1 Re :
  68          !---------------------------------------------------------------
  69
  70          !---------------------------------------------
  71          ! [2.0] get input of da_rttov_ad
  72          !---------------------------------------------
  73
  74          do n = 1, num_rad
  75             con_vars(n) % nlevels = nlevels
  76
  77             allocate (con_vars(n)    % t(nlevels))
  78             allocate (con_vars(n)    % q(nlevels))
  79             allocate (con_vars_ad(n) % t(nlevels))
  80             allocate (con_vars_ad(n) % q(nlevels))
  81
  82             con_vars(n)    % t(:) = 0.0
  83             con_vars(n)    % q(:) = 0.0
  84             con_vars_ad(n) % t(:) = 0.0
  85             con_vars_ad(n) % q(:) = 0.0
  86
  87             con_vars(n) % t(:) = iv%instid(i)%t(:,n)
  88             con_vars(n) % q(:) = iv%instid(i)%mr(:,n)
  89             con_vars(n) % ps   = iv%instid(i)%ps(n)
  90
  91             aux_vars(n) % t2m = iv%instid(i)%t2m(n)
  92             aux_vars(n) % q2m = iv%instid(i)%mr2m(n)
  93             aux_vars(n) % u10 = iv%instid(i)%u10(n)
  94             aux_vars(n) % v10 = iv%instid(i)%v10(n)
  95
  96             aux_vars(n) % surftype = iv%instid(i)%landsea_mask(n)
  97             aux_vars(n) % surft    = iv%instid(i)%ts(n)
  98             ! aux_vars(n) % fastem(:) = 0.0
  99
 100             aux_vars(n) % satzen  = iv%instid(i)%satzen(n)
 101             aux_vars(n) % satazi  = iv%instid(i)%satazi(n)
 102          end do
 103       end if
 104
 105       !-------------------------------------------------
 106       ! [2.9] Distribute evenly across processors
 107       !-------------------------------------------------
 108
 109       d_num_rad=num_tovs_after(i,myproc+1)
 110
 111       if (d_num_rad > 0) then
 112          allocate (d_con_vars(d_num_rad))
 113          allocate (d_con_vars_ad(d_num_rad))
 114          allocate (d_aux_vars(d_num_rad))
 115          allocate (d_tb(nchan,d_num_rad))
 116          d_tb(:,:) = 0.0
 117          d_con_vars(:) % nlevels = nlevels
 118          do n = 1, d_num_rad
 119             allocate (d_con_vars(n) % t(nlevels))
 120             allocate (d_con_vars(n) % q(nlevels))
 121             d_con_vars(n) % t(:) = 0.0
 122             d_con_vars(n) % q(:) = 0.0
 123             allocate (d_con_vars_ad(n) % t(nlevels))
 124             allocate (d_con_vars_ad(n) % q(nlevels))
 125             d_con_vars_ad(n) % t(:) = 0.0
 126             d_con_vars_ad(n) % q(:) = 0.0
 127          end do
 128
 129          ! Fill up with data that stays local
 130
 131          l_num_rad=min(num_rad,d_num_rad)
 132
 133          if (l_num_rad > 0) then
 134             do n = 1, l_num_rad
 135                d_con_vars(n) % t(:) = con_vars(n) % t(:)
 136                d_con_vars(n) % q(:) = con_vars(n) % q(:)
 137                ! d_aux_vars(n) % fastem(:) = 0.0
 138             end do
 139             d_con_vars(1:l_num_rad) % nlevels  = con_vars(1:l_num_rad) % nlevels
 140             d_con_vars(1:l_num_rad) % ps       = con_vars(1:l_num_rad) % ps
 141             d_aux_vars(1:l_num_rad) % t2m      = aux_vars(1:l_num_rad) % t2m
 142             d_aux_vars(1:l_num_rad) % q2m      = aux_vars(1:l_num_rad) % q2m
 143             d_aux_vars(1:l_num_rad) % u10      = aux_vars(1:l_num_rad) % u10
 144             d_aux_vars(1:l_num_rad) % v10      = aux_vars(1:l_num_rad) % v10
 145             d_aux_vars(1:l_num_rad) % surftype = aux_vars(1:l_num_rad) % surftype
 146             d_aux_vars(1:l_num_rad) % surft    = aux_vars(1:l_num_rad) % surft
 147             d_aux_vars(1:l_num_rad) % satzen   = aux_vars(1:l_num_rad) % satzen
 148             d_aux_vars(1:l_num_rad) % satazi   = aux_vars(1:l_num_rad) % satazi
 149
 150             d_tb(:,1:l_num_rad) = jo_grad_y%instid(i)%tb(:,1:l_num_rad)
 151          end if
 152
 153          ! Get data from elsewhere
 154 #ifdef DM_PARALLEL
 155          do c=1,tovs_copy_count(i)
 156             if (tovs_send_pe(i,c)==myproc) then
 157                s=tovs_send_start(i,c)
 158                len=tovs_send_count(i,c)
 159                e=s+len-1
 160                nl=len*nlevels
 161                allocate(temp_t(nlevels,len))
 162                do j=1,len
 163                   temp_t(:,j)=con_vars(s+j-1) % t(:)
 164                end do
 165                call mpi_send( temp_t,nl, true_mpi_real, tovs_recv_pe(i,c), &
 166                   c*8+1, comm, ierr)
 167
 168                allocate(temp_q(nlevels,len))
 169                do j=1,len
 170                   temp_q(:,j)=con_vars(s+j-1) % q(:)
 171                end do
 172                call mpi_send (temp_q,nl, true_mpi_real, tovs_recv_pe(i,c), &
 173                   c*8+2, comm, ierr)
 174
 175                allocate (temp_tb(nchan,len))
 176                temp_tb(:,:)=jo_grad_y%instid(i)%tb(:,s:e)
 177                call mpi_send (temp_tb,len*nchan,  &
 178                   true_mpi_real, tovs_recv_pe(i,c), c*8+3, comm, ierr)
 179
 180                allocate (rtemp(len,8))
 181                rtemp(:,1)= con_vars(s:e) % ps
 182                rtemp(:,2)= aux_vars(s:e) % t2m
 183                rtemp(:,3)= aux_vars(s:e) % q2m
 184                rtemp(:,4)= aux_vars(s:e) % u10
 185                rtemp(:,5)= aux_vars(s:e) % v10
 186                rtemp(:,6)= aux_vars(s:e) % surft
 187                rtemp(:,7)= aux_vars(s:e) % satzen
 188                rtemp(:,8)= aux_vars(s:e) % satazi
 189                call mpi_send(rtemp,len*8, true_mpi_real, tovs_recv_pe(i,c), &
 190                   c*8+4, comm, ierr)
 191
 192                allocate (itemp(len,2))
 193                itemp(:,1)= con_vars(s:e) % nlevels
 194                itemp(:,2)= aux_vars(s:e) % surftype
 195                call mpi_send(itemp,len*2, mpi_integer, tovs_recv_pe(i,c), &
 196                   c*8+5, comm, ierr)
 197             end if
 198
 199             if (tovs_recv_pe(i,c)==myproc) then
 200                s=tovs_recv_start(i,c)
 201                len=tovs_send_count(i,c)
 202                e=s+len-1
 203                nl=len*nlevels
 204                allocate(temp_t(nlevels,len))
 205                call mpi_recv(temp_t,nl, true_mpi_real, tovs_send_pe(i,c), &
 206                   c*8+1, comm, status, ierr)
 207                do j=1,len
 208                   d_con_vars(s+j-1) % t(:)=temp_t(:,j)
 209                end do
 210
 211                allocate(temp_q(nlevels,len))
 212                call mpi_recv(temp_q,nl, true_mpi_real, tovs_send_pe(i,c), &
 213                   c*8+2, comm, status, ierr)
 214                do j=1,len
 215                   d_con_vars(s+j-1) % q(:)=temp_q(:,j)
 216                end do
 217
 218                allocate (temp_tb(nchan,len))
 219                call mpi_recv(temp_tb,len*nchan, true_mpi_real, &
 220                   tovs_send_pe(i,c), c*8+3, comm, status, ierr)
 221                d_tb(:,s:e)=temp_tb(:,:)
 222
 223                allocate (rtemp(len,8))
 224                call mpi_recv(rtemp,len*8, true_mpi_real, &
 225                   tovs_send_pe(i,c), c*8+4, comm, status, ierr)
 226                d_con_vars(s:e) % ps       = rtemp(:,1)
 227                d_aux_vars(s:e) % t2m      = rtemp(:,2)
 228                d_aux_vars(s:e) % q2m      = rtemp(:,3)
 229                d_aux_vars(s:e) % u10      = rtemp(:,4)
 230                d_aux_vars(s:e) % v10      = rtemp(:,5)
 231                d_aux_vars(s:e) % surft    = rtemp(:,6)
 232                d_aux_vars(s:e) % satzen   = rtemp(:,7)
 233                d_aux_vars(s:e) % satazi   = rtemp(:,8)
 234
 235                allocate (itemp(len,2))
 236                call mpi_recv(itemp,len*2, mpi_integer, &
 237                   tovs_send_pe(i,c), c*8+5, comm, status, ierr)
 238                d_con_vars(s:e) % nlevels  = itemp(:,1)
 239                d_aux_vars(s:e) % surftype = itemp(:,2)
 240
 241             end if
 242
 243             if (tovs_send_pe(i,c)==myproc) then
 244                deallocate (temp_tb)
 245                deallocate (temp_t)
 246                deallocate (temp_q)
 247                deallocate (itemp)
 248                deallocate (rtemp)
 249             end if
 250             if (tovs_recv_pe(i,c)==myproc) then
 251                deallocate (temp_tb)
 252                deallocate (temp_t)
 253                deallocate (temp_q)
 254                deallocate (rtemp)
 255                deallocate (itemp)
 256             end if
 257          end do
 258 #endif
 259
 260          if (d_num_rad > 0) then
 261             if (tovs_batch) then
 262                call da_rttov_ad (i, nchan, d_num_rad, d_con_vars, &
 263                d_aux_vars, d_con_vars_ad, d_tb)
 264             else
 265                !$OMP PARALLEL DO  &
 266                !$OMP PRIVATE ( n )
 267                 do n=1,d_num_rad
 268                    call da_rttov_ad (i, nchan, 1, d_con_vars(n:n), &
 269                    d_aux_vars(n:n), d_con_vars_ad(n:n), d_tb(:,n:n))
 270                end do
 271                !$OMP END PARALLEL DO
 272             end if
 273          end if
 274
 275          ! Transfer data back
 276
 277          if (l_num_rad > 0) then
 278             do n = 1, l_num_rad
 279                con_vars_ad(n) % t(:) = d_con_vars_ad(n) % t(:)
 280                con_vars_ad(n) % q(:) = d_con_vars_ad(n) % q(:)
 281             end do
 282
 283             con_vars_ad(1:l_num_rad) % ps = d_con_vars_ad(1:l_num_rad) % ps
 284          end if
 285       end if ! d_num_rad > 0
 286
 287       ! Return the data to other processors. Note the meaning of send_pe and
 288       ! recv_pe is swapped here
 289
 290       nl=nlevels
 291 #ifdef DM_PARALLEL
 292       do c=1,tovs_copy_count(i)
 293          if (tovs_recv_pe(i,c)==myproc) then
 294             s=tovs_recv_start(i,c)
 295             len=tovs_send_count(i,c)
 296             e=s+len-1
 297             nl=len*nlevels
 298             allocate(temp_t(nlevels,len))
 299             do j=1,len
 300                temp_t(:,j)=d_con_vars_ad(s+j-1) % t(:)
 301             end do
 302             call mpi_send (temp_t,nl, true_mpi_real, tovs_send_pe(i,c), &
 303                c*8+6, comm, ierr)
 304
 305             allocate(temp_q(nlevels,len))
 306             do j=1,len
 307                temp_q(:,j)=d_con_vars_ad(s+j-1) % q(:)
 308             end do
 309             call mpi_send (temp_q,nl, true_mpi_real, tovs_send_pe(i,c), &
 310                c*8+7, comm, ierr)
 311
 312             allocate(rtemp(len,1))
 313             rtemp(:,1) = d_con_vars_ad(s:e) % ps
 314             call mpi_send (rtemp,len, true_mpi_real, tovs_send_pe(i,c), &
 315                c*8+8, comm, ierr)
 316          end if
 317
 318          if (tovs_send_pe(i,c)==myproc) then
 319             s=tovs_send_start(i,c)
 320             len=tovs_send_count(i,c)
 321             e=s+len-1
 322             nl=len*nlevels
 323             allocate(temp_t(nlevels,len))
 324             call mpi_recv (temp_t,nl, true_mpi_real, tovs_recv_pe(i,c), &
 325                c*8+6, comm, status, ierr)
 326             do j=1,len
 327                con_vars_ad(s+j-1) % t(:)=temp_t(:,j)
 328             end do
 329
 330             allocate(temp_q(nlevels,len))
 331             call mpi_recv (temp_q,nl, true_mpi_real, tovs_recv_pe(i,c), &
 332                c*8+7, comm, status, ierr)
 333             do j=1,len
 334                con_vars_ad(s+j-1) % q(:)=temp_q(:,j)
 335             end do
 336
 337             allocate(rtemp(len,1))
 338             call mpi_recv (rtemp,len, true_mpi_real, tovs_recv_pe(i,c), &
 339                c*8+8, comm, status, ierr)
 340             con_vars_ad(s:e) % ps=rtemp(:,1)
 341          end if
 342
 343          if (tovs_recv_pe(i,c)==myproc) then
 344             deallocate (temp_t)
 345             deallocate (temp_q)
 346             deallocate (rtemp)
 347          end if
 348          if (tovs_send_pe(i,c)==myproc) then
 349             deallocate (temp_t)
 350             deallocate (temp_q)
 351             deallocate (rtemp)
 352          end if
 353       end do
 354 #endif
 355       if (d_num_rad > 0) then
 356          do n=1,d_num_rad
 357             deallocate (d_con_vars(n) % t)
 358             deallocate (d_con_vars(n) % q)
 359             deallocate (d_con_vars_ad(n) % t)
 360             deallocate (d_con_vars_ad(n) % q)
 361          end do
 362
 363          deallocate (d_tb)
 364          deallocate (d_con_vars)
 365          deallocate (d_aux_vars)
 366          deallocate (d_con_vars_ad)
 367       end if
 368
 369       ! adjoint of convert to hPa
 370
 371       if (num_rad > 0) then
 372          con_vars_ad(:)% ps = con_vars_ad(:)%ps * 0.01
 373       end if
 374
 375       do n=1,num_rad
 376          ! 4.2 scale transform
 377
 378          do k=1, nlevels
 379             model_mr(k) = con_vars_ad(n) % q(k)
 380
 381             if (iv%instid(i)%info%zk(k,n) <= 0.0) then
 382                con_vars_ad(n)%t(k)  = 0.0  !coefs(i) % ref_prfl_t(k,gas_id_watervapour)
 383                con_vars_ad(n)%q(k)  = 0.0  !coefs(i) % ref_prfl_mr(k,gas_id_watervapour)
 384             else
 385                ! adjoint of q(kg/kg) to ppmv
 386
 387                con_vars_ad(n)%q(k) = model_mr(k) * q2ppmv
 388                ! call rttov_v2q (&
 389                !    gas_unit_ppmv,      &! in
 390                !    model_mr(k),        &! in
 391                !    gas_id_watervapour, &! in
 392                !    model_mr(k),        &! in
 393                !    con_vars_ad(n)%q(k)     )    ! inout
 394             end if
 395          end do
 396       end do
 397
 398       if (num_rad > 0) then
 399          allocate(temp_t(nlevels,num_rad))
 400          do n=1,num_rad
 401             temp_t(:,n) = con_vars_ad(n)% t(:)
 402          end do
 403          call da_interp_lin_3d_adj (jo_grad_x%t, iv%instid(i)%info, temp_t)
 404          deallocate(temp_t)
 405
 406          allocate(temp_q(nlevels,num_rad))
 407          do n=1,num_rad
 408             temp_q(:,n) = con_vars_ad(n)% q(:)
 409          end do
 410          call da_interp_lin_3d_adj (jo_grad_x%q, iv%instid(i)%info, temp_q)
 411          deallocate(temp_q)
 412
 413          allocate(temp_ps(num_rad))
 414          do n=1,num_rad
 415             temp_ps(n) = con_vars_ad(n)% ps
 416          end do
 417          call da_interp_lin_2d_adj (jo_grad_x% psfc, iv%instid(i)%info, 1, temp_ps)
 418          deallocate(temp_ps)
 419       end if
 420
 421       if (num_rad > 0) then
 422          do n = 1, num_rad
 423             deallocate (con_vars(n)    % t)
 424             deallocate (con_vars(n)    % q)
 425             deallocate (con_vars_ad(n) % t)
 426             deallocate (con_vars_ad(n) % q)
 427          end do
 428
 429          deallocate (con_vars)
 430          deallocate (aux_vars)
 431          deallocate (con_vars_ad)
 432          deallocate (model_mr)
 433       end if
 434    end do        ! end loop for sensor
 435
 436    if (trace_use) call da_trace_exit("da_transform_xtoy_rttov_adj")
 437 #else
 438     call da_error(__FILE__,__LINE__, &
 439        (/"Must compile with $RTTOV option for radiances"/))
 440 #endif
 441
 442 end subroutine da_transform_xtoy_rttov_adj
 443