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