#include <machine.h>
subroutine read_guess(jiter,lat1,lon1,nsig,sigi,sigl,& 1,8
del2,del21,a3,factsml,factvml,coriolis,slatd,&
mype,nlath,nlon,jcap,nc,nc1,pln,qln,rln,&
#ifdef ibm_sp
trigs,trigs_2,trigs_3,&
#else
trigs,&
#endif
wgts,ltosi,ltosj,iglobal,ltosi_s,ltosj_s,itotsub,&
qsat,esat,desdtv,eges,dedlps,dedq,tropprs,&
deltat,hrdif,fact10,sfct,isli,sno,&
veg_type,veg_frac,soil_type,soil_temp,soil_moi,&
ggrid_g31,ggrid_g32,rdivs,&
rcwm,rtlat,rtlon,rqlat,rqlon,&
dstlast,dstb,ampdivt,npe,&
my1,nsig1,k5,npes,nnjc,nnc1,nnc1l,jnpe,nlast,ncommnc,ns2,ncd2,nclevel,&
nc1m,nc2,nc3,js1,js2,&
spec_6,spec1_3,spec1_6,spec1_9,&
tracers_6,vtid_6,pdryini_6,xncld_6,&
phat_3,phat_6,phat_9,k4)
!$$$ subprogram documentation block
! . . . .
! subprogram: inguessv same as inguess, but add vort,div, del(ps)
! prgmmr: parrish org: w/nmc22 date: 94-02-11
!
! abstract: augment inguess with vort, div, grad (log(psfc))
!
! program history log:
! 94-02-11 parrish
! 98-04-03 weiyu yang
! 99-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
!
! input argument list:
! jiter - outloop iteration number, at here jiter=1
! lat1 - number of gaussian lats in the sub-domain array
! lon1 - number of gaussian longitudes in the sub-domain array
! nsig - number of sigma levels
! nsig1 - number of sigma levels in each pe
! sigi - sigma values at interfaces of sigma layers
! sigl - sigma values at mid-point of each sigma layer
! del2 - n*(n+1)/a**2
! a3 - hydrostatic matrix
! factsml - factor to ensure proper scalar coefficients are zero
! factvml - factor to ensure proper vector coefficients are zero
! coriolis - coriolis coefficients
! mype - pe number
! nlath - number of gaussian lats in one hemisphere
! nlon - number of longitudes
! jcap - triangular truncation
! pln,qln,rln - used by fft
! trigs - used by fft
! wgts - gaussian integration weights
! ltosi,ltosj - array element indices for gathering local subdomain grid field array
! to the global grid field array, not including the boundary grid points
! of the subdomain
! iglobal - number of the global grid points, 2*nalth*nlon
! ltosi_s,ltosj_s - array element indices for scattering global grid field array
! to local subdomain array, including the boundary grid points of
! the subdomain
! itotsub - total number of grid points of all subdomains, including their boundary area
! ampdivt - amplitude for divergence stats
! npe - total number of pes available
!
! output argument list:
! qsat - 1/(saturation) specific humidity
! deltat - time interval of first guess fields
! hrdif - time interval between first guess fields and
! its former and later guess fields
! fact10 - 10 meter wind factor
! sfct - guess skin temperature
! isli - snow-land-ice mask
! sno - snow-ice mask
! rz_ges - guess mountains on grid
! rp_ges - guess log(sfcp) on grid
! ru_ges - guess u on grid
! rv_ges - guess v on grid
! rt_ges - guess t on grid
! rq_ges - guess specific humidity on grid
! roz_ges - guess ozone on grid in each pe subdomain
! rvorts,rplons,rplats - guess vorticity and grad(log(psfc))
! dstlast,dstb - coefs of last and guess div-tend
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use type_kinds, only: fp_kind,single
implicit none
include 'mpi_inc.h'
include 'ijstart.h'
character(12) sigg03,sigg09
logical ice
integer nc,jcap,iglobal,itotsub,lon1,lat1,jiter,nsig,nlon,nlath,mype,ns2
integer j,i,ncommnc,nlast,ncd2,izero,nclevel,npes,nsig1,npe,nnjc,jnpe
integer nnc1l,nnc1,nming_ges,nmings,ingsx,nmings2,k1,k,iwan,nrec,nming2
integer ii,nmype,lon2,lat2,it1,ierr,it,inges,isfc,ns1,nc1,n1,iix,ncy
integer mm1,nsfc,ncyc,nhdif
integer,dimension(2*nlath*nlon):: ltosi,ltosj
integer,dimension(3*nlath*nlon):: ltosi_s,ltosj_s
integer,dimension(nsig1):: k5,my1,k4
integer,dimension(nclevel):: nc2,nc3,js1,js2,nc1m
integer idateg(4),idate5(5)
real(fp_kind) ampdivt,deltat,zero
real(fp_kind),dimension(ncd2,nlath):: pln,qln,rln
#ifdef ibm_sp
real(fp_kind),dimension(50000+nlon*4):: trigs,trigs_2,trigs_3
#else
real(fp_kind),dimension(nlon+15):: trigs
#endif
real(fp_kind) wgts(2*nlath)
real(fp_kind),dimension(lat1+2,lon1+2,nsig,3):: rdivs
real(fp_kind),dimension(lat1+2,lon1+2,nsig,3):: rcwm,rtlat,rtlon,rqlat,rqlon
real(fp_kind),dimension(lat1+2,lon1+2,nsig):: rdivs1
real(fp_kind),dimension(lat1+2):: coriolis,slatd
real(single) hourg4
real(fp_kind) hourg,hrdif(2)
real(fp_kind),dimension(nc):: factsml,factvml,del21
real(fp_kind),dimension(nnc1l):: del2
real(fp_kind),dimension(nnc1l,nsig):: dstlast,dstb
real(fp_kind),dimension(nsig):: sigl
real(fp_kind),dimension(nsig+1):: sigi
real(fp_kind),dimension(nsig,nsig):: a3
real(fp_kind),dimension(nc,nsig1):: spec_6
real(fp_kind),dimension(3,nnc1l,npes):: spec1_3,spec1_6,spec1_9,phat_3,&
phat_6,phat_9
real(single) tracers_3,vtid_3,pdryini_3,xncld_3
real(single) tracers_6,vtid_6,pdryini_6,xncld_6
real(single) tracers_9,vtid_9,pdryini_9,xncld_9
real(fp_kind) buffer_time(3)
#ifdef ibm_sp
integer ireq_time
integer,dimension(MPI_STATUS_SIZE)::istat_time
integer,dimension(9):: ireq_sfc
integer,dimension(MPI_STATUS_SIZE,9)::istat_sfc
#endif
real(fp_kind),dimension(lat1+2,lon1+2,nsig):: qsat,esat,desdtv,eges,dedlps,dedq
integer,dimension(lat1+2,lon1+2):: isli
real(fp_kind),dimension(lat1+2,lon1+2):: fact10,sfct,sno,veg_type,veg_frac,tropprs
real(fp_kind),dimension(lat1+2,lon1+2):: soil_type,soil_temp,soil_moi,sli
real(fp_kind),dimension(3,lat1+2,lon1+2,npes):: ggrid_g31
real(fp_kind),dimension(3,lat1+2,lon1+2,npes,2):: ggrid_g32
real(fp_kind),allocatable,dimension(:,:):: sfc_g
! Start read_guess here.
!
!ibm* prefetch_for_store(zero,izero,n1,ns1,inges)
zero=0.0; izero=0
n1=nsig+1; ns1=2*nsig+1
!
! Set 3 and 9 hour sigma file names based on outer iteration number
inges=11
if (jiter==1) then
!
! Read surface fields needed for near surface winds and radiance data.
!ibm* prefetch_for_store(isfc,nsfc,ncyc,ireq_sfc)
isfc=15
nsfc=9
ncyc = (nsfc-1)/npe+1
#ifdef ibm_sp
ireq_sfc=mpi_request_null
#endif
open(isfc,file='sfcanl',form='unformatted')
allocate (sfc_g(itotsub,ncyc) )
!ibm* prefetch_for_store(mm1,iix)
mm1=mype+1
iix=izero
do ncy=1,ncyc
do ii=1,npe
iix=iix+1
if(iix<nsfc+1)then
!ibm* prefetch_for_store(nmype,nrec)
nmype=npe-iix+(ncy-1)*npe
if(iix == 1)nrec=3 ! surface skin temperature
if(iix == 2)nrec=5 ! snow depth
if(iix == 3)nrec=13 ! sea land ice factor
if(iix == 4)nrec=16 ! 10m wind factor
if(iix == 5)nrec=17 ! veg type
if(iix == 6)nrec=14 ! veg fraction
if(iix == 7)nrec=18 ! soil type
if(iix == 8)nrec=6 ! soil temperature
if(iix == 9)nrec=4 ! soil moisture
if (mype==nmype)call rdbges(sfc_g(1,ncy),nlath,nlon,isfc,nrec, &
ltosi_s,ltosj_s,itotsub)
#ifdef ibm_sp
if(iix == 1)then
call mpe_iscatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,sfct,ijn_s(mm1),&
mpi_rtype,nmype,mpi_comm_world,ireq_sfc(iix),ierror)
else if(iix == 2)then
call mpe_iscatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,sno,ijn_s(mm1),&
mpi_rtype,nmype,mpi_comm_world,ireq_sfc(iix),ierror)
else if(iix == 3)then
call mpe_iscatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,sli,ijn_s(mm1),&
mpi_rtype,nmype,mpi_comm_world,ireq_sfc(iix),ierror)
else if(iix == 4)then
call mpe_iscatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,fact10,ijn_s(mm1),&
mpi_rtype,nmype,mpi_comm_world,ireq_sfc(iix),ierror)
else if(iix == 5)then
call mpe_iscatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,veg_type,ijn_s(mm1),&
mpi_rtype,nmype,mpi_comm_world,ireq_sfc(iix),ierror)
else if(iix == 6)then
call mpe_iscatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,veg_frac,ijn_s(mm1),&
mpi_rtype,nmype,mpi_comm_world,ireq_sfc(iix),ierror)
else if(iix == 7)then
call mpe_iscatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,soil_type,ijn_s(mm1),&
mpi_rtype,nmype,mpi_comm_world,ireq_sfc(iix),ierror)
else if(iix == 8)then
call mpe_iscatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,soil_temp,ijn_s(mm1),&
mpi_rtype,nmype,mpi_comm_world,ireq_sfc(iix),ierror)
else if(iix == 9)then
call mpe_iscatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,soil_moi,ijn_s(mm1),&
mpi_rtype,nmype,mpi_comm_world,ireq_sfc(iix),ierror)
end if
#else
if(iix==1)then
call mpi_scatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,sfct,&
ijn_s(mm1),mpi_rtype,nmype,mpi_comm_world,ierror)
else if(iix==2)then
call mpi_scatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,sno,&
ijn_s(mm1),mpi_rtype,nmype,mpi_comm_world,ierror)
else if(iix==3)then
call mpi_scatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,sli,&
ijn_s(mm1),mpi_rtype,nmype,mpi_comm_world,ierror)
else if(iix==4)then
call mpi_scatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,fact10,&
ijn_s(mm1),mpi_rtype,nmype,mpi_comm_world,ierror)
else if(iix==5)then
call mpi_scatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,veg_type,&
ijn_s(mm1),mpi_rtype,nmype,mpi_comm_world,ierror)
else if(iix==6)then
call mpi_scatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,veg_frac,&
ijn_s(mm1),mpi_rtype,nmype,mpi_comm_world,ierror)
else if(iix==7)then
call mpi_scatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,soil_type,&
ijn_s(mm1),mpi_rtype,nmype,mpi_comm_world,ierror)
else if(iix==8)then
call mpi_scatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,soil_temp,&
ijn_s(mm1),mpi_rtype,nmype,mpi_comm_world,ierror)
else if(iix==9)then
call mpi_scatterv(sfc_g(1,ncy),ijn_s,displs_s,mpi_rtype,soil_moi,&
ijn_s(mm1),mpi_rtype,nmype,mpi_comm_world,ierror)
end if
#endif
end if
end do
end do
close(isfc)
!ibm* prefetch_for_store(sigg03,sigg09)
open(inges,file='sigf06',form='unformatted')
open(12,file='sigf03',form='unformatted')
open(13,file='sigf09',form='unformatted')
endif
!! else
!ibm* prefetch_for_store(sigg03,sigg09)
!! open(inges,file='siganl',form='unformatted')
!! open(12,file='sanl03',form='unformatted')
!! open(13,file='sanl09',form='unformatted')
!! endif
! ggrid_g3--- u, z, q, plat, --- v, empty, oz, plon, ---
! vort, ps, hs, sfcts.
call inguessv(ggrid_g31,inges,lat1,lon1,nsig,hourg,idateg,rdivs,&
rcwm,&
rtlat,rtlon,rqlat,rqlon,&
del21,factsml,factvml,mype,nlath,nlon,jcap,nc,nc1,&
pln,qln,rln,trigs,ltosi_s,ltosj_s,itotsub,npe,&
#ifdef ibm_sp
trigs_2,trigs_3,&
#endif
my1,nsig1,npes,ns2,ncd2,nclevel,nnc1,nnc1l,nnjc,&
nc1m,nc2,nc3,js1,js2,&
spec_6,spec1_6,&
tracers_6,vtid_6,pdryini_6,xncld_6,jiter,phat_6,jnpe,nlast,ncommnc,k4)
close(inges)
! Compute 1/(saturation) specific humidity from 6 hour guess fields.
! This is used in the moisture (limq) constraint in pcgsoi.
!ibm* prefetch_for_store(k1)
k1=n1
do k=1,nsig
k1=k1+1
do j=1,lon1+2
do i=1,lat1+2
!ibm* prefetch_for_store(qsat(i,j,k),rdivs1(i,j,k))
qsat(i,j,k)=ggrid_g31(1,i,j,k1) ! q.
rdivs1(i,j,k)=rdivs(i,j,k,1)
end do
end do
end do
!ibm* prefetch_for_store(ice)
ice=.true.
call genqsat(ggrid_g31,qsat,esat,lat1,lon1,nsig,sigl,desdtv,ice,npes)
call vapres(ggrid_g31,sigl,lat1,lon1,nsig,eges,dedlps,dedq,npes)
!
! Compute full divergence tendency from 6 hour guess fields.
! This is used for the divergence tendency constraint in pcgsoi
!
! ***NOTE: rdivs contains divergence on input; vertical velocity on output
!
if(ampdivt>zero) then
call fulldivt(ggrid_g31,rdivs1,nsig,lat1,lon1,del2,&
a3,sigl,sigi,dstlast,coriolis,mype,nlath,nlon,jcap,nc,nc1,&
pln,qln,rln,trigs,wgts,ltosi,ltosj,iglobal,npe,nsig1,npes,my1,&
nnc1,nnc1l,jnpe,nlast,ncommnc,ns2,k5,ncd2)
if(jiter==1.and.nnjc/=izero) then
do k=1,nsig
do i=1,nnc1l
!ibm* prefetch_for_store(dstb(i,k))
dstb(i,k)=dstlast(i,k)
end do
end do
end if
end if
! Compute tropopause level (in pressure, hPa)
call tpause(lat1,lon1,nsig,npes,ggrid_g31,sigl,slatd,tropprs)
!
! Check for consistency of times from guess files. Do all this on pe 0
if(jiter == 1)then
if(mype==0) then
!ibm* prefetch_for_store(iwan,deltat,hrdif)
iwan=izero; deltat=-9999.0; hrdif=zero
!ibm* prefetch_for_store(idate5(1),idate5(2))
idate5(1)=idateg(4); idate5(2)=idateg(2)
!ibm* prefetch_for_store(idate5(3),idate5(4),idate5(5))
idate5(3)=idateg(3); idate5(4)=idateg(1); idate5(5)=0
call w3fs21(idate5,nmings)
!ibm* prefetch_for_store(nming_ges)
nming_ges=nmings+60*hourg
write(6,*) 'READ_GUESS: guess file hourg,idateg,nming_ges=',hourg,idateg,&
nming_ges
if(abs(hourg)>6.001) iwan=iwan+1
do i=1,2
!ibm* prefetch_for_store(ingsx)
ingsx=11+i
rewind ingsx
read(ingsx,end=1112)
read(ingsx,end=1112) hourg4,idateg
rewind ingsx
!ibm* prefetch_for_store(hourg)
hourg = hourg4
!ibm* prefetch_for_store(idate5(1),idate5(2))
idate5(1)=idateg(4); idate5(2)=idateg(2)
!ibm* prefetch_for_store(idate5(3),idate5(4),idate5(5))
idate5(3)=idateg(3); idate5(4)=idateg(1); idate5(5)=izero
call w3fs21(idate5,nmings2)
!ibm* prefetch_for_store(nming2)
nming2=nmings2+60*hourg
write(6,*)'READ_GUESS: supplementary guess file ',hourg,idateg
write(6,*)'READ_GUESS: for suppl. guess file, nming2=',nming2
if(nmings2/=nmings.and.(nmings-nmings2)/=360) iwan=iwan+1
!ibm* prefetch_for_store(nhdif, hrdif(i))
nhdif=(nming2-nming_ges)/60.0
hrdif(i)=nhdif
if(abs(nhdif)/=3) iwan=iwan+1
go to 1113
1112 continue
iwan=iwan+1
1113 continue
end do
if(iwan==0) then
!ibm* prefetch_for_store(deltat)
deltat=180.0
end if
write(6,*)'READ_GUESS: deltat = ',deltat
!ibm* prefetch_for_store(buffer_time(1),buffer_time(2),buffer_time(3))
buffer_time(1) = deltat
buffer_time(2) = hrdif(1)
buffer_time(3) = hrdif(2)
end if
!ibm* prefetch_for_store(ireq_time)
#ifdef ibm_sp
ireq_time=mpi_request_null
call mpe_ibcast(buffer_time,3,mpi_rtype,0,mpi_comm_world,ireq_time,ierr)
#else
call mpi_bcast(buffer_time,3,mpi_rtype,0,mpi_comm_world,ierr)
#endif
!
! Wait to get times and surface fields. Load snow-land-ice mask
! into integer array, isli. Then, deallocate work arrays.
#ifdef ibm_sp
call mpi_wait(ireq_time,istat_time,ierror)
#endif
!ibm* prefetch_for_store(deltat,hrdif(1),hrdif(2))
deltat = buffer_time(1)
hrdif(1) = buffer_time(2)
hrdif(2) = buffer_time(3)
#ifdef ibm_sp
call mpi_waitall(nsfc,ireq_sfc,istat_sfc,ierror)
#endif
do j=1,lon1+2
do i=1,lat1+2
!ibm* prefetch_for_store(isli(i,j))
isli(i,j)=nint(sli(i,j)+0.0000001)
end do
end do
deallocate(sfc_g)
end if
!
! If 3 and/or 9 hour guess fields are available, read them in.
if(deltat>zero) then
do it=1,2
!ibm* prefetch_for_store(it1,ingsx)
ingsx=11+it
if(it==1) then
call inguessv(ggrid_g32(1,1,1,1,it),ingsx,lat1,lon1,nsig,&
hourg,idateg,rdivs(1,1,1,it+1),&
rcwm(1,1,1,it+1),&
rtlat(1,1,1,it+1),rtlon(1,1,1,it+1),&
rqlat(1,1,1,it+1),rqlon(1,1,1,it+1),&
del21,factsml,factvml,mype,nlath,nlon,jcap,nc,nc1,&
pln,qln,rln,trigs,ltosi_s,ltosj_s,itotsub,npe,&
#ifdef ibm_sp
trigs_2,trigs_3,&
#endif
my1,nsig1,npes,ns2,ncd2,nclevel,nnc1,nnc1l,nnjc,&
nc1m,nc2,nc3,js1,js2,&
spec_6,spec1_3,&
tracers_3,vtid_3,pdryini_3,xncld_3,jiter,phat_3,jnpe,nlast,ncommnc,k4)
else
call inguessv(ggrid_g32(1,1,1,1,it),ingsx,lat1,lon1,nsig,&
hourg,idateg,rdivs(1,1,1,it+1),&
rcwm(1,1,1,it+1),&
rtlat(1,1,1,it+1),rtlon(1,1,1,it+1),&
rqlat(1,1,1,it+1),rqlon(1,1,1,it+1),&
del21,factsml,factvml,mype,nlath,nlon,jcap,nc,nc1,&
pln,qln,rln,trigs,ltosi_s,ltosj_s,itotsub,npe,&
#ifdef ibm_sp
trigs_2,trigs_3,&
#endif
my1,nsig1,npes,ns2,ncd2,nclevel,nnc1,nnc1l,nnjc,&
nc1m,nc2,nc3,js1,js2,&
spec_6,spec1_9,&
tracers_9,vtid_9,pdryini_9,xncld_9,jiter,phat_9,jnpe,nlast,ncommnc,k4)
end if
close(ingsx)
end do
! No 3 and/or 9 hour guess fields. Set equal to 6 hour forecast
else
do it=1,2
do k=1,npes
do j=1,lon1+2
do i=1,lat1+2
!ibm* prefetch_for_store(ggrid_g32(1,i,j,k,it))
ggrid_g32(1,i,j,k,it)=ggrid_g31(1,i,j,k)
!ibm* prefetch_for_store(ggrid_g32(2,i,j,k,it))
ggrid_g32(2,i,j,k,it)=ggrid_g31(2,i,j,k)
!ibm* prefetch_for_store(ggrid_g32(3,i,j,k,it))
ggrid_g32(3,i,j,k,it)=ggrid_g31(3,i,j,k)
end do
end do
end do
end do
lon2=lon1+2
lat2=lat1+2
do it=2,3
do k=1,nsig
do j=1,lon2
do i=1,lat2
rdivs(i,j,k,it) = rdivs(i,j,k,1)
rcwm(i,j,k,it) = rcwm(i,j,k,1)
rtlat(i,j,k,it) = rtlat(i,j,k,1)
rtlon(i,j,k,it) = rtlon(i,j,k,1)
rqlat(i,j,k,it) = rqlat(i,j,k,1)
rqlon(i,j,k,it) = rqlon(i,j,k,1)
end do
end do
end do
end do
endif
!
! End of routine
return
end