var/da/da_transfer_model/da_transfer_xatowrf_nmm_regional.inc

   1 subroutine da_transfer_xatowrf_nmm_regional(grid)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Convert analysis increments into WRF-NMM increments
   5    !           Writes analysis increments
   6    !
   7    !  Author :  Syed RH Rizvi,     MMM/NCAR
   8    !            06/07/2008
   9    !---------------------------------------------------------------------------
  10
  11    implicit none
  12
  13    type(domain), intent(inout)        :: grid
  14
  15    integer :: i, j, k
  16
  17    real    :: sdmd, s1md
  18
  19    ! arrays to hold wrf increments on the c-grid
  20
  21    real, dimension(ims:ime,jms:jme, &
  22       kms:kme) :: &
  23       u_cgrid, v_cgrid, q_cgrid, ph_cgrid
  24
  25    real, dimension(ims:ime,jms:jme) :: mu_cgrid
  26
  27
  28    if (trace_use) call da_trace_entry("da_transfer_xatowrf_nmm_regional")
  29
  30
  31    ! Adjust grid%xa%q to makte grid%xb%q + grid%xa%q > 0.0
  32
  33    if (check_rh == check_rh_tpw) then
  34       call da_check_rh(grid)
  35    else if (check_rh == check_rh_simple) then
  36       call da_check_rh_simple(grid)
  37    end if
  38
  39
  40    if (print_detail_xa) then
  41       write(unit=stdout, fmt='(a, e24.12)') &
  42          'sum(abs(grid%xa%u(its:ite,jts:jte,kts:kte)))=', &
  43          sum(abs(grid%xa%u(its:ite,jts:jte,kts:kte))), &
  44          'sum(abs(grid%xa%v(its:ite,jts:jte,kts:kte)))=', &
  45          sum(abs(grid%xa%v(its:ite,jts:jte,kts:kte))), &
  46          'sum(abs(grid%xa%t(its:ite,jts:jte,kts:kte)))=', &
  47          sum(abs(grid%xa%t(its:ite,jts:jte,kts:kte))), &
  48          'sum(abs(grid%xa%p(its:ite,jts:jte,kts:kte)))=', &
  49          sum(abs(grid%xa%p(its:ite,jts:jte,kts:kte))), &
  50          'sum(abs(grid%xb%t(its:ite,jts:jte,kts:kte)))=', &
  51          sum(abs(grid%xb%t(its:ite,jts:jte,kts:kte))), &
  52          'sum(abs(grid%xb%p(its:ite,jts:jte,kts:kte)))=', &
  53          sum(abs(grid%xb%p(its:ite,jts:jte,kts:kte))), &
  54          'sum(abs(grid%t_2 (its:ite,jts:jte,kts:kte)))=', &
  55          sum(abs(grid%t_2 (its:ite,jts:jte,kts:kte)))
  56
  57
  58        write(unit=stdout, fmt='(2(2x, a, e20.12))') &
  59           'maxval(grid%xa%u(its:ite,jts:jte,kts:kte))=', &
  60           maxval(grid%xa%u(its:ite,jts:jte,kts:kte)), &
  61           'minval(grid%xa%u(its:ite,jts:jte,kts:kte))=', &
  62           minval(grid%xa%u(its:ite,jts:jte,kts:kte)), &
  63           'maxval(grid%xa%v(its:ite,jts:jte,kts:kte))=', &
  64           maxval(grid%xa%v(its:ite,jts:jte,kts:kte)), &
  65           'minval(grid%xa%v(its:ite,jts:jte,kts:kte))=', &
  66           minval(grid%xa%v(its:ite,jts:jte,kts:kte)), &
  67           'maxval(grid%xa%t(its:ite,jts:jte,kts:kte))=', &
  68           maxval(grid%xa%t(its:ite,jts:jte,kts:kte)), &
  69           'minval(grid%xa%t(its:ite,jts:jte,kts:kte))=', &
  70           minval(grid%xa%t(its:ite,jts:jte,kts:kte)), &
  71           'maxval(grid%xa%q(its:ite,jts:jte,kts:kte))=', &
  72           maxval(grid%xa%q(its:ite,jts:jte,kts:kte)), &
  73           'minval(grid%xa%q(its:ite,jts:jte,kts:kte))=', &
  74           minval(grid%xa%q(its:ite,jts:jte,kts:kte)), &
  75           'maxval(grid%xa%p(its:ite,jts:jte,kts:kte))=', &
  76           maxval(grid%xa%p(its:ite,jts:jte,kts:kte)), &
  77           'minval(grid%xa%p(its:ite,jts:jte,kts:kte))=', &
  78           minval(grid%xa%p(its:ite,jts:jte,kts:kte)), &
  79           'maxval(grid%xa%psfc(its:ite,jts:jte))   =', &
  80           maxval(grid%xa%psfc(its:ite,jts:jte)), &
  81           'minval(grid%xa%psfc(its:ite,jts:jte))   =', &
  82           minval(grid%xa%psfc(its:ite,jts:jte))
  83    end if
  84
  85
  86    ! CONVERT FROM A-GRID TO C-GRID
  87
  88    !TBH:  NOTE that grid%xp%halo_id3 = HALO_PSICHI_UV_ADJ which its currently defined
  89    !TBH:  in the Regitstry as "dyn_em 24:grid%xa%u,grid%xa%v,grid%xa%psfc".  Clearly it its not
  90    !TBH:  necessary to update halos in grid%xa%psfc here!  Also, 24-point stencil its
  91    !TBH:  too thick, 9-point should suffice.  Apparently, grid%xp%halo_id3 its used
  92    !TBH:  in many places!  Thits needs to be fixed.
  93
  94 #ifdef DM_PARALLEL
  95 #include "HALO_PSICHI_UV_ADJ.inc"
  96 #endif
  97
  98    ! Fill the boundary
  99
 100    ! The southern boundary
 101    if (jts == jds) then
 102
 103
 104
 105       grid%xa%v(its:ite,jts-1,kts:kte)=2.0*grid%xa%v(its:ite,jts  ,kts:kte) &
 106                             -    grid%xa%v(its:ite,jts+1,kts:kte)
 107
 108    end if
 109
 110    ! The northern boundary
 111    if (jte == jde) then
 112
 113       grid%xa%v(its:ite,jte+1,kts:kte)=2.0*grid%xa%v(its:ite,jte  ,kts:kte) &
 114                             -    grid%xa%v(its:ite,jte-1,kts:kte)
 115
 116    end if
 117
 118    ! The western boundary
 119    if (its == ids) then
 120       grid%xa%u(its-1,jts:jte,kts:kte)=2.0*grid%xa%u(its  ,jts:jte,kts:kte) &
 121                             -    grid%xa%u(its+1,jts:jte,kts:kte)
 122    end if
 123
 124    ! The eastern boundary
 125    if (ite == ide) then
 126       grid%xa%u(ite+1,jts:jte,kts:kte)=2.0*grid%xa%u(ite  ,jts:jte,kts:kte) &
 127                             -    grid%xa%u(ite-1,jts:jte,kts:kte)
 128    end if
 129
 130    ! ========================
 131    ! Write out the increment:
 132    ! ========================
 133
 134       call da_write_increments_for_wrf_nmm_regional (grid)
 135
 136    do k=kts,kte
 137       do j=jts,jte+1
 138          do i=its,ite+1
 139             u_cgrid(i,j,k)=0.5*(grid%xa%u(i-1,j  ,k)+grid%xa%u(i,j,k))
 140             v_cgrid(i,j,k)=0.5*(grid%xa%v(i  ,j-1,k)+grid%xa%v(i,j,k))
 141          end do
 142       end do
 143    end do
 144
 145    !------------------------------------------------------------------------
 146    ! Zero out the unused info
 147    !------------------------------------------------------------------------
 148
 149    ! The northern boundary
 150    if (jte == jde) u_cgrid(its:ite+1,jte+1,kts:kte)=0.0
 151
 152    ! The eastern boundary
 153    if (ite == ide) v_cgrid(ite+1,jts:jte+1,kts:kte)=0.0
 154
 155    !---------------------------------------------------------------------------
 156    ! Add increment to the original guess and update "grid"
 157    !---------------------------------------------------------------------------
 158
 159    do j=jts,jte
 160
 161       do k=kts,kte
 162          do i=its,ite
 163             grid%u_2(i,j,k) = grid%u_2(i,j,k) + u_cgrid(i,j,k)
 164             grid%v_2(i,j,k) = grid%v_2(i,j,k) + v_cgrid(i,j,k)
 165             if( k == kts) &
 166             grid%mu_2(i,j)= grid%mu_2(i,j) + grid%xa%psfc(i,j)
 167             grid%w_2(i,j,k) = grid%w_2(i,j,k) + grid%xa%w(i,j,k)
 168             grid%moist(i,j,k,P_QV) = grid%moist(i,j,k,P_QV) + q_cgrid(i,j,k)
 169             ! makte sure qv its positive.
 170             if (num_pseudo ==0 .and. grid%moist(i,j,k,P_QV) < 1.0e-6) &
 171                 grid%moist(i,j,k,P_QV) = 1.0e-6
 172
 173             if (size(grid%moist,dim=4) >= 4) then
 174                grid%moist(i,j,k,p_qc) = grid%moist(i,j,k,p_qc) + grid%xa%qcw(i,j,k)
 175                grid%moist(i,j,k,p_qr) = grid%moist(i,j,k,p_qr) + grid%xa%qrn(i,j,k)
 176                if (grid%moist(i,j,k,p_qc) < 0.0) grid%moist(i,j,k,p_qc) = 0.0
 177                if (grid%moist(i,j,k,p_qr) < 0.0) grid%moist(i,j,k,p_qr) = 0.0
 178             end if
 179
 180             if (size(grid%moist,dim=4) >= 6) then
 181                grid%moist(i,j,k,p_qi) = grid%moist(i,j,k,p_qi) + grid%xa%qci(i,j,k)
 182                grid%moist(i,j,k,p_qs) = grid%moist(i,j,k,p_qs) + grid%xa%qsn(i,j,k)
 183                if (grid%moist(i,j,k,p_qi) < 0.0) grid%moist(i,j,k,p_qi) = 0.0
 184                if (grid%moist(i,j,k,p_qs) < 0.0) grid%moist(i,j,k,p_qs) = 0.0
 185             end if
 186
 187             if (size(grid%moist,dim=4) >= 7) then
 188                grid%moist(i,j,k,p_qg) = grid%moist(i,j,k,p_qg) + grid%xa%qgr(i,j,k)
 189                if (grid%moist(i,j,k,p_qg) < 0.0) grid%moist(i,j,k,p_qg) = 0.0
 190             end if
 191          end do
 192       end do
 193    end do
 194
 195    ! The northern boundary
 196    if (jte == jde) then
 197       j=jte+1
 198       do k=kts,kte
 199          do i=its,ite
 200             grid%v_2(i,j,k) = grid%v_2(i,j,k) + v_cgrid(i,j,k)
 201          end do
 202       end do
 203    end if
 204
 205    ! The eastern boundary
 206    if (ite == ide) then
 207       i=ite+1
 208       do k=kts,kte
 209          do j=jts,jte
 210             grid%u_2(i,j,k) = grid%u_2(i,j,k) + u_cgrid(i,j,k)
 211          end do
 212       end do
 213    end if
 214
 215    if (print_detail_xa) then
 216       write(unit=stdout, fmt=*) 'simple variables:'
 217
 218       if (ite == ide) then
 219          write (unit=stdout,fmt=*)  ' '
 220
 221          do k=kts+5,kte,10
 222             do j=jts,jte,10
 223                write(unit=stdout, fmt=*) &
 224                     '  grid%u_2(', ide+1, ',', j, ',', k, ')=', &
 225                        grid%u_2(ide+1,j,k)
 226             end do
 227             write(unit=stdout, fmt=*) ' '
 228          end do
 229       end if
 230
 231       if (jte == jde) then
 232          write(unit=stdout, fmt=*) ' '
 233
 234          do k=kts+5,kte,10
 235             do i=its,ite,10
 236                write(unit=stdout, fmt=*) &
 237                     '  grid%v_2(', i, ',', jde+1, ',', k, ')=', &
 238                        grid%v_2(i, jde+1,k)
 239             end do
 240             write(unit=stdout, fmt=*) ' '
 241          end do
 242       end if
 243
 244       write(unit=stdout, fmt='(2(2x, a, e20.12))') &
 245          'maxval(mu_cgrid(its:ite,jts:jte))       =', &
 246          maxval(mu_cgrid(its:ite,jts:jte)), &
 247          'minval(mu_cgrid(its:ite,jts:jte))       =', &
 248          minval(mu_cgrid(its:ite,jts:jte)), &
 249          'maxval(u_cgrid(its:ite,jts:jte,kts:kte))  =', &
 250          maxval(u_cgrid(its:ite,jts:jte,kts:kte)), &
 251          'minval(u_cgrid(its:ite,jts:jte,kts:kte))  =', &
 252          minval(u_cgrid(its:ite,jts:jte,kts:kte)), &
 253          'maxval(v_cgrid(its:ite,jts:jte,kts:kte))  =', &
 254          maxval(v_cgrid(its:ite,jts:jte,kts:kte)), &
 255          'minval(v_cgrid(its:ite,jts:jte,kts:kte))  =', &
 256          minval(v_cgrid(its:ite,jts:jte,kts:kte)), &
 257          'maxval(q_cgrid(its:ite,jts:jte,kts:kte))  =', &
 258          maxval(q_cgrid(its:ite,jts:jte,kts:kte)), &
 259          'minval(q_cgrid(its:ite,jts:jte,kts:kte))  =', &
 260          minval(q_cgrid(its:ite,jts:jte,kts:kte))
 261
 262       do k=kts,kte
 263          write(unit=stdout, fmt='(a, i3)') 'k=', k
 264
 265          write(unit=stdout, fmt='(2(2x, a, e20.12))') &
 266             'maxval(u_cgrid(its:ite,jts:jte,k))  =', maxval(u_cgrid(its:ite,jts:jte,k)), &
 267             'minval(u_cgrid(its:ite,jts:jte,k))  =', minval(u_cgrid(its:ite,jts:jte,k)), &
 268             'maxval(v_cgrid(its:ite,jts:jte,k))  =', maxval(v_cgrid(its:ite,jts:jte,k)), &
 269             'minval(v_cgrid(its:ite,jts:jte,k))  =', minval(v_cgrid(its:ite,jts:jte,k)), &
 270             'maxval(q_cgrid(its:ite,jts:jte,k))  =', maxval(q_cgrid(its:ite,jts:jte,k)), &
 271             'minval(q_cgrid(its:ite,jts:jte,k))  =', minval(q_cgrid(its:ite,jts:jte,k))
 272       end do
 273    end if
 274
 275    if (trace_use) call da_trace_exit("da_transfer_xatowrf_nmm_regional")
 276
 277 end subroutine da_transfer_xatowrf_nmm_regional