module cressanl_common
!$$$ module documentation block
! . . . .
! module: cressanl_common
! prgmmr: pondeca org: np23 date: 2012-10-15
!
! abstract:
!
!
! program history log:
! 2012-10-15 pondeca - run a cressman analysis
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use kinds, only: i_kind,r_single
implicit none
integer(i_kind) nobsmax !the number of obs of all ob types (including ob types
!such as ugrel-ob ,w2-ob, wd-ob, etc.) on each memory
!tile
integer(i_kind) nmax !the largest number of obs of a given ob type on any one memory tile
!this number is therefore constant across processors
integer(i_kind) nflds,kflds
integer(i_kind) kps,kts,kqs,kus,kvs,kugrds,kvgrds,kws,kw2s,kwds, &
ktds,kgusts,kvis,kpblhs,kdists, &
kwspd10ms,ktd2ms,kmxtms,kmitms,kpmsls,khowvs,ktcamts,klcbas, &
kuwnd10ms,kvwnd10ms,kaltws,kaltw2s,kaltwds !alt stands for "alternate"
integer(i_kind),allocatable,dimension(:,:) :: ipointer
integer(i_kind),allocatable,dimension(:,:) :: jpointer
integer(i_kind),allocatable,dimension(:,:) :: uvpointer
real(r_single),allocatable,dimension(:) :: xlocs
real(r_single),allocatable,dimension(:) :: ylocs
real(r_single),allocatable,dimension(:) :: hgt0s
real(r_single),allocatable,dimension(:) :: hgts
real(r_single),allocatable,dimension(:) :: hobs
real(r_single),allocatable,dimension(:) :: rmuses
real(r_single),allocatable,dimension(:) :: oberrs
real(r_single),allocatable,dimension(:) :: dtimes
real(r_single),allocatable,dimension(:) :: rtvts
real(r_single),allocatable,dimension(:,:,:) :: bckgs
real(r_single),allocatable,dimension(:,:,:) :: xberrs !sqrt of bckg error variance
real(r_single),allocatable,dimension(:,:) :: terrain
character(8),allocatable,dimension(:) :: cstations
character(12),allocatable,dimension(:) :: obstypes
logical(4),allocatable,dimension(:) :: insubdoms
logical(4),allocatable,dimension(:) :: dups
contains
!*****************************************************************************************
!*****************************************************************************************
subroutine alloc_cressanl_common(ista,iend,jsta,jend,nx,ny,mype)
implicit none
integer(i_kind),intent(in)::ista,iend,jsta,jend
integer(i_kind),intent(in)::nx,ny
integer(i_kind),intent(in)::mype
nflds=26 !order is ps,t,q,u,v,w,w2,wdir,td,gust,vis,pblh,dist,wspd10m,td2m,mxtm,mitm,pmsl,howv,tcamt,lcbas,uwnd10m,vwnd10m,altw,altw2,altwdir
kflds=9 !perform cressman analysis for kflds fields only
!(currently: u,v,t,ps,q,gust,vis,pblh,dist)
if (mype==0) print*,'in alloc_cressanl_common: nmax=',nmax
allocate(ipointer(nmax,nflds))
allocate(jpointer(nmax,kflds))
allocate(uvpointer(nmax,2))
allocate(xlocs(nobsmax))
allocate(ylocs(nobsmax))
allocate(hgt0s(nobsmax))
allocate(hgts(nobsmax))
allocate(hobs(nobsmax))
allocate(rmuses(nobsmax))
allocate(oberrs(nobsmax))
allocate(dtimes(nobsmax))
allocate(rtvts(nobsmax))
allocate(cstations(nobsmax))
allocate(obstypes(nobsmax))
allocate(insubdoms(nobsmax))
allocate(dups(nobsmax))
allocate(bckgs(ista:iend,jsta:jend,kflds))
allocate(xberrs(ista:iend,jsta:jend,kflds))
allocate(terrain(nx,ny))
end subroutine alloc_cressanl_common
!*****************************************************************************************
!*****************************************************************************************
!*****************************************************************************************
!*****************************************************************************************
subroutine destroy_cressanl_common
implicit none
deallocate(ipointer)
deallocate(jpointer)
deallocate(uvpointer)
deallocate(xlocs)
deallocate(ylocs)
deallocate(hgt0s)
deallocate(hgts)
deallocate(hobs)
deallocate(rmuses)
deallocate(oberrs)
deallocate(dtimes)
deallocate(rtvts)
deallocate(cstations)
deallocate(obstypes)
deallocate(insubdoms)
deallocate(dups)
deallocate(bckgs)
deallocate(xberrs)
deallocate(terrain)
end subroutine destroy_cressanl_common
!*****************************************************************************************
!*****************************************************************************************
!*****************************************************************************************
!*****************************************************************************************
subroutine make_cressanl_common(cgrid, &
ista,iend,jsta,jend,ribuffer_km,mype,npe)
!
!*******************************************************************************
use mpi
implicit none
!Declare passed variables
integer(4),intent(in):: mype,npe
integer(4),intent(in):: ista,iend,jsta,jend
real(4),intent(in):: ribuffer_km
character(60),intent(in):: cgrid
!Declare local parameters
integer(4),parameter::ntypes=19! number of "otype"s
!Declare local variables
real(8) alat1,elon1,dx,elonv,alatan
real(4) dg2rad,elonv4,alatan4,xn
integer(4) idate,nchar,nreal,i,ii,mypegsi
character(8),allocatable,dimension(:):: cdiagbuf
character(8),allocatable,dimension(:):: cprvstg
character(8),allocatable,dimension(:):: csprvstg
real(4),allocatable,dimension(:,:)::rdiagbuf
real(4) rlat,rlon,oberr,ob, &
rmuse,xx,yy
real(4) uob,vob,qtflg
real(4) uobe,vobe
real(4) wobe,wdire
real(4) t0,td0
real(4) p0
real(4) hgt0 !observation elevation
real(4) hgt !model terrain at ob location
real(4) dtime
real(4) rtvts
real(4) bb
character(3) otype
character(12) obstype
real(8) rlat8,rlon8,xx8,yy8
character(8) cstation
logical fexist, r15887_diagfile_fmt
logical lboundtd
logical lambconform
logical polarstereo
integer(4) nn(ntypes)
integer(4) k,n,n0
integer(4) nmax0
real(4) xsta,xend,ysta,yend
real(4) xsta2,xend2
integer(4) ibuffer
logical insubdom
logical insubdom2
integer(4) ierror
integer(4) nx,ny
real(4) da
namelist/diagfile_fmt/r15887_diagfile_fmt
!
!*******************************************************************************
data r15887_diagfile_fmt/.true./
inquire(file='diagfile_fmt_input',exist=fexist)
if (fexist) then
open (55,file='diagfile_fmt_input',form='formatted')
read(55,*) r15887_diagfile_fmt
close(55)
endif
if (mype==0) print*,'in make_cressanl_common: r15887_diagfile_fmt=', &
r15887_diagfile_fmt
call proj_info(cgrid,dx,alat1,elon1,elonv,alatan)
if (mype==0) then
print*,'in make_cressanl_common: cgrid=',trim(cgrid)
print*,'in make_cressanl_common: alat1=',alat1
print*,'in make_cressanl_common: elon1=',elon1
print*,'in make_cressanl_common: dx=',dx
print*,'in make_cressanl_common: elonv=',elonv
print*,'in make_cressanl_common: alatan=',alatan
endif
call domain_dims(cgrid,nx,ny,da)
if (mype==0) then
print*,'in make_cressanl_common: cgrid=',trim(cgrid)
print*,'in make_cressanl_common: nx,ny,da=',nx,ny,da
endif
lambconform=trim(cgrid)=='conus'.or.trim(cgrid)=='cohres'.or. &
trim(cgrid)=='hrrr'.or.trim(cgrid)=='cohresext'.or. &
trim(cgrid)=='cohreswexp'
polarstereo=trim(cgrid)=='alaska'.or.trim(cgrid)=='akhres'.or. &
trim(cgrid)=='juneau'
dg2rad=atan(1.)/45.
elonv4=elonv
alatan4=alatan
xn=1.
if (lambconform) xn=sin(alatan4*dg2rad)
if (mype==0) print*,'in make_cressanl_common: lambconform,polarstereo,xn=', &
lambconform,polarstereo,xn
print*,'in make_cressanl_common: mype,npe,ista,iend,jsta,jend=', &
mype,npe,ista,iend,jsta,jend
if (mype==0) print*,'in make_cressanl_common: ribuffer_km=',ribuffer_km
xsta=float(ista)
xend=float(iend)
ysta=float(jsta)
yend=float(jend)
!====================================================================================
!==>
!====================================================================================
ibuffer=1000.*ribuffer_km/da
if (mype==0) print*, &
'in make_cressanl_common: number of ibuffer points in x: ', ibuffer
xsta2=float(max(1,ista-ibuffer))
xend2=float(min(nx,iend+ibuffer))
print*,'in make_cressanl_common: mype,xsta,xend,xsta2,xend2=', &
mype,xsta,xend,xsta2,xend2
!====================================================================================
!==> determine nobsmax and nmax in diagnostic file
!====================================================================================
open (7,file='diag_conv_ges.dat',form='unformatted')
read(7) idate
nobsmax=0 ; nn(:)=0
read_step_one: do
read(7,end=75 )otype,nchar,nreal,ii,mypegsi
if (ii==0) then
read(7,end=75)
if (r15887_diagfile_fmt) read(7,end=75)
cycle read_step_one
endif
allocate(cdiagbuf(ii))
allocate(cprvstg(ii))
allocate(csprvstg(ii))
allocate(rdiagbuf(nreal,ii))
if (r15887_diagfile_fmt) then
read(7) cdiagbuf,rdiagbuf
read(7) cprvstg,csprvstg
else
read(7) cdiagbuf,rdiagbuf,cprvstg,csprvstg
endif
if (otype(1:3)==' t') then ; n0=1 ; k=1 ; endif
if (otype(1:3)==' q') then ; n0=2 ; k=2 ; endif !q, td
if (otype(2:3)=='ps') then ; n0=1 ; k=3 ; endif
if (otype(2:3)=='uv') then
n0=5 !u, v, w, w2, wd
if (lambconform .or. polarstereo) n0=7 !plus ugrel, vgrel
k=4
endif
if (otype(1:3)=='spd') then ; n0=1 ; k=5 ; endif
if (otype(1:3)=='gst') then ; n0=1 ; k=6 ; endif
if (otype(1:3)=='vis') then ; n0=1 ; k=7 ; endif
if (otype(1:3)=='pbl') then ; n0=1 ; k=8 ; endif
if (otype(1:3)=='cei') then ; n0=1 ; k=9 ; endif
if (otype(1:3)=='wst') then ; n0=1 ; k=10 ; endif
if (otype(1:3)=='td2') then ; n0=1 ; k=11 ; endif
if (otype(1:3)=='mxt') then ; n0=1 ; k=12 ; endif
if (otype(1:3)=='mit') then ; n0=1 ; k=13 ; endif
if (otype(1:3)=='psl') then ; n0=1 ; k=14 ; endif
if (otype(1:3)=='hwv') then ; n0=1 ; k=15 ; endif
if (otype(1:3)=='tca') then ; n0=1 ; k=16 ; endif
if (otype(1:3)=='lcb') then ; n0=1 ; k=17 ; endif
if (otype(1:3)=='uwn') then ; n0=1 ; k=18 ; endif
if (otype(1:3)=='vwn') then ; n0=4 ; k=19 ; endif !vwn, altw, altw2, altwd
do i=1,ii
rlat=rdiagbuf(3,i) ; rlat8=rlat*1._8
rlon=rdiagbuf(4,i) ; rlon8=rlon*1._8
oberr=rdiagbuf(16,i)
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
nobsmax=nobsmax+n0
nn(k)=nn(k)+1
endif
enddo
deallocate(cdiagbuf,rdiagbuf)
deallocate(cprvstg,csprvstg)
enddo read_step_one
75 continue
print*,'in make_cressanl_common: nobsmax=',nobsmax
nmax0=maxval(nn(:))
call mpi_allreduce(nmax0,nmax,1,mpi_integer,mpi_max,mpi_comm_world,ierror)
if (mype==0) print*,'in make_cressanl_common: nmax=',nmax
!
!====================================================================================
!==> allocate cressman fields
!====================================================================================
!
call alloc_cressanl_common(ista,iend,jsta,jend,nx,ny,mype)
!
!====================================================================================
!==> read in from diagnostic file and load cressman arrays
!====================================================================================
rewind(7)
read(7) idate
print*,'in make_cressanl_common, idate=',idate
kps=0; kts=0; kqs=0; kus=0; kvs=0; kws=0; kw2s=0; kwds=0
ktds=0; kgusts=0; kvis=0; kpblhs=0; kdists=0
kugrds=0; kvgrds=0
kwspd10ms=0; ktd2ms=0; kmxtms=0; kmitms=0
kpmsls=0; khowvs=0; ktcamts=0; klcbas=0
kuwnd10ms=0; kvwnd10ms=0
kaltws=0; kaltw2s=0; kaltwds=0
n=0
read_step_two: do
read(7,end=200)otype,nchar,nreal,ii,mypegsi
! print*,'in make_cressanl_common, otype,nchar,nreal,ii,mypegsi=', &
! otype,nchar,nreal,ii,mypegsi
if (ii==0) then
read(7,end=200)
if (r15887_diagfile_fmt) read(7,end=200)
cycle read_step_two
endif
allocate(cdiagbuf(ii))
allocate(cprvstg(ii))
allocate(csprvstg(ii))
allocate(rdiagbuf(nreal,ii))
if (r15887_diagfile_fmt) then
read(7) cdiagbuf,rdiagbuf
read(7) cprvstg,csprvstg
else
read(7) cdiagbuf,rdiagbuf,cprvstg,csprvstg
endif
do i=1,ii
if (cprvstg(i)(1:4)=='B7Hv') cprvstg(i)(5:8)=' ' !this provider name comes with strange
if (csprvstg(i)(1:4)=='B7Hv') csprvstg(i)(5:8)=' ' !characters
enddo
bb=0. !dummy variable
if (otype(2:3)=='ps') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='p-ob'
n=n+1
kps=kps+1 ; ipointer(kps,1)=n
hgt=rdiagbuf(21,i)
oberr=oberr*100. !convert to Pa
ob=ob*100. !convert to Pa
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)==' t') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='t-ob'
n=n+1
kts=kts+1 ; ipointer(kts,2)=n
hgt=rdiagbuf(21,i)
rtvts=qtflg
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,rtvts,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)==' q') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='q-ob'
n=n+1
kqs=kqs+1 ; ipointer(kqs,3)=n
hgt=rdiagbuf(22,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
obstype='td-ob'
n=n+1
ktds=ktds+1 ; ipointer(ktds,9)=n
p0=rdiagbuf(6,i)
t0=-9999.
lboundtd=.false.
call get_stndewpt_r4(100.*p0,ob,t0,td0,lboundtd)
call load_cressanl(n,xx,yy,hgt0,hgt,td0,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(2:3)=='uv') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='u-ob'
n=n+1
kus=kus+1 ; ipointer(kus,4)=n
uobe=rdiagbuf(17,i)
hgt=rdiagbuf(25,i)
call load_cressanl(n,xx,yy,hgt0,hgt,uobe,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
obstype='v-ob'
n=n+1
kvs=kvs+1 ; ipointer(kvs,5)=n
vobe=rdiagbuf(20,i)
call load_cressanl(n,xx,yy,hgt0,hgt,vobe,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
obstype='w-ob'
n=n+1
kws=kws+1 ; ipointer(kws,6)=n
call speeddir(uobe,vobe,wdire,wobe)
call load_cressanl(n,xx,yy,hgt0,hgt,wobe,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
obstype='w2-ob'
n=n+1
kw2s=kw2s+1 ; ipointer(kw2s,7)=n
call load_cressanl(n,xx,yy,hgt0,hgt,wobe,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
obstype='wd-ob'
n=n+1
kwds=kwds+1 ; ipointer(kwds,8)=n
call load_cressanl(n,xx,yy,hgt0,hgt,wdire,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
if (lambconform .or. polarstereo) then
call get_grid_winds(uobe,vobe,uob,vob,xn,elonv4,rlon,dg2rad)
obstype='ugrel-ob'
n=n+1
kugrds=kugrds+1 ; uvpointer(kugrds,1)=n
call load_cressanl(n,xx,yy,hgt0,hgt,uob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
obstype='vgrel-ob'
n=n+1
kvgrds=kvgrds+1 ; uvpointer(kvgrds,2)=n
call load_cressanl(n,xx,yy,hgt0,hgt,vob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
endif
enddo
elseif (otype(1:3)=='spd') then !not used in cressman analysis
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='spd-ob'
n=n+1
! no pointer
hgt=rdiagbuf(22,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='gst') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='gust-ob'
n=n+1
kgusts=kgusts+1 ; ipointer(kgusts,10)=n
hgt=rdiagbuf(22,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='vis') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='vis-ob'
n=n+1
kvis=kvis+1 ; ipointer(kvis,11)=n
hgt=rdiagbuf(22,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='pbl') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='pblh-ob'
n=n+1
kpblhs=kpblhs+1 ; ipointer(kpblhs,12)=n
hgt=-9999.
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='cei') then !"dist" is actually cloud ceiling
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='dist-ob'
n=n+1
kdists=kdists+1 ; ipointer(kdists,13)=n
hgt=-9999.
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='wst') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='wspd10m-ob'
n=n+1
kwspd10ms=kwspd10ms+1 ; ipointer(kwspd10ms,14)=n
hgt=rdiagbuf(22,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='td2') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='td2m-ob'
n=n+1
ktd2ms=ktd2ms+1 ; ipointer(ktd2ms,15)=n
hgt=rdiagbuf(21,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='mxt') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='mxtm-ob'
n=n+1
kmxtms=kmxtms+1 ; ipointer(kmxtms,16)=n
hgt=rdiagbuf(21,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='mit') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='mitm-ob'
n=n+1
kmitms=kmitms+1 ; ipointer(kmitms,17)=n
hgt=rdiagbuf(21,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='psl') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='pmsl-ob'
n=n+1
kpmsls=kpmsls+1 ; ipointer(kpmsls,18)=n
hgt=rdiagbuf(21,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='hwv') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='howv-ob'
n=n+1
khowvs=khowvs+1 ; ipointer(khowvs,19)=n
hgt=rdiagbuf(21,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='tca') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='tcamt-ob'
n=n+1
ktcamts=ktcamts+1 ; ipointer(ktcamts,20)=n
hgt=rdiagbuf(22,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='lcb') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='lcbas-ob'
n=n+1
klcbas=klcbas+1 ; ipointer(klcbas,21)=n
hgt=rdiagbuf(22,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='uwn') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='uwnd10m-ob'
n=n+1
kuwnd10ms=kuwnd10ms+1 ; ipointer(kuwnd10ms,22)=n
hgt=rdiagbuf(22,i)
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
elseif (otype(1:3)=='vwn') then
do i=1,ii
cstation=cdiagbuf(i)
rlat=rdiagbuf(3,i)
rlon=rdiagbuf(4,i)
dtime=rdiagbuf(8,i)
hgt0=rdiagbuf(7,i)
qtflg=rdiagbuf(10,i)
if (rdiagbuf(11,i) .gt. 1.) then
rmuse=rdiagbuf(11,i)*rdiagbuf(12,i)
else
rmuse=rdiagbuf(12,i)
endif
oberr=rdiagbuf(16,i)
ob=rdiagbuf(17,i)
uobe=rdiagbuf(20,i)
rlat8=rlat*1._8
rlon8=rlon*1._8
call latlon_to_grid(dx,alat1,elon1,elonv,alatan, &
rlat8,rlon8,cgrid,xx8,yy8)
xx=xx8
yy=yy8
insubdom=xx.ge.xsta.and.xx.lt.(xend+1.).and.yy.ge.ysta.and.yy.lt.(yend+1.)
insubdom2=xx.ge.xsta2.and.xx.le.xend2.and.yy.ge.ysta.and.yy.le.yend
if (insubdom2 .and. oberr.gt.1.e-05) then
oberr=1./oberr
obstype='vwnd10m-ob'
n=n+1
kvwnd10ms=kvwnd10ms+1 ; ipointer(kvwnd10ms,23)=n
hgt=rdiagbuf(22,i)
bb=uobe !note this , note this , note this , note this , note this , note this!!!
call load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
obstype='altw-ob'
n=n+1
kaltws=kaltws+1 ; ipointer(kaltws,24)=n
call speeddir(uobe,ob,wdire,wobe)
call load_cressanl(n,xx,yy,hgt0,hgt,wobe,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
obstype='altw2-ob'
n=n+1
kaltw2s=kaltw2s+1 ; ipointer(kaltw2s,25)=n
call load_cressanl(n,xx,yy,hgt0,hgt,wobe,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
obstype='altwd-ob'
n=n+1
kaltwds=kaltwds+1 ; ipointer(kaltwds,26)=n
call load_cressanl(n,xx,yy,hgt0,hgt,wdire,rmuse, &
oberr,dtime,bb,cstation,obstype,insubdom)
endif
enddo
endif
deallocate(cdiagbuf,rdiagbuf)
deallocate(cprvstg,csprvstg)
enddo read_step_two
200 continue
close(7)
print*,' in make_cressanl_common: kps,kts,kqs,kus,kvs,kugrds,kvgrds,kws,kw2s,kwds,&
&ktds,kgusts,kvis,kpblhs,kdists=',&
kps,kts,kqs,kus,kvs,kugrds,kvgrds,kws,kw2s,kwds,&
ktds,kgusts,kvis,kpblhs,kdists
print*,' in make_cressanl_common: kwspd10ms,ktd2ms,kmxtms,kmitms,kpmsls,khowvs,ktcamts,klcbas=',&
&kwspd10ms,ktd2ms,kmxtms,kmitms,kpmsls,khowvs,ktcamts,klcbas
print*,' in make_cressanl_common: kuwnd10ms,kvwnd10ms=',kuwnd10ms,kvwnd10ms
print*,' in make_cressanl_common: kaltws,kaltw2s,kaltwds=',kaltws,kaltw2s,kaltwds
print*,' in make_cressanl_common: mype, number of obs,nobsmax=',mype,n,nobsmax
if (n /= nobsmax) then
print*,'in make_cressanl_common: n should be equal to nobsmax ... aborting'
call mpi_finalize(ierror)
stop
endif
jpointer(:,1)=uvpointer(:,1)
jpointer(:,2)=uvpointer(:,2)
jpointer(:,3)=ipointer(:,2)
jpointer(:,4)=ipointer(:,1)
jpointer(:,5)=ipointer(:,3)
jpointer(:,6)=ipointer(:,10)
jpointer(:,7)=ipointer(:,11)
jpointer(:,8)=ipointer(:,12)
jpointer(:,9)=ipointer(:,13)
end subroutine make_cressanl_common
!*****************************************************************************************
!*****************************************************************************************
!=========================================================================================
subroutine load_cressanl(n,xx,yy,hgt0,hgt,ob,rmuse,&
oberr,dtime,bb,cstation,obstype,insubdom)
implicit none
integer(4),intent(in)::n
real(4),intent(in):: xx,yy,hgt0,hgt,ob,rmuse,&
oberr,dtime,bb
character(8),intent(in):: cstation
character(12),intent(in):: obstype
logical,intent(in):: insubdom
!
xlocs(n)=xx
ylocs(n)=yy
hgt0s(n)=hgt0
hgts(n)=hgt
hobs(n)=ob
rmuses(n)=rmuse
oberrs(n)=oberr
dtimes(n)=dtime
rtvts(n)=bb
cstations(n)(1:8)=cstation(1:8)
obstypes(n)(1:12)=obstype(1:12)
insubdoms(n)=insubdom
end subroutine load_cressanl
!=========================================================================================
end module cressanl_common
!*****************************************************************************************
!*****************************************************************************************