      PROGRAM glassim
      include 'params.h'
      include 'files.h'

      
C***********************************************************************
C************       NASA/GODDARD SPACE FLIGHT CENTER        ************
C********              LABORATORY FOR ATMOSPHERES               ********
C******                                                           ******
C******          Data Assimilation Office  -  Code 910.3          ******
C********         Optimum Interpolation Analysis System         ********
C************     Frozen Version 1.5  -- October, 1993      ************
C***********************************************************************

      COMMON / contrl / im0, jnp0, nb0, ndalt0, nhms0, nlay0, nlaym10, n
     &layp10, nymd0, strat, idimen0, slp, fdcos, fexp, saveae, sigdat, m
     &axrep
      LOGICAL strat, slp, fdcos, fexp, saveae, sigdat, first, qbudget_sw
     &itch
      COMMON / obs / maxdex, nfgge, lun, hh (13, 13, 20),
     &vcoroe (18, 18, 4), obserr (mxqobs, 4, 19), ww (13, 13, 20),
     &qobs (mxqobs, 5), offlev (mxqobs)
      common / obschar / dataname (5), obsname(mxqobs)

      COMMON / scan / maxobs, rscan, onestage
      REAL rscan, data (MaxObservPoints * 10), sat_wind (MaxObservPoints
     &), pmand (nlevoi), pmidm (nlevoi - 1), phisz (imr, jmr + 1)
      LOGICAL onestage, repack, suspect (MaxObservPoints, 2), skippuv, s
     &kipmix, skipzuv, qobs
      INTEGER bq (MaxObservPoints)
      DATA first / .true. /
c  Mandatory pressure levels.
      DATA pmand / 20., 30., 50., 70., 100., 150., 200., 250., 300., 400
     &., 500., 700., 850., 1000. /
c  Mid-mandatory pressure levels.
      DATA pmidm / 24.9, 39.7, 59.8, 84.7, 124.4, 174.6, 224.7, 274.7, 3
     &49.1, 449.3, 598.0, 774.1, 924.3 /
      EXTERNAL mixcov, puvcov, zuvcov
C***********************************************************************
C*       FIELDS ON ANALYSIS LEVELS
C***********************************************************************
C
C     lwflag_oi (imr,jmr+1)              LAND/WATER FLAG
C     psfc_oi   (imr,jmr+1)              SURFACE PRESSURE
C     tsfc_oi   (imr,jmr+1)              SURFACE TEMPERATURE (EXTRAPOLAT
C     tvsfc_oi  (imr,jmr+1)              SURFACE VIRTUAL TEMPERATURE
C     psl_oi    (imr,jmr+1)              SEA LEVEL PRESSURE
C     usl_oi    (imr,jmr+1)              SEA LEVEL U WIND
C     vsl_oi    (imr,jmr+1)              SEA LEVEL V WIND
C     u_oi      (imr,jmr+1,nlevoi)       U WIND
C     v_oi      (imr,jmr+1,nlevoi)       V WIND
C     tv_oi     (imr,jmr+1,nlevoi)       VIRTUAL TEMPERATURE
C     t_oi      (imr,jmr+1,nlevoi)       DRY BULB TEMPERATURE
C     tgnd     (imr,jmr+1)              GROUND TEMPERATURE
C     rh_oi     (imr,jmr+1,nlevoi)       RELATIVE HUMIDITY
C     rmix_oi   (imr,jmr+1,nlevoi)       MIXING RATIO
C     z_oi      (imr,jmr+1,nlevoi)       GEOPOTENTIAL HEIGHT
      REAL lwflag_oi (imr, jmr + 1)
      REAL psfc_oi (imr, jmr + 1)
      REAL tsfc_oi (imr, jmr + 1)
      REAL tvsfc_oi (imr, jmr + 1)
      REAL psl_oi (imr, jmr + 1)
      REAL phis (imr, jmr)
      REAL usl_oi (imr, jmr + 1)
      REAL vsl_oi (imr, jmr + 1)
      REAL u_oi (imr, jmr + 1, nlevoi)
      REAL v_oi (imr, jmr + 1, nlevoi)
      REAL tv_oi (imr, jmr + 1, nlevoi)
      REAL tgnd (imr, jmr + 1)
      REAL t_oi (imr, jmr + 1, nlevoi)
      REAL rh_oi (imr, jmr + 1, nlevoi)
      REAL rmix_oi (imr, jmr + 1, nlevoi)
      REAL z_oi (imr, jmr + 1, nlevoi)
      REAL zlwflag (imr, jmr + 1)
C***********************************************************************
C*       ANALYSIS INCREMENTS
C***********************************************************************
      REAL usl_inc (imr, jmr + 1)
      REAL vsl_inc (imr, jmr + 1)
      REAL u_inc (imr, jmr + 1, nlevoi)
      REAL v_inc (imr, jmr + 1, nlevoi)
      REAL rmix_inc (imr, jmr + 1, nlevoi)
      REAL tv_inc (imr, jmr + 1, nlevoi)
C
C***********************************************************************
C*       FIRST GUESS/ANALYSIS ERROR STANDARD DEVIATIONS                *
C***********************************************************************
C     psl_err  (imr,jmr+1)              SEA LEVEL PRESSURE
C     tsl_err  (imr,jmr+1)              SEA LEVEL TEMPERATURE
C     rmix_err (imr,jmr+1,nlevoi)       MIXING RATIO
C     z_err    (imr,jmr+1,nlevoi)       GEOPOTENTIAL HEIGHT
C     u_err    (imr,jmr+1,nlevoi)       U-WIND                    <1.3>
C     v_err    (imr,jmr+1,nlevoi)       V_WIND                    <1.3>
C     usl_err  (imr,jmr+1)              SEA-LEVEL U-WIND          <1.3>
C     vsl_err  (imr,jmr+1)              SEA-LEVEL V-WIND          <1.3>
C**********************************************************************
      REAL psl_err (imr, jmr + 1)
      REAL tsl_err (imr, jmr + 1)
      REAL rmix_err (imr, jmr + 1, nlevoi)
      REAL z_err (imr, jmr + 1, nlevoi)
      REAL u_err (imr, jmr + 1, nlevoi)
      REAL v_err (imr, jmr + 1, nlevoi)
      REAL usl_err (imr, jmr + 1)
      REAL vsl_err (imr, jmr + 1)
C***********************************************************************
C*       LOCAL VARIABLES                                               *
C***********************************************************************
C
      DATA fill / 1.e10 /
      CHARACTER * 8 atime, adate, xlabel (10), err_exp_id
      CHARACTER * 4 job (2), exp_id (2)
      CHARACTER * 60 message
C  GRID VOLUME SETUP FOR OI VERSIONS 1.2 AND LATER
C
C     Grid volume information is contained in dataset 'mvtable', unit 88
C     This dataset contains the following:
C
C     NAME            TYPE                DESCRIPTION
C
C     ngv             simple int          Number of grid volumes
C     volcentre_x     array real          Longitude volume centres
C     volcentre_y     array real          Latitude volume centres
C     glatup          array real          Volume reference latitude
C     glonup          array real          Volume reference longitude
C     nh              array int           No. horiz gridpoints in volume
C
      COMMON / sparam / lsolve, ns, nr, nv, nh (13104), nq, kq (3), volc
     &entre_x (13104), ngv, volcentre_y (13104), glatup (13104), glonup 
     &(13104)

C***********************************************************************
C*       WORKSPACE                                                     *
C***********************************************************************
C
      real WKSPACE (IMR*(JMR+1)*(NLEVOI+1)*4 + 60000 + IMR*(JMR+1))
      common /space/ wkspace

      real,allocatable:: data_mm(:,:)
      integer,allocatable:: sspct_mm(:,:)

C***********************************************
C     Data layout
C*******************************************************
!hpf$ distribute data_mm(block,block)
!hpf$ distribute sspct_mm(block,block)
!hpf$ distribute bq(block)
!hpf$ distribute data(block)
!hpf$ distribute dataname(block)
!hpf$ distribute glatup(block)
!hpf$ distribute glonup(block)
!hpf$ distribute hh(block,block,block)
!hpf$ distribute kq(block)
!hpf$ distribute lwflag_oi(block,block)
!hpf$ distribute nh(block)
!hpf$ distribute obserr(block,block,block)
!hpf$ distribute obsname(block)
!hpf$ distribute offlev(block)
!hpf$ distribute phis(block,block)
!hpf$ distribute phisz(block,block)
!hpf$ distribute pmand(block)
!hpf$ distribute pmidm(block)
!hpf$ distribute psfc_oi(block,block)
!hpf$ distribute psl_err(block,block)
!hpf$ distribute psl_oi(block,block)
!hpf$ distribute qobs(block,block)
!hpf$ distribute rh_oi(block,block,block)
!hpf$ distribute rmix_err(block,block,block)
!hpf$ distribute rmix_inc(block,block,block)
!hpf$ distribute rmix_oi(block,block,block)
!hpf$ distribute sat_wind(block)
!hpf$ distribute suspect(block,block)
!hpf$ distribute tgnd(block,block)
!hpf$ distribute tsfc_oi(block,block)
!hpf$ distribute tsl_err(block,block)
!hpf$ distribute tvsfc_oi(block,block)
!hpf$ distribute tv_inc(block,block,block)
!hpf$ distribute tv_oi(block,block,block)
!hpf$ distribute t_oi(block,block,block)
!hpf$ distribute usl_err(block,block)
!hpf$ distribute usl_inc(block,block)
!hpf$ distribute usl_oi(block,block)
!hpf$ distribute u_err(block,block,block)
!hpf$ distribute u_inc(block,block,block)
!hpf$ distribute u_oi(block,block,block)
!hpf$ distribute vcoroe(block,block,block)
!hpf$ distribute volcentre_x(block)
!hpf$ distribute volcentre_y(block)
!hpf$ distribute vsl_err(block,block)
!hpf$ distribute vsl_inc(block,block)
!hpf$ distribute vsl_oi(block,block)
!hpf$ distribute v_err(block,block,block)
!hpf$ distribute v_inc(block,block,block)
!hpf$ distribute v_oi(block,block,block)
!hpf$ distribute ww(block,block,block)
!hpf$ distribute zlwflag(block,block)
!hpf$ distribute z_err(block,block,block)
!hpf$ distribute z_oi(block,block,block)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c   Workspace arrays
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      real,dimension(imr,jmr+1,NLEVOI+1,3):: wk_aquan
      real,dimension(imr,jmr+1,NLEVOI+1):: wk_asdae
      real, dimension(20000):: wk_ze,wk_zelat,wk_zelon
!hpf$ distribute wk_zelat(block)
!hpf$ distribute wk_zelon(block)
!hpf$ distribute wk_ze(block)
!hpf$ distribute wk_aquan(block,block,block,block)
!hpf$ distribute wk_asdae(block,block,block)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c    Files
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      OPEN( UNIT=89,FILE=ferrordata,FORM='unformatted',STATUS='old' ) 
      OPEN( UNIT=97,FILE=fanalysisdata,STATUS='old' ) 
      OPEN( UNIT=88,FILE=fminivolumetab,STATUS='old' ) 
c      OPEN( UNIT=79,FILE=fhuvdata,STATUS='old' ) 
c      OPEN( UNIT=78,FILE=fmixdata,STATUS='old' ) 
      OPEN( UNIT=77,FILE=fslpdata,STATUS='old' ) 
      OPEN( UNIT=65,FILE=falputznames,STATUS='old' ) 
      OPEN( UNIT=66,FILE=fobservation,STATUS='old' ) 
c     Optinal output files
      OPEN( UNIT=91,FILE=fdel,STATUS='replace' ) 
      OPEN( UNIT=95,FILE=fprs,STATUS='replace' ) 
      OPEN( UNIT=99,FILE=fcanary,STATUS='replace' ) 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C  Read the pressure level analysis restart file.
C  This file replaces the model and was created from a model restart
C  following sigma-to-pressure conversion.
cCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      PRINT *, 'Read previous cycle analysis data'
      READ( 97,* ) nymd,nhms
      PRINT *,nymd,nhms
      DO 7002 ii = 1 , imr
        DO 7000 jj = 1 , jmr + 1
          READ( 97,* ) lwflag_oi(ii,jj),tsfc_oi(ii,jj),psfc_oi(ii,jj),
     &        psl_oi(ii,jj),usl_oi(ii,jj),vsl_oi(ii,jj),phisz(ii,jj),
     &        tvsfc_oi(ii,jj)
 7000   CONTINUE
 7002 CONTINUE
      nymd0 = nymd
      nhms0 = nhms
c**********************************************************************
c
c       Set up analysis.
c
c**********************************************************************
      PRINT *, 'Initailization'
      repack = .TRUE.
      CALL initial(pmand,repack)
      CALL cortable(pmand,nlevoi)
      CALL trigset()
c*********************************************************************
c     For development phase do PUV analysis only
c********************************************************************
      skippuv = .FALSE.
      skipmix = .TRUE.
      skipzuv = .TRUE.
c
c**********************************************************************
c
c        Read dataset containing minivolume information.
c
c**********************************************************************
      write(6,*) 'Reading minivolume info'
      READ( 88,* ) ngv
      write(6,*) ngv
      DO 7010 ii = 1 , ngv
        READ( 88,* ) volcentre_x(ii),volcentre_y(ii),glatup(ii),
     &      glonup(ii),nh(ii)
 7010 CONTINUE
C
C***********************************************************************
C        WORKSPACE POINTERS AND ANALYSIS PARAMETERS
C***********************************************************************
C
      CALL getconh('DLAT    ',dlat)
C  Read analysis error from input file.
      print *,'  Read analysis error from input file'
      CALL rderror(nymd,nhms,psl_err,tsl_err,rmix_err,z_err,istart,
     &    exp_id)
      IF (istart .NE. 0) THEN
        print *,'  error growth'
        CALL growth(psl_err,tsl_err,rmix_err,z_err)
      ENDIF
      istart = 1
C Save pressure level first guess
      usl_inc=usl_oi
      vsl_inc=vsl_oi
c************************!!!!!!!!!!!!!!!!!!!!!*******************
c     These guys are  not initialized, so I commented them out, since 
c     they not used in PUV analysis anyway. 
c     This means tdrybulb does nothing useful here.
c      tv_inc=tv_oi
c      u_inc=u_oi
c      v_inc=v_oi
c      rmix_inc=rmix_oi
c************************!!!!!!!!!!!!!!!!!!!!!*******************
C
C***********************************************************************
C*       COMPUTATION OF PRESSURE LEVEL DRY BULB TEMPERATURES           *
C*       BEFORE ANALYSIS PRESSURE LEVEL HISTORY WRITE                  *
C***********************************************************************
C
      CALL tdrybulb(tv_oi,rmix_oi,t_oi,pmand,pmidm,nlevoi)
      print *, 'Finished tdrybulb'
C**********************************************************************
C  Compute the first guess wind standard deviations.
      print *, 'before werrinit'
      CALL werrinit(z_err,imr,jmr + 1,nlevoi,u_err,v_err)
      print *, 'Finished werrinit'
C**********************************************************************
C
C        Sea level pressure and wind analysis.
C
C*********************************************************************
      IF (.NOT.skippuv) THEN
c************************!!!!!!!!!!!!!!!!!!!!!*******************
c   This added ad hoc, since tgnd is not defined in the code we received
c   Again, dragon is completely useless here.
c************************!!!!!!!!!!!!!!!!!!!!!*******************
        tgnd=tsfc_oi
        CALL dragon(phisz,lwflag_oi,usl_oi,vsl_oi,tsfc_oi,tgnd)
        write(6,*) 'Finished dragon'
        nv = 1
        nq = 3
        CALL rdslp(data,nobs,suspect)
        allocate(data_mm(nobs,10))
        allocate(sspct_mm(nobs,2))
        data_mm(:,1) = data(1:nobs)
        data_mm(:,2) = data(1 + nobs:2*nobs)
        data_mm(:,3)=data(1 + 2 * nobs:3*nobs)
        data_mm(:,4)=data(1 + 3 * nobs:4*nobs)
        data_mm(:,5)=data(1 + 4 * nobs:5*nobs)
        data_mm(:,6)=data(1 + 5 * nobs:6*nobs)
        data_mm(:,7)=data(1 + 6 * nobs:7*nobs)
        data_mm(:,8)=data(1 + 7 * nobs:8*nobs)
        where(suspect(1:nobs,:))
            sspct_mm(1:nobs,:)=1
        elsewhere
            sspct_mm(1:nobs,:)=0
        endwhere

        IF (nobs .GT. 0) THEN
          CALL qcslp(phisz,psl_oi,usl_oi,vsl_oi,data_mm,nobs,t_oi,wk_ze,
     &        wk_zelat,wk_zelon,nze,psfc_oi,sspct_mm,psl_err)
c          CALL puvanl(psl_err,phisz,psl_oi,usl_oi,vsl_oi,tsfc_oi,
c     &        wk_aquan,wk_asdae,t_oi,wk_ze,wk_zelat,wk_zelon,nze,
c     &        psfc_oi,data_mm,nobs,sspct_mm)

        ELSE
          CALL flash(3,'GLASSIM: No slp observations.',1)
          nze = 0
        ENDIF
C
      ENDIF
C***********************************************************************
c     Compute analysis increments
C**********************************************************************
      WHERE (usl_inc .NE. fill) 
          usl_inc = usl_oi - usl_inc
      ENDWHERE
      WHERE (vsl_inc .NE. fill) 
          vsl_inc = vsl_oi - vsl_inc
      ENDWHERE

      STOP
      END
