var/da/da_transfer_model/da_transfer_xatowrftl_adj.inc

   1 subroutine da_transfer_xatowrftl_adj(grid, config_flags, filnam)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Convert WRFTL variables to analysis increments
   5    !           (inverse of the incremental part of xatowrf)
   6    !---------------------------------------------------------------------------
   7
   8    implicit none
   9
  10    type(domain), intent(inout)               :: grid
  11    type(grid_config_rec_type), intent(inout) :: config_flags
  12
  13    character*4, intent(in) :: filnam
  14
  15    ! Local variables
  16
  17    integer :: i, j, k, bid, ierr, ii, jj, spec_bdy_width
  18    real    :: sdmd, s1md
  19    real    :: rho_cgrid
  20
  21    real, dimension(ims:ime,jms:jme,kms:kme) :: a_press
  22 #ifndef A2C
  23    real, dimension(ims:ime,jms:jme,kms:kme) :: utmp
  24    real, dimension(ims:ime,jms:jme,kms:kme) :: vtmp
  25    real, dimension(ims:ime,jms:jme) :: mut
  26    real, dimension(ims:ime,jms:jme) :: muu
  27    real, dimension(ims:ime,jms:jme) :: muv
  28 #endif
  29
  30    integer ndynopt
  31
  32    if (trace_use) call da_trace_entry("da_transfer_xatowrftl_adj")
  33
  34    !---------------------------------------------------------------------------
  35    ! [7.0] Adjoint of outPUT (inPUT)
  36    !---------------------------------------------------------------------------
  37
  38    if ( .not. adj_sens ) then
  39
  40 #ifdef VAR4D
  41
  42    call kj_swap_reverse (model_grid%a_u_2, grid%a_u_2, &
  43                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
  44    call kj_swap_reverse (model_grid%a_v_2, grid%a_v_2, &
  45                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
  46    call kj_swap_reverse (model_grid%a_w_2, grid%a_w_2, &
  47                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
  48    call kj_swap_reverse (model_grid%a_t_2, grid%a_t_2, &
  49                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
  50    call kj_swap_reverse (model_grid%a_ph_2, grid%a_ph_2, &
  51                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
  52    call kj_swap_reverse (model_grid%a_p, grid%a_p, &
  53                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
  54
  55    grid%a_mu_2 = model_grid%a_mu_2
  56    grid%a_rainnc  = model_grid%a_rainnc
  57    grid%a_rainncv = model_grid%a_rainncv
  58    grid%a_rainc  = model_grid%a_rainc
  59    grid%a_raincv = model_grid%a_raincv
  60
  61    do i = PARAM_FIRST_SCALAR, num_a_moist
  62       call kj_swap_reverse (model_grid%a_moist(:,:,:,i), grid%a_moist(:,:,:,i), &
  63                     grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
  64    enddo
  65
  66    call da_transfer_wrftl_lbc_t0_adj ( grid )
  67
  68    ! output gradient here, without LBC control
  69    if ( .not. trajectory_io .or. var4d_detail_out ) then
  70       config_flags%auxhist7_outname = "gradient_d<domain>_<date>"
  71       call  med_hist_out ( grid , AUXHIST7_ALARM , config_flags )
  72       WRITE(message(1), FMT='(A)') 'Output gradient at 0h (including LBC gradient)'
  73       CALL da_message ( message(1:1) )
  74       config_flags%auxhist7_outname = "auxhist7_d<domain>_<date>"
  75    endif
  76
  77 #endif
  78
  79    else
  80
  81    if ( grid%auxinput17_oid .NE. 0 ) then
  82       call close_dataset ( grid%auxinput17_oid , config_flags , "DATASET=AUXINPUT17" )
  83    endif
  84
  85    call open_r_dataset ( grid%auxinput17_oid, TRIM(filnam) , grid , config_flags, &
  86                               "DATASET=AUXINPUT17", ierr )
  87    if ( ierr .NE. 0 ) then
  88       WRITE( message(1) , * ) 'Error opening ', TRIM( filnam )
  89       CALL wrf_error_fatal( TRIM( message(1) ) )
  90    endif
  91
  92    call wrf_debug (00 , 'Calling input_auxinput17' )
  93    call input_auxinput17 ( grid%auxinput17_oid, grid , config_flags , ierr )
  94
  95    call close_dataset ( grid%auxinput17_oid , config_flags , "DATASET=AUXINPUT17" )
  96
  97    endif
  98
  99    !---------------------------------------------------------------------------
 100    ! [6.0] Adjoint of save OTHERinCREMENT
 101    !---------------------------------------------------------------------------
 102
 103    do k=kts,kte+1
 104       do j=jts,jte
 105          do i=its,ite
 106             grid%xa%w(i,j,k)=grid%xa%w(i,j,k)+grid%a_w_2(i,j,k)
 107          end do
 108       end do
 109    end do
 110
 111    grid%a_w_2 = 0.0
 112
 113 #ifndef A2C
 114    !---------------------------------------------------------------------------
 115    ! [5.0] Adjoint of CONVERT FROM A-GRID TO C-GRID
 116    !---------------------------------------------------------------------------
 117
 118    ! Fill the halo region for a_u and a_v.
 119    utmp=grid%xa%u
 120    vtmp=grid%xa%v
 121 #endif
 122    grid%xa%u=grid%a_u_2
 123    grid%xa%v=grid%a_v_2
 124
 125 #ifdef A2C
 126   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 127        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
 128      ipe = ipe + 1
 129      ide = ide + 1
 130   end if
 131
 132   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 133        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
 134      jpe = jpe + 1
 135      jde = jde + 1
 136   end if
 137 #endif
 138 #ifdef DM_PARALLEL
 139 #include "HALO_XA_A.inc"
 140 #endif
 141
 142 #ifdef A2C
 143   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 144        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
 145      ipe = ipe - 1
 146      ide = ide - 1
 147   end if
 148
 149   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 150        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
 151      jpe = jpe - 1
 152      jde = jde - 1
 153   end if
 154 #else
 155    grid%a_u_2=grid%xa%u
 156    grid%a_v_2=grid%xa%v
 157    grid%xa%u=utmp
 158    grid%xa%v=vtmp
 159    utmp=0.0
 160    vtmp=0.0
 161
 162    do k=kts,kte
 163       do j=jts,jte
 164          do i=its,ite
 165             utmp(i,j,k)=utmp(i,j,k)+0.5*(grid%a_u_2(i+1,j  ,k)+grid%a_u_2(i,j,k))
 166             vtmp(i,j,k)=vtmp(i,j,k)+0.5*(grid%a_v_2(i  ,j+1,k)+grid%a_v_2(i,j,k))
 167          end do
 168       end do
 169    end do
 170
 171    utmp(its-1,jts:jte,kts:kte)=utmp(its-1,jts:jte,kts:kte)+0.5*grid%a_u_2(its,jts:jte,kts:kte)
 172    utmp(ite+1,jts:jte,kts:kte)=utmp(ite+1,jts:jte,kts:kte)+0.5*grid%a_u_2(ite+1,jts:jte,kts:kte)
 173    vtmp(its:ite,jts-1,kts:kte)=vtmp(its:ite,jts-1,kts:kte)+0.5*grid%a_v_2(its:ite,jts,kts:kte)
 174    vtmp(its:ite,jte+1,kts:kte)=vtmp(its:ite,jte+1,kts:kte)+0.5*grid%a_v_2(its:ite,jte+1,kts:kte)
 175
 176    ! The western boundary
 177    if (its == grid%xp%ids) then
 178       grid%xa%u(its  ,jts:jte,kts:kte)=grid%xa%u(its  ,jts:jte,kts:kte)+2.0*utmp(its-1,jts:jte,kts:kte)
 179       grid%xa%u(its+1,jts:jte,kts:kte)=grid%xa%u(its+1,jts:jte,kts:kte)-utmp(its-1,jts:jte,kts:kte)
 180    end if
 181
 182    ! The eastern boundary
 183    if (ite == grid%xp%ide) then
 184       grid%xa%u(ite  ,jts:jte,kts:kte)=grid%xa%u(ite  ,jts:jte,kts:kte)+2.0*utmp(ite+1,jts:jte,kts:kte)
 185       grid%xa%u(ite-1,jts:jte,kts:kte)=grid%xa%u(ite-1,jts:jte,kts:kte)-utmp(ite+1,jts:jte,kts:kte)
 186    end if
 187
 188    grid%xa%u=grid%xa%u+utmp
 189
 190    ! The southern boundary
 191    if (jts == grid%xp%jds) then
 192       grid%xa%v(its:ite,jts  ,kts:kte)=grid%xa%v(its:ite,jts  ,kts:kte)+2.0*vtmp(its:ite,jts-1,kts:kte)
 193       grid%xa%v(its:ite,jts+1,kts:kte)=grid%xa%v(its:ite,jts+1,kts:kte)-vtmp(its:ite,jts-1,kts:kte)
 194    end if
 195
 196    ! The northern boundary
 197    if (jte == grid%xp%jde) then
 198       grid%xa%v(its:ite,jte  ,kts:kte)=grid%xa%v(its:ite,jte  ,kts:kte)+2.0*vtmp(its:ite,jte+1,kts:kte)
 199       grid%xa%v(its:ite,jte-1,kts:kte)=grid%xa%v(its:ite,jte-1,kts:kte)-vtmp(its:ite,jte+1,kts:kte)
 200    end if
 201
 202    grid%xa%v=grid%xa%v+vtmp
 203 #endif
 204
 205    grid%a_u_2 = 0.0
 206    grid%a_v_2 = 0.0
 207
 208    !---------------------------------------------------------------------------
 209    ! [4.0] Adjoint of CONVERT TEMPERATURE inCREMENTS inTO THETA inCREMENTS
 210    !       EVALUATE ALSO THE inCREMENTS OF (1/rho) AND GEOPOTENTIAL
 211    !---------------------------------------------------------------------------
 212
 213    a_press(its:ite,jts:jte,kts:kte+1)=0.0
 214
 215    do k=kte,kts,-1
 216       do j=jts,jte
 217          do i=its,ite
 218             grid%xa%p(i,j,k)= grid%xa%p(i,j,k)+grid%a_p(i,j,k)
 219          end do
 220       end do
 221    end do
 222
 223    grid%a_p = 0.0
 224
 225 !  do k=kte,kts,-1
 226 !     do j=jts,jte
 227 !        do i=its,ite
 228 !           rho_cgrid=-(grid%ph_2(i,j,k+1)-grid%ph_2(i,j,k))*grid%a_ph_2(i,j,k+1)/grid%xb%rho(i,j,k)
 229 !           a_press(i,j,k )=a_press(i,j,k )+grid%a_ph_2(i,j,k+1)/grid%xb%rho(i,j,k)
 230 !           a_press(i,j,k+1)=a_press(i,j,k+1)-grid%a_ph_2(i,j,k+1)/grid%xb%rho(i,j,k)
 231 !           grid%a_ph_2(i,j,k ) =grid%a_ph_2(i,j,k)   +grid%a_ph_2(i,j,k+1)
 232 !           grid%xa%q(i,j,k)=grid%xa%q(i,j,k)-grid%xb%rho(i,j,k)*0.61*rho_cgrid/(1.0+0.61*grid%xb%q(i,j,k))
 233 !           grid%xa%t(i,j,k)=grid%xa%t(i,j,k)-grid%xb%rho(i,j,k)*rho_cgrid/grid%xb%t(i,j,k)
 234 !           grid%xa%p(i,j,k)= grid%xa%p(i,j,k)+grid%xb%rho(i,j,k)*rho_cgrid/grid%xb%p(i,j,k)
 235 !        end do
 236 !     end do
 237 !  end do
 238
 239    do k=kts,kte
 240       do j=jts,jte
 241          do i=its,ite
 242             grid%xa%p(i,j,k)=grid%xa%p(i,j,k)-(t0+grid%t_2(i,j,k))*kappa*grid%a_t_2(i,j,k)/grid%xb%p(i,j,k)
 243             grid%xa%t(i,j,k)=grid%xa%t(i,j,k)+(t0+grid%t_2(i,j,k))*grid%a_t_2(i,j,k)/grid%xb%t(i,j,k)
 244          end do
 245       end do
 246    end do
 247
 248    grid%a_t_2 = 0.0
 249    grid%a_ph_2 = 0.0
 250
 251    !---------------------------------------------------------------------------
 252    ! [3.0] Adjoint of COMPUTE pressure increments (for computing theta increments)
 253    !---------------------------------------------------------------------------
 254
 255    do k=kts,kte
 256       do j=jts,jte
 257          do i=its,ite
 258             a_press(i,j,k+1)=a_press(i,j,k+1)+0.5*grid%xa%p(i,j,k)
 259             a_press(i,j,k )=a_press(i,j,k )+0.5*grid%xa%p(i,j,k)
 260             grid%xa%p(i,j,k)=0.0
 261             grid%a_moist(i,j,k,P_A_QV)=grid%a_moist(i,j,k,P_A_QV)-(grid%mu_2(i,j)+grid%mub(i,j))*a_press(i,j,k)*grid%dn(k)
 262             grid%a_mu_2(i,j)=grid%a_mu_2(i,j)-a_press(i,j,k)*(1.0+grid%moist(i,j,k,P_QV))*grid%dn(k)
 263             a_press(i,j,k+1)=a_press(i,j,k+1)+a_press(i,j,k)
 264          end do
 265       end do
 266    end do
 267
 268    !---------------------------------------------------------------------------
 269    ! [2.0] Adjoint of COMPUTE increments of dry-column air mass per unit area
 270    !---------------------------------------------------------------------------
 271
 272    do j=jts,jte
 273       do i=its,ite
 274          sdmd=0.0
 275          s1md=0.0
 276          do k=kts,kte
 277             s1md=s1md+(1.0+grid%moist(i,j,k,P_QV))*grid%dnw(k)
 278          end do
 279          sdmd=sdmd-grid%xb%psac(i,j)*grid%a_mu_2(i,j)/s1md
 280          grid%xa%psfc(i,j)=grid%xa%psfc(i,j)-grid%a_mu_2(i,j)/s1md
 281          do k=kts,kte
 282             grid%a_moist(i,j,k,P_A_QV)=grid%a_moist(i,j,k,P_A_QV)+sdmd*grid%dnw(k)
 283          end do
 284       end do
 285    end do
 286
 287    grid%a_mu_2 = 0.0
 288    !---------------------------------------------------------------------------
 289    ! [1.0] Adjoint of Get the mixing ratio of moisture
 290    !---------------------------------------------------------------------------
 291    do k=kts,kte
 292       do j=jts,jte
 293          do i=its,ite
 294             grid%xa%q(i,j,k)=grid%xa%q(i,j,k)+grid%a_moist(i,j,k,P_A_QV)/(1.0-grid%xb%q(i,j,k))**2
 295          end do
 296       end do
 297    end do
 298
 299    !if ( config_flags%mp_physics_ad == warmrain_ad ) then
 300       if ( f_a_qc .and. cloud_cv_options >= 1 ) then
 301          grid%xa%qcw(its:ite,jts:jte,kts:kte)=grid%a_moist(its:ite,jts:jte,kts:kte,p_a_qc)
 302       end if
 303       if ( f_a_qr .and. cloud_cv_options >= 1 ) then
 304          grid%xa%qrn(its:ite,jts:jte,kts:kte)=grid%a_moist(its:ite,jts:jte,kts:kte,p_a_qr)
 305       end if
 306       !placeholder
 307       !if ( config_flags%mp_physics_ad == icecld_ad ) then
 308          if ( f_a_qi .and. cloud_cv_options >= 2 ) then
 309             grid%xa%qci(its:ite,jts:jte,kts:kte)=grid%a_moist(its:ite,jts:jte,kts:kte,p_a_qi)
 310          end if
 311          if ( f_a_qs .and. cloud_cv_options >= 2 ) then
 312             grid%xa%qsn(its:ite,jts:jte,kts:kte)=grid%a_moist(its:ite,jts:jte,kts:kte,p_a_qs)
 313          end if
 314          if ( f_a_qg .and. cloud_cv_options >= 2 ) then
 315             grid%xa%qgr(its:ite,jts:jte,kts:kte)=grid%a_moist(its:ite,jts:jte,kts:kte,p_a_qg)
 316          end if
 317       !end if
 318    !end if
 319
 320    grid%a_moist = 0.0
 321
 322    grid%a_rainnc = 0.0
 323    grid%a_rainncv = 0.0
 324    grid%a_rainc = 0.0
 325    grid%a_raincv = 0.0
 326
 327    if (trace_use) call da_trace_exit("da_transfer_xatowrftl_adj")
 328
 329 end subroutine da_transfer_xatowrftl_adj
 330
 331