c
      program flxproc
c
c Coupled processor -- plots tgcm/t21 or tgcm/maccm2
c
      include "flxproc.h"
      include "ccm.h"
      include "pltopt.h"
      character*16 modelstr
      real fglb(imx,jmx)
      data ihpabort/1/  ! if 0, hpdeallc will not abort on error
c
      call assign('assign -f77 u:5')	! for f77 style namelist read
c
c Initialization:
c
      call preset
c
c User proc input:
c
      call getinp 
c
c Open ncarg workstation(s) and get corresponding gks workstation id's:
c (wkstation id's are in common in pltopt.h)
c
      if (mkplt.gt.0) then
        if (lenstr(outplt(1)).eq.0.and.lenstr(outplt(2)).eq.0) 
     +    iwk_cgm = 1
        if (outplt(1).eq.'cgm     '.or.outplt(2).eq.'cgm     ') 
     +    iwk_cgm = 1
        if (outplt(1).eq.'ps      '.or.outplt(2).eq.'ps      ')
     +    iwk_ps = 1
c
c Initialize HLU package, create application object, and needed
c workstations for cgm and/or ps: 
c
        call opnwrk(appid,flnm_cgm,iwk_cgm,igks_cgm,
     +    flnm_ps,iwk_ps,igks_ps,psmode,icolor)
      endif
c
c Set global plot parameters:
c
      call setplt(spval)
c
c Time (history) loop:
c
      istart = 1
      iprint = 1
      do it=1,ntms
c       write(6,"(/'flxproc: time ',i3,' of ',i3,': tgcm time ',
c    +    i3,':',i2,':',i2)") it,ntms,(mtimes(i,it),i=1,3)
c
c Get selected tgcm fields:
c (ixfhist(ixz) has to be reset according to modelhts after chflds)
c
        ier = 0
        call chflds(ixfhist,mxfhist,ifproc,mxfproc,ifdep,
     +    mxfproc-mxfhist,1,flab8,1)
        ixfhist(ixz) = 1			! always get z
        if (modelhts.gt.0) ixfhist(ixz) = 2
        call rdtgcm(histvols,nhvols,tmpdir,istart,mtimes(1,it),pfhist,
     +    ixfhist,mxfhist,imx,kmx,jmx,nfhist,zp1,dzp,iden,ionvel,luhist,
     +    iprint,ivol,isdyn,istimes,issech,iier)
        if (ier.eq.1) goto 100 ! history not found
        if (ier.ne.0) then
          write(6,"('>>> flxproc: Error ',i3,' from rdtgcm')") ier
          stop 'rdtgcm'
        endif
        call chflds(ixfhist,mxfhist,ifproc,mxfproc,ifdep,
     +    mxfproc-mxfhist,0,flab8,1)
        nfproc = 0
        do ip=1,mxfproc
          lcfield = lenstr(cfields(ip))
          if (lcfield.gt.0) then
            ifld = ixfindc(flab8,mxfproc,cfields(ip))
            if (ifproc(ifld).gt.0) then
              nfproc = nfproc+1
              ifproc(ifld) = nfproc
            endif
          endif
        enddo
        if (nfproc.eq.0) then
          write(6,"(/'>>> No valid field names to process <<<')")
          write(6,"('Field names given by cfields input parameter',
     +      ' were invalid or not available on the histories'/)")
          stop 'fields'
        endif
        if (dzp.gt..5-.0001.and.dzp.lt..5+.0001) dzp = .50
        model = 'TIGCM           '
        if (.not.istimes.and.isdyn) model = 'TIEGCM          '
        if (istimes) model = 'TIME-GCM        '	! should always be true
        ut = float(mtimes(2,it)) + float(mtimes(3,it))/60.
c       call rptf(fhist,imx,kmx,jmx,nfhist,ixfhist,flab8,mxfhist)
c
c Get all ccm fields at tgcm horizontal grid and ccm2 zp levels:
c fccm(ngcmlon,nccmlev,ngcmlat,nccmfld) (pointer pfccm is in common)
c
        call rdlsd(imx,gcmlon,jmx,gcmlat,p0)
        if (it.eq.1)
     +    write(6,"(/'flxproc: nccmlev=',i2,' ccmlev_mb=',/(5e13.5))")
     +      nccmlev,(ccmlev_mb(kk),kk=1,nccmlev)
c
c Now have p0, so set up gcmpmb (also defind gcmzp):
c
        call alloc(pgcmzp,kmx)
        call alloc(pgcmpmb,kmx)
        do k=1,kmx
          gcmzp(k) = zp1+(k-1)*dzp
          gcmpmb(k) = p0*exp(-gcmzp(k))*1.e-3
        enddo
c
c Now we know t21 or maccm (ccmres). Use kmx to confirm that the tgcm 
c   history read at least has the vertical dimension of a history that
c   was coupled with the determined ccm version:
c tgcm history from a tgcm/t21   coupled run: gcmzp = -16 to +5 (kmx=43)
c tgcm history from a tgcm/maccm coupled run: gcmzp = -11 to +5 (kmx=33)
c
        if (ccmres.eq.'CCM2-T21        ') then
          if (kmx.ne.43) then
            write(6,"(/'>>> ccm model = ',a,' but tgcm kmx=',i2,
     +        ' --> wrong type of tgcm history? (kmx should = 43)')") 
     +        ccmres(1:lenstr(ccmres)),kmx
            write(6,"('nccmlev=',i2,' ccmlev_zp=',/(8f8.2))")
     +        nccmlev,(ccmlev_zp(k),k=1,nccmlev)
            write(6,"('kmx=',i2,' gcmzp=',/(8f8.2))")
     +        kmx,(gcmzp(k),k=1,kmx)
            stop 'tgcmhist'
          endif 
        elseif (ccmres.eq.'CCM2-MACCM      ') then
          if (kmx.ne.33) then
            write(6,"(/'>>> ccm model = ',a,' but tgcm kmx=',i2,
     +        ' --> wrong type of tgcm history? (kmx should = 33)')") 
     +        ccmres(1:lenstr(ccmres)),kmx
            write(6,"('nccmlev=',i2,' ccmlev_zp=',/(8f8.2))")
     +        nccmlev,(ccmlev_zp(k),k=1,nccmlev)
            write(6,"('kmx=',i2,' gcmzp=',/(8f8.2))")
     +        kmx,(gcmzp(k),k=1,kmx)
            stop 'tgcmhist'
          endif 
        elseif (ccmres.eq.'CCM3-T42        ') then
          if (kmx.ne.43) then
            write(6,"(/'>>> ccm model = ',a,' but tgcm kmx=',i2,
     +        ' --> wrong type of tgcm history? (kmx should = 43)')") 
     +        ccmres(1:lenstr(ccmres)),kmx
            write(6,"('nccmlev=',i2,' ccmlev_zp=',/(8f8.2))")
     +        nccmlev,(ccmlev_zp(k),k=1,nccmlev)
            write(6,"('kmx=',i2,' gcmzp=',/(8f8.2))")
     +        kmx,(gcmzp(k),k=1,kmx)
            stop 'tgcmhist'
          endif 
        else
          write(6,"('>>> Unknown ccmres = ',a)") ccmres
          stop 'ccmres'
        endif
c
c Determine indices to coupled fields:
c ifcpl(mxfproc): if ifcpl(i) > 0, then field i is coupled
c   (All ccm fields are coupled fields)
c ixfccm(mxfproc) = indices to ccm fields in fccm (0 if not available)
c tgcm TN = ccm T
c tgcm UN = ccm U
c tgcm VN = ccm V
c tgcm H2O = ccm Q
c tgcm W = ccm OMEGA
c tgcm Z = ccm Z2
c tgcm CH4 = ccm CH4 (tgcm/maccm only)
c
        do ip=1,mxfproc
          do iip=1,mxfccm
            if ((flab8(ip).eq.'TN      '.and.
     +           ccmflab8(iip).eq.'T       ').or. 
     +          (flab8(ip).eq.'UN      '.and.
     +           ccmflab8(iip).eq.'U       ').or. 
     +          (flab8(ip).eq.'VN      '.and.
     +           ccmflab8(iip).eq.'V       ').or. 
     +          (flab8(ip).eq.'H2O     '.and.
     +           ccmflab8(iip).eq.'Q       ').or. 
     +          (flab8(ip).eq.'Z       '.and.
     +           ccmflab8(iip).eq.'Z2      ').or. 
     +          (flab8(ip).eq.'W       '.and.
     +           ccmflab8(iip).eq.'OMEGA   ').or.
     +          (ccmres.eq.'CCM2-MACCM'.and.flab8(ip).eq.'CH4     '.and.
     +           ccmflab8(iip).eq.'CH4     ')) then
              ifcpl(ip) = 1
              ixfccm(ip) = iip
              if (it.eq.1)
     +          write(6,"('Coupled field ip=',i2,' proc field=',a,
     +            ' iip=',i2,' ccm field=',a)") 
     +            ip,flab8(ip),iip,ccmflab8(iip)
            endif
          enddo	! iip=1,mxfccm
        enddo	! ip=1,mxfproc
        ifcpl(ixunvn) = 1
c
c Report field indices:
c
c       write(6,"(/'Field     ip  ifcpl  ifproc  ixfhist  ixfccm')")
c       do ip=1,mxfproc
c         write(6,"(a,i4,i7,i8,i9,i8)") flab8(ip),ip,ifcpl(ip),
c    +      ifproc(ip),ixfhist(ip),ixfccm(ip)
c         if (ip.eq.mxfhist) write(6,"(50('-'),' mxfhist=',i2)")
c    +      mxfhist
c         if (ip.eq.mxfproc) write(6,"(50('-'),' mxfproc=',i2)")
c    +      mxfproc
c       enddo
c       write(6,"(' ')")
c
c Define coupled zp levels:
c
        ncplzp = nccmlev+(kmx-1)
        call alloc(pcplzp,ncplzp)
        do k=1,nccmlev
          cplzp(k) = ccmlev_zp(k)
        enddo
        do k=nccmlev+1,ncplzp
          cplzp(k) = gcmzp(k-nccmlev+1)
        enddo 
        if (it.eq.1) then
          write(6,"(/'Interface between ',a,' and ',a,' is ',
     +      'at zp = ',f7.2)") ccmres(1:lenstr(ccmres)),
     +      model(1:lenstr(model)),cplzp(nccmlev)
          write(6,"('Number of vertical levels=',i2,
     +      ' (nccmlev=',i2,' kmx=',i2,' p0=',1pe13.6,')')") 
     +      ncplzp,nccmlev,kmx,p0
          write(6,"(/'LEVEL  ZP(LN(P0/P))    PRESS(Mb)     AVE HT (KM)',
     +      '   MODEL     MODEL TIME')")
          write(6,"(71('-'))")
        endif
c
c fhist(imx,kmx,jmx,nfhist) was returned by rdtgcm 
c fccm(imx,nccmlev,jmx,nccmfld) was returned by rdlsd
c
        do k=ncplzp,1,-1
          if (k.le.nccmlev) then
            do j=1,jmx
              do i=1,imx
                fglb(i,j) = fccm((ixfccm(ixz)-1)*jmx*nccmlev*imx+
     +            (j-1)*nccmlev*imx+(k-1)*imx+i)
              enddo
            enddo
            aveht = fglbm(fglb,imx,jmx,gcmlat,dlat,dlon,spval)
            modelstr = ccmres
          else
            kk = k-nccmlev+1
            do j=1,jmx
              do i=1,imx
                fglb(i,j) = fhist((ixfhist(ixz)-1)*jmx*kmx*imx+
     +            (j-1)*kmx*imx+(kk-1)*imx+i)
              enddo
            enddo
            aveht = fglbm(fglb,imx,jmx,gcmlat,dlat,dlon,spval)
            modelstr = model
          endif
          pmb = p0*exp(-cplzp(k))*1.e-3
          if (it.eq.1) then
            if (modelstr.eq.model) then
              write(6,"(i3,f11.2,4x,1pe15.6,0pf10.2,8x,a,
     +          ' (',i3,':',i2,':',i2,')')")
     +          k,cplzp(k),pmb,aveht,modelstr(1:lenstr(modelstr)),
     +          (mtimes(i,it),i=1,3)
            else
              write(6,"(i3,f11.2,4x,1pe15.6,0pf10.2,8x,a,
     +          ' (DAY ',f10.5,')')")
     +          k,cplzp(k),pmb,aveht,modelstr(1:lenstr(modelstr)),
     +          ccmday
            endif
          endif
        enddo
c
c Activate workstations for LLU calls:
c
        if (igks_cgm.gt.0) call gacwk(igks_cgm)
        if (igks_ps.gt.0) call gacwk(igks_ps)
c
c Make maps:
c
        if (ipltmaps.gt.0) call mkmaps(mtimes(1,it))
c
c Make latitude slices:
c
        if (ipltlat.gt.0) call mklats(mtimes(1,it))
c
c Make longitude slices:
c
        if (ipltlon.gt.0) then

c         call hpdump(1,"hpdump.dat")
c         len_pgcmzp=ihplen(pgcmzp,ier,1) 
c         len_pgcmpmb=ihplen(pgcmpmb,ier,1) 
c         len_pfhist=ihplen(pfhist,ier,1) 
c         len_pcplzp=ihplen(pcplzp,ier,1) 
c         len_pccmlon=ihplen(pccmlon,ier,1) 
c         len_pccmlat=ihplen(pccmlat,ier,1) 
c         len_pccmlev_mb=ihplen(pccmlev_mb,ier,1) 
c         len_pccmlev_zp=ihplen(pccmlev_zp,ier,1) 
c         len_pfccm=ihplen(pfccm,ier,1) 
c         write(6,"(/'flxproc before mklons: it=',i2,/
c    +      '  len_pgcmzp    =',i6,/,
c    +      '  len_pgcmpmb   =',i6,/,
c    +      '  len_pfhist    =',i6,/,
c    +      '  len_pcplzp    =',i6,/,
c    +      '  len_pccmlon   =',i6,/,
c    +      '  len_pccmlat   =',i6,/,
c    +      '  len_pccmlev_mb=',i6,/,
c    +      '  len_pccmlev_zp=',i6,/,
c    +      '  len_pfccm     =',i6)")
c    +      it,len_pgcmzp,len_pgcmpmb,len_pfhist,len_pcplzp,len_pccmlon,
c    +      len_pccmlat,len_pccmlev_mb,len_pccmlev_zp,len_pfccm

          call mklons(mtimes(1,it))

c         call hpdump(1,"hpdump.dat")
c         len_pgcmzp=ihplen(pgcmzp,ier,1) 
c         len_pgcmpmb=ihplen(pgcmpmb,ier,1) 
c         len_pfhist=ihplen(pfhist,ier,1) 
c         len_pcplzp=ihplen(pcplzp,ier,1) 
c         len_pccmlon=ihplen(pccmlon,ier,1) 
c         len_pccmlat=ihplen(pccmlat,ier,1) 
c         len_pccmlev_mb=ihplen(pccmlev_mb,ier,1) 
c         len_pccmlev_zp=ihplen(pccmlev_zp,ier,1) 
c         len_pfccm=ihplen(pfccm,ier,1) 
c         write(6,"(/'flxproc after mklons: it=',i2,/
c    +      '  len_pgcmzp    =',i6,/,
c    +      '  len_pgcmpmb   =',i6,/,
c    +      '  len_pfhist    =',i6,/,
c    +      '  len_pcplzp    =',i6,/,
c    +      '  len_pccmlon   =',i6,/,
c    +      '  len_pccmlat   =',i6,/,
c    +      '  len_pccmlev_mb=',i6,/,
c    +      '  len_pccmlev_zp=',i6,/,
c    +      '  len_pfccm     =',i6)")
c    +      it,len_pgcmzp,len_pgcmpmb,len_pfhist,len_pcplzp,len_pccmlon,
c    +      len_pccmlat,len_pccmlev_mb,len_pccmlev_zp,len_pfccm

        endif
c
c Load utvert array for current ut:
c
        if (ipltutvert.gt.0) call ldutvert(it)
c
c Deactivate gks workstations:
c
        if (igks_cgm.gt.0) call gdawk(igks_cgm)
        if (igks_ps.gt.0) call gdawk(igks_ps)
c
c Release ccm space:
c
        call hpdeallc(pfccm,ier,ihpabort)
        call hpdeallc(pccmlon,ier,ihpabort)
        call hpdeallc(pccmlat,ier,ihpabort)
        call hpdeallc(pccmlev_mb,ier,ihpabort)
        call hpdeallc(pccmlev_zp,ier,ihpabort)
c
c Release tgcm space:
c
        call hpdeallc(pfhist,ier,ihpabort)
        if (it.lt.ntms) call hpdeallc(pgcmzp,ier,ihpabort)
        call hpdeallc(pgcmpmb,ier,ihpabort)
c
c Release coupled space:
c
        if (it.lt.ntms) call hpdeallc(pcplzp,ier,ihpabort)
c
c End time loop:
c
        iprint = 0	! print flag for rdtgcm, chflds, etc.
        istart = 0
 100    continue
c
c Run check on integrity of memory arena:
c
        call hpcheck(ier)
        if (ier.ne.0) then
          write(6,"('>>> WARNING flxproc: error from hpcheck: ',
     +      'heap has been corrupted: ier=',i3)") ier
        endif
      enddo ! it=1,ntms
c
c Plot ut vs vert:
c
      if (ipltutvert.gt.0) then
        if (igks_cgm.gt.0) call gacwk(igks_cgm)
        if (igks_ps.gt.0) call gacwk(igks_ps)
        call mkutvert
        if (igks_cgm.gt.0) call gdawk(igks_cgm)
        if (igks_ps.gt.0) call gdawk(igks_ps)
        nflds = 0
        do ip=1,mxfproc
          if (ifproc(ip).gt.0) nflds = nflds+1
        enddo
        if (nflds.gt.0) call hpdeallc(pfutvert,ier,ihpabort)
        call hpdeallc(pfutvertz,ier,ihpabort)
      endif
c
c Release space:
c
      call hpdeallc(pgcmzp,ier,ihpabort)
      call hpdeallc(pcplzp,ier,ihpabort)
c
c Retire HLU package:
c
      if (mkplt.gt.0) then
        if (iwk_cgm.gt.0) call NhlFDestroy(iwk_cgm,ier)
        if (iwk_ps.gt.0) call NhlFDestroy(iwk_ps,ier)
        call NhlFDestroy(appid,ier)
        call NhlFClose
      endif
c
c Send cgm file to remote:
c
      write(6,"(' ')")
      if (lenstr(sendcgm).gt.0.and.iwk_cgm.gt.0)
     +  call rcpfile(0,flnm_cgm,sendcgm)
c
c Send ps files to remote:
c
      if (lenstr(sendps).gt.0.and.iwk_ps.gt.0)
     +  call rcpfile(0,flnm_ps,sendps)
c
c Send dat file to remote:
c
      if (lenstr(senddat).gt.0)
     +  call rcpfile(ludat,'flxproc.dat',senddat)
c
c Send xdr file to remote and/or mss:
c
      if (lenstr(sendxdr).gt.0.or.lenstr(sendxdrms).gt.0) then
        iclose = 1
        istat = wrxdrc(trim(flnm_xdr),
     +    dum,dum,dum,dum,dum,idum,idum,dum,dum,dum,dum,idum,iclose)
        if (lenstr(sendxdr).gt.0) call rcpfile(0,flnm_xdr,sendxdr)
        if (lenstr(sendxdrms).gt.0) call sendms(flnm_xdr,sendxdrms)
      endif
c
      stop
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine rptf(fhist,nlon,nzp,nlat,nflds,ixhist,flab,mxhist)
      real fhist(nlon,nzp,nlat,nflds)
      integer ixhist(mxhist)
      character*(*) flab(mxhist)
c
      do ip=1,mxhist
        if (ixhist(ip).gt.0) then
          call fminmax(fhist(1,1,1,ixhist(ip)),nlon*nzp*nlat,
     +      fmin,fmax,1.e36)
          write(6,"('rptf: field ',a,' ip=',i2,' ixhist(ip)=',i2,
     +      ' fmin,max=',2e12.4)") flab(ip)(1:lenstr(flab(ip))),
     +      ip,ixhist(ip),fmin,fmax
        endif
      enddo
      return
      end