      PROGRAM glassim
      include '/usr/include/cm/CMF_defs.h'
      
C  Set maximum allowable number of observations:
      INTEGER MaxObservPoints
      PARAMETER (MaxObservPoints = 150000)
      
      CHARACTER * 80 FErrorData, FAnalysisData, FMiniVolumeTab, FHuvData
      CHARACTER * 80 FMixData, FSlpData, FAlputzNames, FObservation
      CHARACTER * 80 FDel, FPrs, FCanary
      PARAMETER (FErrorData = 'emss005.err.data')
      PARAMETER (FAnalysisData = 'emss005.analysis.ascii')
      real cos
c
      PARAMETER (FMiniVolumeTab = 'mvtable.ascii')
c
      PARAMETER (FHuvData = 'huv.data.920101')
C
      PARAMETER (FMixData = 'mix.data.920101')
C***********************************************************************
      PARAMETER (FSlpData = 'slp.data.920101')
C************       NASA/GODDARD SPACE FLIGHT CENTER        ************
C********              LABORATORY FOR ATMOSPHERES               ********
C******                                                           ******
C******          Data Assimilation Office  -  Code 910.3          ******
C********         Optimum Interpolation Analysis System         ********
      PARAMETER (FAlputzNames = 'alputz.namelist')
C************     Frozen Version 1.5  -- October, 1993      ************
      PARAMETER (FObservation = 'obs')
C***********************************************************************
C
      PARAMETER (FDel = 'emss005.del')
      PARAMETER (FPrs = 'emss005.prs')
c     Optinal output files
      PARAMETER (FCanary = 'canary')
C  Set analysis resolution:
C   imr:  Grid points East-West
C   jmr:  Grid points North-South (minus 1)
C   nlevoi: Number of analysis levels
      PARAMETER (imr = 144, jmr = 90, nlevoi = 14, nlayr = 20)
      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
      PARAMETER (mxqobs = 200)
      COMMON / obs / maxdex, nfgge, lun, hh (13, 13, 20),
     &vcoroe (18, 18, 4), obserr (200, 4, 19), ww (13, 13, 20),
     &qobs (200, 5), offlev (200)
      common / obschar / dataname (5), obsname(200)

      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)
!hpf$ distribute bq(block)
!hpf$ distribute wk_aquan1(block,block,block,block)
!hpf$ distribute wk_aquan3(block,block,block,block)
!hpf$ distribute wk_asdae(block,block,block)
!hpf$ distribute wk_ze(block)
!hpf$ distribute wk_zelat(block)
!hpf$ distribute wk_zelon(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)
c!hpf$ distribute wkspace(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)
      DATA first / .true. /
      SAVE first, istart
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 ANLAYSIS 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     tgrnd     (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>
      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***********************************************************************
C*           WORKSPACE                                                 *
C***********************************************************************
C
      real wk_aquan1(imr , jmr, nlevoi + 1, 1),
     &     wk_aquan3(imr , jmr, nlevoi + 1, 3),
     &     wk_asdae(imr , jmr, nlevoi + 1),
     &     wk_ze(20000),wk_zelat(20000),wk_zelon(20000)
      COMMON / space / wk_aquan1, wk_aquan3,wk_asdae,
     &                 wk_ze,wk_zelat,wk_zelon
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)

      real,allocatable:: data_mm(:,:)
!hpf$ distribute data_mm(block,block)

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c    Files
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      INTEGER Ti1, Ti2, Ti3
      OPEN( UNIT=89,FILE=ferrordata,FORM='unformatted',STATUS='new' ) 
      OPEN( UNIT=97,FILE=fanalysisdata,STATUS='new' ) 
      OPEN( UNIT=88,FILE=fminivolumetab,STATUS='new' ) 
      OPEN( UNIT=79,FILE=fhuvdata,STATUS='new' ) 
      OPEN( UNIT=78,FILE=fmixdata,STATUS='new' ) 
      OPEN( UNIT=77,FILE=fslpdata,STATUS='new' ) 
      OPEN( UNIT=65,FILE=falputznames,STATUS='new' ) 
      OPEN( UNIT=66,FILE=fobservation,STATUS='new' ) 
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.
      PRINT *, 'GVL: Read the Analysis data'
      READ( 97,* ) nymd,nhms
      PRINT *, 'GVL: ',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
      DO 7008 kk = 1 , nlevoi
        DO 7006 ii = 1 , imr
          DO 7004 jj = 1 , jmr + 1
            READ( 97,* ) u_oi(ii,jj,kk),v_oi(ii,jj,kk),tv_oi(ii,jj,kk),
     &          z_oi(ii,jj,kk),rmix_oi(ii,jj,kk),rh_oi(ii,jj,kk)
 7004     CONTINUE
 7006   CONTINUE
 7008 CONTINUE
      nymd0 = nymd
      nhms0 = nhms
c**********************************************************************
c
c       Set up analysis.
c
c**********************************************************************
      PRINT *, 'GVL: Do the initailization'
      repack = .TRUE.
      CALL initial(imr,jmr + 1,nlayr,nlevoi,pmand,repack)
      CALL cortable(pmand,nlevoi)
      CALL trigset()
      skippuv = .FALSE.
      skipmix = .TRUE.
      skipzuv = .TRUE.
c
c**********************************************************************
c
c        Read dataset containing minivolume information.
c
c**********************************************************************
      PRINT *, 'GVL: Read the Minivolume data set'
      READ( 88,* ) 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
      iaquan = 1
      iasdae = iaquan + imr * (jmr + 1) * (nlevoi + 1) * 3
      ize = iasdae + imr * (jmr + 1) * (nlevoi + 1)
      izelon = ize + 20000
      izelat = izelon + 20000
      iz1000 = izelat + 20000
CCC   NEXT   = IZ1000 + imr*(jmr+1)
      CALL getconh('DLAT    ',dlat)
C  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
        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
      u_inc=u_oi
      v_inc=v_oi
      rmix_inc=rmix_oi
      tv_inc=tv_oi
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)
C**********************************************************************
C  Compute the first guess wind standard deviations.
      CALL werrinit(z_err,imr,jmr + 1,nlevoi,u_err,v_err)
C**********************************************************************
C
C        Sea level pressure and wind analysis.
C
C*********************************************************************
      IF (.NOT.skippuv) THEN
        CALL dragon(phisz,lwflag_oi,usl_oi,vsl_oi,tsfc_oi,tgz)
        nv = 1
        nq = 3
        CALL rdslp(data,nobs,suspect)
        allocate(data_mm(nobs,10))
      DO i = 1 , nobs
        data_mm(i,1) = data(i)
        data_mm(i,2) = data(i + nobs)
        data_mm(i,3)=data(i + 2 * nobs)
        data_mm(i,4)=data(i + 3 * nobs)
        data_mm(i,5)=data(i + 4 * nobs)
        data_mm(i,6)=data(i + 5 * nobs)
        data_mm(i,7)=data(i + 6 * nobs)
        data_mm(i,8)=data(i + 7 * nobs)
      ENDDO

        IF (nobs .GT. 0) THEN
          CALL qcslp(phisz,psl_oi,usl_oi,vsl_oi,data,nobs,t_oi,wk_ze,
     &        wk_zelat,wk_zelon,nze,psfc_oi,suspect,psl_err)
          CALL puvanl(psl_err,phisz,psl_oi,usl_oi,vsl_oi,tsfc_oi,puvcov,
     &        wk_aquan3,wk_asdae,t_oi,wk_ze,wk_zelat,wk_zelon,nze,
     &        psfc_oi,data,nobs,suspect)
         ELSE
          CALL flash(3,'GLASSIM: No slp observations.',1)
          nze = 0
        ENDIF
      ELSE
        IF (repack) THEN
          DO 7012 logical_units = 1 , nfgge
            CALL rdthru(2,lun - 1 + logical_units)
 7012     CONTINUE
        ENDIF
      ENDIF
C
C**********************************************************************
C
C      Upper air water vapour mixing ratio analysis
C
C**********************************************************************
      IF (.NOT.skipmix) THEN
        nv = 2
        nq = 1
        mmix = 6
        CALL rdmix(nymd,nhms,data,bq,nobs,pmand,psfc_oi,imr,jmr + 1,
     &      nlevoi,mmix)
        IF (nobs .GT. 0) THEN
          CALL delmix(nymd,nhms,rh_oi,rmix_oi,pmand,psfc_oi,data,nobs,
     &        imr,jmr + 1,nlevoi,mmix,sigdat,bq,t_oi)
          CALL mixanl(rmix_err,rh_oi,rmix_oi,pmand,psfc_oi,mixcov,
     &        wkspace(iaquan),wkspace(iasdae),imr,jmr + 1,nlevoi,mmix,
     &        nobs,data,bq)
        ELSE
          CALL flash(3,'GLASSIM: No mix observations.',1)
        ENDIF
      ELSE
        IF (repack) THEN
          DO 7014 logical_units = 1 , nfgge
            CALL rdthru(3,lun - 1 + logical_units)
 7014     CONTINUE
        ENDIF
      ENDIF
C
C**********************************************************************
C
C      Upper air height and wind analysis
C
C**********************************************************************
C
      IF (.NOT.skipzuv) THEN
C
C  Recalculation of 1000mb heights using analyzed sea level analysis.
C
        CALL chgz1000(psfc_oi,tvsfc_oi,psl_oi,phisz,wkspace(iz1000),
     &      pmand,imr,jmr + 1,nlevoi)
C       print *,pmand
        nv = 2
        nq = 3
        CALL rdhuv(z_oi,tv_oi,wkspace(iz1000),data,nobs,suspect,wkspace(
     &      ize),wkspace(izelat),wkspace(izelon),nze)
        IF (nobs .GT. 0) THEN
          CALL delhuv(z_oi,u_oi,v_oi,tv_oi,psfc_oi,wkspace(iz1000),z_err
     &        ,u_err,v_err,data,nobs,wkspace(ize),wkspace(izelat),
     &        wkspace(izelon),nze,nymd,nhms,suspect)
          CALL zuvanl(z_err,z_oi,u_oi,v_oi,t_oi,psfc_oi,wkspace(iz1000),
     &        zuvcov,wkspace(iaquan),wkspace(iasdae),wkspace(ize),
     &        wkspace(izelat),wkspace(izelon),nze,data,nobs)
        ELSE
          CALL flash(3,'GLASSIM: No zuv observations.',1)
        ENDIF
      ENDIF
      CALL prptly()
C
C***********************************************************************
C*       COMPUTE ANALYSIS INCREMENTS                                   *
C***********************************************************************
      DO i = 1 , imr * (jmr + 1)
        IF (usl_inc(i,1) .NE. fill) THEN
          usl_inc(i,1) = usl_oi(i,1) - usl_inc(i,1)
        ENDIF
      ENDDO
      DO i = 1 , imr * (jmr + 1)
        IF (vsl_inc(i,1) .NE. fill) THEN
          vsl_inc(i,1) = vsl_oi(i,1) - vsl_inc(i,1)
        ENDIF
      ENDDO
      DO i = 1 , imr * (jmr + 1) * nlevoi
        IF (u_inc(i,1,1) .NE. fill) THEN
          u_inc(i,1,1) = u_oi(i,1,1) - u_inc(i,1,1)
        ENDIF
      ENDDO
      DO i = 1 , imr * (jmr + 1) * nlevoi
        IF (v_inc(i,1,1) .NE. fill) THEN
          v_inc(i,1,1) = v_oi(i,1,1) - v_inc(i,1,1)
        ENDIF
      ENDDO
      DO i = 1 , imr * (jmr + 1) * nlevoi
        IF (rmix_inc(i,1,1) .NE. fill) THEN
          rmix_inc(i,1,1) = rmix_oi(i,1,1) - rmix_inc(i,1,1)
        ENDIF
      ENDDO
      CALL ztot(z_oi,tv_oi,pmand,pmidm,imr,jmr + 1,nlevoi)
      DO i = 1 , imr * (jmr + 1) * nlevoi
        tv_inc(i,1,1) = tv_oi(i,1,1) - tv_inc(i,1,1)
      ENDDO
C
C***********************************************************************
C*       COMPUTATION OF PRESSURE LEVEL DRY BULB TEMPERATURES           *
C*       AFTER ANALYSIS PRESSURE LEVEL HISTORY WRITE                   *
C***********************************************************************
C
      CALL tdrybulb_v1(tv_oi,rmix_oi,t_oi,pmand,1,nlevoi,pmidm,1,nlevoi 
     &    - 1)
      CALL werrinit(z_err,imr,jmr + 1,nlevoi,u_err,v_err)
      CALL mixtorh(rmix_oi,t_oi,pmand,imr,jmr + 1,nlevoi,rh_oi,fill)
      DO 7032 i = 1 , imr
        DO 7030 j = 1 , jmr + 1
          DO 7028 k = 1 , nlevoi
            IF (rh_oi(i,j,k) .NE. fill) THEN
              rh_oi(i,j,k) = rh_oi(i,j,k) * 100.
            ENDIF
 7028     CONTINUE
 7030   CONTINUE
 7032 CONTINUE
C***********************************************************************
C*       WRITE ANALYSIS ERROR STANDARD DEVIATIONS                      *
C***********************************************************************
C
      REWIND 89
      CALL oierrwrt(89,exp_id,nymd,nhms,psl_err,tsl_err,rmix_err,z_err,
     &    imr,jmr + 1,nlevoi,wkspace,istat)
      IF (istat .NE. 0) THEN
        PRINT *, ' FATAL: GLASSIM> Error in analysis error data conver',
     &      'sion.'
        STOP 'Fatal.'
      ENDIF
      END
