subroutine da_update_firstguess(grid, out_filename) 4,34
!---------------------------------------------------------------------------
! Purpose: Only replace the fields touched by WRFDA, rather than re-generate
! whole wrfvar_output from the scratch.
!
! Update : jliu@ucar.edu Apr 1, 2013
! Minor change for Purpose
! Added copy file mods instead of copying file outside
! Requires Fortran 2003 standard compiler
! Creator : Junmei Ban, Mar 14, 2013
!
! The following WRF fields are modified:
! grid%u_2
! grid%v_2
! grid%w_2
! grid%mu_2
! grid%ph_2
! grid%t_2
! grid%moist
! grid%p
! grid%psfc
! grid%t2, grid%q2, grid%u10, grid%v10, grid%th2
!
! Note about dry and mosit theta perturbation:
!-----------------------------------------------
! In the WRF model, grid%t_2 is either moist or dry theta perturbation,
! depending on the use_theta_m switch. The grid%th_phy_m_t0 is always
! the dry theta perturbation. The T/THM name is a character string of
! the input/output field for the grid%t_2 variable.
!
! Inside of the DA Registry, grid%t_2 is always the dry theta perturbation.
! We should be careful about referring to the internal field (such as grid%t_2)
! and the string that is used to name the field in the input/output file.
!---------------------------------------------------------------------------
use module_domain
, only : get_ijk_from_grid, program_name
use da_control, only : use_radarobs, use_rad, crtm_cloud, &
use_radar_rhv, use_radar_rqv
use module_state_description, only : p_qv, p_qc, p_qr, p_qi, &
p_qs, p_qg, &
f_qc, f_qr, f_qi, f_qs, f_qg
use module_model_constants, only: R_d, R_v, T0
#if (WRF_CHEM == 1)
use module_domain_type, only : fieldlist
use da_control, only : use_chemic_surfobs, stdout
use module_state_description, only : PARAM_FIRST_SCALAR, num_chem
#endif
implicit none
INTERFACE
integer(c_int32_t) function copyfile(ifile, ofile) bind(c)
import :: c_int32_t, C_CHAR
CHARACTER(KIND=C_CHAR), DIMENSION(*), intent(in) :: ifile, ofile
END function copyfile
END INTERFACE
include 'netcdf.inc'
type(domain), intent(in) :: grid
character(*), intent(in), optional :: out_filename
! Declare local parameters
character(len=120) :: file_name
character(len=19) :: DateStr1
character(len=4) :: staggering=' N/A'
character(len=3) :: ordering
character(len=80), dimension(3) :: dimnames
character(len=80) :: rmse_var
integer :: dh1
integer :: i,j,k
integer :: ndim1
integer :: WrfType
integer :: it, ierr, Status, Status_next_time
integer :: wrf_real
integer :: nlon_regional,nlat_regional,nsig_regional
integer :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe
integer, dimension(4) :: start_index, end_index1
real, dimension(:), allocatable :: globbuf
real*4,allocatable :: field3(:,:,:),field2(:,:)
real*4,allocatable :: field3u(:,:,:),field3v(:,:,:),field3ph(:,:,:)
character(len=4) :: fgname
integer :: julyr, julday
real :: gmt
real*4 :: gmt4
real :: qvf
#if (WRF_CHEM == 1)
integer :: ic
character(len=200) :: varname
type( fieldlist ), pointer :: p
#endif
wrf_real=104
end_index1=0
call get_ijk_from_grid ( grid , &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe )
!
! update wrfvar_output file with analysis variables from 3dvar
!
if (present(out_filename)) then
file_name = trim(out_filename)
if (file_name == 'ana02') then
fgname = 'fg02'
else
fgname = 'fg'
endif
else
file_name = 'wrfvar_output'
fgname = 'fg'
endif
if (rootproc) then
ierr = copyfile(trim(fgname)//C_NULL_CHAR, trim(file_name)//C_NULL_CHAR)
if ( ierr /= 0 ) then
write(unit=message(1),fmt='(a)') "Failed to create "//trim(file_name)//" from "//trim(fgname)
call da_error(__FILE__,__LINE__,message(1:1))
endif
call ext_ncd_open_for_update( trim(file_name), 0, 0, "", dh1, Status)
if ( Status /= 0 )then
write(unit=message(1),fmt='(a)') "Failed to open "//trim(file_name)
call da_error(__FILE__,__LINE__,message(1:1))
endif
!------------- get date info
call ext_ncd_get_next_time(dh1, DateStr1, Status_next_time)
if ( var4d .or. num_fgat_time == 1 ) then ! Don't do it for FGAT
if ( DateStr1 /= start_date )then
! impossible scenario
! start_date is set to be equal to file date in da_med_initialdata_input.inc
write(unit=message(1),fmt='(a)') 'date info mismatch '//trim(DateStr1)//" != "//trim(start_date)
call da_error(__FILE__,__LINE__,message(1:1))
endif
endif
! update analysis time info in global attributes
! needs to be done before the ext_ncd_write_field calls
if ( var4d .or. num_fgat_time == 1 ) then ! For 4dvar or 3dvar
call get_julgmt(start_date, julyr, julday, gmt)
CALL ext_ncd_put_dom_ti_char (dh1, 'TITLE', ' OUTPUT FROM '//trim(program_name), ierr)
CALL ext_ncd_put_dom_ti_char (dh1, 'START_DATE', trim(start_date), ierr)
CALL ext_ncd_put_dom_ti_char (dh1, 'SIMULATION_START_DATE', trim(start_date), ierr)
else ! For FGAT
call get_julgmt(DateStr1//' ', julyr, julday, gmt)
CALL ext_ncd_put_dom_ti_char (dh1, 'TITLE', ' OUTPUT FROM '//trim(program_name), ierr)
CALL ext_ncd_put_dom_ti_char (dh1, 'START_DATE', trim(DateStr1), ierr)
CALL ext_ncd_put_dom_ti_char (dh1, 'SIMULATION_START_DATE', trim(DateStr1), ierr)
endif
gmt4 = real(gmt)
CALL ext_ncd_put_dom_ti_real (dh1, 'GMT', gmt4, 1, ierr)
CALL ext_ncd_put_dom_ti_integer (dh1, 'JULYR', julyr, 1, ierr)
CALL ext_ncd_put_dom_ti_integer (dh1, 'JULDAY', julday, 1, ierr)
!------------- get grid info
rmse_var='T'
call ext_ncd_get_var_info (dh1,rmse_var,ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
nlon_regional=end_index1(1)
nlat_regional=end_index1(2)
nsig_regional=end_index1(3)
!write(6,*)' nlon,nlat,nsig_regional=',nlon_regional,nlat_regional,nsig_regional
allocate(field2(nlon_regional,nlat_regional),field3(nlon_regional,nlat_regional,nsig_regional))
allocate(field3u(nlon_regional+1,nlat_regional,nsig_regional))
allocate(field3v(nlon_regional,nlat_regional+1,nsig_regional))
allocate(field3ph(nlon_regional,nlat_regional,nsig_regional+1))
end if ! end of rootproc
!
! update MU
!
#ifdef DM_PARALLEL
if (rootproc) then
ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))
else
ALLOCATE(globbuf(1))
end if
globbuf=0.0
call wrf_patch_to_global_double(grid%mu_2,globbuf,1,' ','xy', &
ids, ide-1, jds, jde-1, 1, 1, &
ims, ime, jms, jme, 1, 1, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
#endif
if (rootproc) then
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
#else
field2(i,j)=grid%mu_2(i,j)
#endif
end do
end do
rmse_var='MU'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field2,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update PSFC
!
#ifdef DM_PARALLEL
if (rootproc) then
ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3))) ! global_mu_2(ids:ide-1,jds:jde-1) )
else
ALLOCATE(globbuf(1))
end if
globbuf=0.0
call wrf_patch_to_global_double(grid%psfc,globbuf,1,' ','xy', &
ids, ide-1, jds, jde-1, 1, 1, &
ims, ime, jms, jme, 1, 1, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
#endif
if (rootproc) then
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
#else
field2(i,j)=grid%psfc(i,j)
#endif
end do
end do
rmse_var='PSFC'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field2,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update T2
!
#ifdef DM_PARALLEL
if (rootproc) then
ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))
else
ALLOCATE(globbuf(1))
end if
globbuf=0.0
call wrf_patch_to_global_double(grid%t2,globbuf,1,' ','xy', &
ids, ide-1, jds, jde-1, 1, 1, &
ims, ime, jms, jme, 1, 1, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
#endif
if (rootproc) then
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
#else
field2(i,j)=grid%t2(i,j)
#endif
end do
end do
rmse_var='T2'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field2,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update TH2
!
#ifdef DM_PARALLEL
if (rootproc) then
ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))
else
ALLOCATE(globbuf(1))
end if
globbuf=0.0
call wrf_patch_to_global_double(grid%th2,globbuf,1,' ','xy', &
ids, ide-1, jds, jde-1, 1, 1, &
ims, ime, jms, jme, 1, 1, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
#endif
if (rootproc) then
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
#else
field2(i,j)=grid%th2(i,j)
#endif
end do
end do
rmse_var='TH2'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field2,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update Q2
!
#ifdef DM_PARALLEL
if (rootproc) then
ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))
else
ALLOCATE(globbuf(1))
end if
globbuf=0.0
call wrf_patch_to_global_double(grid%q2,globbuf,1,' ','xy', &
ids, ide-1, jds, jde-1, 1, 1, &
ims, ime, jms, jme, 1, 1, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
#endif
if (rootproc) then
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
#else
field2(i,j)=grid%q2(i,j)
#endif
end do
end do
rmse_var='Q2'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field2,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update U10
!
#ifdef DM_PARALLEL
if (rootproc) then
ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3))) ! global_mu_2(ids:ide-1,jds:jde-1) )
else
ALLOCATE(globbuf(1))
end if
globbuf=0.0
call wrf_patch_to_global_double(grid%u10,globbuf,1,' ','xy', &
ids, ide-1, jds, jde-1, 1, 1, &
ims, ime, jms, jme, 1, 1, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
#endif
if (rootproc) then
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
#else
field2(i,j)=grid%u10(i,j)
#endif
end do
end do
rmse_var='U10'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field2,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update V10
!
#ifdef DM_PARALLEL
if (rootproc) then
ALLOCATE(globbuf((ide-1-ids+3)*3*(jde-1-jds+3)))
else
ALLOCATE(globbuf(1))
end if
globbuf=0.0
call wrf_patch_to_global_double(grid%v10,globbuf,1,' ','xy', &
ids, ide-1, jds, jde-1, 1, 1, &
ims, ime, jms, jme, 1, 1, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1),1, 1)
#endif
if (rootproc) then
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field2(i,j)=globbuf(i+(j-1)*(nlon_regional+1))
#else
field2(i,j)=grid%v10(i,j)
#endif
end do
end do
rmse_var='V10'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field2,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update P
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%p, globbuf, 1, '', 'xyz' , &
ids, ide-1, jds, jde-1, kds, kde-1, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1 )
#endif
if (rootproc) then
do k= 1,nsig_regional
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3(i,j,k)=grid%p(i,j,k)
#endif
end do
end do
end do
rmse_var='P'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if ! end of rootproc
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update T
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%t_2, globbuf, 1, '', 'xyz' , &
ids, ide-1, jds, jde-1, kds, kde-1, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1 )
#endif
if (rootproc) then
do k= 1,nsig_regional
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3(i,j,k)=grid%t_2(i,j,k)
#endif
end do
end do
end do
rmse_var='T'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! collect QVAPOR from patches to a global array
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qv), globbuf, 1, '', 'xyz' , &
ids, ide-1, jds, jde-1, kds, kde-1, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1 )
#endif
if (rootproc) then
!
! Update THM (moist theta perturbation) before updating QVAPOR
!
if (grid%use_theta_m == 1) then ! convert dry theta perturbation to moist one
write(unit=message(1),fmt='(A, I2)') 'convert T to THM when use_theta_m = ', grid%use_theta_m
call da_message(message(1:1))
do k= 1,nsig_regional
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
qvf=1.+(R_v/R_d)*globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
qvf=1.+(R_v/R_d)*grid%moist(i,j,k,p_qv)
#endif
! field3 here is dry theta perturbation generated earlier
field3(i,j,k)=(field3(i,j,k)+T0)*qvf - T0
end do
end do
end do
else ! THM = T when use_theta_m = 0
write(unit=message(1),fmt='(A, I2)') 'THM = T when use_theta_m = ', grid%use_theta_m
call da_message(message(1:1))
end if
rmse_var='THM'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
!
! Update QVAPOR
!
do k= 1,nsig_regional
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3(i,j,k)=grid%moist(i,j,k,p_qv)
#endif
end do
end do
end do
rmse_var='QVAPOR'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if ! end of rootproc
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update PH
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%ph_2, globbuf, 1, 'Z', 'xyz' , &
ids, ide-1, jds, jde-1, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe )
#endif
if (rootproc) then
do k= 1,nsig_regional+1
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3ph(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3ph(i,j,k)=grid%ph_2(i,j,k)
#endif
end do
end do
end do
rmse_var='PH'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3ph,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if ! end of rootproc
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update U
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%u_2, globbuf, 1, 'X', 'xyz' , &
ids, ide, jds, jde-1, kds, kde-1, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide), jps, min(jpe,jde-1), kps, kpe-1 )
#endif
if (rootproc) then
do k= 1,nsig_regional
do j= 1,nlat_regional
do i= 1,nlon_regional+1
#ifdef DM_PARALLEL
field3u(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3u(i,j,k)=grid%u_2(i,j,k)
#endif
end do
end do
end do
rmse_var='U'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3u,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if ! end of rootproc
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update V
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%v_2, globbuf, 1, 'Y', 'xyz' , &
ids, ide-1, jds, jde, kds, kde-1, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde), kps, kpe-1 )
#endif
if (rootproc) then
do k= 1,nsig_regional
do j= 1,nlat_regional+1
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3v(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3v(i,j,k)=grid%v_2(i,j,k)
#endif
end do
end do
end do
rmse_var='V'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3v,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if ! end of rootproc
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!
! update W
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%w_2, globbuf, 1, 'Z', 'xyz' , &
ids, ide-1, jds, jde-1, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe )
#endif
if (rootproc) then
do k= 1,nsig_regional+1
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3ph(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3ph(i,j,k)=grid%w_2(i,j,k)
#endif
end do
end do
end do
rmse_var='W'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3ph,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
!-------------Update QCLOUD, QRAIN, QICE, QSNOW & QGROUP
if ( cloud_cv_options >= 1 ) then ! update qcw and qrn
if ( f_qc ) then
!
! update QCLOUD
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qc), globbuf, 1, '', 'xyz' , &
ids, ide-1, jds, jde-1, kds, kde-1, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1 )
#endif
if (rootproc) then
do k= 1,nsig_regional
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3(i,j,k)=grid%moist(i,j,k,p_qc)
#endif
end do
end do
end do
rmse_var='QCLOUD'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if ! end of rootproc
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
end if ! f_qc
if ( f_qr ) then
!
! update QRAIN
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qr), globbuf, 1, '', 'xyz' , &
ids, ide-1, jds, jde-1, kds, kde-1, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1 )
#endif
if (rootproc) then
do k= 1,nsig_regional
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3(i,j,k)=grid%moist(i,j,k,p_qr)
#endif
end do
end do
end do
rmse_var='QRAIN'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if ! end of rootproc
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
end if ! f_qr
end if ! cloud_cv_options >= 1
if ( cloud_cv_options >= 2 ) then ! update qci, qsn, qgr
if ( f_qi ) then
!
! update QICE
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qi), globbuf, 1, '', 'xyz' , &
ids, ide-1, jds, jde-1, kds, kde-1, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1 )
#endif
if (rootproc) then
do k= 1,nsig_regional
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3(i,j,k)=grid%moist(i,j,k,p_qi)
#endif
end do
end do
end do
rmse_var='QICE'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if ! end of rootproc
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
end if ! f_qi
if ( f_qs ) then
!
! update QSNOW
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qs), globbuf, 1, '', 'xyz' , &
ids, ide-1, jds, jde-1, kds, kde-1, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1 )
#endif
if (rootproc) then
do k= 1,nsig_regional
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3(i,j,k)=grid%moist(i,j,k,p_qs)
#endif
end do
end do
end do
rmse_var='QSNOW'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if ! end of rootproc
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
end if ! f_qs
if ( f_qg ) then
!
! update QGRAUP
!
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%moist(:,:,:,p_qg), globbuf, 1, '', 'xyz' , &
ids, ide-1, jds, jde-1, kds, kde-1, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1 )
#endif
if (rootproc) then
do k= 1,nsig_regional
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3(i,j,k)=grid%moist(i,j,k,p_qg)
#endif
end do
end do
end do
rmse_var='QGRAUP'
call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, &
start_index,end_index1, WrfType, ierr )
! write(6,*)' rmse_var=',trim(rmse_var)
! write(6,*)' ordering=',ordering
! write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(6,*)' ndim1=',ndim1
! write(6,*)' staggering=',staggering
! write(6,*)' start_index=',start_index
! write(6,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if ! end of rootproc
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
end if ! f_qg
end if ! cloud_cv_options >= 2
!-------------End of update QCLOUD, QRAIN, QICE, QSNOW & QGROUP
#if (WRF_CHEM == 1)
!The code can be modified if chem IC for gocart are already present in the original wrfinput file,
! but should not need to be. -Wei Sun 02/2019
if ( use_chemic_surfobs ) then
p => grid%head_statevars
varname = ''
do while (trim(varname) .ne. 'CHEM' .and. ASSOCIATED( p ))
p => p%next
varname = trim( p%DataName )
end do
if ( trim(varname) .eq. 'CHEM' ) then
ndim1 = p%Ndim - 1
ordering = p%MemoryOrder
staggering = p%Stagger
start_index(1) = p%sd1
start_index(2) = p%sd2
start_index(3) = p%sd3
end_index1(1) = p%ed1
end_index1(2) = p%ed2
end_index1(3) = p%ed3
!! do ic = 1, num_chem-1
do ic = PARAM_FIRST_SCALAR, num_chem
!!! do ic = 2,2
#ifdef DM_PARALLEL
if (rootproc) then
allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3)))
else
allocate( globbuf( 1 ) )
endif
globbuf=0.0
CALL wrf_patch_to_global_double ( grid%chem(:,:,:,ic), globbuf, 1, '', 'xyz' , &
ids, ide-1, jds, jde-1, kds, kde-1, &
ims, ime, jms, jme, kms, kme, &
ips, min(ipe,ide-1), jps, min(jpe,jde-1), kps, kpe-1 )
#endif
if (rootproc) then
do k= 1,nsig_regional
do j= 1,nlat_regional
do i= 1,nlon_regional
#ifdef DM_PARALLEL
field3(i,j,k)=globbuf(i+(j-1)*(nlon_regional+1)+(k-1)*(nlon_regional+1)*(nlat_regional+1))
#else
field3(i,j,k)=grid%chem(i,j,k,ic)
#endif
end do
end do
end do
rmse_var = TRIM( p%dname_table(grid%id,ic) )
dimnames(1) = TRIM(p%dimname1)
dimnames(2) = TRIM(p%dimname2)
dimnames(3) = TRIM(p%dimname3)
! write(*,*)' rmse_var=',trim(rmse_var)
! write(*,*)' nlon, nlat, nsig= ',nlon_regional,nlat_regional,nsig_regional
! write(*,*)' ordering=',ordering
! write(*,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL
! write(*,*)' ndim1=',ndim1
! write(*,*)' staggering=',staggering
! write(*,*)' start_index=',start_index
! write(*,*)' end_index1=',end_index1
call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), &
field3,WRF_REAL,0,0,0,ordering, &
staggering, dimnames , &
start_index,end_index1, & !dom
start_index,end_index1, & !mem
start_index,end_index1, & !pat
ierr )
end if ! end of rootproc
#ifdef DM_PARALLEL
deallocate(globbuf)
#endif
end do
else
write(unit=message(1),fmt='(a)') "Failed to find CHEM in statevar list!"
call da_error(__FILE__,__LINE__,message(1:1))
end if
end if
#endif
if (rootproc) then
deallocate(field2,field3,field3u,field3v,field3ph)
call ext_ncd_ioclose(dh1, Status)
end if ! end of rootproc
end subroutine da_update_firstguess