var/da/da_main/da_update_firstguess.inc

   1 subroutine da_update_firstguess(grid, out_filename)
   2
   3 !---------------------------------------------------------------------------
   4 !  Purpose: Only replace the fields touched by WRFDA, rather than re-generate
   5 !           whole wrfvar_output from the scratch.
   6 !
   7 !  Update  :   jliu@ucar.edu Apr 1, 2013
   8 !              Minor change for Purpose
   9 !              Added copy file mods instead of copying file outside
  10 !              Requires Fortran 2003 standard compiler
  11 !  Creator :   Junmei Ban, Mar 14, 2013
  12 !
  13 !  The following WRF fields are modified:
  14 !    grid%u_2
  15 !    grid%v_2
  16 !    grid%w_2
  17 !    grid%mu_2
  18 !    grid%ph_2
  19 !    grid%t_2
  20 !    grid%moist
  21 !    grid%p
  22 !    grid%psfc
  23 !    grid%t2, grid%q2, grid%u10, grid%v10, grid%th2
  24 !
  25 ! Note about dry and mosit theta perturbation:
  26 !-----------------------------------------------
  27 ! In the WRF model, grid%t_2 is either moist or dry theta perturbation,
  28 ! depending on the use_theta_m switch. The grid%th_phy_m_t0 is always
  29 ! the dry theta perturbation. The T/THM name is a character string of
  30 ! the input/output field for the grid%t_2 variable.
  31 !
  32 ! Inside of the DA Registry, grid%t_2 is always the dry theta perturbation.
  33 ! We should be careful about referring to the internal field (such as grid%t_2)
  34 ! and the string that is used to name the field in the input/output file.
  35 !---------------------------------------------------------------------------
  36
  37   use module_domain, only : get_ijk_from_grid, program_name
  38   use da_control,    only : use_radarobs, use_rad, crtm_cloud, &
  39                             use_radar_rhv, use_radar_rqv
  40   use module_state_description, only : p_qv, p_qc, p_qr, p_qi, &
  41                                        p_qs, p_qg, &
  42                                        f_qc, f_qr, f_qi, f_qs, f_qg
  43   use module_model_constants, only: R_d, R_v, T0
  44 #if (WRF_CHEM == 1)
  45   use module_domain_type, only : fieldlist
  46   use da_control,    only : use_chemic_surfobs, stdout
  47   use module_state_description, only : PARAM_FIRST_SCALAR, num_chem
  48 #endif
  49
  50   implicit none
  51
  52   INTERFACE
  53     integer(c_int32_t) function copyfile(ifile, ofile) bind(c)
  54       import :: c_int32_t, C_CHAR
  55       CHARACTER(KIND=C_CHAR), DIMENSION(*), intent(in) :: ifile, ofile
  56     END function copyfile
  57   END INTERFACE
  58
  59   include 'netcdf.inc'
  60
  61   type(domain), intent(in)             :: grid
  62   character(*), intent(in), optional   :: out_filename
  63
  64 ! Declare local parameters
  65   character(len=120) :: file_name
  66   character(len=19)  :: DateStr1
  67   character(len=4)   :: staggering=' N/A'
  68   character(len=3)   :: ordering
  69   character(len=80), dimension(3)  ::  dimnames
  70   character(len=80)  :: rmse_var
  71   integer            :: dh1
  72   integer            :: i,j,k
  73   integer            :: ndim1
  74   integer            :: WrfType
  75   integer            :: it, ierr, Status, Status_next_time
  76   integer            :: wrf_real
  77   integer            :: nlon_regional,nlat_regional,nsig_regional
  78   integer            :: ids, ide, jds, jde, kds, kde,    &
  79                         ims, ime, jms, jme, kms, kme,    &
  80                         ips, ipe, jps, jpe, kps, kpe
  81   integer, dimension(4)           :: start_index, end_index1
  82   real, dimension(:), allocatable :: globbuf
  83   real*4,allocatable :: field3(:,:,:),field2(:,:)
  84   real*4,allocatable :: field3u(:,:,:),field3v(:,:,:),field3ph(:,:,:)
  85   character(len=4) ::  fgname
  86   integer :: julyr, julday
  87   real    :: gmt
  88   real*4  :: gmt4
  89   real    :: qvf
  90
  91 #if (WRF_CHEM == 1)
  92   integer :: ic
  93   character(len=200)  :: varname
  94   type( fieldlist ), pointer :: p
  95 #endif
  96
  97   wrf_real=104
  98   end_index1=0
  99
 100   call get_ijk_from_grid (  grid ,                         &
 101                           ids, ide, jds, jde, kds, kde,    &
 102                           ims, ime, jms, jme, kms, kme,    &
 103                           ips, ipe, jps, jpe, kps, kpe    )
 104
 105 !
 106 ! update wrfvar_output file with analysis variables from 3dvar
 107 !
 108   if (present(out_filename)) then
 109      file_name = trim(out_filename)
 110      if (file_name == 'ana02') then
 111         fgname = 'fg02'
 112      else
 113         fgname = 'fg'
 114      endif
 115   else
 116      file_name = 'wrfvar_output'
 117      fgname = 'fg'
 118   endif
 119
 120   if (rootproc) then
 121      ierr = copyfile(trim(fgname)//C_NULL_CHAR, trim(file_name)//C_NULL_CHAR)
 122      if ( ierr /= 0 ) then
 123         write(unit=message(1),fmt='(a)') "Failed to create "//trim(file_name)//" from "//trim(fgname)
 124         call da_error(__FILE__,__LINE__,message(1:1))
 125      endif
 126
 127      call ext_ncd_open_for_update( trim(file_name), 0, 0, "", dh1, Status)
 128      if ( Status /= 0 )then
 129         write(unit=message(1),fmt='(a)') "Failed to open "//trim(file_name)
 130         call da_error(__FILE__,__LINE__,message(1:1))
 131      endif
 132
 133 !-------------  get date info
 134
 135      call ext_ncd_get_next_time(dh1, DateStr1, Status_next_time)
 136      if ( var4d .or. num_fgat_time == 1 )  then ! Don't do it for FGAT
 137         if ( DateStr1 /= start_date )then
 138            ! impossible scenario
 139            ! start_date is set to be equal to file date in da_med_initialdata_input.inc
 140            write(unit=message(1),fmt='(a)') 'date info mismatch '//trim(DateStr1)//" != "//trim(start_date)
 141            call da_error(__FILE__,__LINE__,message(1:1))
 142         endif
 143      endif
 144
 145      ! update analysis time info in global attributes
 146      ! needs to be done before the ext_ncd_write_field calls
 147      if ( var4d .or. num_fgat_time == 1 ) then  ! For 4dvar or 3dvar
 148         call get_julgmt(start_date, julyr, julday, gmt)
 149         CALL ext_ncd_put_dom_ti_char (dh1, 'TITLE', ' OUTPUT FROM '//trim(program_name), ierr)
 150         CALL ext_ncd_put_dom_ti_char (dh1, 'START_DATE', trim(start_date), ierr)
 151         CALL ext_ncd_put_dom_ti_char (dh1, 'SIMULATION_START_DATE', trim(start_date), ierr)
 152      else                                  ! For FGAT
 153         call get_julgmt(DateStr1//'     ', julyr, julday, gmt)
 154         CALL ext_ncd_put_dom_ti_char (dh1, 'TITLE', ' OUTPUT FROM '//trim(program_name), ierr)
 155         CALL ext_ncd_put_dom_ti_char (dh1, 'START_DATE', trim(DateStr1), ierr)
 156         CALL ext_ncd_put_dom_ti_char (dh1, 'SIMULATION_START_DATE', trim(DateStr1), ierr)
 157      endif
 158      gmt4 = real(gmt)
 159      CALL ext_ncd_put_dom_ti_real    (dh1, 'GMT',    gmt4,   1, ierr)
 160      CALL ext_ncd_put_dom_ti_integer (dh1, 'JULYR',  julyr,  1, ierr)
 161      CALL ext_ncd_put_dom_ti_integer (dh1, 'JULDAY', julday, 1, ierr)
 162
 163 !-------------  get grid info
 164      rmse_var='T'
 165      call ext_ncd_get_var_info (dh1,rmse_var,ndim1,ordering,staggering, &
 166                                start_index,end_index1, WrfType, ierr    )
 167      nlon_regional=end_index1(1)
 168      nlat_regional=end_index1(2)
 169      nsig_regional=end_index1(3)
 170
 171      !write(6,*)' nlon,nlat,nsig_regional=',nlon_regional,nlat_regional,nsig_regional
 172      allocate(field2(nlon_regional,nlat_regional),field3(nlon_regional,nlat_regional,nsig_regional))
 173      allocate(field3u(nlon_regional+1,nlat_regional,nsig_regional))
 174      allocate(field3v(nlon_regional,nlat_regional+1,nsig_regional))
 175      allocate(field3ph(nlon_regional,nlat_regional,nsig_regional+1))
 176
 177   end if ! end of rootproc
 178
 179 !
 180 ! update MU
 181 !
 182 #ifdef DM_PARALLEL
 183   if (rootproc) then
 184      ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))
 185   else
 186      ALLOCATE(globbuf(1))
 187   end if
 188   globbuf=0.0
 189   call wrf_patch_to_global_double(grid%mu_2,globbuf,1,' ','xy', &
 190                                   ids, ide-1,          jds, jde-1,         1, 1,  &
 191                                   ims, ime,            jms, jme,           1, 1,  &
 192                                   ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
 193 #endif
 194   if (rootproc) then
 195      do j= 1,nlat_regional
 196         do i= 1,nlon_regional
 197 #ifdef DM_PARALLEL
 198            field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
 199 #else
 200            field2(i,j)=grid%mu_2(i,j)
 201 #endif
 202         end do
 203      end do
 204      rmse_var='MU'
 205      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 206           start_index,end_index1, WrfType, ierr    )
 207 !    write(6,*)' rmse_var=',trim(rmse_var)
 208 !    write(6,*)' ordering=',ordering
 209 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 210 !    write(6,*)' ndim1=',ndim1
 211 !    write(6,*)' staggering=',staggering
 212 !    write(6,*)' start_index=',start_index
 213 !    write(6,*)' end_index1=',end_index1
 214      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),       &
 215                              field2,WRF_REAL,0,0,0,ordering,   &
 216                              staggering, dimnames ,              &
 217                              start_index,end_index1,             & !dom
 218                              start_index,end_index1,             & !mem
 219                              start_index,end_index1,             & !pat
 220                              ierr                               )
 221   end if
 222
 223 #ifdef DM_PARALLEL
 224   deallocate(globbuf)
 225 #endif
 226
 227 !
 228 ! update PSFC
 229 !
 230 #ifdef DM_PARALLEL
 231   if (rootproc) then
 232      ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))    ! global_mu_2(ids:ide-1,jds:jde-1) )
 233   else
 234      ALLOCATE(globbuf(1))
 235   end if
 236   globbuf=0.0
 237   call wrf_patch_to_global_double(grid%psfc,globbuf,1,' ','xy', &
 238                                   ids, ide-1,          jds, jde-1,         1, 1,  &
 239                                   ims, ime,            jms, jme,           1, 1,  &
 240                                   ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
 241 #endif
 242   if (rootproc) then
 243      do j= 1,nlat_regional
 244         do i= 1,nlon_regional
 245 #ifdef DM_PARALLEL
 246            field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
 247 #else
 248            field2(i,j)=grid%psfc(i,j)
 249 #endif
 250         end do
 251      end do
 252      rmse_var='PSFC'
 253      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 254                                start_index,end_index1, WrfType, ierr    )
 255 !    write(6,*)' rmse_var=',trim(rmse_var)
 256 !    write(6,*)' ordering=',ordering
 257 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 258 !    write(6,*)' ndim1=',ndim1
 259 !    write(6,*)' staggering=',staggering
 260 !    write(6,*)' start_index=',start_index
 261 !    write(6,*)' end_index1=',end_index1
 262      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),        &
 263                              field2,WRF_REAL,0,0,0,ordering,    &
 264                              staggering, dimnames ,               &
 265                              start_index,end_index1,              & !dom
 266                              start_index,end_index1,              & !mem
 267                              start_index,end_index1,              & !pat
 268                              ierr                   )
 269   end if
 270
 271 #ifdef DM_PARALLEL
 272   deallocate(globbuf)
 273 #endif
 274
 275 !
 276 ! update T2
 277 !
 278 #ifdef DM_PARALLEL
 279   if (rootproc) then
 280      ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))
 281   else
 282      ALLOCATE(globbuf(1))
 283   end if
 284   globbuf=0.0
 285   call wrf_patch_to_global_double(grid%t2,globbuf,1,' ','xy', &
 286                                   ids, ide-1,          jds, jde-1,         1, 1,  &
 287                                   ims, ime,            jms, jme,           1, 1,  &
 288                                   ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
 289 #endif
 290   if (rootproc) then
 291      do j= 1,nlat_regional
 292         do i= 1,nlon_regional
 293 #ifdef DM_PARALLEL
 294            field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
 295 #else
 296            field2(i,j)=grid%t2(i,j)
 297 #endif
 298         end do
 299      end do
 300      rmse_var='T2'
 301      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 302           start_index,end_index1, WrfType, ierr    )
 303 !    write(6,*)' rmse_var=',trim(rmse_var)
 304 !    write(6,*)' ordering=',ordering
 305 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 306 !    write(6,*)' ndim1=',ndim1
 307 !    write(6,*)' staggering=',staggering
 308 !    write(6,*)' start_index=',start_index
 309 !    write(6,*)' end_index1=',end_index1
 310      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 311                              field2,WRF_REAL,0,0,0,ordering,     &
 312                              staggering, dimnames ,                &
 313                              start_index,end_index1,               & !dom
 314                              start_index,end_index1,               & !mem
 315                              start_index,end_index1,               & !pat
 316                              ierr                                 )
 317   end if
 318
 319 #ifdef DM_PARALLEL
 320   deallocate(globbuf)
 321 #endif
 322
 323 !
 324 ! update TH2
 325 !
 326 #ifdef DM_PARALLEL
 327   if (rootproc) then
 328      ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))
 329   else
 330      ALLOCATE(globbuf(1))
 331   end if
 332   globbuf=0.0
 333   call wrf_patch_to_global_double(grid%th2,globbuf,1,' ','xy', &
 334                                   ids, ide-1,          jds, jde-1,         1, 1,  &
 335                                   ims, ime,            jms, jme,           1, 1,  &
 336                                   ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
 337 #endif
 338   if (rootproc) then
 339      do j= 1,nlat_regional
 340         do i= 1,nlon_regional
 341 #ifdef DM_PARALLEL
 342            field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
 343 #else
 344            field2(i,j)=grid%th2(i,j)
 345 #endif
 346         end do
 347      end do
 348      rmse_var='TH2'
 349      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 350           start_index,end_index1, WrfType, ierr    )
 351 !    write(6,*)' rmse_var=',trim(rmse_var)
 352 !    write(6,*)' ordering=',ordering
 353 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 354 !    write(6,*)' ndim1=',ndim1
 355 !    write(6,*)' staggering=',staggering
 356 !    write(6,*)' start_index=',start_index
 357 !    write(6,*)' end_index1=',end_index1
 358      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 359                              field2,WRF_REAL,0,0,0,ordering,     &
 360                              staggering, dimnames ,                &
 361                              start_index,end_index1,               & !dom
 362                              start_index,end_index1,               & !mem
 363                              start_index,end_index1,               & !pat
 364                              ierr                                 )
 365   end if
 366
 367 #ifdef DM_PARALLEL
 368   deallocate(globbuf)
 369 #endif
 370
 371 !
 372 ! update Q2
 373 !
 374 #ifdef DM_PARALLEL
 375   if (rootproc) then
 376      ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))
 377   else
 378      ALLOCATE(globbuf(1))
 379   end if
 380   globbuf=0.0
 381   call wrf_patch_to_global_double(grid%q2,globbuf,1,' ','xy', &
 382                                   ids, ide-1,          jds, jde-1,         1, 1,  &
 383                                   ims, ime,            jms, jme,           1, 1,  &
 384                                   ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
 385 #endif
 386   if (rootproc) then
 387      do j= 1,nlat_regional
 388         do i= 1,nlon_regional
 389 #ifdef DM_PARALLEL
 390            field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
 391 #else
 392            field2(i,j)=grid%q2(i,j)
 393 #endif
 394         end do
 395      end do
 396      rmse_var='Q2'
 397      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 398           start_index,end_index1, WrfType, ierr    )
 399 !    write(6,*)' rmse_var=',trim(rmse_var)
 400 !    write(6,*)' ordering=',ordering
 401 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 402 !    write(6,*)' ndim1=',ndim1
 403 !    write(6,*)' staggering=',staggering
 404 !    write(6,*)' start_index=',start_index
 405 !    write(6,*)' end_index1=',end_index1
 406      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),        &
 407                              field2,WRF_REAL,0,0,0,ordering,    &
 408                              staggering, dimnames ,               &
 409                              start_index,end_index1,              & !dom
 410                              start_index,end_index1,              & !mem
 411                              start_index,end_index1,              & !pat
 412                              ierr                                 )
 413   end if
 414
 415 #ifdef DM_PARALLEL
 416   deallocate(globbuf)
 417 #endif
 418
 419 !
 420 ! update U10
 421 !
 422 #ifdef DM_PARALLEL
 423   if (rootproc) then
 424      ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))    ! global_mu_2(ids:ide-1,jds:jde-1) )
 425   else
 426      ALLOCATE(globbuf(1))
 427   end if
 428   globbuf=0.0
 429   call wrf_patch_to_global_double(grid%u10,globbuf,1,' ','xy', &
 430                                   ids, ide-1,          jds, jde-1,         1, 1,  &
 431                                   ims, ime,            jms, jme,           1, 1,  &
 432                                   ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
 433 #endif
 434   if (rootproc) then
 435      do j= 1,nlat_regional
 436         do i= 1,nlon_regional
 437 #ifdef DM_PARALLEL
 438            field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
 439 #else
 440            field2(i,j)=grid%u10(i,j)
 441 #endif
 442         end do
 443      end do
 444      rmse_var='U10'
 445      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 446           start_index,end_index1, WrfType, ierr    )
 447 !    write(6,*)' rmse_var=',trim(rmse_var)
 448 !    write(6,*)' ordering=',ordering
 449 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 450 !    write(6,*)' ndim1=',ndim1
 451 !    write(6,*)' staggering=',staggering
 452 !    write(6,*)' start_index=',start_index
 453 !    write(6,*)' end_index1=',end_index1
 454      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 455                              field2,WRF_REAL,0,0,0,ordering,     &
 456                              staggering, dimnames ,                &
 457                              start_index,end_index1,               & !dom
 458                              start_index,end_index1,               & !mem
 459                              start_index,end_index1,               & !pat
 460                              ierr                                 )
 461   end if
 462
 463 #ifdef DM_PARALLEL
 464   deallocate(globbuf)
 465 #endif
 466
 467 !
 468 ! update V10
 469 !
 470 #ifdef DM_PARALLEL
 471   if (rootproc) then
 472      ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))
 473   else
 474      ALLOCATE(globbuf(1))
 475   end if
 476   globbuf=0.0
 477   call wrf_patch_to_global_double(grid%v10,globbuf,1,' ','xy', &
 478                                   ids, ide-1,          jds, jde-1,         1, 1,  &
 479                                   ims, ime,            jms, jme,           1, 1,  &
 480                                   ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
 481 #endif
 482   if (rootproc) then
 483      do j= 1,nlat_regional
 484         do i= 1,nlon_regional
 485 #ifdef DM_PARALLEL
 486            field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
 487 #else
 488            field2(i,j)=grid%v10(i,j)
 489 #endif
 490         end do
 491      end do
 492      rmse_var='V10'
 493      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 494           start_index,end_index1, WrfType, ierr    )
 495 !    write(6,*)' rmse_var=',trim(rmse_var)
 496 !    write(6,*)' ordering=',ordering
 497 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 498 !    write(6,*)' ndim1=',ndim1
 499 !    write(6,*)' staggering=',staggering
 500 !    write(6,*)' start_index=',start_index
 501 !    write(6,*)' end_index1=',end_index1
 502      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 503                              field2,WRF_REAL,0,0,0,ordering,     &
 504                              staggering, dimnames ,                &
 505                              start_index,end_index1,               & !dom
 506                              start_index,end_index1,               & !mem
 507                              start_index,end_index1,               & !pat
 508                              ierr                                 )
 509   end if
 510
 511 #ifdef DM_PARALLEL
 512   deallocate(globbuf)
 513 #endif
 514
 515 !
 516 !  update P
 517 !
 518 #ifdef DM_PARALLEL
 519   if (rootproc) then
 520     allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
 521   else
 522     allocate( globbuf( 1 ) )
 523   endif
 524   globbuf=0.0
 525
 526   CALL wrf_patch_to_global_double ( grid%p, globbuf, 1, '', 'xyz' ,             &
 527                                     ids, ide-1,          jds, jde-1,          kds, kde-1, &
 528                                     ims, ime,            jms, jme,            kms, kme,   &
 529                                     ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1  )
 530 #endif
 531   if (rootproc) then
 532      do k= 1,nsig_regional
 533         do j= 1,nlat_regional
 534            do i= 1,nlon_regional
 535 #ifdef DM_PARALLEL
 536               field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
 537 #else
 538               field3(i,j,k)=grid%p(i,j,k)
 539 #endif
 540            end do
 541         end do
 542      end do
 543      rmse_var='P'
 544      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 545           start_index,end_index1, WrfType, ierr    )
 546 !    write(6,*)' rmse_var=',trim(rmse_var)
 547 !    write(6,*)' ordering=',ordering
 548 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 549 !    write(6,*)' ndim1=',ndim1
 550 !    write(6,*)' staggering=',staggering
 551 !    write(6,*)' start_index=',start_index
 552 !    write(6,*)' end_index1=',end_index1
 553      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 554                              field3,WRF_REAL,0,0,0,ordering,     &
 555                              staggering, dimnames ,                &
 556                              start_index,end_index1,               & !dom
 557                              start_index,end_index1,               & !mem
 558                              start_index,end_index1,               & !pat
 559                              ierr                                 )
 560   end if  ! end of rootproc
 561
 562 #ifdef DM_PARALLEL
 563   deallocate(globbuf)
 564 #endif
 565
 566 !
 567 !  update T
 568 !
 569 #ifdef DM_PARALLEL
 570   if (rootproc) then
 571     allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
 572   else
 573     allocate( globbuf( 1 ) )
 574   endif
 575   globbuf=0.0
 576   CALL wrf_patch_to_global_double ( grid%t_2, globbuf, 1, '', 'xyz' ,                       &
 577                                     ids, ide-1,          jds, jde-1,          kds, kde-1, &
 578                                     ims, ime,            jms, jme,            kms, kme,   &
 579                                     ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1  )
 580 #endif
 581   if (rootproc) then
 582      do k= 1,nsig_regional
 583         do j= 1,nlat_regional
 584            do i= 1,nlon_regional
 585 #ifdef DM_PARALLEL
 586               field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
 587 #else
 588               field3(i,j,k)=grid%t_2(i,j,k)
 589 #endif
 590            end do
 591         end do
 592      end do
 593      rmse_var='T'
 594      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 595           start_index,end_index1, WrfType, ierr    )
 596 !    write(6,*)' rmse_var=',trim(rmse_var)
 597 !    write(6,*)' ordering=',ordering
 598 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 599 !    write(6,*)' ndim1=',ndim1
 600 !    write(6,*)' staggering=',staggering
 601 !    write(6,*)' start_index=',start_index
 602 !    write(6,*)' end_index1=',end_index1
 603      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 604                              field3,WRF_REAL,0,0,0,ordering,       &
 605                              staggering, dimnames ,                &
 606                              start_index,end_index1,               & !dom
 607                              start_index,end_index1,               & !mem
 608                              start_index,end_index1,               & !pat
 609                              ierr                                 )
 610
 611   end if
 612
 613 #ifdef DM_PARALLEL
 614   deallocate(globbuf)
 615 #endif
 616
 617 !
 618 !  collect QVAPOR from patches to a global array
 619 !
 620 #ifdef DM_PARALLEL
 621   if (rootproc) then
 622     allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
 623   else
 624     allocate( globbuf( 1 ) )
 625   endif
 626   globbuf=0.0
 627   CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qv), globbuf, 1, '', 'xyz' ,         &
 628                                     ids, ide-1,          jds, jde-1,          kds, kde-1, &
 629                                     ims, ime,            jms, jme,            kms, kme,   &
 630                                     ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1  )
 631 #endif
 632
 633   if (rootproc) then
 634 !
 635 ! Update THM (moist theta perturbation) before updating QVAPOR
 636 !
 637       if (grid%use_theta_m == 1) then ! convert dry theta perturbation to moist one
 638          write(unit=message(1),fmt='(A, I2)') 'convert T to THM when use_theta_m = ', grid%use_theta_m
 639          call da_message(message(1:1))
 640          do k= 1,nsig_regional
 641             do j= 1,nlat_regional
 642                do i= 1,nlon_regional
 643 #ifdef DM_PARALLEL
 644                  qvf=1.+(R_v/R_d)*globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
 645 #else
 646                  qvf=1.+(R_v/R_d)*grid%moist(i,j,k,p_qv)
 647 #endif
 648                ! field3 here is dry theta perturbation generated earlier
 649                  field3(i,j,k)=(field3(i,j,k)+T0)*qvf - T0
 650                end do
 651             end do
 652          end do
 653       else  ! THM = T when use_theta_m = 0
 654          write(unit=message(1),fmt='(A, I2)') 'THM = T when use_theta_m = ', grid%use_theta_m
 655          call da_message(message(1:1))
 656       end if
 657
 658       rmse_var='THM'
 659       call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 660            start_index,end_index1, WrfType, ierr    )
 661
 662       call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 663                               field3,WRF_REAL,0,0,0,ordering,       &
 664                               staggering, dimnames ,                &
 665                               start_index,end_index1,               & !dom
 666                               start_index,end_index1,               & !mem
 667                               start_index,end_index1,               & !pat
 668                               ierr                                 )
 669
 670
 671 !
 672 ! Update QVAPOR
 673 !
 674      do k= 1,nsig_regional
 675         do j= 1,nlat_regional
 676            do i= 1,nlon_regional
 677 #ifdef DM_PARALLEL
 678               field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
 679 #else
 680               field3(i,j,k)=grid%moist(i,j,k,p_qv)
 681 #endif
 682            end do
 683         end do
 684      end do
 685      rmse_var='QVAPOR'
 686      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 687           start_index,end_index1, WrfType, ierr    )
 688 !    write(6,*)' rmse_var=',trim(rmse_var)
 689 !    write(6,*)' ordering=',ordering
 690 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 691 !    write(6,*)' ndim1=',ndim1
 692 !    write(6,*)' staggering=',staggering
 693 !    write(6,*)' start_index=',start_index
 694 !    write(6,*)' end_index1=',end_index1
 695      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 696                              field3,WRF_REAL,0,0,0,ordering,       &
 697                              staggering, dimnames ,                &
 698                              start_index,end_index1,               & !dom
 699                              start_index,end_index1,               & !mem
 700                              start_index,end_index1,               & !pat
 701                              ierr                                 )
 702   end if  ! end of rootproc
 703
 704 #ifdef DM_PARALLEL
 705   deallocate(globbuf)
 706 #endif
 707
 708 !
 709 !  update PH
 710 !
 711 #ifdef DM_PARALLEL
 712   if (rootproc) then
 713     allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
 714   else
 715     allocate( globbuf( 1 ) )
 716   endif
 717   globbuf=0.0
 718   CALL wrf_patch_to_global_double ( grid%ph_2, globbuf, 1, 'Z', 'xyz' ,                   &
 719                                     ids, ide-1,          jds, jde-1,          kds, kde, &
 720                                     ims, ime,            jms, jme,            kms, kme, &
 721                                     ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe  )
 722 #endif
 723   if (rootproc) then
 724      do k= 1,nsig_regional+1
 725         do j= 1,nlat_regional
 726            do i= 1,nlon_regional
 727 #ifdef DM_PARALLEL
 728               field3ph(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
 729 #else
 730               field3ph(i,j,k)=grid%ph_2(i,j,k)
 731 #endif
 732            end do
 733         end do
 734      end do
 735      rmse_var='PH'
 736      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 737           start_index,end_index1, WrfType, ierr    )
 738 !    write(6,*)' rmse_var=',trim(rmse_var)
 739 !    write(6,*)' ordering=',ordering
 740 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 741 !    write(6,*)' ndim1=',ndim1
 742 !    write(6,*)' staggering=',staggering
 743 !    write(6,*)' start_index=',start_index
 744 !    write(6,*)' end_index1=',end_index1
 745      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 746                              field3ph,WRF_REAL,0,0,0,ordering,     &
 747                              staggering, dimnames ,                &
 748                              start_index,end_index1,               & !dom
 749                              start_index,end_index1,               & !mem
 750                              start_index,end_index1,               & !pat
 751                              ierr                                 )
 752   end if  ! end of rootproc
 753
 754 #ifdef DM_PARALLEL
 755   deallocate(globbuf)
 756 #endif
 757
 758 !
 759 !  update U
 760 !
 761 #ifdef DM_PARALLEL
 762   if (rootproc) then
 763     allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
 764   else
 765     allocate( globbuf( 1 ) )
 766   endif
 767   globbuf=0.0
 768   CALL wrf_patch_to_global_double ( grid%u_2, globbuf, 1, 'X', 'xyz' ,                    &
 769                                     ids, ide,          jds, jde-1,          kds, kde-1, &
 770                                     ims, ime,          jms, jme,            kms, kme,   &
 771                                     ips, min(ipe,ide), jps, min(jpe,jde-1), kps, kpe-1  )
 772 #endif
 773   if (rootproc) then
 774      do k= 1,nsig_regional
 775         do j= 1,nlat_regional
 776            do i= 1,nlon_regional+1
 777 #ifdef DM_PARALLEL
 778               field3u(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
 779 #else
 780               field3u(i,j,k)=grid%u_2(i,j,k)
 781 #endif
 782            end do
 783         end do
 784      end do
 785      rmse_var='U'
 786      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 787           start_index,end_index1, WrfType, ierr    )
 788 !    write(6,*)' rmse_var=',trim(rmse_var)
 789 !    write(6,*)' ordering=',ordering
 790 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 791 !    write(6,*)' ndim1=',ndim1
 792 !    write(6,*)' staggering=',staggering
 793 !    write(6,*)' start_index=',start_index
 794 !    write(6,*)' end_index1=',end_index1
 795      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 796                              field3u,WRF_REAL,0,0,0,ordering,      &
 797                              staggering, dimnames ,                &
 798                              start_index,end_index1,               & !dom
 799                              start_index,end_index1,               & !mem
 800                              start_index,end_index1,               & !pat
 801                              ierr                                 )
 802   end if  ! end of rootproc
 803
 804 #ifdef DM_PARALLEL
 805   deallocate(globbuf)
 806 #endif
 807
 808 !
 809 !  update V
 810 !
 811 #ifdef DM_PARALLEL
 812   if (rootproc) then
 813     allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
 814   else
 815     allocate( globbuf( 1 ) )
 816   endif
 817   globbuf=0.0
 818   CALL wrf_patch_to_global_double ( grid%v_2, globbuf, 1, 'Y', 'xyz' ,                    &
 819                                     ids, ide-1,          jds, jde,          kds, kde-1, &
 820                                     ims, ime,            jms, jme,          kms, kme,   &
 821                                     ips, min(ipe,ide-1), jps, min(jpe,jde), kps, kpe-1  )
 822 #endif
 823   if (rootproc) then
 824      do k= 1,nsig_regional
 825         do j= 1,nlat_regional+1
 826            do i= 1,nlon_regional
 827 #ifdef DM_PARALLEL
 828               field3v(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
 829 #else
 830               field3v(i,j,k)=grid%v_2(i,j,k)
 831 #endif
 832            end do
 833         end do
 834      end do
 835      rmse_var='V'
 836      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 837           start_index,end_index1, WrfType, ierr    )
 838 !    write(6,*)' rmse_var=',trim(rmse_var)
 839 !    write(6,*)' ordering=',ordering
 840 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 841 !    write(6,*)' ndim1=',ndim1
 842 !    write(6,*)' staggering=',staggering
 843 !    write(6,*)' start_index=',start_index
 844 !    write(6,*)' end_index1=',end_index1
 845      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 846                              field3v,WRF_REAL,0,0,0,ordering,      &
 847                              staggering, dimnames ,                &
 848                              start_index,end_index1,               & !dom
 849                              start_index,end_index1,               & !mem
 850                              start_index,end_index1,               & !pat
 851           ierr                                 )
 852   end if  ! end of rootproc
 853
 854 #ifdef DM_PARALLEL
 855   deallocate(globbuf)
 856 #endif
 857
 858 !
 859 !  update W
 860 !
 861 #ifdef DM_PARALLEL
 862   if (rootproc) then
 863     allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
 864   else
 865     allocate( globbuf( 1 ) )
 866   endif
 867   globbuf=0.0
 868   CALL wrf_patch_to_global_double ( grid%w_2, globbuf, 1, 'Z', 'xyz' ,                    &
 869                                     ids, ide-1,          jds, jde-1,          kds, kde, &
 870                                     ims, ime,            jms, jme,            kms, kme, &
 871                                     ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe  )
 872 #endif
 873   if (rootproc) then
 874      do k= 1,nsig_regional+1
 875         do j= 1,nlat_regional
 876            do i= 1,nlon_regional
 877 #ifdef DM_PARALLEL
 878               field3ph(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
 879 #else
 880               field3ph(i,j,k)=grid%w_2(i,j,k)
 881 #endif
 882            end do
 883         end do
 884      end do
 885      rmse_var='W'
 886      call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 887           start_index,end_index1, WrfType, ierr    )
 888 !    write(6,*)' rmse_var=',trim(rmse_var)
 889 !    write(6,*)' ordering=',ordering
 890 !    write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 891 !    write(6,*)' ndim1=',ndim1
 892 !    write(6,*)' staggering=',staggering
 893 !    write(6,*)' start_index=',start_index
 894 !    write(6,*)' end_index1=',end_index1
 895      call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 896                              field3ph,WRF_REAL,0,0,0,ordering,     &
 897                              staggering, dimnames ,                &
 898                              start_index,end_index1,               & !dom
 899                              start_index,end_index1,               & !mem
 900                              start_index,end_index1,               & !pat
 901                              ierr                                  )
 902   end if
 903
 904 #ifdef DM_PARALLEL
 905   deallocate(globbuf)
 906 #endif
 907
 908 !-------------Update QCLOUD, QRAIN, QICE, QSNOW & QGROUP
 909 if ( cloud_cv_options >= 1 ) then ! update qcw and qrn
 910    if ( f_qc ) then
 911 !
 912 !  update QCLOUD
 913 !
 914 #ifdef DM_PARALLEL
 915       if (rootproc) then
 916          allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
 917       else
 918          allocate( globbuf( 1 ) )
 919       endif
 920       globbuf=0.0
 921       CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qc), globbuf, 1, '', 'xyz' ,     &
 922                                     ids, ide-1,          jds, jde-1,          kds, kde-1, &
 923                                     ims, ime,            jms, jme,            kms, kme,   &
 924                                     ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1  )
 925 #endif
 926       if (rootproc) then
 927          do k= 1,nsig_regional
 928             do j= 1,nlat_regional
 929                do i= 1,nlon_regional
 930 #ifdef DM_PARALLEL
 931                   field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
 932 #else
 933                   field3(i,j,k)=grid%moist(i,j,k,p_qc)
 934 #endif
 935                end do
 936             end do
 937          end do
 938          rmse_var='QCLOUD'
 939          call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 940               start_index,end_index1, WrfType, ierr    )
 941 !        write(6,*)' rmse_var=',trim(rmse_var)
 942 !        write(6,*)' ordering=',ordering
 943 !        write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 944 !        write(6,*)' ndim1=',ndim1
 945 !        write(6,*)' staggering=',staggering
 946 !        write(6,*)' start_index=',start_index
 947 !        write(6,*)' end_index1=',end_index1
 948          call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
 949                                  field3,WRF_REAL,0,0,0,ordering,       &
 950                                  staggering, dimnames ,                &
 951                                  start_index,end_index1,               & !dom
 952                                  start_index,end_index1,               & !mem
 953                                  start_index,end_index1,               & !pat
 954                                  ierr                                  )
 955       end if  ! end of rootproc
 956
 957 #ifdef DM_PARALLEL
 958       deallocate(globbuf)
 959 #endif
 960    end if ! f_qc
 961
 962    if ( f_qr ) then
 963 !
 964 !  update QRAIN
 965 !
 966 #ifdef DM_PARALLEL
 967       if (rootproc) then
 968          allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
 969       else
 970          allocate( globbuf( 1 ) )
 971       endif
 972       globbuf=0.0
 973       CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qr), globbuf, 1, '', 'xyz' ,     &
 974                                     ids, ide-1,          jds, jde-1,          kds, kde-1, &
 975                                     ims, ime,            jms, jme,            kms, kme,   &
 976                                     ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1  )
 977 #endif
 978       if (rootproc) then
 979          do k= 1,nsig_regional
 980             do j= 1,nlat_regional
 981                do i= 1,nlon_regional
 982 #ifdef DM_PARALLEL
 983                   field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
 984 #else
 985                   field3(i,j,k)=grid%moist(i,j,k,p_qr)
 986 #endif
 987                end do
 988             end do
 989          end do
 990          rmse_var='QRAIN'
 991          call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
 992               start_index,end_index1, WrfType, ierr    )
 993 !        write(6,*)' rmse_var=',trim(rmse_var)
 994 !        write(6,*)' ordering=',ordering
 995 !        write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
 996 !        write(6,*)' ndim1=',ndim1
 997 !        write(6,*)' staggering=',staggering
 998 !        write(6,*)' start_index=',start_index
 999 !        write(6,*)' end_index1=',end_index1
1000          call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
1001                                  field3,WRF_REAL,0,0,0,ordering,       &
1002                                  staggering, dimnames ,                &
1003                                  start_index,end_index1,               & !dom
1004                                  start_index,end_index1,               & !mem
1005                                  start_index,end_index1,               & !pat
1006                                  ierr                                  )
1007       end if  ! end of rootproc
1008
1009 #ifdef DM_PARALLEL
1010       deallocate(globbuf)
1011 #endif
1012
1013    end if ! f_qr
1014 end if ! cloud_cv_options >= 1
1015
1016 if ( cloud_cv_options >= 2 ) then ! update qci, qsn, qgr
1017    if ( f_qi ) then
1018 !
1019 !  update QICE
1020 !
1021 #ifdef DM_PARALLEL
1022       if (rootproc) then
1023          allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
1024       else
1025          allocate( globbuf( 1 ) )
1026       endif
1027       globbuf=0.0
1028       CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qi), globbuf, 1, '', 'xyz' ,     &
1029                                     ids, ide-1,          jds, jde-1,          kds, kde-1, &
1030                                     ims, ime,            jms, jme,            kms, kme,   &
1031                                     ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1  )
1032 #endif
1033       if (rootproc) then
1034          do k= 1,nsig_regional
1035             do j= 1,nlat_regional
1036                do i= 1,nlon_regional
1037 #ifdef DM_PARALLEL
1038                   field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
1039 #else
1040                   field3(i,j,k)=grid%moist(i,j,k,p_qi)
1041 #endif
1042                end do
1043             end do
1044          end do
1045          rmse_var='QICE'
1046          call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
1047               start_index,end_index1, WrfType, ierr    )
1048 !        write(6,*)' rmse_var=',trim(rmse_var)
1049 !        write(6,*)' ordering=',ordering
1050 !        write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
1051 !        write(6,*)' ndim1=',ndim1
1052 !        write(6,*)' staggering=',staggering
1053 !        write(6,*)' start_index=',start_index
1054 !        write(6,*)' end_index1=',end_index1
1055          call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
1056                                  field3,WRF_REAL,0,0,0,ordering,       &
1057                                  staggering, dimnames ,                &
1058                                  start_index,end_index1,               & !dom
1059                                  start_index,end_index1,               & !mem
1060                                  start_index,end_index1,               & !pat
1061                                  ierr                                 )
1062       end if  ! end of rootproc
1063
1064 #ifdef DM_PARALLEL
1065       deallocate(globbuf)
1066 #endif
1067
1068    end if ! f_qi
1069
1070    if ( f_qs ) then
1071 !
1072 !  update QSNOW
1073 !
1074 #ifdef DM_PARALLEL
1075       if (rootproc) then
1076          allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
1077       else
1078          allocate( globbuf( 1 ) )
1079       endif
1080       globbuf=0.0
1081       CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qs), globbuf, 1, '', 'xyz' ,     &
1082                                     ids, ide-1,          jds, jde-1,          kds, kde-1, &
1083                                     ims, ime,            jms, jme,            kms, kme,   &
1084                                     ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1  )
1085 #endif
1086       if (rootproc) then
1087          do k= 1,nsig_regional
1088             do j= 1,nlat_regional
1089                do i= 1,nlon_regional
1090 #ifdef DM_PARALLEL
1091                   field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
1092 #else
1093                   field3(i,j,k)=grid%moist(i,j,k,p_qs)
1094 #endif
1095                end do
1096             end do
1097          end do
1098          rmse_var='QSNOW'
1099          call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
1100               start_index,end_index1, WrfType, ierr    )
1101 !        write(6,*)' rmse_var=',trim(rmse_var)
1102 !        write(6,*)' ordering=',ordering
1103 !        write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
1104 !        write(6,*)' ndim1=',ndim1
1105 !        write(6,*)' staggering=',staggering
1106 !        write(6,*)' start_index=',start_index
1107 !        write(6,*)' end_index1=',end_index1
1108          call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
1109                                  field3,WRF_REAL,0,0,0,ordering,       &
1110                                  staggering, dimnames ,                &
1111                                  start_index,end_index1,               & !dom
1112                                  start_index,end_index1,               & !mem
1113                                  start_index,end_index1,               & !pat
1114                                  ierr                                 )
1115       end if  ! end of rootproc
1116
1117 #ifdef DM_PARALLEL
1118       deallocate(globbuf)
1119 #endif
1120
1121    end if ! f_qs
1122
1123    if ( f_qg ) then
1124 !
1125 !  update QGRAUP
1126 !
1127 #ifdef DM_PARALLEL
1128       if (rootproc) then
1129          allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
1130       else
1131          allocate( globbuf( 1 ) )
1132       endif
1133       globbuf=0.0
1134       CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qg), globbuf, 1, '', 'xyz' ,     &
1135                                     ids, ide-1,          jds, jde-1,          kds, kde-1, &
1136                                     ims, ime,            jms, jme,            kms, kme,   &
1137                                     ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1  )
1138 #endif
1139       if (rootproc) then
1140          do k= 1,nsig_regional
1141             do j= 1,nlat_regional
1142                do i= 1,nlon_regional
1143 #ifdef DM_PARALLEL
1144                   field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
1145 #else
1146                   field3(i,j,k)=grid%moist(i,j,k,p_qg)
1147 #endif
1148                end do
1149             end do
1150          end do
1151          rmse_var='QGRAUP'
1152          call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
1153               start_index,end_index1, WrfType, ierr    )
1154 !        write(6,*)' rmse_var=',trim(rmse_var)
1155 !        write(6,*)' ordering=',ordering
1156 !        write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
1157 !        write(6,*)' ndim1=',ndim1
1158 !        write(6,*)' staggering=',staggering
1159 !        write(6,*)' start_index=',start_index
1160 !        write(6,*)' end_index1=',end_index1
1161          call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
1162                                  field3,WRF_REAL,0,0,0,ordering,       &
1163                                  staggering, dimnames ,                &
1164                                  start_index,end_index1,               & !dom
1165                                  start_index,end_index1,               & !mem
1166                                  start_index,end_index1,               & !pat
1167                                  ierr                                 )
1168       end if  ! end of rootproc
1169
1170 #ifdef DM_PARALLEL
1171       deallocate(globbuf)
1172 #endif
1173
1174    end if ! f_qg
1175 end if ! cloud_cv_options >= 2
1176 !-------------End of update QCLOUD, QRAIN, QICE, QSNOW & QGROUP
1177
1178 #if (WRF_CHEM == 1)
1179 !The code can be modified if chem IC for gocart are already present in the original wrfinput file,
1180 ! but should not need to be. -Wei Sun 02/2019
1181 if ( use_chemic_surfobs ) then
1182
1183   p => grid%head_statevars
1184   varname = ''
1185   do while (trim(varname) .ne. 'CHEM' .and. ASSOCIATED( p ))
1186      p => p%next
1187      varname = trim( p%DataName )
1188   end do
1189
1190   if ( trim(varname) .eq. 'CHEM' ) then
1191      ndim1 = p%Ndim - 1
1192      ordering = p%MemoryOrder
1193      staggering = p%Stagger
1194      start_index(1) = p%sd1
1195      start_index(2) = p%sd2
1196      start_index(3) = p%sd3
1197      end_index1(1) = p%ed1
1198      end_index1(2) = p%ed2
1199      end_index1(3) = p%ed3
1200
1201      do ic = PARAM_FIRST_SCALAR, num_chem
1202 #ifdef DM_PARALLEL
1203         if (rootproc) then
1204           allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
1205         else
1206           allocate( globbuf( 1 ) )
1207         endif
1208         globbuf=0.0
1209         CALL wrf_patch_to_global_double ( grid%chem(:,:,:,ic), globbuf, 1, '', 'xyz' ,      &
1210                                           ids, ide-1,          jds, jde-1,          kds, kde-1, &
1211                                           ims, ime,            jms, jme,            kms, kme,   &
1212                                           ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1  )
1213 #endif
1214         if (rootproc) then
1215            do k= 1,nsig_regional
1216               do j= 1,nlat_regional
1217                  do i= 1,nlon_regional
1218 #ifdef DM_PARALLEL
1219                     field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
1220 #else
1221                     field3(i,j,k)=grid%chem(i,j,k,ic)
1222 #endif
1223                  end do
1224               end do
1225            end do
1226            rmse_var = TRIM( p%dname_table(grid%id,ic) )
1227            dimnames(1) = TRIM(p%dimname1)
1228            dimnames(2) = TRIM(p%dimname2)
1229            dimnames(3) = TRIM(p%dimname3)
1230 !          write(*,*)' rmse_var=',trim(rmse_var)
1231 !          write(*,*)' nlon, nlat, nsig= ',nlon_regional,nlat_regional,nsig_regional
1232 !          write(*,*)' ordering=',ordering
1233 !          write(*,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
1234 !          write(*,*)' ndim1=',ndim1
1235 !          write(*,*)' staggering=',staggering
1236 !          write(*,*)' start_index=',start_index
1237 !          write(*,*)' end_index1=',end_index1
1238            call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var),         &
1239                                    field3,WRF_REAL,0,0,0,ordering,       &
1240                                    staggering, dimnames ,                &
1241                                    start_index,end_index1,               & !dom
1242                                    start_index,end_index1,               & !mem
1243                                    start_index,end_index1,               & !pat
1244                                    ierr                                 )
1245         end if  ! end of rootproc
1246 #ifdef DM_PARALLEL
1247         deallocate(globbuf)
1248 #endif
1249      end do
1250   else
1251      write(unit=message(1),fmt='(a)') "Failed to find CHEM in statevar list!"
1252      call da_error(__FILE__,__LINE__,message(1:1))
1253   end if
1254
1255 end if
1256 #endif
1257
1258   if (rootproc) then
1259      deallocate(field2,field3,field3u,field3v,field3ph)
1260      call ext_ncd_ioclose(dh1, Status)
1261   end if   ! end of rootproc
1262
1263 end subroutine da_update_firstguess