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
!*****************************************************************************************
!*****************************************************************************************