c
c------------------------------------------------------------------
c Begin file /home/sting/foster/timesgcm/addch4/addch4.f
c------------------------------------------------------------------
c
      program addch4
c
c Read existing timesgcm history, and write a new one with ch4 and o21d
c inserted after o+ (ch4 from glbmean and o21d = spval 1.e-20)
c
      include "glbm.h"
      parameter(jmx=36,imx=73,kmx=kmxglbm,imxm1=imx-1,imxp1=imx+1)
c
c nfrd = number of fields to read
c nfwr = number of fields to write (nfrd + ch4 and o21d)
c nwlatrd = number of words to read for a latitude slice
c nwlatwr = number of words to write for a latitude slice
c
      parameter(nfrd=47,nwlatrd=1+imxp1*kmx*nfrd)
      parameter(nfwr=50,nwlatwr=1+imxp1*kmx*nfwr)
      dimension frd(nwlatrd),fr(imxp1,kmx,nfrd)
      dimension fwr(nwlatwr),fw(imxp1,kmx,nfwr)
      equivalence(frd(2),fr), (fwr(2),fw)
      dimension gcmlat(jmx),summary(100),iheader(500)
c
      data luhist/8/, luglbm/20/, lunew/9/
c     data md/21/,mh/0/,mm/0/	! model time for TISS04 (solstice)
      data md/26/,mh/0/,mm/0/	! model time for TISE01 (equinox)
      data glbmvol/'/ROBLE/GLBMEANS/GLBMS01 '/ ! global means history
      data igmch4/43/, igmn2/10/
c
      do j=1,jmx
        gcmlat(j) = -87.5+(j-1)*5.
      enddo
c
c Read glbmean history (will be assigned and acquired in getglbm):
c (result is glbmfld(kmxglbm,nglbmfld)
c
      call getglbm(luglbm)
      write(6,"('main: ch4 from glbmean=',/(6e12.4))")
     +  (glbmfld(k,igmch4),k=1,kmxglbm)
c
c Read timegcm history (must already by acquired and assigned in script
c to unit luhist):
c
      it = 0
 100  continue
      read(luhist,end=900) iheader
      it = it+1
      if (iheader(2).eq.md.and.iheader(3).eq.mh.and.iheader(4).eq.mm) 
     +    then
        read(luhist) summary	! summary
        goto 101
      else
c
c Read past unwanted history:
        write(6,"('Searching for history ',i2,':',i2,':',i2,
     +    ' Found history ',i2,':',i2,':',i2)") md,mh,mm,
     +    iheader(2),iheader(3),iheader(4)
        read(luhist) dummy	! summary
        do j=1,jmx
          read(luhist) dummy
        enddo
        goto 100
      endif
c
c Found wanted history -- write header and summary to new file:
 101  continue
      write(6,"('Found history ',i2,':',i2,':',i2,' at it=',i3)")
     +  md,mh,mm,it
      write(lunew) iheader
      write(lunew) summary
c
c Read latitude slices:
      do j=1,jmx
        buffer in(luhist,1)(frd(1),frd(nwlatrd))
        if (unit(luhist)) 10,11,11
 11     write(6,"(' >>> Gettime: io problem buffering in j=',
     +    i3,' lu=',i3,' unit(luhist)=',e12.4)")
     +    j,luhist,unit(luhist)
        stop 'bufin'
 10     continue
c
c Write out some samples:        
c       if (mod(j,12).eq.0) call samples(fr,imxp1,kmx,nfrd,j)
c 
c Insert ch4 and o21d into f array after o+ :
c (ch4 is converted to mass mixing ratios):
c
c First 13 fields are same as read (NT->NOP):
        do ip=1,13
          fw(:,:,ip) = fr(:,:,ip)
        enddo        
c
c Field 14 is ch4 from global mean, converted to mass mixing ratios, 
c and field 43 is "time-1" for ch4 (NCH4NM) (since it is a prognostic):
c (o1 and o2 are fields 17 and 4 in frd, and n2 is field igmn2 in glbmfld)
c
        do k=1,kmx
          fw(:,k,14) = glbmfld(k,igmch4)*16. / (fr(:,k,17)*16. + 
     +                  fr(:,k,4)*32. + glbmfld(k,igmn2)*28.)
          fw(:,k,43) = fw(:,k,14)
        enddo
        if (mod(j,12).eq.0) then
          write(6,"(' ')")
          write(6,"('j=',i3,' ch4 from glbmfld=',/(5e14.6))")
     +      j,(glbmfld(k,igmch4),k=1,kmx)
          write(6,"('j=',i3,' ch4 converted for hist(i=1)=',
     +      /(5e14.6))") j,(fw(1,k,14),k=1,kmx)
        endif
c
c Special value for o21d:
        fw(:,:,15) = 1.e-20
c
c New fields 16-42 are same as fields 14-40 that were read (NPNO2->NOPNM):
c New fields 44-50 are same as fields 41-47 that were read (NMS->NKLDP2):
c (field 43 is ch4 at time-1, i.e., same as field 14 -- see above)
c
        do ip=16,42
          fw(:,:,ip) = fr(:,:,ip-2)
        enddo
        do ip=44,50
          fw(:,:,ip) = fr(:,:,ip-3)
        enddo
c
c Write new latitude slice:
        buffer out(lunew,1)(fwr(1),fwr(nwlatwr))   
      enddo	! j loop
c
      stop 'done'
 900  write(6,"('EOF encountered on unit ',i3)") luhist
      stop 'eof'
      end