#ifdef WRF
subroutine read_wrf_mass_binary_guess(mype)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_wrf_mass_guess read wrf_mass interface file
! prgmmr: parrish org: np22 date: 2003-09-05
!
! abstract: in place of read_guess for global application, read guess
! from regional model, in this case the wrf mass core model.
! This version reads a binary file created
! in a previous step that interfaces with the wrf infrastructure.
! A later version will read directly from the wrf restart file.
! The guess is read in by complete horizontal fields, one field
! per processor, in parallel. Each horizontal input field is
! converted from the staggered c-grid to an unstaggered a-grid.
! On the c-grid, u is shifted 1/2 point in the negative x direction
! and v 1/2 point in the negative y direction, but otherwise the
! three grids are regular. When the fields are read in, nothing
! is done to mass variables, but wind variables are interpolated to
! mass points.
!
! program history log:
! 2004--7-15 parrish
! 2004-08-02 treadon - add only to module use, add intent in/out
! 2004-09-10 parrish - correct error in land-sea mask interpretation
! 2004-11-08 parrish - change to mpi-io for binary file format
! 2004-12-15 treadon - remove variable mype from call load_geop_hgt,
! rename variable wrf_ges_filename as wrfges
! 2005-02-17 todling - ifdef'ed wrf code out
! 2005-02-23 wu - setup for qoption=2 and output stats of RH
! 2005-04-01 treadon - add initialization of ges_oz, prsi_oz; comestic format changes
! 2005-05-27 parrish - add call get_derivatives
! 2005-11-21 kleist - new call to genqsat
! 2005-11-21 derber - make qoption=1 work same as qoption=2
! 2005-11-29 derber - remove external iteration dependent calculations
! 2005-12-15 treadon - remove initialization of certain guess arrays (done elsewhere)
! 2006-02-02 treadon - remove load_prsges,load_geop_hgt,prsl,prslk (not used)
! 2006-02-15 treadon - convert moisture mixing ratio to specific humidity
! 2006-03-07 treadon - convert guess potential temperature to vritual temperature
! 2006-04-06 middlecoff - changed nfcst from 11 to lendian_in
! 2006-07-28 derber - include sensible temperatures
! 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-04-16 safford - rm unused uses
!
! input argument list:
! mype - pe number
!
! NOTES: need to pay special attention to various surface fields to make sure
! they are correct (check units). fact10 needs to be computed from
! 10m wind, which is included in wrf mass interface file.
!
! Cloud water currently set to zero. Need to fix before doing precip
! assimilation--wait until global adopts Brad Ferrier's multi-species
! scheme??
!
! Ozone currently set to zero. No ozone variable in wrf mass core--climatology
! instead. Do we use climatology?
!
! No background bias yet. (biascor ignored)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,r_single,i_long,i_llong,i_kind
use mpimod, only: mpi_sum,mpi_integer,mpi_comm_world,npe,ierror, &
mpi_offset_kind,mpi_info_null,mpi_mode_rdonly,mpi_status_size
use guess_grids, only: ges_z,ges_ps,ges_tv,ges_q,ges_u,ges_v,&
fact10,soil_type,veg_frac,veg_type,sfct,sno,soil_temp,soil_moi,&
isli,nfldsig,ifilesig,ges_tsen
use gridmod, only: lat2,lon2,nlat_regional,nlon_regional,&
nsig,eta1_ll,pt_ll,itotsub,aeta1_ll
use constants, only: izero,ione,zero,one,grav,fv,zero_single,rd_over_cp_mass,one_tenth,h300
use gsi_io, only: lendian_in
implicit none
! Declare passed variables
integer(i_kind),intent(in):: mype
! Declare local parameters
real(r_kind),parameter:: r0_01 = 0.01_r_kind
real(r_kind),parameter:: r10 = 10.0_r_kind
real(r_kind),parameter:: r100 = 100.0_r_kind
! Declare local variables
integer(i_kind) kt,kq,ku,kv
! MASS variable names stuck in here
integer(i_kind) mfcst
! other internal variables
real(r_single),allocatable::tempa(:,:)
real(r_single),allocatable::temp1(:,:),temp1u(:,:),temp1v(:,:)
real(r_single),allocatable::all_loc(:,:,:)
integer(i_kind),allocatable::igtype(:),kdim(:),kord(:)
integer(kind=mpi_offset_kind),allocatable::offset(:)
integer(kind=mpi_offset_kind) this_offset
integer(i_kind),allocatable::length(:)
integer(i_kind) this_length
character(6) filename
character(9) wrfges
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_mub,i_mu,i_fis,i_t,i_q,i_u,i_v,i_sno,i_u10,i_v10,i_smois,i_tslb
integer(i_kind) i_sm,i_xice,i_sst,i_tsk,i_ivgtyp,i_isltyp,i_vegfrac
integer(i_kind) isli_this
real(r_kind) psfc_this,psfc_this_dry,sm_this,xice_this
real(r_kind),dimension(lat2,lon2):: q_integral
real(r_kind),dimension(lat2,lon2,nsig):: ges_pot
integer(i_kind) num_doubtful_sfct,num_doubtful_sfct_all
real(r_kind) deltasigma
integer(i_llong) n_position
integer(i_kind) iskip,ksize,jextra,nextra
integer(i_kind) status(mpi_status_size)
integer(i_kind) jbegin(0:npe),jend(0:npe-ione),jend2(0:npe-ione)
integer(i_kind) kbegin(0:npe),kend(0:npe-ione)
integer(i_long),allocatable:: ibuf(:,:)
integer(i_long),allocatable:: jbuf(:,:,:)
integer(i_long) dummy9(9)
real(r_single) pt_regional_single
real(r_kind):: work_prsl,work_prslk
integer(i_kind) iadd
character(132) memoryorder
! WRF MASS input grid dimensions in module gridmod
! These are the following:
! im -- number of x-points on C-grid
! jm -- number of y-points on C-grid
! lm -- number of vertical levels ( = nsig for now)
num_doubtful_sfct=izero
im=nlon_regional
jm=nlat_regional
lm=nsig
if(jm<=npe)then
write(6,*)' in read_wrf_mass_binary_guess, jm <= npe, ',&
'so program will end.'
call stop2(1)
endif
if(mype==izero) write(6,*)' in read_wrf_mass_binary_guess, im,jm,lm=',im,jm,lm
! Following is for convenient WRF MASS input
num_mass_fields=21_i_kind+5_i_kind*lm+2_i_kind
num_loc_groups=num_mass_fields/npe
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, lm =",i6)')lm
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, num_mass_fields=",i6)')num_mass_fields
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, nfldsig =",i6)')nfldsig
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, npe =",i6)')npe
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, num_loc_groups=",i6)')num_loc_groups
allocate(offset(num_mass_fields))
allocate(igtype(num_mass_fields),kdim(num_mass_fields),kord(num_mass_fields))
allocate(length(num_mass_fields))
! 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
! igtype = 0 for dummy (nothing to read)
! offset is the byte count preceding each record to be read from the wrf binary file.
! used as individual file pointers by mpi_file_read_at
do it=1,nfldsig
write(filename,'("sigf",i2.2)')ifilesig(it)
open(lendian_in,file=filename,form='unformatted') ; rewind lendian_in
write(6,*)'READ_WRF_MASS_BINARY_GUESS: open lendian_in=',lendian_in,&
' to filename=',filename,' on it=',it
if(mype == izero) write(6,*)'READ_WRF_MASS_OFFSET_FILE: open lendian_in=',lendian_in,' to file=',filename
read(lendian_in) dummy9,pt_regional_single
write(6,*)'READ_WRF_MASS_BINARY_GUESS: dummy9=',dummy9
do iskip=2,5
read(lendian_in)
end do
read(lendian_in) wrfges
read(lendian_in) ! n_position ! offset for START_DATE record
write(6,*)'READ_WRF_MASS_BINARY_GUESS: read wrfges,n_position= ',wrfges,' ',n_position
i=izero
i=i+ione ; i_mub=i ! mub
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ione
if(mype == izero) write(6,*)' mub, i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
i=i+ione ; i_mu =i ! mu
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ione
if(mype == izero) write(6,*)' mu, i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
i_fis=i+ione ! sfc geopotential
read(lendian_in) n_position,memoryorder
do k=1,lm+ione
i=i+ione
if(trim(memoryorder)=='XZY') then
iadd=izero
kord(i)=lm+ione
else
iadd=(k-ione)*im*jm*4
kord(i)=ione
end if
offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=lm+ione
if(mype == izero.and.k==ione) write(6,*)' sfc geopot i,igtype(i),offset(i),kord(i) = ', &
i,igtype(i),offset(i),kord(i)
end do
i_t=i+ione
read(lendian_in) n_position,memoryorder
do k=1,lm
i=i+ione ! theta(k) (pot temp)
if(trim(memoryorder)=='XZY') then
iadd=izero
kord(i)=lm
else
iadd=(k-ione)*im*jm*4
kord(i)=ione
end if
offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=lm
if(mype == izero.and.k==ione) write(6,*)' temp i,igtype(i),offset(i),kord(i) = ', &
i,igtype(i),offset(i),kord(i)
end do
i_q=i+ione
read(lendian_in) n_position,memoryorder
do k=1,lm
i=i+ione ! q(k)
if(trim(memoryorder)=='XZY') then
iadd=izero
kord(i)=lm
else
iadd=(k-ione)*im*jm*4
kord(i)=ione
end if
offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=lm
if(mype == izero.and.k==ione) write(6,*)' q i,igtype(i),offset(i),kord(i) = ', &
i,igtype(i),offset(i),kord(i)
end do
i_u=i+ione
read(lendian_in) n_position,memoryorder
do k=1,lm
i=i+ione ! u(k)
if(trim(memoryorder)=='XZY') then
iadd=izero
kord(i)=lm
else
iadd=(k-ione)*(im+ione)*jm*4
kord(i)=ione
end if
offset(i)=n_position+iadd
igtype(i)=2_i_kind ; kdim(i)=lm
length(i)=(im+ione)*jm
if(mype == izero.and.k==ione) write(6,*)' u i,igtype(i),offset(i),kord(i) = ', &
i,igtype(i),offset(i),kord(i)
end do
i_v=i+ione
read(lendian_in) n_position,memoryorder
do k=1,lm
i=i+ione ! v(k)
if(trim(memoryorder)=='XZY') then
iadd=izero
kord(i)=lm
else
iadd=(k-ione)*im*(jm+ione)*4
kord(i)=ione
end if
offset(i)=n_position+iadd ; length(i)=im*(jm+ione) ; igtype(i)=3_i_kind ; kdim(i)=lm
if(mype == izero.and.k==ione) write(6,*)' v i,igtype(i),offset(i),kord(i) = ', &
i,igtype(i),offset(i),kord(i)
end do
i=i+ione ; i_sm=i ! landmask
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ione
if(mype == izero) write(6,*)' landmask i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
i=i+ione ; i_xice=i ! xice
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ione
if(mype == izero) write(6,*)' xice i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
i=i+ione ; i_sst=i ! sst
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ione
if(mype == izero) write(6,*)' sst i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
i=i+ione ; i_ivgtyp=i ! ivgtyp
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=-ione ; kdim(i)=ione
if(mype == izero) write(6,*)' ivgtyp i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
i=i+ione ; i_isltyp=i ! isltyp
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=-ione ; kdim(i)=ione
if(mype == izero) write(6,*)' isltyp i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
i=i+ione ; i_vegfrac=i ! vegfrac
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ione
if(mype == izero) write(6,*)' vegfrac i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
i=i+ione ; i_sno=i ! sno
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ione
if(mype == izero) write(6,*)' sno i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
i=i+ione ; i_u10=i ! u10
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ione
if(mype == izero) write(6,*)' u10 i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
i=i+ione ; i_v10=i ! v10
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ione
if(mype == izero) write(6,*)' v10 i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
i=i+ione ; i_smois=i ! smois
read(lendian_in) n_position,ksize,memoryorder
if(trim(memoryorder)=='XZY') then
kord(i)=ksize
else
kord(i)=ione
end if
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ksize
if(mype == izero) write(6,*)' smois i,igtype(i),offset(i),kord(i) = ', &
i,igtype(i),offset(i),kord(i)
do k=2,ksize
i=i+ione
if(trim(memoryorder)=='XZY') then
iadd=izero
kord(i)=ksize
else
iadd=(k-ione)*im*jm*4
kord(i)=ione
end if
offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ksize
end do
i=i+ione ; i_tslb=i ! tslb
read(lendian_in) n_position,ksize,memoryorder
if(trim(memoryorder)=='XZY') then
kord(i)=ksize
else
kord(i)=ione
end if
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ksize
if(mype == izero) write(6,*)' tslb i,igtype(i),offset(i),kord(i) = ', &
i,igtype(i),offset(i),kord(i)
do k=2,ksize
i=i+ione
if(trim(memoryorder)=='XZY') then
iadd=izero
kord(i)=ksize
else
iadd=(k-ione)*im*jm*4
kord(i)=ione
end if
offset(i)=n_position+iadd ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ksize
if(mype == izero) write(6,*)' i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
end do
i=i+ione ; i_tsk=i ! tsk
read(lendian_in) n_position
offset(i)=n_position ; length(i)=im*jm ; igtype(i)=ione ; kdim(i)=ione
if(mype == izero) write(6,*)' tsk i,igtype(i),offset(i) = ',i,igtype(i),offset(i)
close(lendian_in)
! End of stuff from MASS restart file
! 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)=ione
if(nextra > izero) then
do k=1,nextra
kbegin(k)=kbegin(k-ione)+ione+num_loc_groups
end do
end if
do k=nextra+ione,npe
kbegin(k)=kbegin(k-ione)+num_loc_groups
end do
do k=0,npe-ione
kend(k)=kbegin(k+ione)-ione
end do
if(mype == izero) then
write(6,*)' kbegin=',kbegin
write(6,*)' kend= ',kend
end if
num_j_groups=jm/npe
jextra=jm-num_j_groups*npe
jbegin(0)=ione
if(jextra > izero) then
do j=1,jextra
jbegin(j)=jbegin(j-ione)+ione+num_j_groups
end do
end if
do j=jextra+ione,npe
jbegin(j)=jbegin(j-ione)+num_j_groups
end do
do j=0,npe-ione
jend(j)=min(jbegin(j+ione)-ione,jm)
end do
if(mype == izero) then
write(6,*)' jbegin=',jbegin
write(6,*)' jend= ',jend
end if
allocate(ibuf((im+ione)*(jm+ione),kbegin(mype):kend(mype)))
call mpi_file_open(mpi_comm_world,trim(wrfges),mpi_mode_rdonly,mpi_info_null,mfcst,ierror)
! read geopotential
if(kord(i_fis)/=ione) then
allocate(jbuf(im,lm+ione,jbegin(mype):jend(mype)))
this_offset=offset(i_fis)+(jbegin(mype)-ione)*4*im*(lm+ione)
this_length=(jend(mype)-jbegin(mype)+ione)*im*(lm+ione)
call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length, &
mpi_integer,status,ierror)
call transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
jbegin,jend,kbegin,kend,mype,npe,im,jm,lm+ione,im+ione,jm+ione,i_fis,i_fis)
deallocate(jbuf)
end if
! read temps
if(kord(i_t)/=ione) then
allocate(jbuf(im,lm,jbegin(mype):jend(mype)))
this_offset=offset(i_t)+(jbegin(mype)-ione)*4*im*lm
this_length=(jend(mype)-jbegin(mype)+ione)*im*lm
call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length, &
mpi_integer,status,ierror)
call transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+ione,jm+ione,i_t,i_t+lm-ione)
deallocate(jbuf)
end if
! read q
if(kord(i_q)/=ione) then
allocate(jbuf(im,lm,jbegin(mype):jend(mype)))
this_offset=offset(i_q)+(jbegin(mype)-ione)*4*im*lm
this_length=(jend(mype)-jbegin(mype)+ione)*im*lm
call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length, &
mpi_integer,status,ierror)
call transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
jbegin,jend,kbegin,kend,mype,npe,im,jm,lm,im+ione,jm+ione,i_q,i_q+lm-ione)
deallocate(jbuf)
end if
! read u
if(kord(i_u)/=ione) then
allocate(jbuf(im+ione,lm,jbegin(mype):jend(mype)))
this_offset=offset(i_u)+(jbegin(mype)-ione)*4*(im+ione)*lm
this_length=(jend(mype)-jbegin(mype)+ione)*(im+ione)*lm
call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length, &
mpi_integer,status,ierror)
call transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
jbegin,jend,kbegin,kend,mype,npe,im+ione,jm,lm,im+ione,jm+ione,i_u,i_u+lm-ione)
deallocate(jbuf)
end if
! read v
if(kord(i_v)/=ione) then
jend2=jend
! Account for extra lat for v
jend2(npe-ione)=jend2(npe-ione)+ione
allocate(jbuf(im,lm,jbegin(mype):jend2(mype)))
this_offset=offset(i_v)+(jbegin(mype)-ione)*4*im*lm
this_length=(jend2(mype)-jbegin(mype)+ione)*im*lm
call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length, &
mpi_integer,status,ierror)
call transfer_jbuf2ibuf(jbuf,jbegin(mype),jend2(mype),ibuf,kbegin(mype),kend(mype), &
jbegin,jend2,kbegin,kend,mype,npe,im,jm+ione,lm,im+ione,jm+ione,i_v,i_v+lm-ione)
deallocate(jbuf)
end if
! read smois
if(kord(i_smois)/=ione) then
allocate(jbuf(im,ksize,jbegin(mype):jend(mype)))
this_offset=offset(i_smois)+(jbegin(mype)-ione)*4*im*ksize
this_length=(jend(mype)-jbegin(mype)+ione)*im*ksize
call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length, &
mpi_integer,status,ierror)
call transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
jbegin,jend,kbegin,kend,mype,npe,im,jm,ksize,im+ione,jm+ione,i_smois,i_smois)
deallocate(jbuf)
end if
! read tslb
if(kord(i_tslb)/=ione) then
allocate(jbuf(im,ksize,jbegin(mype):jend(mype)))
this_offset=offset(i_tslb)+(jbegin(mype)-ione)*4*im*ksize
this_length=(jend(mype)-jbegin(mype)+ione)*im*ksize
call mpi_file_read_at(mfcst,this_offset,jbuf(1,1,jbegin(mype)),this_length, &
mpi_integer,status,ierror)
call transfer_jbuf2ibuf(jbuf,jbegin(mype),jend(mype),ibuf,kbegin(mype),kend(mype), &
jbegin,jend,kbegin,kend,mype,npe,im,jm,ksize,im+ione,jm+ione,i_tslb,i_tslb)
deallocate(jbuf)
end if
!---------------------- read surface files last
do k=kbegin(mype),kend(mype)
if(kdim(k)==ione.or.kord(k)==ione) then
call mpi_file_read_at(mfcst,offset(k),ibuf(1,k),length(k),mpi_integer,status,ierror)
if(igtype(k)==ione) call expand_ibuf(ibuf(1,k),im ,jm ,im+ione,jm+ione)
if(igtype(k)==2_i_kind) call expand_ibuf(ibuf(1,k),im+ione,jm ,im+ione,jm+ione)
if(igtype(k)==3_i_kind) call expand_ibuf(ibuf(1,k),im ,jm+ione,im+ione,jm+ione)
end if
end do
call mpi_file_close(mfcst,ierror)
! next interpolate to analysis grid, then distribute to subdomains
allocate(temp1(im,jm),temp1u(im+ione,jm),temp1v(im,jm+ione))
allocate(tempa(itotsub,kbegin(mype):kend(mype)))
do ifld=kbegin(mype),kend(mype)
if(igtype(ifld) == ione) then
call move_ibuf_hg(ibuf(1,ifld),temp1,im+ione,jm+ione,im,jm)
call fill_mass_grid2t(temp1,im,jm,tempa(1,ifld),ione)
else if(igtype(ifld) == -ione) then
call move_ibuf_ihg(ibuf(1,ifld),temp1,im+ione,jm+ione,im,jm)
call fill_mass_grid2t(temp1,im,jm,tempa(1,ifld),ione)
else if(igtype(ifld) == 2_i_kind) then
call move_ibuf_hg(ibuf(1,ifld),temp1u,im+ione,jm+ione,im+ione,jm)
call fill_mass_grid2u(temp1u,im,jm,tempa(1,ifld),ione)
else if(igtype(ifld) == 3_i_kind) then
call move_ibuf_hg(ibuf(1,ifld),temp1v,im+ione,jm+ione,im,jm+ione)
call fill_mass_grid2v(temp1v,im,jm,tempa(1,ifld),ione)
end if
end do
deallocate(ibuf)
deallocate(temp1,temp1u,temp1v)
allocate(all_loc(lat2,lon2,num_mass_fields))
call generic_grid2sub(tempa,all_loc,kbegin(mype),kend(mype),kbegin,kend,mype,num_mass_fields)
! Next do conversion of units as necessary and
! reorganize into WeiYu's format--
kt=i_t-ione
kq=i_q-ione
ku=i_u-ione
kv=i_v-ione
! wrf pressure variable is dry air partial pressure--need to add water vapor contribution
! so accumulate 1 + total water vapor to use as correction factor
q_integral=one
do k=1,nsig
deltasigma=eta1_ll(k)-eta1_ll(k+ione)
kt=kt+ione
kq=kq+ione
ku=ku+ione
kv=kv+ione
do i=1,lon2
do j=1,lat2
ges_u(j,i,k,it) = all_loc(j,i,ku)
ges_v(j,i,k,it) = all_loc(j,i,kv)
ges_q(j,i,k,it) = all_loc(j,i,kq)
q_integral(j,i) = q_integral(j,i)+deltasigma*ges_q(j,i,k,it)
! Convert guess mixing ratio to specific humidity
ges_q(j,i,k,it) = ges_q(j,i,k,it)/(one+ges_q(j,i,k,it))
! Add offset to get guess potential temperature
ges_pot(j,i,k) = all_loc(j,i,kt) + h300
end do
end do
end do
do i=1,lon2
do j=1,lat2
! NOTE: MASS surface elevation is multiplied by g, so divide by g below
ges_z(j,i,it) = all_loc(j,i,i_fis)/grav
! Convert psfc units of mb and then convert to log(psfc) in cb
psfc_this_dry=r0_01*(all_loc(j,i,i_mub)+all_loc(j,i,i_mu)+pt_regional_single)
psfc_this=(psfc_this_dry-pt_ll)*q_integral(j,i)+pt_ll
ges_ps(j,i,it)=one_tenth*psfc_this ! convert from mb to cb
sno(j,i,it)=all_loc(j,i,i_sno)
soil_moi(j,i,it)=all_loc(j,i,i_smois)
soil_temp(j,i,it)=all_loc(j,i,i_tslb)
end do
end do
! Convert potenital temperature to temperature. Then convert
! sensible to virtual temperature
do k=1,nsig
do i=1,lon2
do j=1,lat2
work_prsl = one_tenth*(aeta1_ll(k)*(r10*ges_ps(j,i,it)-pt_ll)+pt_ll)
work_prslk = (work_prsl/r100)**rd_over_cp_mass
ges_tsen(j,i,k,it)= ges_pot(j,i,k)*work_prslk
ges_tv(j,i,k,it) = ges_tsen(j,i,k,it) * (one+fv*ges_q(j,i,k,it))
end do
end do
end do
! Transfer surface fields
do i=1,lon2
do j=1,lat2
fact10(j,i,it)=one ! later fix this by using correct w10/w(1)
veg_type(j,i,it)=all_loc(j,i,i_ivgtyp)
veg_frac(j,i,it)=r0_01*all_loc(j,i,i_vegfrac)
soil_type(j,i,it)=all_loc(j,i,i_isltyp)
sm_this=zero
if(all_loc(j,i,i_sm) /= zero_single) sm_this=one
xice_this=zero
if(all_loc(j,i,i_xice) /= zero_single) xice_this=one
isli_this=izero
if(xice_this==one) isli_this=2_i_kind
if(xice_this==zero.and.sm_this==one) isli_this=ione
isli(j,i,it)=isli_this
!?????????????????????????????????check to see if land skin temp is pot temp--if so, need to convert
sfct(j,i,it)=all_loc(j,i,i_sst)
if(isli(j,i,it) /= izero) sfct(j,i,it)=all_loc(j,i,i_tsk)
if(sfct(j,i,it) < one) then
! For now, replace missing skin temps with 1st sigma level temp
sfct(j,i,it)=ges_tsen(j,i,1,it)
num_doubtful_sfct=num_doubtful_sfct+ione
if(num_doubtful_sfct <= 100) &
write(6,*)'READ_WRF_MASS_BINARY_GUESS: ',&
'doubtful skint replaced with 1st sigma level t, j,i,mype,sfct=',&
j,i,mype,sfct(j,i,it)
end if
end do
end do
end do
call mpi_reduce(num_doubtful_sfct,num_doubtful_sfct_all,ione,mpi_integer,mpi_sum,&
izero,mpi_comm_world,ierror)
if(mype==izero) write(6,*)' in read_wrf_mass_guess, num_doubtful_sfct_all = ',num_doubtful_sfct_all
if(mype==izero) write(6,*)' in read_wrf_mass_guess, num_doubtful_sfct_all = ',num_doubtful_sfct_all
if(mype==10_i_kind) write(6,*)' in read_wrf_mass_guess, min,max(sfct)=', &
minval(sfct),maxval(sfct)
deallocate(all_loc,igtype,kdim,kord)
return
end subroutine read_wrf_mass_binary_guess
subroutine read_wrf_mass_netcdf_guess(mype)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_wrf_mass_guess read wrf_mass interface file
! prgmmr: parrish org: np22 date: 2003-09-05
!
! abstract: in place of read_guess for global application, read guess
! from regional model, in this case the wrf mass core model.
! This version reads a binary file created
! in a previous step that interfaces with the wrf infrastructure.
! A later version will read directly from the wrf restart file.
! The guess is read in by complete horizontal fields, one field
! per processor, in parallel. Each horizontal input field is
! converted from the staggered c-grid to an unstaggered a-grid.
! On the c-grid, u is shifted 1/2 point in the negative x direction
! and v 1/2 point in the negative y direction, but otherwise the
! three grids are regular. When the fields are read in, nothing
! is done to mass variables, but wind variables are interpolated to
! mass points.
!
! program history log:
! 2004--7-15 parrish
! 2004-08-02 treadon - add only to module use, add intent in/out
! 2004-09-10 parrish - correct error in land-sea mask interpretation
! 2005-02-23 wu - setup for qoption=2 and output stats of RH
! 2005-04-01 treadon - add initialization of ges_oz, prsi_oz; comestic format changes
! 2005-05-27 parrish - add call get_derivatives
! 2005-11-21 derber - make qoption=1 work same as qoption=2
! 2005-11-29 derber - remove external iteration dependent calculations
! 2005-12-15 treadon - remove initialization of certain guess arrays (done elsewhere)
! 2006-02-02 treadon - remove load_prsges,load_geop_hgt,prsl,prslk (not used)
! 2006-02-15 treadon - convert moisture mixing ratio to specific humidity
! 2006-03-07 treadon - convert guess potential temperature to vritual temperature
! 2006-04-06 middlecoff - changed nfcst from 11 to lendian_in
! 2006-07-28 derber - include sensible temperatures
! 2006-07-31 kleist - change to use ges_ps instead of lnps
! 2007-03-13 derber - remove unused qsinv2 from jfunc use list
! 2008-04-16 safford - rm unused uses
! 2010-03-29 hu - add code to read in cloud/hydrometeor fields
! and distributed them to all processors
!
! input argument list:
! mype - pe number
!
! NOTES: need to pay special attention to various surface fields to make sure
! they are correct (check units). fact10 needs to be computed from
! 10m wind, which is included in wrf mass interface file.
!
! Cloud water currently set to zero. Need to fix before doing precip
! assimilation--wait until global adopts Brad Ferrier's multi-species
! scheme??
!
! Ozone currently set to zero. No ozone variable in wrf mass core--climatology
! instead. Do we use climatology?
!
! No background bias yet. (biascor ignored)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,r_single,i_kind
use mpimod, only: mpi_sum,mpi_integer,mpi_real4,mpi_comm_world,npe,ierror
use guess_grids, only: ges_z,ges_ps,ges_tv,ges_q,ges_u,ges_v,&
fact10,soil_type,veg_frac,veg_type,sfct,sno,soil_temp,soil_moi,&
isli,nfldsig,ifilesig,ges_tsen
use guess_grids, only: ges_qc,ges_qi,ges_qr,ges_qs,ges_qg, &
ges_xlon,ges_xlat,soil_temp_cld,isli_cld,ges_tten
use gridmod, only: lat2,lon2,nlat_regional,nlon_regional,&
nsig,ijn_s,displs_s,eta1_ll,pt_ll,itotsub,aeta1_ll
use constants, only: izero,ione,zero,one,grav,fv,zero_single,rd_over_cp_mass,one_tenth
use gsi_io, only: lendian_in
use rapidrefresh_cldsurf_mod, only: l_cloud_analysis
implicit none
! Declare passed variables
integer(i_kind),intent(in):: mype
! Declare local parameters
real(r_kind),parameter:: r0_01=0.01_r_kind
real(r_kind),parameter:: r10=10.0_r_kind
real(r_kind),parameter:: r100=100.0_r_kind
! Declare local variables
integer(i_kind) kt,kq,ku,kv
! MASS variable names stuck in here
! other internal variables
real(r_single) tempa(itotsub)
real(r_single),allocatable::temp1(:,:),temp1u(:,:),temp1v(:,:)
real(r_single),allocatable::all_loc(:,:,:)
integer(i_kind),allocatable::itemp1(:,:)
integer(i_kind),allocatable::igtype(:),jsig_skip(:)
character(60),allocatable::identity(:)
character(6) filename
integer(i_kind) irc_s_reg(npe),ird_s_reg(npe)
integer(i_kind) ifld,im,jm,lm,num_mass_fields
integer(i_kind) num_all_fields,num_loc_groups,num_all_pad
integer(i_kind) i,icount,icount_prev,it,j,k
integer(i_kind) i_0,i_psfc,i_fis,i_t,i_q,i_u,i_v,i_sno,i_u10,i_v10,i_smois,i_tslb
integer(i_kind) i_sm,i_xice,i_sst,i_tsk,i_ivgtyp,i_isltyp,i_vegfrac
integer(i_kind) isli_this
real(r_kind) psfc_this,psfc_this_dry,sm_this,xice_this
real(r_kind),dimension(lat2,lon2):: q_integral
real(r_kind),dimension(lat2,lon2,nsig):: ges_pot
integer(i_kind) num_doubtful_sfct,num_doubtful_sfct_all
real(r_kind) deltasigma
real(r_kind):: work_prsl,work_prslk
integer(i_kind) i_qc,i_qi,i_qr,i_qs,i_qg,kqc,kqi,kqr,kqs,kqg,i_xlon,i_xlat,i_tt,ktt
! WRF MASS input grid dimensions in module gridmod
! These are the following:
! im -- number of x-points on C-grid
! jm -- number of y-points on C-grid
! lm -- number of vertical levels ( = nsig for now)
num_doubtful_sfct=izero
if(mype==izero) write(6,*)' at 0 in read_wrf_mass_guess'
! Big section of operations done only on first outer iteration
if(mype==izero) write(6,*)' at 0.1 in read_wrf_mass_guess'
im=nlon_regional
jm=nlat_regional
lm=nsig
! Following is for convenient WRF MASS input
num_mass_fields=14_i_kind+4_i_kind*lm
if(l_cloud_analysis) num_mass_fields=14_i_kind+4_i_kind*lm+6_i_kind*lm+2_i_kind
num_all_fields=num_mass_fields*nfldsig
num_loc_groups=num_all_fields/npe
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, lm =",i6)')lm
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, num_mass_fields=",i6)')num_mass_fields
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, nfldsig =",i6)')nfldsig
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, num_all_fields=",i6)')num_all_fields
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, npe =",i6)')npe
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, num_loc_groups=",i6)')num_loc_groups
do
num_all_pad=num_loc_groups*npe
if(num_all_pad >= num_all_fields) exit
num_loc_groups=num_loc_groups+ione
end do
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, num_all_pad =",i6)')num_all_pad
if(mype==izero) write(6,'(" at 1 in read_wrf_mass_guess, num_loc_groups=",i6)')num_loc_groups
allocate(all_loc(lat2,lon2,num_all_pad))
allocate(jsig_skip(num_mass_fields))
allocate(igtype(num_mass_fields))
allocate(identity(num_mass_fields))
! 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 u-grid
! igtype < 0 for integer field
i=izero
! for cloud analysis
if(l_cloud_analysis) then
i=i+ione ; i_xlat=i ! xlat
write(identity(i),'("record ",i3,"--xlat")')i
jsig_skip(i)=3_i_kind ! number of files to skip before getting to xlat
igtype(i)=ione
i=i+ione ; i_xlon=i ! xlon
write(identity(i),'("record ",i3,"--xlon")')i
jsig_skip(i)=izero !
igtype(i)=ione
endif
i=i+ione ; i_psfc=i ! psfc
write(identity(i),'("record ",i3,"--psfc")')i
jsig_skip(i)=5_i_kind ! number of files to skip before getting to psfc
if(l_cloud_analysis) jsig_skip(i)=0_i_kind ! number of files to skip before getting to psfc
igtype(i)=ione
i=i+ione ; i_fis=i ! sfc geopotential
write(identity(i),'("record ",i3,"--fis")')i
jsig_skip(i)=izero
igtype(i)=ione
i_t=i+ione
do k=1,lm
i=i+ione ! theta(k) (pot temp)
write(identity(i),'("record ",i3,"--t(",i2,")")')i,k
jsig_skip(i)=izero
igtype(i)=ione
end do
i_q=i+ione
do k=1,lm
i=i+ione ! q(k)
write(identity(i),'("record ",i3,"--q(",i2,")")')i,k
jsig_skip(i)=izero ; igtype(i)=ione
end do
i_u=i+ione
do k=1,lm
i=i+ione ! u(k)
write(identity(i),'("record ",i3,"--u(",i2,")")')i,k
jsig_skip(i)=izero ; igtype(i)=2_i_kind
end do
i_v=i+ione
do k=1,lm
i=i+ione ! v(k)
write(identity(i),'("record ",i3,"--v(",i2,")")')i,k
jsig_skip(i)=izero ; igtype(i)=3_i_kind
end do
i=i+ione ; i_sm=i ! landmask
write(identity(i),'("record ",i3,"--sm")')i
jsig_skip(i)=izero ; igtype(i)=ione
i=i+ione ; i_xice=i ! xice
write(identity(i),'("record ",i3,"--xice")')i
jsig_skip(i)=izero ; igtype(i)=ione
i=i+ione ; i_sst=i ! sst
write(identity(i),'("record ",i3,"--sst")')i
jsig_skip(i)=izero ; igtype(i)=ione
i=i+ione ; i_ivgtyp=i ! ivgtyp
write(identity(i),'("record ",i3,"--ivgtyp")')i
jsig_skip(i)=izero ; igtype(i)=-ione
i=i+ione ; i_isltyp=i ! isltyp
write(identity(i),'("record ",i3,"--isltyp")')i
jsig_skip(i)=izero ; igtype(i)=-ione
i=i+ione ; i_vegfrac=i ! vegfrac
write(identity(i),'("record ",i3,"--vegfrac")')i
jsig_skip(i)=izero ; igtype(i)=ione
i=i+ione ; i_sno=i ! sno
write(identity(i),'("record ",i3,"--sno")')i
jsig_skip(i)=izero ; igtype(i)=ione
i=i+ione ; i_u10=i ! u10
write(identity(i),'("record ",i3,"--u10")')i
jsig_skip(i)=izero ; igtype(i)=ione
i=i+ione ; i_v10=i ! v10
write(identity(i),'("record ",i3,"--v10")')i
jsig_skip(i)=izero ; igtype(i)=ione
i=i+ione ; i_smois=i ! smois
write(identity(i),'("record ",i3,"--smois(",i2,")")')i,k
jsig_skip(i)=izero ; igtype(i)=ione
i=i+ione ; i_tslb=i ! tslb
write(identity(i),'("record ",i3,"--tslb(",i2,")")')i,k
jsig_skip(i)=izero ; igtype(i)=ione
i=i+ione ; i_tsk=i ! tsk
write(identity(i),'("record ",i3,"--sst")')i
jsig_skip(i)=izero ; igtype(i)=ione
! for cloud array
if(l_cloud_analysis) then
i_qc=i+ione
do k=1,lm
i=i+ione ! qc(k)
write(identity(i),'("record ",i3,"--qc(",i2,")")')i,k
jsig_skip(i)=izero ; igtype(i)=ione
end do
i_qr=i+ione
do k=1,lm
i=i+ione ! qi(k)
write(identity(i),'("record ",i3,"--qr(",i2,")")')i,k
jsig_skip(i)=izero ; igtype(i)=ione
end do
i_qs=i+ione
do k=1,lm
i=i+ione ! qr(k)
write(identity(i),'("record ",i3,"--qs(",i2,")")')i,k
jsig_skip(i)=izero ; igtype(i)=ione
end do
i_qi=i+ione
do k=1,lm
i=i+ione ! qs(k)
write(identity(i),'("record ",i3,"--qi(",i2,")")')i,k
jsig_skip(i)=izero ; igtype(i)=ione
end do
i_qg=i+ione
do k=1,lm
i=i+ione ! qg(k)
write(identity(i),'("record ",i3,"--qg(",i2,")")')i,k
jsig_skip(i)=izero ; igtype(i)=ione
end do
i_tt=i+ione
do k=1,lm
i=i+ione ! tt(k)
write(identity(i),'("record ",i3,"--tt(",i2,")")')i,k
jsig_skip(i)=izero ; igtype(i)=ione
end do
endif
! End of stuff from MASS restart file
allocate(temp1(im,jm),itemp1(im,jm),temp1u(im+ione,jm),temp1v(im,jm+ione))
do i=1,npe
irc_s_reg(i)=ijn_s(mype+ione)
end do
ird_s_reg(1)=izero
do i=1,npe
if(i /= ione) ird_s_reg(i)=ird_s_reg(i-ione)+irc_s_reg(i-ione)
end do
! Read wrf MASS fixed format file created from external interface
! This is done by reading in parallel from every pe, and redistributing
! to local domains once for every npe fields read in, using
! mpi_all_to_allv
icount=izero
icount_prev=ione
do it=1,nfldsig
write(filename,'("sigf",i2.2)')ifilesig(it)
open(lendian_in,file=filename,form='unformatted') ; rewind lendian_in
write(6,*)'READ_WRF_MASS_GUESS: open lendian_in=',lendian_in,' to file=',filename
! Read, interpolate, and distribute MASS restart fields
do ifld=1,num_mass_fields
icount=icount+ione
if(jsig_skip(ifld) > izero) then
do i=1,jsig_skip(ifld)
read(lendian_in)
end do
end if
if(mype==mod(icount-ione,npe)) then
if(igtype(ifld)==ione) then
read(lendian_in)((temp1(i,j),i=1,im),j=1,jm)
write(6,'(" ifld, temp1(im/2,jm/2)=",i6,e15.5)')ifld,temp1(im/2,jm/2)
call fill_mass_grid2t(temp1,im,jm,tempa,ione)
end if
if(igtype(ifld)==2_i_kind) then
read(lendian_in)((temp1u(i,j),i=1,im+ione),j=1,jm)
write(6,'(" ifld, temp1u(im/2,jm/2)=",i6,e15.5)')ifld,temp1u(im/2,jm/2)
call fill_mass_grid2u(temp1u,im,jm,tempa,ione)
end if
if(igtype(ifld)==3_i_kind) then
read(lendian_in)((temp1v(i,j),i=1,im),j=1,jm+ione)
write(6,'(" ifld, temp1v(im/2,jm/2)=",i6,e15.5)')ifld,temp1v(im/2,jm/2)
call fill_mass_grid2v(temp1v,im,jm,tempa,ione)
end if
if(igtype(ifld) < izero) then
read(lendian_in)((itemp1(i,j),i=1,im),j=1,jm)
do j=1,jm
do i=1,im
temp1(i,j)=itemp1(i,j)
end do
end do
write(6,'(" ifld, temp1(im/2,jm/2)=",i6,e15.5)')ifld,temp1(im/2,jm/2)
call fill_mass_grid2t(temp1,im,jm,tempa,ione)
end if
else
read(lendian_in)
end if
! Distribute to local domains everytime we have npe fields
if(mod(icount,npe) == izero.or.icount==num_all_fields) then
call mpi_alltoallv(tempa,ijn_s,displs_s,mpi_real4, &
all_loc(1,1,icount_prev),irc_s_reg,ird_s_reg,mpi_real4,mpi_comm_world,ierror)
icount_prev=icount+ione
end if
end do
close(lendian_in)
end do
! do kv=i_v,i_v+nsig-ione
! if(mype==izero) write(6,*)' at 1.15, kv,mype,j,i,v=', &
! kv,mype,2,1,all_loc(2,1,kv)
! end do
! Next do conversion of units as necessary and
! reorganize into WeiYu's format--
do it=1,nfldsig
i_0=(it-ione)*num_mass_fields
kt=i_0+i_t-ione
kq=i_0+i_q-ione
ku=i_0+i_u-ione
kv=i_0+i_v-ione
! hydrometeors
if(l_cloud_analysis) then
kqc=i_0+i_qc-ione
kqr=i_0+i_qr-ione
kqs=i_0+i_qs-ione
kqi=i_0+i_qi-ione
kqg=i_0+i_qg-ione
ktt=i_0+i_tt-ione
endif
! wrf pressure variable is dry air partial pressure--need to add water vapor contribution
! so accumulate 1 + total water vapor to use as correction factor
q_integral=one
do k=1,nsig
deltasigma=eta1_ll(k)-eta1_ll(k+ione)
kt=kt+ione
kq=kq+ione
ku=ku+ione
kv=kv+ione
! hydrometeors
if(l_cloud_analysis) then
kqc=kqc+ione
kqr=kqr+ione
kqs=kqs+ione
kqi=kqi+ione
kqg=kqg+ione
ktt=ktt+ione
endif
do i=1,lon2
do j=1,lat2
ges_u(j,i,k,it) = all_loc(j,i,ku)
ges_v(j,i,k,it) = all_loc(j,i,kv)
ges_pot(j,i,k) = all_loc(j,i,kt)
ges_q(j,i,k,it) = all_loc(j,i,kq)
q_integral(j,i) = q_integral(j,i)+deltasigma*ges_q(j,i,k,it)
! Convert guess mixing ratio to specific humidity
ges_q(j,i,k,it) = ges_q(j,i,k,it)/(one+ges_q(j,i,k,it))
! hydrometeors
if(l_cloud_analysis) then
ges_qc(j,i,k,it) = all_loc(j,i,kqc)
ges_qi(j,i,k,it) = all_loc(j,i,kqi)
ges_qr(j,i,k,it) = all_loc(j,i,kqr)
ges_qs(j,i,k,it) = all_loc(j,i,kqs)
ges_qg(j,i,k,it) = all_loc(j,i,kqg)
ges_tten(j,i,k,it) = all_loc(j,i,ktt)
endif
end do
end do
end do
do i=1,lon2
do j=1,lat2
! NOTE: MASS surface elevation is multiplied by g, so divide by g below
ges_z(j,i,it) = all_loc(j,i,i_0+i_fis)/grav
! Convert psfc units of mb and then convert to log(psfc) in cb
psfc_this_dry=r0_01*all_loc(j,i,i_0+i_psfc)
psfc_this=(psfc_this_dry-pt_ll)*q_integral(j,i)+pt_ll
ges_ps(j,i,it)=one_tenth*psfc_this ! convert from mb to cb
sno(j,i,it)=all_loc(j,i,i_0+i_sno)
soil_moi(j,i,it)=all_loc(j,i,i_0+i_smois)
soil_temp(j,i,it)=all_loc(j,i,i_0+i_tslb)
! for cloud analysis
if(l_cloud_analysis) then
soil_temp_cld(j,i,it)=soil_temp(j,i,it)
ges_xlon(j,i,it)=all_loc(j,i,i_0+i_xlon)
ges_xlat(j,i,it)=all_loc(j,i,i_0+i_xlat)
endif
end do
end do
if(mype==10_i_kind) write(6,*)' in read_wrf_mass_guess, min,max(soil_moi)=', &
minval(soil_moi),maxval(soil_moi)
if(mype==10_i_kind) write(6,*)' in read_wrf_mass_guess, min,max(soil_temp)=', &
minval(soil_temp),maxval(soil_temp)
! Convert potenital temperature to temperature
do k=1,nsig
do i=1,lon2
do j=1,lat2
work_prsl = one_tenth*(aeta1_ll(k)*(r10*ges_ps(j,i,it)-pt_ll)+pt_ll)
work_prslk = (work_prsl/r100)**rd_over_cp_mass
ges_tsen(j,i,k,it) = ges_pot(j,i,k)*work_prslk
ges_tv(j,i,k,it) = ges_tsen(j,i,k,it) * (one+fv*ges_q(j,i,k,it))
end do
end do
end do
end do
! Transfer surface fields
do it=1,nfldsig
i_0=(it-ione)*num_mass_fields
do i=1,lon2
do j=1,lat2
fact10(j,i,it)=one ! later fix this by using correct w10/w(1)
veg_type(j,i,it)=all_loc(j,i,i_0+i_ivgtyp)
veg_frac(j,i,it)=r0_01*all_loc(j,i,i_0+i_vegfrac)
soil_type(j,i,it)=all_loc(j,i,i_0+i_isltyp)
sm_this=zero
if(all_loc(j,i,i_0+i_sm) /= zero_single) sm_this=one
xice_this=zero
if(all_loc(j,i,i_0+i_xice) /= zero_single) xice_this=one
isli_this=izero
if(xice_this==one) isli_this=2_i_kind
if(xice_this==zero.and.sm_this==one) isli_this=ione
isli(j,i,it)=isli_this
!?????????????????????????????????check to see if land skin temp is pot temp--if so, need to convert
sfct(j,i,it)=all_loc(j,i,i_0+i_sst)
if(isli(j,i,it) /= izero) sfct(j,i,it)=all_loc(j,i,i_0+i_tsk)
if(sfct(j,i,it) < one) then
! For now, replace missing skin temps with 1st sigma level temp
sfct(j,i,it)=all_loc(j,i,i_0+i_t)
write(6,*)' doubtful skint replaced with 1st sigma level t, j,i,mype,sfct=',&
j,i,mype,sfct(j,i,it)
num_doubtful_sfct=num_doubtful_sfct+ione
end if
if(l_cloud_analysis) then
isli_cld(j,i,it)=isli(j,i,it)
endif
end do
end do
end do
call mpi_reduce(num_doubtful_sfct,num_doubtful_sfct_all,ione,mpi_integer,mpi_sum,&
izero,mpi_comm_world,ierror)
if(mype==izero) write(6,*)' in read_wrf_mass_guess, num_doubtful_sfct_all = ',num_doubtful_sfct_all
if(mype==izero) write(6,*)' in read_wrf_mass_guess, num_doubtful_sfct_all = ',num_doubtful_sfct_all
if(mype==10_i_kind) write(6,*)' in read_wrf_mass_guess, min,max(sfct)=', &
minval(sfct),maxval(sfct)
if(mype==10_i_kind) write(6,*)' in read_wrf_mass_guess, min,max(veg_type)=', &
minval(veg_type),maxval(veg_type)
if(mype==10_i_kind) write(6,*)' in read_wrf_mass_guess, min,max(veg_frac)=', &
minval(veg_frac),maxval(veg_frac)
if(mype==10_i_kind) write(6,*)' in read_wrf_mass_guess, min,max(soil_type)=', &
minval(soil_type),maxval(soil_type)
if(mype==10_i_kind) write(6,*)' in read_wrf_mass_guess, min,max(isli)=', &
minval(isli),maxval(isli)
deallocate(all_loc,jsig_skip,igtype,identity)
deallocate(temp1,itemp1,temp1u,temp1v)
return
end subroutine read_wrf_mass_netcdf_guess
#else /* Start no WRF-library block */
subroutine read_wrf_mass_binary_guess()
!$$$ subprogram documentation block
! . . . .
! subprogram: read_wrf_mass_binary_guess
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-12-07 lueken - added subprogram doc block and implicit none
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
write(6,*)'READ_WRF_MASS_BINARY_GUESS: dummy routine, does nothing!'
end subroutine read_wrf_mass_binary_guess
subroutine read_wrf_mass_netcdf_guess()
!$$$ subprogram documentation block
! . . . .
! subprogram: read_wrf_mass_netcdf_guess
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-12-07 lueken - added subprogram doc block and implicit none
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
write(6,*)'READ_WRF_MASS_NETCDF_GUESS: dummy routine, does nothing!'
end subroutine read_wrf_mass_netcdf_guess
#endif /* End no WRF-library block */
subroutine generic_grid2sub(tempa,all_loc,kbegin_loc,kend_loc,kbegin,kend,mype,num_fields)
!$$$ subprogram documentation block
! . . . .
! subprogram: generic_grid2sub converts from full horizontal grid to subdomains
! prgmmr: parrish org: np22 date: 2004-11-29
!
! abstract: variation on subroutine grid2sub, 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 grid2sub for related use with mpi io.
!
! input argument list:
! tempa - input grid values in horizontal slab 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:
! all_loc - output grid values in vertical subdomain mode
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use mpimod, only: ierror,mpi_comm_world,mpi_real4,npe
use gridmod, only: ijn_s,itotsub,lat2,lon2
use kinds, only: r_single,i_kind
use constants, only: izero,ione
implicit none
integer(i_kind),intent(in ) :: kbegin_loc,kend_loc,mype,num_fields
integer(i_kind),intent(in ) :: kbegin(0:npe),kend(0:npe-ione)
real(r_single) ,intent(in ) :: tempa(itotsub,kbegin_loc:kend_loc)
real(r_single) ,intent( out) :: all_loc(lat2*lon2*num_fields)
integer(i_kind) k
integer(i_kind) sendcounts(0:npe-ione),sdispls(0:npe),recvcounts(0:npe-ione),rdispls(0:npe)
! first get alltoallv indices
sdispls(0)=izero
do k=0,npe-ione
sendcounts(k)=ijn_s(k+ione)*(kend_loc-kbegin_loc+ione)
sdispls(k+ione)=sdispls(k)+sendcounts(k)
end do
rdispls(0)=izero
do k=0,npe-ione
recvcounts(k)=ijn_s(mype+ione)*(kend(k)-kbegin(k)+ione)
rdispls(k+ione)=rdispls(k)+recvcounts(k)
end do
! then call reorder2
call reorder2_s(tempa,kend_loc-kbegin_loc+ione)
! then alltoallv and i think we are done??
call mpi_alltoallv(tempa,sendcounts,sdispls,mpi_real4, &
all_loc,recvcounts,rdispls,mpi_real4,mpi_comm_world,ierror)
end subroutine generic_grid2sub
subroutine reorder2_s(work,k_in)
!$$$ subprogram documentation block
! . . .
! subprogram: reorder2_s
!
! prgrmmr: kleist org: np20 date: 2004-01-25
!
! abstract: adapt reorder2 to single precision
!
! program history log:
! 2004-01-25 kleist
! 2004-05-14 kleist, documentation
! 2004-07-15 todling, protex-complaint prologue
! 2004-11-29 parrish, adapt reorder2 to single precision
! 2008-04-16 safford -- add subprogram doc block
!
! input argument list:
! k_in ! number of levs in work array
! work
!
! output argument list:
! work
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp; sgi origin 2000; compaq/hp
!
!$$$
! !USES:
use constants, only: izero,ione,zero_single
use mpimod, only: npe
use gridmod, only: ijn_s,itotsub
use kinds, only: r_single,i_kind
implicit none
! !INPUT PARAMETERS:
integer(i_kind) ,intent(in ) :: k_in ! number of levs in work array
! !INPUT/OUTPUT PARAMETERS:
real(r_single),dimension(itotsub,k_in),intent(inout) :: work
integer(i_kind) iloc,iskip,i,k,n
real(r_single),dimension(itotsub*k_in):: temp
! Zero out temp array
do k=1,itotsub*k_in
temp(k)=zero_single
end do
! Load temp array in order of subdomains
iloc=izero
iskip=izero
do n=1,npe
if (n/=ione) then
iskip=iskip+ijn_s(n-ione)
end if
do k=1,k_in
do i=1,ijn_s(n)
iloc=iloc+ione
temp(iloc)=work(iskip+i,k)
end do
end do
end do
! Now load the tmp array back into work
iloc=izero
do k=1,k_in
do i=1,itotsub
iloc=iloc+ione
work(i,k)=temp(iloc)
end do
end do
return
end subroutine reorder2_s
subroutine expand_ibuf(ibuf,im,jm,imp,jmp)
!$$$ subprogram documentation block
! . . . .
! subprogram: expand_ibuf expand array in place
! prgmmr: parrish org: np22 date: 2004-11-29
!
! abstract: expand array in place from im,jm to imp,jmp
!
! program history log:
! 2004-11-29 parrish
! 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 im,jm
! im - first grid index
! jm - second grid index
!
! output argument list:
! ibuf - output grid values in im+1,jm+1
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
! field of dim im*jm read into array of dim imp*jmp--need to readjust
use kinds, only: i_long,i_kind
use constants, only: izero,ione
implicit none
integer(i_kind),intent(in ) :: im,jm,imp,jmp
integer(i_long),intent(inout) :: ibuf(imp*jmp)
integer(i_long) itemp(imp,jmp)
integer(i_kind) i,ii,j
do j=1,jmp
do i=1,imp
itemp(i,j)=0_i_long
end do
end do
ii=izero
do j=1,jm
do i=1,im
ii=ii+ione
itemp(i,j)=ibuf(ii)
end do
end do
ii=izero
do j=1,jmp
do i=1,imp
ii=ii+ione
ibuf(ii)=itemp(i,j)
end do
end do
end subroutine expand_ibuf
subroutine transfer_jbuf2ibuf(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 ikj to ijk
! prgmmr: parrish org: np22 date: 2004-11-29
!
! abstract: redistribute 3-d field from ikj order to ijk order across processors
!
! program history log:
! 2004-11-29 parrish
! 2005-02-16 todling - replaced "use mpi" by use of mpimod
!
! input argument list:
! jbuf - input grid values distributed as all ik and a range of j on each processor
! 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:
! ibuf - output grid redistributed to ijk order (full i, full j, k_start <= k <= k_end)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
! flip around from ikj to ijk, moving result from jbuf to ibuf
use mpimod, only: mpi_comm_world,mpi_integer
use kinds, only: i_long,i_kind
use constants, only: izero,ione
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(in ) :: jbuf(im_jbuf,lm_jbuf,jbegin_loc:jend_loc)
integer(i_long),intent( out) :: ibuf(im_ibuf,jm_ibuf,kbegin_loc:kend_loc)
integer(i_kind),intent(in ) :: jbegin(0:npe),jend(0:npe-ione)
integer(i_kind),intent(in ) :: kbegin(0:npe),kend(0:npe-ione)
integer(i_long) sendbuf(im_jbuf*lm_jbuf*(jend_loc-jbegin_loc+2_i_kind))
integer(i_long) recvbuf(im_jbuf*jm_jbuf*(kend_loc-kbegin_loc+ione))
integer(i_long) recvcounts(0:npe-ione),displs(0:npe)
integer(i_kind) i,ipe,j,ierror,k,n,ii,k_t_start,k_t_end,sendcount
do ipe=0,npe-ione
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-ione
recvcounts(i)=im_jbuf*(k_t_end-k_t_start+ione)*(jend(i)-jbegin(i)+ione)
displs(i+ione)=displs(i)+recvcounts(i)
end do
! gather everything to ipe
ii=izero
do k=k_t_start,k_t_end
do j=jbegin_loc,jend_loc
do i=1,im_jbuf
ii=ii+ione
sendbuf(ii)=jbuf(i,k-k_start+ione,j)
end do
end do
end do
sendcount=ii
call mpi_gatherv(sendbuf,sendcount,mpi_integer,recvbuf,recvcounts, &
displs,mpi_integer,ipe,mpi_comm_world,ierror)
if(ipe==mype) then
ii=izero
do n=0,npe-ione
do k=k_t_start,k_t_end
do j=jbegin(n),jend(n)
do i=1,im_jbuf
ii=ii+ione
ibuf(i,j,k)=recvbuf(ii)
end do
end do
end do
end do
end if
end do
end subroutine transfer_jbuf2ibuf
subroutine move_ibuf_hg(ibuf,temp1,im_buf,jm_buf,im_out,jm_out)
!$$$ subprogram documentation block
! . . . .
! subprogram: move_ibuf_hg 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:
! ibuf - 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:
! temp1 - output grid values
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
! cp buf to temp1
use kinds, only: r_single,i_kind,i_long
use constants, only: zero_single
implicit none
integer(i_kind),intent(in ) :: im_buf,jm_buf,im_out,jm_out
integer(i_long),intent(in ) :: ibuf(im_buf,jm_buf)
real(r_single) ,intent( out) :: temp1(im_out,jm_out)
integer(i_kind) i,j
do j=1,jm_out
do i=1,im_out
temp1(i,j)=transfer(ibuf(i,j),zero_single)
end do
end do
end subroutine move_ibuf_hg
subroutine move_ibuf_ihg(ibuf,temp1,im_buf,jm_buf,im_out,jm_out)
!$$$ subprogram documentation block
! . . . .
! subprogram: move_ibuf_hg copy from one array to another
! prgmmr: parrish org: np22 date: 2004-11-29
!
! abstract: copy from one array to another, converting from int to real
!
! program history log:
! 2004-11-29 parrish
!
! input argument list:
! ibuf - 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:
! temp1 - output grid values
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
! cp buf to temp1
use kinds, only: i_long,r_single,i_kind
implicit none
integer(i_kind),intent(in ) :: im_buf,jm_buf,im_out,jm_out
integer(i_long),intent(in ) :: ibuf(im_buf,jm_buf)
real(r_single) ,intent( out) :: temp1(im_out,jm_out)
integer(i_kind) i,j
do j=1,jm_out
do i=1,im_out
temp1(i,j)=ibuf(i,j)
end do
end do
end subroutine move_ibuf_ihg