var/da/da_transfer_model/da_transfer_xatowrftl.inc

   1 subroutine da_transfer_xatowrftl(grid, config_flags, filnam, timestr)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Convert analysis increments into WRFTL increments
   5    !           (following xatowrf, but only keep the increments)
   6    !---------------------------------------------------------------------------
   7
   8    implicit none
   9
  10    type(domain), intent(inout)               :: grid
  11    type(grid_config_rec_type), intent(inout) :: config_flags
  12    character*4, intent(in) :: filnam
  13    character(len=256), intent(in) :: timestr
  14
  15 #ifdef VAR4D
  16
  17    integer :: i, j, k
  18    integer :: is, ie, js, je, ks, ke
  19    real    :: sdmd, s1md
  20    real    :: rho_cgrid
  21
  22 #ifdef A2C
  23    real, allocatable, dimension(:,:,:) :: g_press
  24 #else
  25    real, allocatable, dimension(:,:,:) :: utmp, vtmp, g_press
  26 #endif
  27
  28    integer ndynopt
  29
  30    if (trace_use) call da_trace_entry("da_transfer_xatowrftl")
  31
  32    is=grid%xp%its
  33    ie=grid%xp%ite
  34    js=grid%xp%jts
  35    je=grid%xp%jte
  36    ks=grid%xp%kts
  37    ke=grid%xp%kte
  38
  39    allocate (g_press(grid%xp%ims:grid%xp%ime,grid%xp%jms:grid%xp%jme,grid%xp%kms:grid%xp%kme))
  40
  41    !---------------------------------------------------------------------------
  42    ! [1.0] Get the mixing ratio of moisture first, as it is easy.
  43    !---------------------------------------------------------------------------
  44
  45    do k=ks,ke
  46       do j=js,je
  47          do i=is,ie
  48             grid%g_moist(i,j,k,P_G_QV)=grid%xa%q(i,j,k)/(1.0-grid%xb%q(i,j,k))**2
  49          end do
  50       end do
  51    end do
  52
  53    !---------------------------------------------------------------------------
  54    ! [2.0] COMPUTE increments of dry-column air mass per unit area
  55    !---------------------------------------------------------------------------
  56
  57    do j=js,je
  58       do i=is,ie
  59          sdmd=0.0
  60          s1md=0.0
  61          do k=ks,ke
  62             sdmd=sdmd+grid%g_moist(i,j,k,P_G_QV)*grid%dnw(k)
  63             s1md=s1md+(1.0+grid%moist(i,j,k,P_QV))*grid%dnw(k)
  64          end do
  65          grid%g_mu_2(i,j)=-(grid%xa%psfc(i,j)+grid%xb%psac(i,j)*sdmd)/s1md
  66       end do
  67    end do
  68
  69    !---------------------------------------------------------------------------
  70    ! [3.0] compute pressure increments (for computing theta increments)
  71    !---------------------------------------------------------------------------
  72
  73    do j=js,je
  74       do i=is,ie
  75          g_press(i,j,ke+1)=0.0
  76          do k=ke,ks,-1
  77             g_press(i,j,k)=g_press(i,j,k+1) &
  78                -(grid%g_mu_2(i,j)*(1.0+grid%moist(i,j,k,P_QV)) &
  79                +(grid%mu_2(i,j)+grid%mub(i,j))*grid%g_moist(i,j,k,P_G_QV))* &
  80                grid%dn(k)
  81             grid%xa%p(i,j,k)=0.5*(g_press(i,j,k)+g_press(i,j,k+1))
  82          end do
  83       end do
  84    end do
  85
  86    !---------------------------------------------------------------------------
  87    ! [4.0] convert temperature increments into theta increments
  88    !       evaluate also the increments of (1/rho) and geopotential
  89    !---------------------------------------------------------------------------
  90
  91    do k=ks,ke
  92       do j=js,je
  93          do i=is,ie
  94             grid%g_t_2(i,j,k)=(t0+grid%t_2(i,j,k))* &
  95                (grid%xa%t(i,j,k)/grid%xb%t(i,j,k) &
  96                         -kappa*grid%xa%p(i,j,k)/grid%xb%p(i,j,k))
  97          end do
  98       end do
  99    end do
 100 !  do j=js,je
 101 !     do i=is,ie
 102 !        grid%g_ph_2(i,j,ks)=0.0
 103 !        do k=ks,ke
 104 !           rho_cgrid=grid%xb%rho(i,j,k) &
 105 !                     *(grid%xa%p(i,j,k)/grid%xb%p(i,j,k) &
 106 !                     -grid%xa%t(i,j,k)/grid%xb%t(i,j,k) &
 107 !                     -0.61*grid%xa%q(i,j,k)/(1.0+0.61*grid%xb%q(i,j,k)))
 108 !           grid%g_ph_2(i,j,k+1)=grid%g_ph_2(i,j,k) &
 109 !              -(g_press(i,j,k+1)-g_press(i,j,k) &
 110 !              +(grid%ph_2(i,j,k+1)-grid%ph_2(i,j,k))*rho_cgrid) &
 111 !              /grid%xb%rho(i,j,k)
 112 !        end do
 113 !     end do
 114 !  end do
 115
 116    grid%g_ph_2 = 0.0
 117
 118 !  grid%xa%p = 0.0
 119    grid%g_p = grid%xa%p
 120
 121    deallocate (g_press)
 122
 123    !---------------------------------------------------------------------------
 124    ! [5.0] convert from a-grid to c-grid
 125    !---------------------------------------------------------------------------
 126
 127 #ifdef A2C
 128   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 129        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
 130      ipe = ipe + 1
 131      ide = ide + 1
 132   end if
 133
 134   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 135        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
 136      jpe = jpe + 1
 137      jde = jde + 1
 138   end if
 139
 140 #ifdef DM_PARALLEL
 141 #include "HALO_XA_A.inc"
 142 #endif
 143
 144 !rizvi's update
 145   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 146        fg_format==fg_format_wrf_arw_global  ) .and. ide == ipe ) then
 147      ipe = ipe - 1
 148      ide = ide - 1
 149   end if
 150
 151   if ((fg_format==fg_format_wrf_arw_regional  .or. &
 152        fg_format==fg_format_wrf_arw_global  ) .and. jde == jpe ) then
 153      jpe = jpe - 1
 154      jde = jde - 1
 155   end if
 156 !rizvi's update
 157
 158    grid%g_u_2 = grid%xa%u
 159    grid%g_v_2 = grid%xa%v
 160
 161 #else
 162 #ifdef DM_PARALLEL
 163 #include "HALO_XA_A.inc"
 164
 165    allocate ( utmp(grid%xp%ims:grid%xp%ime,grid%xp%jms:grid%xp%jme, grid%xp%kms:grid%xp%kme) )
 166    allocate ( vtmp(grid%xp%ims:grid%xp%ime,grid%xp%jms:grid%xp%jme, grid%xp%kms:grid%xp%kme) )
 167
 168    utmp = grid%xa%u
 169    vtmp = grid%xa%v
 170
 171    ! The southern boundary (fill A-GRID boundaries)
 172    ! To keep the gradient, A(0) = 2A(1)-A(2)
 173    if (js == grid%xp%jds) then
 174       vtmp(is:ie,js-1,ks:ke)=2.0*grid%xa%v(is:ie,js  ,ks:ke) &
 175                             -    grid%xa%v(is:ie,js+1,ks:ke)
 176    end if
 177
 178    ! The northern boundary
 179    if (je == grid%xp%jde) then
 180       vtmp(is:ie,je+1,ks:ke)=2.0*grid%xa%v(is:ie,je  ,ks:ke) &
 181                             -    grid%xa%v(is:ie,je-1,ks:ke)
 182    end if
 183
 184    ! The western boundary (fill A-GRID boundaries)
 185    ! To keep the gradient, A(0) = 2A(1)-A(2)
 186    if (is == grid%xp%ids) then
 187       utmp(is-1,js:je,ks:ke)=2.0*grid%xa%u(is  ,js:je,ks:ke) &
 188                             -    grid%xa%u(is+1,js:je,ks:ke)
 189    end if
 190
 191    ! The eastern boundary
 192    if (ie == grid%xp%ide) then
 193       utmp(ie+1,js:je,ks:ke)=2.0*grid%xa%u(ie  ,js:je,ks:ke) &
 194                             -    grid%xa%u(ie-1,js:je,ks:ke)
 195    end if
 196
 197    do k=ks,ke
 198       do j=js,je
 199          do i=is,ie+1
 200             grid%g_u_2(i,j,k)=0.5*(utmp(i-1,j  ,k)+utmp(i,j,k))
 201          end do
 202       end do
 203       do j=js,je+1
 204          do i=is,ie
 205             grid%g_v_2(i,j,k)=0.5*(vtmp(i  ,j-1,k)+vtmp(i,j,k))
 206          end do
 207       end do
 208    end do
 209
 210
 211    deallocate (utmp)
 212    deallocate (vtmp)
 213 #else
 214
 215    do k=ks,ke
 216       do j=js,je
 217          do i=is+1,ie
 218             grid%g_u_2(i,j,k)=0.5*(grid%xa%u(i-1,j,k)+grid%xa%u(i,j,k))
 219          end do
 220       end do
 221       do j=js+1,je
 222          do i=is,ie
 223             grid%g_v_2(i,j,k)=0.5*(grid%xa%v(i,j-1,k)+grid%xa%v(i,j,k))
 224          end do
 225       end do
 226    end do
 227
 228    ! To keep the gradient, A(N+1) = 2A(N)-A(N-1)
 229    ! and on C-Grid, this will lead to C(N+1)=(A(N)+A(N+1))/2=(3A(N)-A(N-1))/2
 230
 231    ! The eastern boundary
 232    grid%g_u_2(ie+1,js:je,ks:ke)=(3.0*grid%xa%u(ie,js:je,ks:ke)-grid%xa%u(ie-1,js:je,ks:ke))/2.0
 233
 234    ! The northern boundary
 235    grid%g_v_2(is:ie,je+1,ks:ke)=(3.0*grid%xa%v(is:ie,je,ks:ke)-grid%xa%v(is:ie,je-1,ks:ke))/2.0
 236
 237    ! To keep the gradient, A(0) = 2A(1)-A(2)
 238    ! and on C-Grid, this will lead to C(1)=(A(0)+A(1))/2=(3A(1)-A(2))/2
 239
 240    ! The western boundary
 241    grid%g_u_2(is,js:je,ks:ke)=(3.0*grid%xa%u(is,js:je,ks:ke)-grid%xa%u(is+1,js:je,ks:ke))/2.0
 242
 243    ! The southern boundary
 244    grid%g_v_2(is:ie,js,ks:ke)=(3.0*grid%xa%v(is:ie,js,ks:ke)-grid%xa%v(is:ie,js+1,ks:ke))/2.0
 245
 246 #endif
 247
 248 #endif
 249    !---------------------------------------------------------------------------
 250    ! [6.0] save OTHERinCREMENT
 251    !---------------------------------------------------------------------------
 252
 253    do j=js,je
 254       do k=ks,ke+1
 255          do i=is,ie
 256             grid%g_w_2(i,j,k)=grid%xa%w(i,j,k)
 257          end do
 258       end do
 259    end do
 260
 261    !if ( config_flags%mp_physics_ad == warmrain_ad ) then
 262       if ( f_g_qc .and. cloud_cv_options >= 1 ) then
 263          grid%g_moist(is:ie,js:je,ks:ke,p_g_qc) = grid%xa%qcw(is:ie,js:je,ks:ke)
 264       end if
 265       if ( f_g_qr .and. cloud_cv_options >= 1 ) then
 266          grid%g_moist(is:ie,js:je,ks:ke,p_g_qr) = grid%xa%qrn(is:ie,js:je,ks:ke)
 267       end if
 268       !placeholder
 269       !if ( config_flags%mp_physics_ad == icecld_ad ) then
 270          if ( f_g_qi .and. cloud_cv_options >= 2 ) then
 271             grid%g_moist(is:ie,js:je,ks:ke,p_g_qi) = grid%xa%qci(is:ie,js:je,ks:ke)
 272          end if
 273          if ( f_g_qs .and. cloud_cv_options >= 2 ) then
 274             grid%g_moist(is:ie,js:je,ks:ke,p_g_qs) = grid%xa%qsn(is:ie,js:je,ks:ke)
 275          end if
 276          if ( f_g_qg .and. cloud_cv_options >= 2 ) then
 277             grid%g_moist(is:ie,js:je,ks:ke,p_g_qg) = grid%xa%qgr(is:ie,js:je,ks:ke)
 278          end if
 279       !end if
 280    !end if
 281
 282    call da_transfer_wrftl_lbc_t0 ( grid )
 283
 284    !---------------------------------------------------------------------------
 285    ! [7.0] output
 286    !---------------------------------------------------------------------------
 287
 288    call kj_swap (grid%g_u_2, model_grid%g_u_2, &
 289                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 290    call kj_swap (grid%g_v_2, model_grid%g_v_2, &
 291                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 292    call kj_swap (grid%g_w_2, model_grid%g_w_2, &
 293                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 294    call kj_swap (grid%g_t_2, model_grid%g_t_2, &
 295                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 296    call kj_swap (grid%g_ph_2, model_grid%g_ph_2, &
 297                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 298    call kj_swap (grid%g_p, model_grid%g_p, &
 299                  grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 300    model_grid%g_mu_2 = grid%g_mu_2
 301    model_grid%g_rainnc  = 0.0
 302    model_grid%g_rainncv = 0.0
 303    model_grid%g_rainc  = 0.0
 304    model_grid%g_raincv = 0.0
 305
 306    do i = PARAM_FIRST_SCALAR, num_g_moist
 307       call kj_swap (grid%g_moist(:,:,:,i), model_grid%g_moist(:,:,:,i), &
 308                     grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, grid%xp%kms, grid%xp%kme)
 309    enddo
 310
 311    if ( .not. trajectory_io .or. var4d_detail_out ) &
 312       call  med_hist_out ( grid , AUXHIST8_ALARM , config_flags )
 313
 314    if (trace_use) call da_trace_exit("da_transfer_xatowrftl")
 315
 316 #endif
 317 end subroutine da_transfer_xatowrftl
 318
 319