module wrwrfmassa_mod
use abstract_wrwrfmassa_mod
  type, extends(abstract_wrwrfmassa_class) :: wrwrfmassa_class
  contains
    procedure, pass(this) :: wrwrfmassa_binary => wrwrfmassa_binary_wrf  
    procedure, pass(this) :: wrwrfmassa_netcdf => wrwrfmassa_netcdf_wrf  
    procedure, nopass :: update_start_date
    procedure, pass(this) :: generic_sub2grid
    procedure, nopass :: reorder_s 
    procedure, nopass :: contract_ibuf
    procedure, nopass :: transfer_ibuf2jbuf
    procedure, nopass :: move_hg_ibuf
  end type wrwrfmassa_class
contains
  subroutine wrwrfmassa_binary_wrf(this,mype)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    wrwrfmassa              write out wrf MASS restart file
  !   prgmmr: parrish          org: np22                date: 2004-06-23
  !
  ! abstract:  read wrf MASS guess restart interface file, add analysis
  !            increment, and write out wrf MASS analysis restart
  !            interface file.
  !
  ! program history log:
  !   2004-06-23  parrish, document
  !   2004-08-03  treadon - add only to module use, add intent in/out
  !   2004-11-18  parrish, rewrite for mpi-io
  !   2004-12-15  treadon - write analysis to file "wrf_inout"
  !   2005-02-17  todling - ifdef'ed wrf code out
  !   2005-11-29  derber - remove qsat_fix array
  !   2006-02-02  treadon - rename prslk as ges_prslk
  !   2006-02-15  treadon - convert specific humidity to moisture mixing ratio
  !   2006-03-07  treadon - convert virtual temperature to potential temperature
  !   2006-04-06  middlecoff - changed nfcst from 11 to lendian_in
  !   2006-07-28  derber  - include sensible temperature
  !   2006-07-31  kleist - change to use ges_ps instead of lnps
  !   2007-03-13  derber - remove unused qsinv2 from jfunc use list
  !   2007-04-12  parrish - add modifications to allow any combination of ikj or ijk
  !                          grid ordering for input 3D fields
  !   2008-03-31  safford - rm unused uses
  !   2008-12-05  todling - adjustment for dsfct time dimension addition
  !   2010-06-24  hu     - add code to write cloud/hydrometeor analysis fields to "wrf_inout"
  !   2011-04-29  todling - introduce MetGuess and wrf_mass_guess_mod
  !   2012-10-11  parrish - add option to swap bytes immediately after every call to mpi_file_read_at and
  !                           before every call to mpi_file_write_at (to handle cases of big-endian
  !                           file/little-endian machine and vice-versa)
  !   2012-11-26  hu     - add code to write updated soil fields to "wrf_inout"
  !   2013-01-26  parrish - WCOSS debug compile type mismatch error -- 
  !                            change to_native_endianness_i4 to to_native_endianness_r4
  !   2013-10-19  todling - metguess now holds background
  !   2014-03-12  hu     - add code to read ges_q2 (2m Q), 
  !                               Qnr(rain number concentration), 
  !                               and nsoil (number of soil levels)
  !   2015-01-13  ladwig - add code to read Qni and Qnc (cloud ice and water
  !                        number concentration)
  !   2015-01-15  hu     - add i_snowT_check to control temperature adjustment
  !                               over snow  
  !   2017-03-23  Hu      - add code to use hybrid vertical coodinate in WRF MASS
  !                               core
  !
  !   input argument list:
  !     mype     - pe number
  !
  !   output argument list:
  !     no output arguments
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm RS/6000 SP
  !
  !$$$
    use kinds, only: r_kind,r_single,i_long,i_llong,i_kind
    use mpimod, only: mpi_byte,mpi_integer4,mpi_real4,mpi_comm_world,npe,ierror, &
         mpi_offset_kind,mpi_info_null,mpi_mode_rdwr,mpi_status_size
    use guess_grids, only: dsfct,&
         ntguessfc,ntguessig,ifilesig,ges_tsen
    use wrf_mass_guess_mod, only: ges_tten
    use gridmod, only: lon1,lat1,nlat_regional,nlon_regional,&
         nsig,nsig_soil,eta1_ll,pt_ll,itotsub,iglobal,update_regsfc,&
         aeta1_ll,eta2_ll,aeta2_ll
    use constants, only: one,zero_single,rd_over_cp_mass,one_tenth,h300,r10,r100
    use constants, only: soilmoistmin
    use gsi_io, only: lendian_in,verbose
    use rapidrefresh_cldsurf_mod, only: l_hydrometeor_bkio,l_gsd_soilTQ_nudge,&
         i_use_2mq4b
    use wrf_mass_guess_mod, only: destroy_cld_grids
    use gsi_bundlemod, only: GSI_BundleGetPointer
    use gsi_metguess_mod, only: gsi_metguess_get,GSI_MetGuess_Bundle
    use native_endianness, only: byte_swap
    use mpeu_util, only: die
    use read_wrf_mass_guess_mod, only: read_wrf_mass_guess_class
  
    implicit none
  
  ! Declare passed variables
    class(wrwrfmassa_class), intent(inout) :: this 
    integer(i_kind),intent(in   ) :: mype
  
  ! Declare local parameters
    type(read_wrf_mass_guess_class) :: read_wrf
    character(len=*),parameter::myname='wrwrfmassa_binary'
    real(r_kind),parameter:: r225=225.0_r_kind
    type(read_wrf_mass_guess_class) :: wrf_mass_guess
  
  ! Declare local variables
    real(r_single),allocatable::tempa(:,:),tempb(:,:)
    real(r_single),allocatable::temp1(:,:),temp1u(:,:),temp1v(:,:)
    real(r_single),allocatable::all_loc(:,:,:)
    integer(i_kind),allocatable::itemp1(:,:)
    integer(i_kind),allocatable::igtype(:),kdim(:),kord(:)
    integer(kind=mpi_offset_kind),allocatable::offset(:)
    integer(kind=mpi_offset_kind) this_offset,offset_mub,offset_start_date
    integer(i_kind),allocatable::length(:)
    integer(i_kind) this_length,length_mub,igtype_mub,length_start_date
    integer(i_llong) num_swap
    character(6) filename
    character(9) wrfanl
    integer(i_kind) ifld,im,jm,lm,num_mass_fields
    integer(i_kind) num_loc_groups,num_j_groups
    integer(i_kind) i,it,j,k
    integer(i_kind) i_mu,i_t,i_q,i_u,i_v
    integer(i_kind) i_qc,i_qi,i_qr,i_qs,i_qg,i_qnr,i_qni,i_qnc
    integer(i_kind) kqc,kqi,kqr,kqs,kqg,kqnr,kqni,kqnc,i_tt,ktt
    integer(i_kind) i_th2,i_q2,i_soilt1,i_tslb,i_smois,ktslb,ksmois,ksize
    integer(i_kind) i_sst,i_tsk
    real(r_kind) psfc_this,psfc_this_dry
    real(r_kind):: work_prsl,work_prslk
    real(r_kind),dimension(lat1+2,lon1+2):: q_integral,q_integralc4h
    integer(i_llong) n_position
    integer(i_kind) iskip,jextra,nextra
    integer(i_kind) status(mpi_status_size)
    integer(i_kind) request,request_ldmk,request_ldmk_snow,request_ldmk_seaice
    integer(i_kind) jbegin(0:npe),jend(0:npe-1),jend2(0:npe-1)
    integer(i_kind) kbegin(0:npe),kend(0:npe-1)
    integer(i_long),allocatable:: ibuf(:,:)
    integer(i_long),allocatable:: jbuf(:,:,:)
    real(r_single),allocatable::mub(:,:), landmask(:,:),snow(:,:),seaice(:,:)
    real(r_single),allocatable::t1st2d(:,:)
    integer(i_kind) kdim_mub,i_snowT_check
    integer(i_kind) kt,kq,ku,kv
    integer(i_kind) mfcst
    integer(i_long) iyear,imonth,iday,ihour,iminute,isecond,dummy3(3)
    real(r_single) pt_regional_single
    real(r_kind) deltasigma,deltasigmac4h
    integer(i_kind) ip1,jp1
    character(1) chdrbuf(2048)
    integer(i_kind) iadd,ier,istatus,n_actual_clouds
    character(132) memoryorder
  
    real(r_kind), pointer :: ges_ps(:,:)=>NULL()
    real(r_kind), pointer :: ges_tsk(:,:)=>NULL()
    real(r_kind), pointer :: ges_th2(:,:)=>NULL()
    real(r_kind), pointer :: ges_q2(:,:)=>NULL()
    real(r_kind), pointer :: ges_soilt1(:,:)=>NULL()
    real(r_kind), pointer :: ges_tslb_it(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_smois_it(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_u (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_v (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_q (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qc(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qi(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qr(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qs(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qg(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qnr(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qni(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qnc(:,:,:)=>NULL()
    logical print_verbose
  
    print_verbose=.false.
    if(verbose .and. mype == 0)print_verbose=.true.
    it=ntguessig
  
    ksize=0
  
  ! Inquire about cloud guess fields
    call gsi_metguess_get('clouds::3d',n_actual_clouds,istatus)
    if (n_actual_clouds>0) then
  !    get pointer to relevant instance of cloud-related backgroud
       ier=0
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql', ges_qc, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi', ges_qi, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr', ges_qr, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qs', ges_qs, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qg', ges_qg, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr,istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qni',ges_qni,istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnc',ges_qnc,istatus );ier=ier+istatus
       if (ier/=0) n_actual_clouds=0
    end if
  
  !   1.  get offsets etc only for records to be updated
  
  !        they are MU, T, Q, U, V,  skint/sst       (for MU, need MU_B and pt_regional)
  
    im=nlon_regional
    jm=nlat_regional
    lm=nsig
  
    num_mass_fields=4*lm+4
    if(l_hydrometeor_bkio .and. n_actual_clouds>0) num_mass_fields=4*lm+4+9*lm
    if(l_gsd_soilTQ_nudge) num_mass_fields=4*lm+4+2*nsig_soil+2
    if(l_hydrometeor_bkio .and. l_gsd_soilTQ_nudge) num_mass_fields=4*lm+4+9*lm+2*nsig_soil+2
    allocate(offset(num_mass_fields))
    allocate(igtype(num_mass_fields),kdim(num_mass_fields),kord(num_mass_fields))
    allocate(length(num_mass_fields))
    allocate(t1st2d(im,jm))
  
  !    igtype is a flag indicating whether each input MASS field is h-, u-, or v-grid
  !    and whether integer or real
  !     abs(igtype)=1 for h-grid
  !                =2 for u-grid
  !                =3 for v-grid
  !
  !     igtype < 0 for integer field
  
  !    offset is the byte count preceding each record to be read/written from/to the wrf binary file.
  !       used as individual file pointers by mpi_file_read/mpi_file_write
  
    write(filename,'("sigf",i2.2)')ifilesig(it)
    open(lendian_in,file=filename,form='unformatted') ; rewind lendian_in
    if(mype == 0) write(6,*)'READ_WRF_MASS_OFFSET_FILE:  open lendian_in=',lendian_in,' to file=',filename
    read(lendian_in) iyear,imonth,iday,ihour,iminute,isecond,dummy3,pt_regional_single
    do iskip=2,5
       read(lendian_in)
    end do
    read(lendian_in) 
    read(lendian_in) n_position          !  offset for START_DATE record
    offset_start_date=n_position
    length_start_date=2048
  
    i=0
    read(lendian_in) n_position                                         ! mub
    offset_mub=n_position
    length_mub=im*jm
    igtype_mub=1
    kdim_mub=1
  
  
  !     open wrf file for mpi-io reading and writing
    wrfanl = 'wrf_inout'
    call mpi_file_open(mpi_comm_world,trim(wrfanl),mpi_mode_rdwr,mpi_info_null,mfcst,ierror)
  
  !     update START_DATE record so it contains new analysis time in place of old starting time
    call mpi_file_read_at(mfcst,offset_start_date,chdrbuf,length_start_date,mpi_byte,status,ierror)
    if(mype==0)  then
       call update_start_date(chdrbuf,iyear,imonth,iday,ihour,iminute,isecond)
       call mpi_file_write_at(mfcst,offset_start_date,chdrbuf,length_start_date,mpi_byte,status,ierror)
    end if
  
  !          begin read of mub while doing rest of offset setup
    allocate(mub(im,jm))
    call mpi_file_iread_at(mfcst,offset_mub,mub,length_mub,mpi_real4,request,ierror)
  
  
    i=i+1 ; i_mu =i                                                ! mu
    read(lendian_in) n_position
    offset(i)=n_position ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=1
    if(print_verbose) write(6,*)' mu, i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
  
    read(lendian_in) n_position                                    !  geopotential  (should this be updated??)
  
    i_t=i+1
    read(lendian_in) n_position,memoryorder
    do k=1,lm
       i=i+1                                                       ! theta(k)  (pot temp)
       if(trim(memoryorder)=='XZY') then
          iadd=0
          kord(i)=lm
       else
          iadd=(k-1)*im*jm*4
          kord(i)=1
       end if
       offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=lm
       if(print_verbose.and.k==1) write(6,*)' temp i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
    end do
  
    i_q=i+1
    read(lendian_in) n_position,memoryorder
    do k=1,lm
       i=i+1                                                       ! q(k)
       if(trim(memoryorder)=='XZY') then
          iadd=0
          kord(i)=lm
       else
          iadd=(k-1)*im*jm*4
          kord(i)=1
       end if
       offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=lm
       if(print_verbose .and.k==1) write(6,*)' q i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
    end do
  
    i_u=i+1
    read(lendian_in) n_position,memoryorder
    do k=1,lm
       i=i+1                                                       ! u(k)
       if(trim(memoryorder)=='XZY') then
          iadd=0
          kord(i)=lm
       else
          iadd=(k-1)*(im+1)*jm*4
          kord(i)=1
       end if
       offset(i)=n_position+iadd
       igtype(i)=2 ; kdim(i)=lm
       length(i)=(im+1)*jm
       if(print_verbose.and.k==1) write(6,*)' u i,igtype,offset,kdim(i),kord(i) = ', &
                                                             i,igtype(i),offset(i),kdim(i),kord(i)
    end do
  
    i_v=i+1
    read(lendian_in) n_position,memoryorder
    do k=1,lm
       i=i+1                                                       ! v(k)
       if(trim(memoryorder)=='XZY') then
          iadd=0
          kord(i)=lm
       else
          iadd=(k-1)*im*(jm+1)*4
          kord(i)=1
       end if
       offset(i)=n_position+iadd ; length(i)=im*(jm+1) ; igtype(i)=3 ; kdim(i)=lm
       if(print_verbose.and.k==1) write(6,*)' v i,igtype,offset,kdim(i),kord(i) = ', &
                                                              i,igtype(i),offset(i),kdim(i),kord(i)
    end do
  
    read(lendian_in) n_position                                         ! landmask
    allocate(landmask(im,jm))
    call mpi_file_iread_at(mfcst,n_position,landmask,im*jm,mpi_real4,request_ldmk,ierror)
  
    read(lendian_in) n_position                                         ! xice
    allocate(seaice(im,jm))
    call mpi_file_iread_at(mfcst,n_position,seaice,im*jm,mpi_real4,request_ldmk_seaice,ierror)
  
    i=i+1 ; i_sst=i                                                ! sst
    read(lendian_in) n_position
    offset(i)=n_position ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=1
    if(print_verbose) write(6,*)' sst i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
  
    read(lendian_in)                                                    ! ivgtyp
    read(lendian_in)                                                    ! isltyp
    read(lendian_in)                                                    ! vegfrac
  
    read(lendian_in) n_position                                         ! sno
    allocate(snow(im,jm))
    call mpi_file_iread_at(mfcst,n_position,snow,im*jm,mpi_real4,request_ldmk_snow,ierror)
  
    read(lendian_in)                                                    ! u10
    read(lendian_in)                                                    ! v10
    if(l_gsd_soilTQ_nudge) then
       i_smois=i+1
       read(lendian_in) n_position,ksize,memoryorder
       do k=1,ksize
          i=i+1                                                         ! smois(k)
          if(trim(memoryorder)=='XZY') then
             iadd=0
             kord(i)=ksize
          else
             iadd=(k-1)*im*jm*4
             kord(i)=1
          end if
          offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=ksize
          if(print_verbose.and.k==1) write(6,*)' smois i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
       end do
  
       i_tslb=i+1
       read(lendian_in) n_position,ksize,memoryorder
       do k=1,ksize
          i=i+1                                                      ! tslb(k)
          if(trim(memoryorder)=='XZY') then
             iadd=0
             kord(i)=ksize
          else
             iadd=(k-1)*im*jm*4
             kord(i)=1
          end if
          offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=ksize
          if(print_verbose.and.k==1) write(6,*)' tslb i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
       end do
    else
       i_tslb=0
       i_smois=0
       read(lendian_in)                                                 ! smois
       read(lendian_in)                                                 ! tslb
    endif
  
    i=i+1 ; i_tsk=i                                                ! tsk
    read(lendian_in) n_position
    offset(i)=n_position ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=1
    if(print_verbose) write(6,*)' tsk i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
  
  
    i=i+1 ; i_q2=i                                                ! q2
    read(lendian_in) n_position
    offset(i)=n_position ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=1
    if(print_verbose) write(6,*)' q2 i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
  
    if(l_gsd_soilTQ_nudge) then
       i=i+1 ; i_soilt1=i                                            ! soilt1
       read(lendian_in) n_position
       offset(i)=n_position ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=1
       if(print_verbose) write(6,*)' soilt1 i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
  
       i=i+1 ; i_th2=i                                               ! th2
       read(lendian_in) n_position
       offset(i)=n_position ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=1
       if(print_verbose) write(6,*)' th2, i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
    else
       i_soilt1=0
       i_th2=0
    endif
  
  ! for cloud/hydrometeor analysis fields
    if(l_hydrometeor_bkio .and. n_actual_clouds>0) then
  
       i_qc=i+1
       read(lendian_in) n_position,memoryorder
       do k=1,lm
          i=i+1                                                       ! qc(k)
          if(trim(memoryorder)=='XZY') then
             iadd=0
             kord(i)=lm
          else
             iadd=(k-1)*im*jm*4
             kord(i)=1
          end if
          offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=lm
          if(print_verbose.and.k==1) write(6,*)' qc i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
       end do
  
       i_qr=i+1
       read(lendian_in) n_position,memoryorder
       do k=1,lm
          i=i+1                                                       ! qr(k)
          if(trim(memoryorder)=='XZY') then
             iadd=0
             kord(i)=lm
          else
             iadd=(k-1)*im*jm*4
             kord(i)=1
           end if
          offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=lm
          if(print_verbose.and.k==1) write(6,*)' qr i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
       end do
  
       i_qi=i+1
       read(lendian_in) n_position,memoryorder
       do k=1,lm
          i=i+1                                                       ! qi(k)
          if(trim(memoryorder)=='XZY') then
             iadd=0
             kord(i)=lm
          else
             iadd=(k-1)*im*jm*4
             kord(i)=1
          end if
          offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=lm
          if(print_verbose.and.k==1) write(6,*)' qi i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
       end do
  
       i_qs=i+1
       read(lendian_in) n_position,memoryorder
       do k=1,lm
          i=i+1                                                       ! qs(k)
          if(trim(memoryorder)=='XZY') then
             iadd=0
             kord(i)=lm
          else
             iadd=(k-1)*im*jm*4
             kord(i)=1
          end if
          offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=lm
          if(print_verbose.and.k==1) write(6,*)' qs i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
       end do
   
       i_qg=i+1
       read(lendian_in) n_position,memoryorder
       do k=1,lm
          i=i+1                                                       ! qg(k)
          if(trim(memoryorder)=='XZY') then
             iadd=0
             kord(i)=lm
          else
             iadd=(k-1)*im*jm*4
             kord(i)=1
          end if
          offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=lm
          if(print_verbose.and.k==1) write(6,*)' qg i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
       end do
  
       i_qnr=i+1
       read(lendian_in) n_position,memoryorder
       do k=1,lm
          i=i+1                                                       !  qnr(k)
          if(trim(memoryorder)=='XZY') then
             iadd=0
             kord(i)=lm
          else
             iadd=(k-1)*im*jm*4
             kord(i)=1
          end if
          offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=lm
          if(print_verbose.and.k==1) write(6,*)' qnr i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
       end do

       i_qni=i+1
       read(lendian_in) n_position,memoryorder
       do k=1,lm
          i=i+1                                                       !  qni(k)
          if(trim(memoryorder)=='XZY') then
             iadd=0
             kord(i)=lm
          else
             iadd=(k-1)*im*jm*4
             kord(i)=1
          end if
          offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=lm
          if(print_verbose.and.k==1) write(6,*)' qni i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
       end do

       i_qnc=i+1
       read(lendian_in) n_position,memoryorder
       do k=1,lm
          i=i+1                                                       !  qnc(k)
          if(trim(memoryorder)=='XZY') then
             iadd=0
             kord(i)=lm
          else
             iadd=(k-1)*im*jm*4
             kord(i)=1
          end if
          offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=lm
          if(print_verbose.and.k==1) write(6,*)' qnc i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
       end do
  
       i_tt=i+1
       read(lendian_in) n_position,memoryorder
       do k=1,lm
          i=i+1                                                       ! tt(k)
          if(trim(memoryorder)=='XZY') then
             iadd=0
             kord(i)=lm
          else
             iadd=(k-1)*im*jm*4
             kord(i)=1
          end if
          offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=1 ; kdim(i)=lm
          if(print_verbose.and.k==1) write(6,*)' tt i,igtype,offset,kdim(i) = ',i,igtype(i),offset(i),kdim(i)
       end do
  
    endif    ! l_hydrometeor_bkio
  
    close(lendian_in)
  
  !          set up evenly distributed index range over all processors for all input fields
  
  
    num_loc_groups=num_mass_fields/npe
    nextra=num_mass_fields-num_loc_groups*npe
    kbegin(0)=1
    if(nextra > 0) then
       do k=1,nextra
          kbegin(k)=kbegin(k-1)+1+num_loc_groups
       end do
    end if
    do k=nextra+1,npe
       kbegin(k)=kbegin(k-1)+num_loc_groups
    end do
    do k=0,npe-1
       kend(k)=kbegin(k+1)-1
    end do
!   if(mype == 0) then
!      write(6,*)' kbegin=',kbegin
!      write(6,*)' kend= ',kend
!   end if
    num_j_groups=jm/npe
    jextra=jm-num_j_groups*npe
    jbegin(0)=1
    if(jextra > 0) then
       do j=1,jextra
          jbegin(j)=jbegin(j-1)+1+num_j_groups
       end do
    end if
    do j=jextra+1,npe
       jbegin(j)=jbegin(j-1)+num_j_groups
    end do
    do j=0,npe-1
     jend(j)=min(jbegin(j+1)-1,jm)
    end do
!   if(mype == 0) then
!      write(6,*)' jbegin=',jbegin
!      write(6,*)' jend= ',jend
!   end if
  
  !     sub2grid to get tempa for fields to be updated
  
  
  ! Create all_loc from ges_*
    allocate(all_loc(lat1,lon1,num_mass_fields))
    all_loc=zero_single
    kt=i_t-1
    kq=i_q-1
    ku=i_u-1
    kv=i_v-1
    q_integral=one
    q_integralc4h=zero_single
  ! for hydrometeors
    if(l_hydrometeor_bkio .and. n_actual_clouds>0) then
  !    get pointer to relevant instance of cloud-related backgroud
       ier=0
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql', ges_qc, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi', ges_qi, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr', ges_qr, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qs', ges_qs, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qg', ges_qg, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr,istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qni',ges_qni,istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnc',ges_qnc,istatus );ier=ier+istatus
       if (ier/=0) then
           write(6,*)'wrwrfmassa_binary: getpointer failed, cannot do cloud analysis'
           call stop2(999)
       endif
       kqc=i_qc-1
       kqi=i_qi-1
       kqr=i_qr-1
       kqs=i_qs-1
       kqg=i_qg-1
       kqnr=i_qnr-1
       kqni=i_qni-1
       kqnc=i_qnc-1
       ktt=i_tt-1
    endif
    ier=0
    call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ps', ges_ps, istatus );ier=ier+istatus
    call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'u' , ges_u , istatus );ier=ier+istatus
    call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'v' , ges_v , istatus );ier=ier+istatus
    call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q' , ges_q , istatus );ier=ier+istatus
    if (ier/=0) then ! doesn't have to die - code can be generalized to bypass missing vars
        write(6,*)'wrwrfmassa_binary: getpointer failed, cannot retrieve ps,u,v,q'
        call stop2(999)
    endif
    do k=1,nsig
       deltasigma=eta1_ll(k)-eta1_ll(k+1)
       deltasigmac4h=eta2_ll(k)-eta2_ll(k+1)
       kt=kt+1
       kq=kq+1
       ku=ku+1
       kv=kv+1
  ! for hydrometeors
       if(l_hydrometeor_bkio .and. n_actual_clouds>0) then
          kqc=kqc+1
          kqi=kqi+1
          kqr=kqr+1
          kqs=kqs+1
          kqg=kqg+1
          kqnr=kqnr+1
          kqni=kqni+1
          kqnc=kqnc+1
          ktt=ktt+1
       endif
       do i=1,lon1
          ip1=i+1
          do j=1,lat1
             jp1=j+1
             all_loc(j,i,ku)=ges_u(jp1,ip1,k)
             all_loc(j,i,kv)=ges_v(jp1,ip1,k)
  
  
  !          Convert sensible temperature to potential temperature
             work_prsl  = one_tenth*(aeta1_ll(k)*(r10*ges_ps(jp1,ip1)-pt_ll)+aeta2_ll(k)+pt_ll)
             work_prslk = (work_prsl/r100)**rd_over_cp_mass
             all_loc(j,i,kt) = ges_tsen(jp1,ip1,k,it)/work_prslk
  
  
  !          Subtract offset to get output potential temperature
             all_loc(j,i,kt)=all_loc(j,i,kt) - h300
  
  
  !          Convert specific humidity to mixing ratio
             all_loc(j,i,kq)= ges_q(jp1,ip1,k)/(one-ges_q(jp1,ip1,k))
  
  ! for hydrometeors
             if(l_hydrometeor_bkio .and. n_actual_clouds>0) then
                all_loc(j,i,kqc)=ges_qc(jp1,ip1,k)
                all_loc(j,i,kqi)=ges_qi(jp1,ip1,k)
                all_loc(j,i,kqr)=ges_qr(jp1,ip1,k)
                all_loc(j,i,kqs)=ges_qs(jp1,ip1,k)
                all_loc(j,i,kqg)=ges_qg(jp1,ip1,k)
                all_loc(j,i,kqnr)=ges_qnr(jp1,ip1,k)
                all_loc(j,i,kqni)=ges_qni(jp1,ip1,k)
                all_loc(j,i,kqnc)=ges_qnc(jp1,ip1,k)
                all_loc(j,i,ktt)=ges_tten(jp1,ip1,k,it)
             endif
  
             q_integral(jp1,ip1)=q_integral(jp1,ip1)+deltasigma*ges_q(jp1,ip1,k)/(one-ges_q(jp1,ip1,k))
             q_integralc4h(jp1,ip1)=q_integralc4h(jp1,ip1)+deltasigmac4h*ges_q(jp1,ip1,k)/(one-ges_q(jp1,ip1,k))
  
          end do
       end do
    end do
  
    do i=1,lon1
       ip1=i+1
       do j=1,lat1
          jp1=j+1
          psfc_this=r10*ges_ps(jp1,ip1)   ! convert from cb to mb
          psfc_this_dry=pt_ll+((psfc_this-pt_ll)-q_integralc4h(jp1,ip1))/q_integral(jp1,ip1)
          all_loc(j,i,i_mu)=r100*psfc_this_dry
       end do
    end do
  
  ! Load updated skin temperature array if writing out to analysis file
    if (update_regsfc) then
       do i=1,lon1
          ip1=i+1
          do j=1,lat1
             jp1=j+1
             all_loc(j,i,i_sst)=dsfct(jp1,ip1,ntguessfc)
             all_loc(j,i,i_tsk)=dsfct(jp1,ip1,ntguessfc)
          end do
       end do
    end if
  
  ! for soil nudging
    if(l_gsd_soilTQ_nudge) then
       ier=0
       call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'smoist',ges_smois_it,istatus)
       ier=ier+istatus
       call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'tslb' ,ges_tslb_it ,istatus)
       ier=ier+istatus
       if (ier/=0) call die(trim(myname),'cannot get pointers for tslb/smois, ier =',ier)
       ktslb=i_tslb-1
       ksmois=i_smois-1
       do k=1,nsig_soil
          ktslb=ktslb+1
          ksmois=ksmois+1
          do i=1,lon1
             ip1=i+1
             do j=1,lat1
                jp1=j+1
                all_loc(j,i,ktslb)=ges_tslb_it(jp1,ip1,k)
                all_loc(j,i,ksmois)=ges_smois_it(jp1,ip1,k)
             end do
          end do
       end do
       ier=0
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'th2m',  ges_th2, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q2m' ,  ges_q2,  istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'tsoil', ges_soilt1, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'tskn' , ges_tsk , istatus );ier=ier+istatus
       if (ier/=0) then ! doesn't have to die - code can be generalized to bypass missing vars
           write(6,*)'wrwrfmassa_binary: getpointer failed, cannot retrieve soil stuff'
           call stop2(999)
       endif
       do i=1,lon1
          ip1=i+1
          do j=1,lat1
             jp1=j+1
             all_loc(j,i,i_tsk)=ges_tsk(jp1,ip1)
             all_loc(j,i,i_th2)=ges_th2(jp1,ip1)
             all_loc(j,i,i_soilt1)=ges_soilt1(jp1,ip1)
          end do
       end do
    endif ! l_gsd_soilTQ_nudge
    if (i_use_2mq4b > 0) then
       do i=1,lon1
          ip1=i+1
          do j=1,lat1
  !  Convert 2m specific humidity to mixing ratio
             jp1=j+1
             all_loc(j,i,i_q2)=ges_q2(jp1,ip1)/(one-ges_q2(jp1,ip1))
          end do
       end do
    endif
  
    allocate(tempa(itotsub,kbegin(mype):kend(mype)))
    call this%generic_sub2grid(all_loc,tempa,kbegin(mype),kend(mype),kbegin,kend,mype,num_mass_fields)
    deallocate(all_loc)
  
    allocate(ibuf((im+1)*(jm+1),kbegin(mype):kend(mype)))
  !      finish reading in mub
    call mpi_wait(request,status,ierror)
    if(byte_swap) then
       num_swap=length_mub
       call to_native_endianness_r4(mub(1,1),num_swap)
    end if
    call mpi_wait(request_ldmk,status,ierror)
    if(byte_swap) then
       num_swap=im*jm
       call to_native_endianness_r4(landmask(1,1),num_swap) !got landmask for
                                                            ! soil nudging in 2X
                                                            ! grid
    end if
    call mpi_wait(request_ldmk_snow,status,ierror)
    if(byte_swap) then
       num_swap=im*jm
       call to_native_endianness_r4(snow(1,1),num_swap)     !got snow for
                                                            ! soil nudging in 2X
                                                            ! grid
    end if
    call mpi_wait(request_ldmk_seaice,status,ierror)
    if(byte_swap) then
       num_swap=im*jm
       call to_native_endianness_r4(seaice(1,1),num_swap)   !got seaice for
                                                            ! soil nudging in 2X
                                                            ! grid
    end if
  
  
  !   2.  create ibuf with records to be updated read in
  
  !                                    read temps
    if(kord(i_t)/=1) then
       allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
       this_offset=offset(i_t)+(jbegin(mype)-1)*4*im*lm
       this_length=(jend(mype)-jbegin(mype)+1)*im*lm
       call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
       if(byte_swap) then
          num_swap=this_length
          call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
       end if
       call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                          jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_t,i_t+lm-1)
       deallocate(jbuf)
    end if
  
  !                                    read q
    if(kord(i_q)/=1) then
       allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
       this_offset=offset(i_q)+(jbegin(mype)-1)*4*im*lm
       this_length=(jend(mype)-jbegin(mype)+1)*im*lm
       call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
       if(byte_swap) then
          num_swap=this_length
          call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
       end if
       call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                          jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_q,i_q+lm-1)
       deallocate(jbuf)
    end if
  
  !                                    read u
    if(kord(i_u)/=1) then
       allocate(jbuf(im+1,lm,jbegin(mype):min(jend(mype),jm)))
       this_offset=offset(i_u)+(jbegin(mype)-1)*4*(im+1)*lm
       this_length=(jend(mype)-jbegin(mype)+1)*(im+1)*lm
       call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
       if(byte_swap) then
          num_swap=this_length
          call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
       end if
       call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                          jbegin,jend,kbegin,kend,mype,npe,im+1,jm,lm,im+1,jm+1,i_u,i_u+lm-1)
       deallocate(jbuf)
    end if
  
  
  !                                    read v
    if(kord(i_v)/=1) then
       jend2=jend
       jend2(npe-1)=jend2(npe-1)+1
       allocate(jbuf(im,lm,jbegin(mype):jend2(mype)))
       this_offset=offset(i_v)+(jbegin(mype)-1)*4*im*lm
       this_length=(jend2(mype)-jbegin(mype)+1)*im*lm
       call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
       if(byte_swap) then
          num_swap=this_length
          call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
       end if
       call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend2(mype),ibuf,kbegin(mype),kend(mype), &
                          jbegin,jend2,kbegin,kend,mype,npe,im,jm+1,lm,im+1,jm+1,i_v,i_v+lm-1)
       deallocate(jbuf)
    end if
  
    if(l_gsd_soilTQ_nudge) then
       lm=nsig_soil
  !                                    read smois
       if(kord(i_smois)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          this_offset=offset(i_smois)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                         jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_smois,i_smois+lm-1)
          deallocate(jbuf)
       end if
  !                                    read tslb
       if(kord(i_tslb)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          this_offset=offset(i_tslb)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                         jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_tslb,i_tslb+lm-1)
          deallocate(jbuf)
       end if
  
       lm=nsig
    endif
  
  ! read hydrometeors
    if(l_hydrometeor_bkio .and. n_actual_clouds>0) then
  !                                    read qc
       if(kord(i_qc)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          this_offset=offset(i_qc)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qc,i_qc+lm-1)
          deallocate(jbuf)
       end if
  
  !                                    read qr
       if(kord(i_qr)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          this_offset=offset(i_qr)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qr,i_qr+lm-1)
          deallocate(jbuf)
       end if
  
  !                                    read qi
       if(kord(i_qi)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          this_offset=offset(i_qi)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qi,i_qi+lm-1)
          deallocate(jbuf)
       end if
  
  !                                    read qs
       if(kord(i_qs)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          this_offset=offset(i_qs)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qs,i_qs+lm-1)
          deallocate(jbuf)
       end if
  
  !                                    read qg
       if(kord(i_qg)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          this_offset=offset(i_qg)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qg,i_qg+lm-1)
          deallocate(jbuf)
       end if
   
  !                                    read qnr
       if(kord(i_qnr)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          this_offset=offset(i_qnr)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qnr,i_qnr+lm-1)
          deallocate(jbuf)
       end if

  !                                    read qni
       if(kord(i_qni)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          this_offset=offset(i_qni)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qni,i_qni+lm-1)
          deallocate(jbuf)
       end if
  
  !                                    read qnc
       if(kord(i_qnc)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          this_offset=offset(i_qnc)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qnc,i_qnc+lm-1)
          deallocate(jbuf)
       end if
  
  !                                    read tt
       if(kord(i_tt)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          this_offset=offset(i_tt)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call read_wrf%transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_tt,i_tt+lm-1)
          deallocate(jbuf)
       end if
  
    endif   ! l_hydrometeor_bkio
  
  !---------------------- read surface files last
    do k=kbegin(mype),kend(mype)
       if(kdim(k)==1.or.kord(k)==1) then
          call mpi_file_read_at(mfcst,offset(k),ibuf(1,k),length(k),mpi_integer4,status,ierror)
          if(byte_swap) then
             num_swap=length(k)
             call to_native_endianness_i4(ibuf(1,k),num_swap)
          end if
          if(igtype(k)==1) call wrf_mass_guess%expand_ibuf(ibuf(1,k),im  ,jm     ,im+1,jm+1)
          if(igtype(k)==2) call wrf_mass_guess%expand_ibuf(ibuf(1,k),im+1,jm     ,im+1,jm+1)
          if(igtype(k)==3) call wrf_mass_guess%expand_ibuf(ibuf(1,k),im  ,jm+1,im+1,jm+1)
       end if
    end do
  ! call mpi_barrier(mpi_comm_world,ierror)
  
  !   5.  tempa --> updated ibuf --> jbuf --> write out
  
    allocate(tempb(itotsub,kbegin(mype):kend(mype)))
    allocate(temp1(im,jm),itemp1(im,jm),temp1u(im+1,jm),temp1v(im,jm+1))
    do ifld=kbegin(mype),kend(mype)
       if((ifld==i_sst).and..not.update_regsfc) cycle
       if((ifld==i_tsk).and..not.(update_regsfc .or. l_gsd_soilTQ_nudge)) cycle
       if(igtype(ifld) == 1) then
          call wrf_mass_guess%move_ibuf_hg(ibuf(1,ifld),temp1,im+1,jm+1,im,jm)
          if(ifld==i_mu) then
             temp1=temp1+mub+pt_regional_single
          end if
          call fill_mass_grid2t(temp1,im,jm,tempb(1,ifld),2)
          if(ifld==i_sst.or.ifld==i_tsk) then
             do i=1,iglobal
                if(tempb(i,ifld) < r225) then
                   tempa(i,ifld)=zero_single
                else
                   tempa(i,ifld)=tempa(i,ifld)-tempb(i,ifld)
                end if
             end do
          else
             do i=1,iglobal
                tempa(i,ifld)=tempa(i,ifld)-tempb(i,ifld)
             end do
          end if
          if( ifld==i_tsk .or. ifld==i_soilt1 .or. & 
             (ifld >=i_smois .and. ifld <=i_smois+ksize-1) .or. &
             (ifld >=i_tslb .and. ifld <=i_tslb+ksize-1) ) then 
  ! for 2X soil nudging
             i_snowT_check=0
             if(ifld==i_tsk) i_snowT_check=1
             if(ifld==i_soilt1) i_snowT_check=4
             if(ifld==i_tslb) i_snowT_check=3
             if(ifld >=i_smois .and. ifld <=i_smois+ksize-1) i_snowT_check=2
             call unfill_mass_grid2t_ldmk(tempa(1,ifld),im,jm,temp1,landmask, &
                                          snow,seaice,t1st2d,i_snowT_check)
  ! make sure soil moisture is larger than soilmoistmin (whih is 0.002 (sand)).
             if(ifld >=i_smois .and. ifld <=i_smois+ksize-1) then
                do i=1,im
                   do j=1,jm
                      temp1(i,j) = max(temp1(i,j),soilmoistmin)
                   end do
                end do
             endif
          else
             call unfill_mass_grid2t(tempa(1,ifld),im,jm,temp1)
          endif
          if(ifld==i_mu) then
             temp1=temp1-mub-pt_regional_single
          end if
          call this%move_hg_ibuf(temp1,ibuf(1,ifld),im+1,jm+1,im,jm)
       else if(igtype(ifld) == 2) then
          call wrf_mass_guess%move_ibuf_hg(ibuf(1,ifld),temp1u,im+1,jm+1,im+1,jm)
          call fill_mass_grid2u(temp1u,im,jm,tempb(1,ifld),2)
          tempa(:,ifld)=tempa(:,ifld)-tempb(:,ifld)
          call unfill_mass_grid2u(tempa(1,ifld),im,jm,temp1u)
          call this%move_hg_ibuf(temp1u,ibuf(1,ifld),im+1,jm+1,im+1,jm)
       else if(igtype(ifld) == 3) then
          call wrf_mass_guess%move_ibuf_hg(ibuf(1,ifld),temp1v,im+1,jm+1,im,jm+1)
          call fill_mass_grid2v(temp1v,im,jm,tempb(1,ifld),2)
          tempa(:,ifld)=tempa(:,ifld)-tempb(:,ifld)
          call unfill_mass_grid2v(tempa(1,ifld),im,jm,temp1v)
          call this%move_hg_ibuf(temp1v,ibuf(1,ifld),im+1,jm+1,im,jm+1)
       end if
    end do
  
  !           finally write ibuf back out ( ibuf --> jbuf -->  mpi_file_write )
  
  !                                    write temps
    if(kord(i_t)/=1) then
       allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
       call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                          jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_t,i_t+lm-1)
       this_offset=offset(i_t)+(jbegin(mype)-1)*4*im*lm
       this_length=(jend(mype)-jbegin(mype)+1)*im*lm
       if(byte_swap) then
          num_swap=this_length
          call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
       end if
       call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
       deallocate(jbuf)
    end if
  
  !                                    write q
    if(kord(i_q)/=1) then
       allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
       call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                          jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_q,i_q+lm-1)
       this_offset=offset(i_q)+(jbegin(mype)-1)*4*im*lm
       this_length=(jend(mype)-jbegin(mype)+1)*im*lm
       if(byte_swap) then
          num_swap=this_length
          call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
       end if
       call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
       deallocate(jbuf)
    end if
  
  !                                    write u
    if(kord(i_u)/=1) then
       allocate(jbuf(im+1,lm,jbegin(mype):min(jend(mype),jm)))
       call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                          jbegin,jend,kbegin,kend,mype,npe,im+1,jm,lm,im+1,jm+1,i_u,i_u+lm-1)
       this_offset=offset(i_u)+(jbegin(mype)-1)*4*(im+1)*lm
       this_length=(jend(mype)-jbegin(mype)+1)*(im+1)*lm
       if(byte_swap) then
          num_swap=this_length
          call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
       end if
       call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
       deallocate(jbuf)
    end if
  
  !                                    write v
    if(kord(i_v)/=1) then
       jend2=jend
       jend2(npe-1)=jend2(npe-1)+1
       allocate(jbuf(im,lm,jbegin(mype):jend2(mype)))
       call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend2(mype),ibuf,kbegin(mype),kend(mype), &
                          jbegin,jend2,kbegin,kend,mype,npe,im,jm+1,lm,im+1,jm+1,i_v,i_v+lm-1)
       this_offset=offset(i_v)+(jbegin(mype)-1)*4*im*lm
       this_length=(jend2(mype)-jbegin(mype)+1)*im*lm
       if(byte_swap) then
          num_swap=this_length
          call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
       end if
       call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
       deallocate(jbuf)
    end if
  
    if(l_gsd_soilTQ_nudge) then
       lm=nsig_soil
  !                                    write smois
       if(kord(i_smois)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_smois,i_smois+lm-1)
          this_offset=offset(i_smois)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          deallocate(jbuf)
       end if
  
  !                                    write tslb
       if(kord(i_tslb)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_tslb,i_tslb+lm-1)
          this_offset=offset(i_tslb)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          deallocate(jbuf)
       end if
  
       lm=nsig
    endif
  
  !  write hydrometeors
    if(l_hydrometeor_bkio .and. n_actual_clouds>0) then
  !                                    write qc
       if(kord(i_qc)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qc,i_qc+lm-1)
          this_offset=offset(i_qc)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          deallocate(jbuf)
       end if
  
  !                                    write qr
       if(kord(i_qr)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qr,i_qr+lm-1)
          this_offset=offset(i_qr)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          deallocate(jbuf)
       end if
  
  !                                    write qi
       if(kord(i_qi)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qi,i_qi+lm-1)
          this_offset=offset(i_qi)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          deallocate(jbuf)
       end if
  
  !                                    write qs
       if(kord(i_qs)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qs,i_qs+lm-1)
          this_offset=offset(i_qs)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          deallocate(jbuf)
       end if
  
  !                                    write qg
       if(kord(i_qg)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qg,i_qg+lm-1)
          this_offset=offset(i_qg)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          deallocate(jbuf)
       end if
  
  !                                    write qnr
       if(kord(i_qnr)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qnr,i_qnr+lm-1)
          this_offset=offset(i_qnr)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          deallocate(jbuf)
       end if
  
  !                                    write qni
       if(kord(i_qni)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qni,i_qni+lm-1)
          this_offset=offset(i_qni)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          deallocate(jbuf)
       end if
  
  !                                    write qnc
       if(kord(i_qnc)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_qnc,i_qnc+lm-1)
          this_offset=offset(i_qnc)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          deallocate(jbuf)
       end if
  
  !                                    write tt
       if(kord(i_tt)/=1) then
          allocate(jbuf(im,lm,jbegin(mype):min(jend(mype),jm)))
          call this%transfer_ibuf2jbuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
                           jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+1,jm+1,i_tt,i_tt+lm-1)
          this_offset=offset(i_tt)+(jbegin(mype)-1)*4*im*lm
          this_length=(jend(mype)-jbegin(mype)+1)*im*lm
          if(byte_swap) then
             num_swap=this_length
             call to_native_endianness_i4(jbuf(1,1,jbegin(mype)),num_swap)
          end if
          call mpi_file_write_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length,mpi_integer4,status,ierror)
          deallocate(jbuf)
       end if
  
    end if ! l_hydrometeor_bkio
  !---------------------- write surface files last
    do k=kbegin(mype),kend(mype)
       if(kdim(k)==1.or.kord(k)==1) then
          if(igtype(k)==1) call this%contract_ibuf(ibuf(1,k),im  ,jm     ,im+1,jm+1)
          if(igtype(k)==2) call this%contract_ibuf(ibuf(1,k),im+1,jm     ,im+1,jm+1)
          if(igtype(k)==3) call this%contract_ibuf(ibuf(1,k),im  ,jm+1,im+1,jm+1)
          if(byte_swap) then
          num_swap=length(k)
             call to_native_endianness_i4(ibuf(1,k),num_swap)
          end if
          call mpi_file_write_at(mfcst,offset(k),ibuf(1,k),length(k),mpi_integer4,status,ierror)
       end if
    end do
  
    deallocate(ibuf)
    deallocate(offset)
    deallocate(igtype)
    deallocate(kdim)
    deallocate(kord)
    deallocate(length)
    deallocate(mub)
    deallocate(landmask)
    deallocate(snow)
    deallocate(seaice)
    deallocate(t1st2d)
    deallocate(tempa)
    deallocate(tempb)
    deallocate(temp1)
    deallocate(itemp1)
    deallocate(temp1u)
    deallocate(temp1v)
    call destroy_cld_grids
    call mpi_file_close(mfcst,ierror)
  
  end subroutine wrwrfmassa_binary_wrf
  
  subroutine generic_sub2grid(this,all_loc,tempa,kbegin_loc,kend_loc,kbegin,kend,mype,num_fields)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    generic_sub2grid   converts from subdomains to full horizontal grid
  !   prgmmr: parrish          org: np22                date: 2004-11-29
  !
  ! abstract: variation on subroutine sub2grid, with more general distribution of variables
  !              along the k index.
  !
  ! program history log:
  !   2004-02-03  kleist, new mpi strategy
  !   2004-05-06  derber
  !   2004-07-15  treadon - handle periodic subdomains
  !   2004-07-28  treadon - add only on use declarations; add intent in/out
  !   2004-10-26  kleist - u,v removed; periodicity accounted for only in
  !               sub2grid routine if necessary
  !   2004-11-29  parrish - adapt sub2grid for related use with mpi io.
  !
  !   input argument list:
  !     all_loc  - input grid values in vertical subdomain mode
  !     kbegin_loc - starting k index for tempa on local processor
  !     kend_loc   - ending k index for tempa on local processor
  !     kbegin     - starting k indices for tempa for all processors
  !     kend       - ending k indices for tempa for all processors
  !     mype       - local processor number
  !     num_fields - total range of k index (1 <= k <= num_fields)
  !
  !   output argument list:
  !     tempa    - output grid values in horizontal slab mode.
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm RS/6000 SP
  !
  !$$$
    use mpimod, only: ierror,mpi_comm_world,mpi_real4,npe
    use gridmod, only: ijn,itotsub,lat1,lon1
    use kinds, only: r_single,i_kind
    implicit none
  
    class(wrwrfmassa_class), intent(inout) :: this 
    integer(i_kind),intent(in   ) :: kbegin_loc,kend_loc,mype,num_fields
    integer(i_kind),intent(in   ) :: kbegin(0:npe),kend(0:npe-1)
    real(r_single) ,intent(  out) :: tempa(itotsub,kbegin_loc:kend_loc)
    real(r_single) ,intent(in   ) :: all_loc(lat1*lon1*num_fields)
  
    integer(i_kind) k
    integer(i_kind) sendcounts(0:npe-1),sdispls(0:npe),recvcounts(0:npe-1),rdispls(0:npe)
  
  ! first get alltoallv indices
  
    sdispls(0)=0
    do k=0,npe-1
       sendcounts(k)=ijn(k+1)*(kend_loc-kbegin_loc+1)
       sdispls(k+1)=sdispls(k)+sendcounts(k)
    end do
    rdispls(0)=0
    do k=0,npe-1
       recvcounts(k)=ijn(mype+1)*(kend(k)-kbegin(k)+1)
       rdispls(k+1)=rdispls(k)+recvcounts(k)
    end do
  
    call mpi_alltoallv(all_loc,recvcounts,rdispls,mpi_real4, &
                  tempa,sendcounts,sdispls,mpi_real4,mpi_comm_world,ierror)
  
    call this%reorder_s(tempa,kend_loc-kbegin_loc+1)
  
  end subroutine generic_sub2grid
  
    subroutine reorder_s(work,k_in)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    reorder_s 
  !   prgmmr: parrish          org: np22                date: 2004-11-29
  !
  ! abstract:  adapt reorder to work with single precision
  !
  ! program history log:
  !   2004-01-25  kleist
  !   2004-05-14  kleist, documentation
  !   2004-07-15  todling, protex-complaint prologue
  !   2004-11-29  adapt reorder to work with single precision
  !   2008-03-27  safford - add standard doc block, rm unused vars
  !
  !   input argument list:
  !     k_in     - number of levs in work array
  !     work     - array to reorder
  !
  !   output argument list:
  !     work     - array to reorder
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm rs/6000 sp; sgi origin 2000; compaq/hp
  !
  !$$$
  
      use mpimod, only: npe
      use kinds, only: r_single,i_kind
      use constants, only: zero
      use gridmod, only: ijn,itotsub
      implicit none
  
      integer(i_kind)                       , intent(in   ) :: k_in    ! number of levs in work array
  
      real(r_single),dimension(itotsub*k_in), intent(inout) :: work ! array to reorder
  
  
      integer(i_kind) iloc,iskip,i,k,n
      real(r_single),dimension(itotsub,k_in):: temp
  
  ! Zero out temp array
      do k=1,k_in
         do i=1,itotsub
            temp(i,k)=zero
         end do
      end do
  
  ! Load temp array in desired order
      do k=1,k_in
         iskip=0
         iloc=0
         do n=1,npe
            if (n/=1) then
               iskip=iskip+ijn(n-1)*k_in
            end if
            do i=1,ijn(n)
               iloc=iloc+1
               temp(iloc,k)=work(i + iskip + &
                        (k-1)*ijn(n))
            end do
         end do
      end do
  
  ! Load the temp array back into work
      iloc=0
      do k=1,k_in
         do i=1,itotsub
            iloc=iloc+1
            work(iloc)=temp(i,k)
         end do
      end do
  
      return
    end subroutine reorder_s
  
  subroutine contract_ibuf(ibuf,im,jm,imp,jmp)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    contract_ibuf    contract array in place
  !   prgmmr: parrish          org: np22                date: 2004-11-29
  !
  ! abstract: contract array in place from imp,jmp to im,jm
  !
  ! program history log:
  !   2004-11-29  parrish
  !   2005-06-10  devenyi/treadon - remove ii=ii+1 from itemp(i,j)=0 loop
  !   2007-04-12  parrish - replace im+1, jm+1 with inputs imp, jmp to allow
  !                           for use with u and v fields, where im=imp or jm=jmp
  !
  !   input argument list:
  !     ibuf     - input grid values in imp,jmp
  !     im       - first grid index
  !     jm       - second grid index
  !     imp
  !     jmp
  !
  !   output argument list:
  !     ibuf     - output grid values in im,jm
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm RS/6000 SP
  !
  !$$$
  
  !   field of dim imp*jmp needs to be converted to field of im*jm
  
    use kinds, only: i_long,i_kind
    implicit none
  
    integer(i_kind),intent(in   ) :: im,jm,imp,jmp
    integer(i_long),intent(inout) :: ibuf(imp*jmp)
  
    integer(i_kind) i,ii,j
    integer(i_long) itemp(imp,jmp)
  
    do j=1,jmp
       do i=1,imp
          itemp(i,j)=0_i_long
       end do
    end do
    ii=0
    do j=1,jmp
       do i=1,imp
          ii=ii+1
          itemp(i,j)=ibuf(ii)
       end do
    end do
  
    ibuf=0_i_long
    ii=0
    do j=1,jm
       do i=1,im
          ii=ii+1
          ibuf(ii)=itemp(i,j)
       end do
    end do
  
  end subroutine contract_ibuf
  
  subroutine transfer_ibuf2jbuf(jbuf,jbegin_loc,jend_loc,ibuf,kbegin_loc,kend_loc, &
                       jbegin,jend,kbegin,kend,mype,npe,im_jbuf,jm_jbuf,lm_jbuf, &
                       im_ibuf,jm_ibuf,k_start,k_end)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    transfer_jbuf2ibuf   flip from ijk to ikj
  !   prgmmr: parrish          org: np22                date: 2004-11-29
  !
  ! abstract: redistribute 3-d field from ijk order to ikj order across processors
  !
  ! program history log:
  !   2004-11-29  parrish
  !   2005-02-16  todling, replaced "use mpi" by mpimod
  !
  !   input argument list:
  !     ibuf     - input grid redistributed to ijk order (full i, full j, k_start <= k <= k_end)
  !     jbegin_loc - local processor starting j index
  !     jend_loc   - local processor ending j index
  !     kbegin_loc - local processor starting k index
  !     kend_loc   - local processor ending k index
  !     jbegin     - starting j indices for all processors
  !     jend       - ending j indices for all processors
  !     kbegin     - starting k indices for all processors
  !     kend       - ending k indices for all processors
  !     mype       - local processor number
  !     npe        - total number of processors
  !     im_jbuf    - full range of i index for jbuf array
  !     jm_jbuf    - full range of j index for jbuf array
  !     lm_jbuf    - full range of k index for jbuf array
  !     im_ibuf    - full range of i index for ibuf array
  !     jm_ibuf    - full range of j index for ibuf array
  !     k_start    - beginning index for range of k in ibuf array
  !     k_end      - ending index for range of k in ibuf array
  !
  !   output argument list:
  !     jbuf     - output grid values distributed as all ik and a range of j on each processor
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm RS/6000 SP
  !
  !$$$
  
  !  flip around from ijk to ikj, moving result from ibuf to jbuf
  
    use mpimod, only: mpi_comm_world,mpi_integer4
    use kinds, only: i_long,i_kind
    implicit none
  
    integer(i_kind),intent(in   ) :: jbegin_loc,jend_loc,kbegin_loc,kend_loc,mype,npe,im_jbuf,jm_jbuf,lm_jbuf
    integer(i_kind),intent(in   ) :: im_ibuf,jm_ibuf,k_start,k_end
  
    integer(i_long),intent(  out) :: jbuf(im_jbuf,lm_jbuf,jbegin_loc:jend_loc)
    integer(i_long),intent(in   ) :: ibuf(im_ibuf,jm_ibuf,kbegin_loc:kend_loc)
    integer(i_kind),intent(in   ) :: jbegin(0:npe),jend(0:npe-1)
    integer(i_kind),intent(in   ) :: kbegin(0:npe),kend(0:npe-1)
  
    integer(i_long) sendbuf(im_jbuf*lm_jbuf*(min(jend_loc,jm_jbuf)-jbegin_loc+1))
    integer(i_long) recvbuf(im_jbuf*jm_jbuf*(kend_loc-kbegin_loc+1))
    integer(i_long) recvcounts(0:npe-1),displs(0:npe)
    integer(i_kind) i,ipe,j,ierror,k,n,ii,k_t_start,k_t_end,sendcount
  
    do ipe=0,npe-1
       k_t_start=max(k_start,kbegin(ipe))
       k_t_end=  min(k_end,kend(ipe))
       if(k_t_end < k_t_start) cycle
  
       displs(0)=0_i_long
       do i=0,npe-1
          recvcounts(i)=im_jbuf*(k_t_end-k_t_start+1_i_long)*(jend(i)-jbegin(i)+1_i_long)
          displs(i+1)=displs(i)+recvcounts(i)
       end do
  
       if(ipe==mype) then
          ii=0
          do n=0,npe-1
             do k=k_t_start,k_t_end
                do j=jbegin(n),jend(n)
                   do i=1,im_jbuf
                      ii=ii+1
                      recvbuf(ii)=ibuf(i,j,k)
                   end do
                end do
             end do
          end do
       end if
       ii=0
       do k=k_t_start,k_t_end
          do j=jbegin_loc,jend_loc
             do i=1,im_jbuf
                ii=ii+1
             end do
          end do
       end do
       sendcount=ii
       call mpi_scatterv(recvbuf,recvcounts,displs,mpi_integer4, &
                         sendbuf,sendcount,mpi_integer4,ipe,mpi_comm_world,ierror)
       ii=0
       do k=k_t_start,k_t_end
          do j=jbegin_loc,jend_loc
             do i=1,im_jbuf
                ii=ii+1
                jbuf(i,k-k_start+1,j)=sendbuf(ii)
             end do
          end do
       end do
  
  
    end do
  
  end subroutine transfer_ibuf2jbuf
  
  subroutine move_hg_ibuf(temp1,ibuf,im_buf,jm_buf,im_out,jm_out)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    move_hg_ibuf  copy from one array to another
  !   prgmmr: parrish          org: np22                date: 2004-11-29
  !
  ! abstract: copy from one array to another
  !
  ! program history log:
  !   2004-11-29  parrish
  !
  !   input argument list:
  !     temp1    - input grid values
  !     im_buf   - first index of input array buf
  !     jm_buf   - second index of input array buf
  !     im_out   - first index of output array temp1
  !     jm_out   - second index of output array temp1
  !
  !   output argument list:
  !     ibuf      - output grid values
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm RS/6000 SP
  !
  !$$$
  
  !        cp temp1 to ibuf
  
    use kinds, only: r_single,i_kind,i_long
    use constants, only: zero_ilong
    implicit none
  
    integer(i_kind),intent(in   ) :: im_buf,jm_buf,im_out,jm_out
    real(r_single) ,intent(in   ) :: temp1(im_out,jm_out)
    integer(i_long),intent(  out) :: ibuf(im_buf,jm_buf)
  
    integer(i_kind) i,j
  
    ibuf=zero_ilong
    do j=1,jm_out
       do i=1,im_out
          ibuf(i,j)=transfer(temp1(i,j),zero_ilong)
       end do
    end do
  
  end subroutine move_hg_ibuf
  
  subroutine wrwrfmassa_netcdf_wrf(this,mype)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    wrwrfmassa              write out wrf MASS restart file
  !   prgmmr: parrish          org: np22                date: 2004-06-23
  !
  ! abstract:  read wrf MASS guess restart interface file, add analysis
  !            increment, and write out wrf MASS analysis restart 
  !            interface file.
  !
  ! program history log:
  !   2004-06-23  parrish, document
  !   2004-08-03  treadon - add only to module use, add intent in/out
  !   2006-02-15  treadon - convert specific humidity to moisture mixing ratio
  !   2006-03-07  treadon - convert virtual temperature to potential temperature
  !   2006-04-06  middlecoff - changed iog  from 11 to lendian_in
  !                            changed ioan from 51 to lendian_out
  !   2006-07-28  derber  - include sensible temperature
  !   2006-07-31  kleist - change to use ges_ps instead of lnps
  !   2008-03-27  safford - rm unused vars and uses
  !   2008-12-05  todling - adjustment for dsfct time dimension addition
  !   2010-03-29  hu     - add code to gether cloud/hydrometeor fields and write out
  !   2010-04-01  treadon - move strip_single to gridmod
  !   2011-04-29  todling - introduce MetGuess and wrf_mass_guess_mod
  !   2011-09-20  hclin   - added 15 wrfchem/gocart fields for aod
  !   2012-04-13  whitaker - don't call GSI_BundleGetPointer if n_actual_clouds = 0
  !   2013-10-19  todling - metguess now holds background
  !   2013-10-24  todling - general interface to strip
  !   2014-03-12  hu     - add code to read ges_q2 (2m Q), 
  !                               Qnr(rain number concentration), 
  !                               and nsoil (number of soil levels)
  !   2015-01-13  ladwig - add code to read Qni (cloud water number
  !                               concentration)
  !   2017-03-23  Hu     - add code to use hybrid vertical coodinate in WRF MASS
  !                        core
  !
  !   input argument list:
  !     mype     - pe number
  !
  !   output argument list:
  !     no output arguments
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm RS/6000 SP
  !
  !$$$
    use kinds, only: r_kind,r_single,i_kind
    use guess_grids, only: ntguessfc,ntguessig,ifilesig,dsfct,&
         ges_tsen, ges_w_btlev
    use wrf_mass_guess_mod, only: ges_tten
    use mpimod, only: mpi_comm_world,ierror,mpi_real4
    use gridmod, only: pt_ll,eta1_ll,lat2,iglobal,itotsub,update_regsfc,&
         lon2,nsig,nsig_soil,lon1,lat1,nlon_regional,nlat_regional,ijn,displs_g,&
         aeta1_ll,strip,eta2_ll,aeta2_ll
    use constants, only: one,zero_single,rd_over_cp_mass,one_tenth,r10,r100
    use constants, only: soilmoistmin
    use gsi_io, only: lendian_in, lendian_out
    use rapidrefresh_cldsurf_mod, only: l_hydrometeor_bkio,l_gsd_soilTQ_nudge,&
         i_use_2mq4b,i_use_2mt4b
    use chemmod, only: laeroana_gocart,wrf_pm2_5
    use gsi_bundlemod, only: GSI_BundleGetPointer
    use gsi_metguess_mod, only: gsi_metguess_get,GSI_MetGuess_Bundle
    use gsi_chemguess_mod, only: GSI_ChemGuess_Bundle, gsi_chemguess_get
    use mpeu_util, only: die
    use control_vectors, only : w_exist, dbz_exist
    use obsmod,only: if_model_dbz
    implicit none
  
  ! Declare passed variables
    class(wrwrfmassa_class), intent(inout) :: this 
    integer(i_kind),intent(in   ) :: mype
  
  ! Declare local parameters
    character(len=*),parameter::myname='wrwrfmassa_netcdf'
    real(r_kind),parameter:: r225=225.0_r_kind
    integer, parameter :: nlvl_qcheck=3
  
  ! Declare local variables
    integer(i_kind) im,jm,lm
    real(r_single),allocatable::temp1(:),temp1u(:),temp1v(:),tempa(:),tempb(:)
    real(r_single),allocatable::all_loc(:,:,:)
    real(r_single),allocatable::strp(:)
    real(r_single),allocatable::landmask(:),snow(:),seaice(:),t1st2d(:),pt2t(:)
    real(r_single),allocatable::qs1st2d(:,:),sfcprs(:)
    real(r_kind),allocatable::tmp2da(:,:),tmp2db(:,:),tmp2dc(:,:)
    character(6) filename
    integer(i_kind) i,j,k,kt,kq,ku,kv,it,i_psfc,i_t,i_q,i_u,i_v,i_w,i_dbz,ij
    integer(i_kind) i_qc,i_qi,i_qr,i_qs,i_qg,i_qnr,i_qni,i_qnc
    integer(i_kind) kqc,kqi,kqr,kqs,kqg,kqnr,kqni,kqnc,i_tt,ktt,kw,kdbz
    integer(i_kind) i_sst,i_skt,i_th2,i_q2,i_soilt1,i_tslb,i_smois,ktslb,ksmois
    integer(i_kind) :: iv, n_gocart_var,i_snowT_check
    integer(i_kind),allocatable :: i_chem(:), kchem(:)
    integer(i_kind) num_mass_fields,num_all_fields,num_all_pad,num_mass_fields_base
    integer(i_kind) regional_time0(6),nlon_regional0,nlat_regional0,nsig0,nsoil
    integer(i_kind) n_actual_clouds,ier,istatus
    real(r_kind) psfc_this,psfc_this_dry
    real(r_kind),dimension(lat2,lon2):: q_integral,q_integralc4h
    real(r_kind) deltasigma,deltasigmac4h
    real(r_kind):: work_prsl,work_prslk
    real(r_single) pt0
    real(r_single) aeta10(nsig),eta10(nsig+1),aeta20(nsig),eta20(nsig+1)
    real(r_single) glon0(nlon_regional,nlat_regional),glat0(nlon_regional,nlat_regional)
    real(r_single) dx_mc0(nlon_regional,nlat_regional),dy_mc0(nlon_regional,nlat_regional)
  
    real(r_kind), pointer :: ges_ps(:,:  )=>NULL()
    real(r_kind), pointer :: ges_tsk(:,:)=>NULL()
    real(r_kind), pointer :: ges_th2(:,:)=>NULL()
    real(r_kind), pointer :: ges_q2(:,:)=>NULL()
    real(r_kind), pointer :: ges_soilt1(:,:)=>NULL()
    real(r_kind), pointer :: ges_tslb_it(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_smois_it(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_u (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_v (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_w (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_q (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qc(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qi(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qr(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qs(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qg(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qnr(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qni(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_qnc(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_dbz(:,:,:)=>NULL()
  
    real(r_kind), pointer :: ges_sulf (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_bc1  (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_bc2  (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_oc1  (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_oc2  (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_dust1(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_dust2(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_dust3(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_dust4(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_dust5(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_seas1(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_seas2(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_seas3(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_seas4(:,:,:)=>NULL()
    real(r_kind), pointer :: ges_p25  (:,:,:)=>NULL()
    real(r_kind), pointer :: ges_pm2_5  (:,:,:)=>NULL()
  
    associate( this => this ) ! eliminates warning for unused dummy argument needed for binding
    end associate
    it=ntguessig
  
  ! Inquire about cloud guess fields
    call gsi_metguess_get('clouds::3d',n_actual_clouds,istatus)
    if (n_actual_clouds>0) then
  !    get pointer to relevant instance of cloud-related backgroud
       ier=0
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql', ges_qc, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi', ges_qi, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr', ges_qr, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qs', ges_qs, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qg', ges_qg, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr,istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qni',ges_qni,istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnc',ges_qnc,istatus );ier=ier+istatus
       if (ier/=0) n_actual_clouds=0
    end if
  
    im=nlon_regional
    jm=nlat_regional
    lm=nsig
  
    num_mass_fields_base=2+4*lm + 1
    num_mass_fields=num_mass_fields_base
!    The 9 3D cloud analysis fields are: ql,qi,qr,qs,qg,qnr,qni,qnc,tt
    if(l_hydrometeor_bkio .and. n_actual_clouds>0) num_mass_fields=num_mass_fields + 9*lm
    if(l_gsd_soilTQ_nudge) num_mass_fields=num_mass_fields+2*nsig_soil+1
    if(i_use_2mt4b > 0 ) num_mass_fields=num_mass_fields+2
    if(i_use_2mt4b <= 0 .and. i_use_2mq4b > 0) num_mass_fields=num_mass_fields+1
    if ( laeroana_gocart ) then
       call gsi_chemguess_get ( 'aerosols::3d', n_gocart_var, ier )
       if ( n_gocart_var > 0 ) then
          num_mass_fields = num_mass_fields + n_gocart_var*lm
          allocate(i_chem(n_gocart_var))
          allocate(kchem(n_gocart_var))
       else
          laeroana_gocart = .false.
       endif
    endif
  
    if ( wrf_pm2_5 ) then
       num_mass_fields = num_mass_fields + lm
       allocate(i_chem(1))
       allocate(kchem(1))
    endif

    if(w_exist) num_mass_fields = num_mass_fields +lm+1
    if(dbz_exist.and.if_model_dbz) num_mass_fields = num_mass_fields +lm
    
  
    num_all_fields=num_mass_fields
    num_all_pad=num_all_fields
    allocate(all_loc(lat2,lon2,num_all_pad))
    allocate(strp(lat1*lon1))
  
    i_psfc=1
    i_t=2
    i_q=i_t+lm
    i_u=i_q+lm
    i_v=i_u+lm
    if(w_exist)then
      i_w=i_v+lm
      i_sst=i_w+lm+1
    else
      i_sst=i_v+lm
    endif
    if(i_use_2mt4b > 0) then
       i_th2=i_sst+1
    else
       i_th2=i_sst
    endif
    if(l_gsd_soilTQ_nudge) then
       i_tslb=i_th2+1
       i_smois=i_tslb+nsig_soil
       i_soilt1=i_smois+nsig_soil
       i_skt=i_soilt1+1
    else
       i_skt=i_th2+1
       i_tslb=0
       i_smois=0
       i_soilt1=0
    endif
    if(i_use_2mt4b > 0 .or. i_use_2mq4b > 0) then
       i_q2=i_skt+1
    else
       i_q2=i_skt
    endif
  
  ! for hydrometeors
    if(l_hydrometeor_bkio .and. n_actual_clouds>0) then
       i_qc=i_q2+1
       i_qr=i_qc+lm
       i_qs=i_qr+lm
       i_qi=i_qs+lm
       i_qg=i_qi+lm
       i_qnr=i_qg+lm
       i_qni=i_qnr+lm
       i_qnc=i_qni+lm
       if(dbz_exist.and.if_model_dbz)then
         i_dbz=i_qnc+lm
         i_tt=i_dbz+lm
       else
         i_tt=i_qnc+lm
       end if
       if ( laeroana_gocart ) then
          do iv = 1, n_gocart_var
             i_chem(iv)=i_tt+(iv-1)*lm+1
          end do
       endif
  
       if ( wrf_pm2_5 ) then
          iv=1
          i_chem(iv)=i_tt+(iv-1)*lm+1
       endif
       
    else
       if ( laeroana_gocart) then
          do iv = 1, n_gocart_var
             i_chem(iv)=i_skt+(iv-1)*lm+1
          end do
       endif
  
       if ( wrf_pm2_5 ) then
          iv=1
          i_chem(iv)=i_skt+(iv-1)*lm+1
       endif
  
  
    endif
    
    allocate(temp1(im*jm),temp1u((im+1)*jm),temp1v(im*(jm+1)))
    if(mype==0) then
       allocate(landmask(im*jm),snow(im*jm),seaice(im*jm))
       allocate(t1st2d(im*jm),pt2t(im*jm))
       allocate(qs1st2d(im*jm,nlvl_qcheck))
       allocate(sfcprs(im*jm))
    endif
  
    if(mype == 0) write(6,*)' at 2 in wrwrfmassa'
  
    if(mype == 0) then
       write(filename,'("sigf",i2.2)')ifilesig(ntguessig)
       print *,'update ',trim(filename)
       open (lendian_in,file=filename,form='unformatted')
       open (lendian_out,file='siganl',form='unformatted')
       rewind lendian_in ; rewind lendian_out
    end if
  
  ! Convert analysis variables to MASS variables
  ! get pointer to relevant instance of cloud-related backgroud
    if (n_actual_clouds>0) then
       ier=0
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql', ges_qc, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi', ges_qi, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr', ges_qr, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qs', ges_qs, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qg', ges_qg, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr,istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qni',ges_qni,istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnc',ges_qnc,istatus );ier=ier+istatus
       if(dbz_exist) &
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'dbz',ges_dbz,istatus );ier=ier+istatus
       if (ier/=0) then
           write(6,*)'READ_WRF_MASS_BINARY_GUESS: getpointer failed, cannot do cloud analysis'
           if (l_hydrometeor_bkio .and. n_actual_clouds>0) call stop2(999)
       endif
    endif
  
    ier=0
    call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ps', ges_ps, istatus );ier=ier+istatus
    call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'u' , ges_u , istatus );ier=ier+istatus
    call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'v' , ges_v , istatus );ier=ier+istatus
    if(w_exist)then
      call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'w' , ges_w , istatus );ier=ier+istatus
    endif
    call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q' , ges_q , istatus );ier=ier+istatus
    if (ier/=0) then ! doesn't have to die - code can be generalized to bypass missing vars
        write(6,*)'wrwrfmassa_binary: getpointer failed, cannot retrieve ps,u,v,q'
        call stop2(999)
    endif
    
  ! Create all_loc from ges_*
    if(mype == 0) write(6,*)' at 3 in wrwrfmassa'
    all_loc=zero_single
    kt=i_t-1
    kq=i_q-1
    ku=i_u-1
    kv=i_v-1
    if(w_exist)kw=i_w-1
  ! for hydrometeors
    if(l_hydrometeor_bkio .and. n_actual_clouds>0) then
       kqc=i_qc-1
       kqi=i_qi-1
       kqr=i_qr-1
       kqs=i_qs-1
       kqg=i_qg-1
       kqnr=i_qnr-1
       kqni=i_qni-1
       kqnc=i_qnc-1
       if(dbz_exist.and.if_model_dbz)kdbz=i_dbz-1
       ktt=i_tt-1
    endif
    if ( laeroana_gocart ) then
       ier = 0
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'sulf',  ges_sulf,  istatus );ier=ier+istatus 
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'bc1',   ges_bc1,   istatus );ier=ier+istatus   
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'bc2',   ges_bc2,   istatus );ier=ier+istatus 
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'oc1',   ges_oc1,   istatus );ier=ier+istatus      
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'oc2',   ges_oc2,   istatus );ier=ier+istatus    
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'dust1', ges_dust1, istatus );ier=ier+istatus  
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'dust2', ges_dust2, istatus );ier=ier+istatus 
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'dust3', ges_dust3, istatus );ier=ier+istatus  
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'dust4', ges_dust4, istatus );ier=ier+istatus  
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'dust5', ges_dust5, istatus );ier=ier+istatus     
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'seas1', ges_seas1, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'seas2', ges_seas2, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'seas3', ges_seas3, istatus );ier=ier+istatus 
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'seas4', ges_seas4, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'p25',   ges_p25,   istatus );ier=ier+istatus
       if (ier/=0 .and. mype == 0) then
           write(6,*)'WRWRFMASSA_NETCDF: getpointer failed for gocart species'
       endif
       do iv = 1, n_gocart_var
          kchem(iv) = i_chem(iv)-1
       end do
    endif
  
    if ( wrf_pm2_5 ) then
       ier = 0
       call GSI_BundleGetPointer ( GSI_ChemGuess_Bundle(it), 'pm2_5',  ges_pm2_5,  istatus )
       ier=ier+istatus
       if (ier/=0 .and. mype == 0) then
           write(6,*)'WRWRFMASSA_NETCDF: getpointer failed for gocart pm2_5'
       endif
       iv=1
       kchem(iv) = i_chem(iv)-1
    endif
  
  
    q_integral=one
    q_integralc4h=zero_single
    do k=1,nsig
       deltasigma=eta1_ll(k)-eta1_ll(k+1)
       deltasigmac4h=eta2_ll(k)-eta2_ll(k+1)
       kt=kt+1
       kq=kq+1
       ku=ku+1
       kv=kv+1
  ! for hydrometeors
       if(l_hydrometeor_bkio .and. n_actual_clouds>0) then
          kqc=kqc+1
          kqi=kqi+1
          kqr=kqr+1
          kqs=kqs+1
          kqg=kqg+1
          kqnr=kqnr+1
          kqni=kqni+1
          kqnc=kqnc+1
          if(dbz_exist.and.if_model_dbz)kdbz=kdbz+1
          ktt=ktt+1
       endif
       if ( laeroana_gocart ) then
          do iv = 1, n_gocart_var
             kchem(iv) = kchem(iv)+1
          end do
       endif
  
       if (wrf_pm2_5) then
          iv=1
          kchem(iv) = kchem(iv)+1
       endif
  
  
       do i=1,lon2
          do j=1,lat2
             all_loc(j,i,ku)=ges_u(j,i,k)
             all_loc(j,i,kv)=ges_v(j,i,k)
  
  !          Convert sensible temperature to potential temperature
             work_prsl  = one_tenth*(aeta1_ll(k)*(r10*ges_ps(j,i)-pt_ll)+aeta2_ll(k)+pt_ll)
             work_prslk = (work_prsl/r100)**rd_over_cp_mass
             all_loc(j,i,kt) = ges_tsen(j,i,k,it)/work_prslk
  
  !          Convert specific humidity to mixing ratio
             all_loc(j,i,kq)= ges_q(j,i,k)/(one-ges_q(j,i,k))
             	
  ! for hydrometeors      
             if(l_hydrometeor_bkio .and. n_actual_clouds>0) then
                all_loc(j,i,kqc)=ges_qc(j,i,k)
                all_loc(j,i,kqi)=ges_qi(j,i,k)
                all_loc(j,i,kqr)=ges_qr(j,i,k)
                all_loc(j,i,kqs)=ges_qs(j,i,k)
                all_loc(j,i,kqg)=ges_qg(j,i,k)
                all_loc(j,i,kqnr)=ges_qnr(j,i,k)
                all_loc(j,i,kqni)=ges_qni(j,i,k)
                all_loc(j,i,kqnc)=ges_qnc(j,i,k)
                if(dbz_exist.and.if_model_dbz)all_loc(j,i,kdbz)=ges_dbz(j,i,k)
                all_loc(j,i,ktt)=ges_tten(j,i,k,it)
             endif
  
             if ( laeroana_gocart ) then
                all_loc(j,i,kchem(1))=ges_sulf(j,i,k)
                all_loc(j,i,kchem(2))=ges_bc1(j,i,k)
                all_loc(j,i,kchem(3))=ges_bc2(j,i,k)
                all_loc(j,i,kchem(4))=ges_oc1(j,i,k)
                all_loc(j,i,kchem(5))=ges_oc2(j,i,k)
                all_loc(j,i,kchem(6))=ges_dust1(j,i,k)
                all_loc(j,i,kchem(7))=ges_dust2(j,i,k)
                all_loc(j,i,kchem(8))=ges_dust3(j,i,k)
                all_loc(j,i,kchem(9))=ges_dust4(j,i,k)
                all_loc(j,i,kchem(10))=ges_dust5(j,i,k)
                all_loc(j,i,kchem(11))=ges_seas1(j,i,k)
                all_loc(j,i,kchem(12))=ges_seas2(j,i,k)
                all_loc(j,i,kchem(13))=ges_seas3(j,i,k)
                all_loc(j,i,kchem(14))=ges_seas4(j,i,k)
                if ( n_gocart_var>=15 ) all_loc(j,i,kchem(15))=ges_p25(j,i,k)
             endif
  
             if ( wrf_pm2_5 ) then
                all_loc(j,i,kchem(1))=ges_pm2_5(j,i,k)
             endif
             
  
             q_integral(j,i)=q_integral(j,i)+deltasigma* &
                  ges_q(j,i,k)/(one-ges_q(j,i,k))
             q_integralc4h(j,i)=q_integralc4h(j,i)+deltasigmac4h* &
                  ges_q(j,i,k)/(one-ges_q(j,i,k))
          end do
       end do
    end do

    if(w_exist) then
     kw=kw+1
     do i=1,lon2
        do j=1,lat2
          all_loc(j,i,kw)=ges_w_btlev(j,i,1,it) ! for w on the bottom not changed
        enddo
     enddo

     do k=1,nsig-1
       kw=kw+1
       do i=1,lon2
          do j=1,lat2
           all_loc(j,i,kw)=0.5*(ges_w(j,i,k)+ges_w(j,i,k+1))
          enddo
       enddo
     enddo

     kw=kw+1
     do i=1,lon2
        do j=1,lat2
          all_loc(j,i,kw)=ges_w_btlev(j,i,2,it) ! for w on the top not changed
        enddo
     enddo
    endif


    do i=1,lon2
       do j=1,lat2
          psfc_this=r10*ges_ps(j,i)   ! convert from cb to mb
          psfc_this_dry=pt_ll+((psfc_this-pt_ll)-q_integralc4h(j,i))/q_integral(j,i)
          all_loc(j,i,i_psfc)=r100*psfc_this_dry
       end do
    end do
    
    if(mype == 0) then
       read(lendian_in) regional_time0,nlon_regional0,nlat_regional0,nsig0,pt0,nsoil
       write(lendian_out) regional_time0,nlon_regional0,nlat_regional0,nsig0,pt0,nsoil
       read(lendian_in) aeta10,aeta20
       write(lendian_out) aeta10,aeta20
       read(lendian_in) eta10,eta20
       write(lendian_out) eta10,eta20
       read(lendian_in) glat0,dx_mc0
       write(lendian_out) glat0,dx_mc0
       read(lendian_in) glon0,dy_mc0
       write(lendian_out) glon0,dy_mc0
    end if
    
  ! Update psfc
    if(mype == 0) write(6,*)' at 6 in wrwrfmassa'
  
    allocate(tempa(itotsub),tempb(itotsub))
    if(mype == 0) read(lendian_in)temp1
    call strip(all_loc(:,:,i_psfc),strp)
    call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
         tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
    if(mype == 0) then
       call fill_mass_grid2t(temp1,im,jm,tempb,2)
       do i=1,iglobal
          tempa(i)=tempa(i)-tempb(i)
       end do
       call unfill_mass_grid2t(tempa,im,jm,temp1)
       write(lendian_out)temp1
       do i=1,im*jm
          sfcprs(i) = temp1(i)
        !  work_prsl = one_tenth*(aeta1_ll(1)*(temp1(i)/r100-pt_ll)+aeta2_ll(1)+pt_ll)
        !  pt2t(i)   = (work_prsl/r100)**rd_over_cp_mass
        !  qs1st2d(i)= work_prsl  ! save surface pressure for qs
       enddo
    end if
  
  !  FIS read/write
    if(mype == 0) then
       read(lendian_in)temp1
       write(lendian_out)temp1
    end if
  
  ! Update t
    kt=i_t-1
    do k=1,nsig
       kt=kt+1
       if(mype == 0) read(lendian_in)temp1
       call strip(all_loc(:,:,kt),strp)
       call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
            tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
       if(mype == 0) then
          call fill_mass_grid2t(temp1,im,jm,tempb,2)
          do i=1,iglobal
             tempa(i)=tempa(i)-tempb(i)
          end do
          call unfill_mass_grid2t(tempa,im,jm,temp1)
          write(lendian_out)temp1
          if(k>=1 .and. k<=nlvl_qcheck) then
             ! calculate moisture saturation at model 1st level  
             allocate(tmp2da(im,jm),tmp2db(im,jm),tmp2dc(im,jm))
             ij=0
             do j=1,jm
                do i=1,im
                   ij=ij+1
                   work_prsl = one_tenth*(aeta1_ll(k)*(sfcprs(ij)/r100-pt_ll)+aeta2_ll(k)+pt_ll)
                   tmp2dc(i,j)=work_prsl                                       ! pressure
                   tmp2db(i,j)=temp1(ij)*(work_prsl/r100)**rd_over_cp_mass     ! T in K
                   if(k==1) t1st2d(ij)  = tmp2db(i,j)
                enddo
             enddo
             call genqsat(tmp2da,tmp2db,tmp2dc,im,jm,1,.false.,0)
             ij=0
             do j=1,jm
                do i=1,im
                   ij=ij+1
                   qs1st2d(ij,k)=tmp2da(i,j) ! moisture saturation 
                enddo
             enddo
             deallocate(tmp2da,tmp2db,tmp2dc)
             if(k==1) pt2t=t1st2d   ! save a copy
          endif
       end if
    end do
  
  ! Update q
    kq=i_q-1
    do k=1,nsig
       kq=kq+1
       if(mype == 0) read(lendian_in)temp1
       call strip(all_loc(:,:,kq),strp)
       call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
            tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
       if(mype == 0) then
          call fill_mass_grid2t(temp1,im,jm,tempb,2)
          do i=1,iglobal
             tempa(i)=tempa(i)-tempb(i)
          end do
          if(k>=1 .and. k<=nlvl_qcheck) then
             call unfill_mass_grid2t_drycheck(tempa,im,jm,temp1,qs1st2d(:,k))
          else
             call unfill_mass_grid2t(tempa,im,jm,temp1)
          endif
          write(lendian_out)temp1
       end if
    end do
  
  ! Update u
    ku=i_u-1
    do k=1,nsig
       ku=ku+1
       if(mype == 0) read(lendian_in)temp1u
       call strip(all_loc(:,:,ku),strp)
       call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
            tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
       if(mype == 0) then
          call fill_mass_grid2u(temp1u,im,jm,tempb,2)
          do i=1,iglobal
             tempa(i)=tempa(i)-tempb(i)
          end do
          call unfill_mass_grid2u(tempa,im,jm,temp1u)
          write(lendian_out)temp1u
       end if
    end do
  
  ! Update v
    kv=i_v-1
    do k=1,nsig
       kv=kv+1
       if(mype == 0) read(lendian_in)temp1v
       call strip(all_loc(:,:,kv),strp)
       call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
            tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
       if(mype == 0) then
          call fill_mass_grid2v(temp1v,im,jm,tempb,2)
          do i=1,iglobal
             tempa(i)=tempa(i)-tempb(i)
          end do
          call unfill_mass_grid2v(tempa,im,jm,temp1v)
          write(lendian_out)temp1v
       end if
    end do

    if(w_exist) then
    ! Update w
      kw=i_w-1
      do k=1,nsig+1
        kw=kw+1
        if(mype == 0) read(lendian_in)temp1
        call strip(all_loc(:,:,kw),strp)
        call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
             tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
        if(mype == 0) then
           call fill_mass_grid2t(temp1,im,jm,tempb,2)
           do i=1,iglobal
              tempa(i)=tempa(i)-tempb(i)
           end do
           call unfill_mass_grid2t(tempa,im,jm,temp1)
           write(lendian_out)temp1
        end if
      end do
    endif ! for w_exist

    
  ! Load updated skin temperature array if writing out to analysis file
    if (update_regsfc) then
       do i=1,lon2
          do j=1,lat2
             all_loc(j,i,i_sst)=dsfct(j,i,ntguessfc)
             all_loc(j,i,i_skt)=dsfct(j,i,ntguessfc)
          end do
       end do
    end if
  
  ! for soil nudging
    if(l_gsd_soilTQ_nudge) then
       ier=0
       call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'smoist',ges_smois_it,istatus)
       ier=ier+istatus
       call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'tslb' ,ges_tslb_it ,istatus)
       ier=ier+istatus
       if (ier/=0) call die(trim(myname),'cannot get pointers for tslb/smois, ier =',ier)
       ktslb=i_tslb-1
       ksmois=i_smois-1
       do k=1,nsig_soil
          ktslb=ktslb+1
          ksmois=ksmois+1
          do i=1,lon2
             do j=1,lat2
                all_loc(j,i,ktslb)=ges_tslb_it(j,i,k)
                all_loc(j,i,ksmois)=ges_smois_it(j,i,k)
             end do
          end do
       end do
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'tsoil', ges_soilt1, istatus );ier=ier+istatus
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'tskn' , ges_tsk , istatus );ier=ier+istatus
       if (ier/=0) then ! doesn't have to die - code can be generalized to bypass missing vars
           write(6,*)'wrwrfmassa_netcdf: getpointer failed, cannot retrieve soil stuff'
           call stop2(999)
       endif
       do i=1,lon2
          do j=1,lat2
             all_loc(j,i,i_skt)=ges_tsk(j,i)
             all_loc(j,i,i_soilt1)=ges_soilt1(j,i)
          end do
       end do
    endif ! l_gsd_soilTQ_nudge

    if(i_use_2mt4b > 0 ) then
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'th2m', ges_th2,istatus );ier=ier+istatus
       if (ier/=0) then ! doesn't have to die - code can be generalized to bypass
           write(6,*)'wrwrfmassa_netcdf: getpointer failed, cannot retrieve th2'
           call stop2(999)
       endif
       do i=1,lon2
          do j=1,lat2
             all_loc(j,i,i_th2)=ges_th2(j,i)
          end do
       end do
    endif
    if(i_use_2mt4b > 0 .or. i_use_2mq4b > 0) then
       call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q2m',   ges_q2, istatus );ier=ier+istatus
       if (ier/=0) then ! doesn't have to die - code can be generalized to bypass missing vars          
           write(6,*)'wrwrfmassa_netcdf: getpointer failed, cannot retrieve q2m'
           call stop2(999)
       endif
       do i=1,lon2
          do j=1,lat2
  ! Convert 2m specific humidity to mixing ratio
             all_loc(j,i,i_q2)=ges_q2(j,i)/(one-ges_q2(j,i))
          end do
       end do
    endif
  
    if(mype == 0) then
  ! SM   This is landmask
       read(lendian_in)temp1
       landmask=temp1
       write(lendian_out)temp1
  ! SICE
       read(lendian_in)temp1
       seaice=temp1
       write(lendian_out)temp1
    end if
  
  ! SST
    if(update_regsfc) then
       if(mype == 0) read(lendian_in)temp1
       if (mype==0)write(6,*)' at 9.1 in wrwrfmassa,max,min(temp1)=',maxval(temp1),minval(temp1)
       call strip(all_loc(:,:,i_sst),strp)
       call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
            tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
       if(mype == 0) then
          if(mype == 0) write(6,*)' at 9.2 in wrwrfmassa,max,min(tempa)=',maxval(tempa),minval(tempa)
          call fill_mass_grid2t(temp1,im,jm,tempb,2)
          do i=1,iglobal
             if(tempb(i) < (r225)) then
                tempa(i)=zero_single
             else
                tempa(i)=tempa(i)-tempb(i)
             end if
          end do
          if(mype == 0) write(6,*)' at 9.4 in wrwrfmassa,max,min(tempa)=',maxval(tempa),minval(tempa)
          call unfill_mass_grid2t(tempa,im,jm,temp1)
          write(6,*)' at 9.6 in wrwrfmassa,max,min(temp1)=',maxval(temp1),minval(temp1)
          write(lendian_out)temp1
       end if     !endif mype==0
    else
       if(mype==0) then
          read(lendian_in)temp1
          write(lendian_out)temp1
       end if
    end if   !end if check updatesfc
    
  ! REST OF FIELDS
    if(l_gsd_soilTQ_nudge) then
       if (mype == 0) then
          do k=4,6
             read(lendian_in)temp1
             write(lendian_out)temp1
          end do
  ! SM   This is snow
          read(lendian_in)temp1
          write(lendian_out)temp1
          snow=temp1
          do k=8,9
             read(lendian_in)temp1
             write(lendian_out)temp1
          end do
       end if
  ! Update smois
       ksmois=i_smois-1
       do k=1,nsig_soil
          ksmois=ksmois+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,ksmois),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             i_snowT_check=2
             call unfill_mass_grid2t_ldmk(tempa,im,jm,temp1,landmask,&
                                          snow,seaice,t1st2d,i_snowT_check)
  ! make sure soil moisture is larger than soilmoistmin (0.002 (sand)).
             do i=1,im*jm
                temp1(i) = max(temp1(i),soilmoistmin)
             enddo
             write(lendian_out)temp1
          end if
       end do
  ! update TSLB
       ktslb=i_tslb-1
       do k=1,nsig_soil
          ktslb=ktslb+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,ktslb),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             i_snowT_check=3
             if(k==1) t1st2d=temp1-t1st2d
             call unfill_mass_grid2t_ldmk(tempa,im,jm,temp1,landmask, &
                                          snow,seaice,t1st2d,i_snowT_check)
             write(lendian_out)temp1
          end if
       end do
    else
       if (mype == 0) then
          do k=4,11
             read(lendian_in)temp1
             write(lendian_out)temp1
          end do
          snow=0
       end if
    endif !  l_gsd_soilTQ_nudge
    
  ! Update SKIN TEMP
    if(update_regsfc .or. l_gsd_soilTQ_nudge) then
       if(mype == 0) read(lendian_in)temp1
       if (mype==0)write(6,*)' at 10.0 in wrwrfmassa,max,min(temp1)=',maxval(temp1),minval(temp1)
       call strip(all_loc(:,:,i_skt),strp)
       call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
            tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
       if(mype == 0) then
          call fill_mass_grid2t(temp1,im,jm,tempb,2)
          do i=1,iglobal
             if(tempb(i) < (r225)) then
                tempa(i)=zero_single
             else 
                tempa(i)=tempa(i)-tempb(i)
             end if
          end do
          i_snowT_check=1
          t1st2d=temp1-pt2t
          call unfill_mass_grid2t_ldmk(tempa,im,jm,temp1,landmask, &
                                       snow,seaice,t1st2d,i_snowT_check)
          write(lendian_out)temp1
       end if
    else
       if (mype == 0) then
          read(lendian_in)temp1
          write(lendian_out)temp1
       end if
    end if
  
  ! Update Q2
    if( i_use_2mq4b>0 .or. i_use_2mt4b > 0) then
       if(mype == 0) read(lendian_in)temp1
       if (mype==0)write(6,*)' at 10.11 in wrwrfmassa,max,min(temp1)=', &
                                           maxval(temp1),minval(temp1)
       call strip(all_loc(:,:,i_q2),strp)
       tempa=zero_single
       call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
            tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
       if(mype == 0) then
          write(6,*)' at 10.12 in wrwrfmassa,max,min(tempa)=', &
                            maxval(tempa),minval(tempa)
          call fill_mass_grid2t(temp1,im,jm,tempb,2)
          do i=1,iglobal
             tempa(i)=tempa(i)-tempb(i)
          end do
          write(6,*)' at 10.13 in wrwrfmassa,max,min(tempa)=', &
                             maxval(tempa),minval(tempa)
          call unfill_mass_grid2t_drycheck(tempa,im,jm,temp1,qs1st2d(:,1))
          write(6,*)' at 10.14 in wrwrfmassa,max,min(temp1)=', &
                             maxval(temp1),minval(temp1)
          write(lendian_out)temp1
       end if     !endif mype==0
    endif
  
    if(l_gsd_soilTQ_nudge) then
  ! update soilt1
       if(mype == 0) read(lendian_in)temp1
       if (mype==0)write(6,*)' at 10.1 in wrwrfmassa,max,min(temp1)=',maxval(temp1),minval(temp1)
       call strip(all_loc(:,:,i_soilt1),strp)
       tempa=zero_single
       call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
            tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
       if(mype == 0) then
           write(6,*)' at 10.2 in wrwrfmassa,max,min(tempa)=',maxval(tempa),minval(tempa)
          call fill_mass_grid2t(temp1,im,jm,tempb,2)
          do i=1,iglobal
             tempa(i)=tempa(i)-tempb(i)
          end do
          write(6,*)' at 10.3 in wrwrfmassa,max,min(tempa)=',maxval(tempa),minval(tempa)
          i_snowT_check=4
          call unfill_mass_grid2t_ldmk(tempa,im,jm,temp1,landmask, &
                                      snow,seaice,t1st2d,i_snowT_check)
          write(6,*)' at 10.4 in wrwrfmassa,max,min(temp1)=',maxval(temp1),minval(temp1)
          write(lendian_out)temp1
       end if     !endif mype==0
    endif  !  l_gsd_soilTQ_nudge
  ! update TH2
    if( i_use_2mt4b>0 ) then
       if(mype == 0) read(lendian_in)temp1
       if (mype==0)write(6,*)' at 10.5 in wrwrfmassa,max,min(temp1)=',maxval(temp1),minval(temp1)
       call strip(all_loc(:,:,i_th2),strp)
       call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
            tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
       if(mype == 0) then
           write(6,*)' at 10.6 in wrwrfmassa,max,min(tempa)=',maxval(tempa),minval(tempa)
          call fill_mass_grid2t(temp1,im,jm,tempb,2)
          do i=1,iglobal
             if(tempb(i) < (r100)) then
                tempa(i)=zero_single
             else
                tempa(i)=tempa(i)-tempb(i)
             end if
          end do
          write(6,*)' at 10.7 in wrwrfmassa,max,min(tempa)=',maxval(tempa),minval(tempa)
          call unfill_mass_grid2t(tempa,im,jm,temp1)
          write(6,*)' at 10.8 in wrwrfmassa,max,min(temp1)=',maxval(temp1),minval(temp1)
          write(lendian_out)temp1
       end if     !endif mype==0
    endif  ! i_use_2mt4b>0
  !
  ! for saving cloud analysis results
    if(l_hydrometeor_bkio .and. n_actual_clouds>0) then
  ! Update qc
       kqc=i_qc-1
       do k=1,nsig
          kqc=kqc+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,kqc),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             call unfill_mass_grid2t(tempa,im,jm,temp1)
             write(lendian_out)temp1
          end if
       end do
  
  ! Update qr     
       kqr=i_qr-1
       do k=1,nsig
          kqr=kqr+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,kqr),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             call unfill_mass_grid2t(tempa,im,jm,temp1)
             write(lendian_out)temp1
          end if
       end do
  
  ! Update qs     
       kqs=i_qs-1
       do k=1,nsig
          kqs=kqs+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,kqs),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             call unfill_mass_grid2t(tempa,im,jm,temp1)
             write(lendian_out)temp1
          end if
       end do
  
  ! Update qi     
       kqi=i_qi-1
       do k=1,nsig
          kqi=kqi+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,kqi),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             call unfill_mass_grid2t(tempa,im,jm,temp1)
             write(lendian_out)temp1
          end if
       end do
  
  ! Update qg     
       kqg=i_qg-1
       do k=1,nsig
          kqg=kqg+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,kqg),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             call unfill_mass_grid2t(tempa,im,jm,temp1)
             write(lendian_out)temp1
          end if
       end do
  
  ! Update qnr     
       kqnr=i_qnr-1
       do k=1,nsig
          kqnr=kqnr+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,kqnr),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             call unfill_mass_grid2t(tempa,im,jm,temp1)
             write(lendian_out)temp1
          end if
       end do
   
  ! Update qni     
       kqni=i_qni-1
       do k=1,nsig
          kqni=kqni+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,kqni),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             call unfill_mass_grid2t(tempa,im,jm,temp1)
             write(lendian_out)temp1
          end if
       end do
   
  ! Update qnc     
       kqnc=i_qnc-1
       do k=1,nsig
          kqnc=kqnc+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,kqnc),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             call unfill_mass_grid2t(tempa,im,jm,temp1)
             write(lendian_out)temp1
          end if
       end do
   
       if(dbz_exist.and.if_model_dbz)then
       ! Update refl_10cm
            kdbz=i_dbz-1
            do k=1,nsig
               kdbz=kdbz+1
               if(mype == 0) read(lendian_in)temp1
               call strip(all_loc(:,:,kdbz),strp)
               call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
                    tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
               if(mype == 0) then
                  call fill_mass_grid2t(temp1,im,jm,tempb,2)
                  do i=1,iglobal
                     tempa(i)=tempa(i)-tempb(i)
                  end do
                  call unfill_mass_grid2t(tempa,im,jm,temp1)
                  write(lendian_out)temp1
               end if
            end do
       end if
   
  ! Update tten     
       ktt=i_tt-1
       do k=1,nsig
          ktt=ktt+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,ktt),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             call unfill_mass_grid2t(tempa,im,jm,temp1)
             write(lendian_out)temp1
          end if
       end do
  
    endif    ! l_hydrometeor_bkio
  
    if ( laeroana_gocart ) then
       do iv = 1, n_gocart_var
          kchem(iv)=i_chem(iv)-1
          do k=1,nsig
             tempa=0.0
             kchem(iv)=kchem(iv)+1
             if(mype == 0) read(lendian_in)temp1
             call strip(all_loc(:,:,kchem(iv)),strp)
             call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
                  tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
             if(mype == 0) then
                call fill_mass_grid2t(temp1,im,jm,tempb,2)
                do i=1,iglobal
                   tempa(i)=tempa(i)-tempb(i)
                end do
                call unfill_mass_grid2t(tempa,im,jm,temp1)
                write(lendian_out)temp1
             end if
          end do
       end do
       deallocate(i_chem)
       deallocate(kchem)
    endif
  
    if ( wrf_pm2_5 ) then
       iv=1
       kchem(iv)=i_chem(iv)-1
       do k=1,nsig
          tempa=0.0
          kchem(iv)=kchem(iv)+1
          if(mype == 0) read(lendian_in)temp1
          call strip(all_loc(:,:,kchem(iv)),strp)
          call mpi_gatherv(strp,ijn(mype+1),mpi_real4, &
               tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror)
          if(mype == 0) then
             call fill_mass_grid2t(temp1,im,jm,tempb,2)
             do i=1,iglobal
                tempa(i)=tempa(i)-tempb(i)
             end do
             call unfill_mass_grid2t(tempa,im,jm,temp1)
             write(lendian_out)temp1
          end if
       end do
       deallocate(i_chem)
       deallocate(kchem)
    endif
  
    if (mype==0) then
       close(lendian_in)
       close(lendian_out)
    endif
  
    deallocate(all_loc)
    deallocate(strp)
    deallocate(temp1)
    deallocate(temp1u)
    deallocate(temp1v)
    deallocate(tempa)
    deallocate(tempb)
    if(mype==0) then
       deallocate(landmask)
       deallocate(snow)
       deallocate(seaice)
       deallocate(t1st2d)
       deallocate(pt2t)
       deallocate(qs1st2d)
       deallocate(sfcprs)
    endif
    
  end subroutine wrwrfmassa_netcdf_wrf

  subroutine update_start_date(chdrbuf,iyear,imonth,iday,ihour,iminute,isecond)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    update_start_date  update start date on wrf file
  !   prgmmr: parrish          org: np22                date: 2004-11-29
  !
  ! abstract: update date record in START_DATE header record
  !
  ! program history log:
  !   2004-11-29  parrish
  !   2006-04-24  middlecoff - Put error handling in search for START-TIME
  !                            and changed search loop increment from 4 to one
  !
  !   input argument list:
  !     chdrbuf  - 2048 byte wrf header record containing "START_DATE"
  !     iyear    - analysis time year
  !     imonth   - analysis time month
  !     iday     - analysis time day
  !     ihour    - analysis time hour
  !     iminute  - analysis time minute
  !     isecond  - analysis time second
  !
  !   output argument list:
  !     chdrbuf  - output header record with new analysis date
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm RS/6000 SP
  !
  !$$$
  
  !    update date record in START_DATE header record
  
    use kinds, only: i_kind
    implicit none
  
    character(1)   ,intent(inout) :: chdrbuf(2048)
    integer(i_kind),intent(in   ) :: iyear,imonth,iday,ihour,iminute,isecond
  
    character(2) c_two
    character(4) c_four
    character(1) c2(2),d2(2),c4(4),d4(4)
    equivalence (c2(1),c_two),(c4(1),c_four)
    integer(i_kind) i,ibegin,j
  
    ibegin=0
    do i=1,1932
       if(chdrbuf(i)=='S'.and.chdrbuf(i+4)=='T'.and.chdrbuf(i+8)=='A' &
         .and.chdrbuf(i+12)=='R'.and.chdrbuf(i+16)=='T'.and.chdrbuf(i+20)=='_' &
         .and.chdrbuf(i+24)=='D'.and.chdrbuf(i+28)=='A'.and.chdrbuf(i+32)=='T' &
         .and.chdrbuf(i+36)=='E') then
          ibegin=i+44
          exit
       end if
    end do
    if(ibegin==0) then
       write(6,*)'UPDATE_START_DATE:  ***ERROR***  ibegin = ',ibegin
       write(6,*)'chdrbuf=',chdrbuf
       call stop2(24)
    endif
  
    i=ibegin-4
    write(c_four,'(i4.4)')iyear
    do j=1,4
       i=i+4
       d4(j)=chdrbuf(i)
       chdrbuf(i)=c4(j)
    end do
    write(6,*) 'UPDATE_START_DATE:  old year, new year =',d4,' , ',c4
  
  !          skip "-"
    i=i+4
  
    write(c_two,'(i2.2)')imonth
    do j=1,2
       i=i+4
       d2(j)=chdrbuf(i)
       chdrbuf(i)=c2(j)
    end do
    write(6,*) 'UPDATE_START_DATE:  old month, new month =',d2,' , ',c2
  
  !          skip "-"
    i=i+4
  
    write(c_two,'(i2.2)')iday
    do j=1,2
       i=i+4
       d2(j)=chdrbuf(i)
       chdrbuf(i)=c2(j)
    end do
    write(6,*) 'UPDATE_START_DATE:  old day, new day =',d2,' , ',c2
  
  !          skip "_"
    i=i+4
  
    write(c_two,'(i2.2)')ihour
    do j=1,2
       i=i+4
       d2(j)=chdrbuf(i)
       chdrbuf(i)=c2(j)
    end do
    write(6,*) 'UPDATE_START_DATE:  old hour, new hour =',d2,' , ',c2
  
  !          skip ":"
    i=i+4
  
    write(c_two,'(i2.2)')iminute
    do j=1,2
       i=i+4
       d2(j)=chdrbuf(i)
       chdrbuf(i)=c2(j)
    end do
    write(6,*) 'UPDATE_START_DATE:  old minute, new minute =',d2,' , ',c2
  
  !          skip ":"
    i=i+4
  
    write(c_two,'(i2.2)')isecond
    do j=1,2
       i=i+4
       d2(j)=chdrbuf(i)
       chdrbuf(i)=c2(j)
    end do
    write(6,*) 'UPDATE_START_DATE:  old second, new second =',d2,' , ',c2
  
  end subroutine update_start_date
end module wrwrfmassa_mod