!
      module mk_derived
      use proc
      use fields,only: field,mxnaf
      use hist,only: history,hdr
      use get_hist,only: vtgcm_n2
      use den_convert,only: denconv,barm,pkt
      use input,only: ie6300,ie5577,ieo200
      use gfactors,only: exrate
      use mk_drifts
      use epflux
!
! qbary module not available on IBM (because fitpack is not available):
#ifndef AIX
      use qbary
#endif
      implicit none
      contains
!-------------------------------------------------------------------
      subroutine mkderived(f,nf,flat,h,imx,kmx,lat,iden,ionvel)
!
! Calculate needed derived fields in f(nf) from flat:
! (note barm and pkt, used in density conversions, were
!  defined by mkdenparms (in den_convert))
!
! Args:
      type(field),intent(inout) :: f(nf)
      type(history),intent(in) :: h
      integer,intent(in) :: nf,lat,imx,kmx,iden,ionvel
! flat is not actually changed, but must be declared inout so it
! can be passed to denconv.
      real,intent(inout) :: flat(imx,kmx,h%nflds)
!
! Locals:
      integer :: ix,k,nfd,i,ifu,ifv,ier,nna
      integer :: idefu=0,idefv=0
      integer :: isavephi,isaveoh,isaveex,isaveepflux,isaveqbary
      integer ::
     +  ixo2,ixo1,ixt,ixz,ixu,ixv,ixw,ixne,ixte,ixo2p,ixo21d,ixo1d,
     +  ixco2,ixn2,ixno,ixh,ixo3,ixho2,ixoh,ixco,ixh2o,ixch4,ixh2,ixhe
      integer :: iepvy,iepvz,iepvdiv
      real :: fo2(imx,kmx), fo1(imx,kmx), fn2(imx,kmx),hmf2,fof2,
     |  nmf2,fo2_cm3(imx,kmx), fo1_cm3(imx,kmx), fn2_cm3(imx,kmx),
     |  fh(imx,kmx),fhe(imx,kmx),fh2(imx,kmx),fik(imx,kmx)
      real :: fo1d_cm3(imx,kmx)
      real,allocatable :: fwk(:,:,:)
      character(len=16) :: dunit
      real :: fmin,fmax
      real :: dzp,zp,press
      real :: reqlocs_tmp
      real :: slt, sza
      real,parameter :: boltz=1.3807e-16, rmh2o=18.
      real :: pzps(kmx)		! scale heights
      real :: zcm(kmx)          ! geopotential in cm
!
! Externals:
      integer,external :: ixfind,ixfindc
      real,external :: getsza, fslt
      external fminmax
      interface
        function vecsum(u,v,id1,id2,spv)
          implicit none
          integer,intent(in) :: id1,id2
          real,intent(in) :: u(id1,id2),v(id1,id2) 
          real:: vecsum(id1,id2)			! result variable
          real,intent(in),optional :: spv
        end function vecsum
      end interface
!
! Check for derived fields:
!
      nfd = 0
      do ix=1,nf
        if (f(ix)%derived) nfd = nfd+1
      enddo
      if (nfd == 0) return
!
! Get fo1, fo2, fn2 in user requested units according to iden:
! (o2,o,tn,z should always be on the histories)
! (denconv calls below do not change flat)
!
      ixt  = ixfindc(h%fnames,h%nflds,'TN      ')
      ixz  = ixfindc(h%fnames,h%nflds,'Z       ')
      ixo1 = ixfindc(h%fnames,h%nflds,'O1      ')
      ixo2 = ixfindc(h%fnames,h%nflds,'O2      ')
!
! Fatal if T or Z are missing:
      if (ixz==0.or.ixt==0) then
        write(6,"('>>> mkderived: need tn, and z on history:',
     +    ' ixt=',i2,' ixz=',i2)") ixt,ixz
        stop 'mkderived'
      endif
!
! o2,o1 are missing: this is fatal if not a jtgcm history:
      if (ixo2==0.or.ixo1==0) then
        if (.not.h%isjtgcm) then
          write(6,"('>>> mkderived: need o2 and o1 on history:',
     +      ' ixo2=',i2,' ixo1=',i2)") ixo2,ixo1
          stop 'mkderived'
        else ! is jtgcm
!         write(6,"('>>> WARNING mkderived: o2 and/or o1 are missing:',
!    |      ' ixo2=',i3,' ixo1=',i3)") ixo2,ixo1
        endif
      endif
!
! N2 is on mtgcm histories, vtgcm N2 mmr has been read into vtgcm_n2,
!   otherwise (tgcm) n2=1-o2-o:
!
      if (h%ismtgcm) then
        ixn2 = ixfindc(h%fnames,h%nflds,'N2      ')
        fn2(:,:) = flat(:,:,ixn2)
      elseif (h%isvtgcm) then
        if (h%isnew) then
          ixn2 = ixfindc(h%fnames,h%nflds,'N2      ')
          fn2(:,:) = flat(:,:,ixn2)
        else
          fn2(:,:) = vtgcm_n2(1:imx,:,lat)
        endif
      else
        fn2(:,:) = max(.00001,(1.-flat(:,:,ixo2)-flat(:,:,ixo1)))
      endif
!
! Conversion to user requested units (do not alter flat, do alter fn2):
      if (.not.h%isjtgcm) then
        if (iden > 0) then
          call denconv(flat(:,:,ixo1),fo1,1,barm,pkt,16.,
     +      imx,kmx,iden,'MMR',' ',dunit)
          call denconv(flat(:,:,ixo2),fo2,1,barm,pkt,32.,
     +      imx,kmx,iden,'MMR',' ',dunit)
          call denconv(fn2,fn2,0,barm,pkt,28.,
     +      imx,kmx,iden,'MMR',' ',dunit)
        else
          fo2(:,:) = flat(:,:,ixo2)
          fo1(:,:) = flat(:,:,ixo1)
        endif
      endif ! not isjtgcm
!
! For jtgcm:
      ixh  = ixfindc(h%fnames,h%nflds,'H       ')
      ixhe = ixfindc(h%fnames,h%nflds,'HE      ')
      if (h%isjtgcm) then
        if (ixh==0.or.ixhe==0) then
          write(6,"('>>> mkderived: need H and HE in flat for jtgcm:',
     |      ' ixh=',i3,' ixhe=',i3)") ixh,ixhe
          stop 'mkderived'
        endif
        fh2(:,:) = max(.00001,(1.-flat(:,:,ixhe)-flat(:,:,ixh)))
        if (iden > 0) then
          call denconv(flat(:,:,ixh),fh,1,barm,pkt,1.,
     |      imx,kmx,iden,'MMR',' ',dunit)
          call denconv(flat(:,:,ixhe),fhe,1,barm,pkt,4.,
     |      imx,kmx,iden,'MMR',' ',dunit)
          call denconv(fh2,fh2,0,barm,pkt,2.2,
     |      imx,kmx,iden,'MMR',' ',dunit)
        else
          fhe(:,:) = flat(:,:,ixhe)
          fh(:,:) = flat(:,:,ixh)
        endif
      endif
!
! fo1_cm3, fo2_cm3, and fn2_cm3 are o2,o,n2 in cm3:
!
      if (.not.h%isjtgcm) then
        call denconv(flat(:,:,ixo2),fo2_cm3(:,:),1,
     +    barm,pkt,32.,imx,kmx,-1,'MMR','CM3',dunit)
        call denconv(flat(:,:,ixo1),fo1_cm3(:,:),1,
     +    barm,pkt,16.,imx,kmx,-1,'MMR','CM3',dunit)
        if (h%ismtgcm) then
          fn2_cm3 = flat(:,:,ixn2)
        elseif (h%isvtgcm) then
          if (h%isnew) then
            fn2_cm3 = flat(:,:,ixn2)
          else
            fn2_cm3(:,:) = vtgcm_n2(1:imx,:,lat)
          endif
        else
          fn2_cm3 = max(.00001,(1.-flat(:,:,ixo2)-flat(:,:,ixo1)))
        endif
        call denconv(fn2_cm3,fn2_cm3,0,barm,pkt,28.,imx,kmx,-1,
     +    'MMR','CM3',dunit)
      endif ! not isjtgcm
!
      isavephi = 0
      isaveoh  = 0
      isaveex  = 0
      isaveepflux = 0
      isaveqbary = 0
      mkepv = .false.	! in module epflux
#ifndef AIX
      mkqbary = .false.	! in module qbary
#endif
      nna = 0		! number of requested na fields
!
! Fields loop:
! (Process only derived fields, requested or not requested)
!
      fields_loop: do ix=1,nf
        if (.not.f(ix)%derived.or..not.associated(f(ix)%data))
     +    cycle fields_loop
!
! bf 7/14/03: n2 is already available in desired units in fn2 from above:
! (denconv call in previous version did not work on IBM because isetout==0)
!
        if (f(ix)%fname8 == 'N2      ') then
          f(ix)%data(:,lat,:) = fn2(:,:)
          allocate(f(ix)%lev(h%nzp),stat=ier)
          if (ier /= 0) 
     +      call allocerr(ier,"mkderived allocating lev for N2")
          f(ix)%lev = h%lev ! midpoints
!
! 3/15/06 btf: add WN as derived field calculated from OMEGA (for new histories):
! (see also known field 'W' from old histories, proclat_old.F)
!
        elseif (f(ix)%fname8 == 'WN      ') then
          ixw = ixfindc(h%fnames,h%nflds,'WN      ')
          if (ixw > 0) then
            write(6,"('Note mkderived: found WN on the history -- ',
     |        'NOT calculating WN as derived field')")
            cycle fields_loop
          endif
          ixw = ixfindc(h%fnames,h%nflds,'OMEGA   ')
          if (ixw==0) then
            write(6,"('>>> WARNING mkderived: need OMEGA to make WN')")
            cycle fields_loop
          endif
          dlev = (h%zpit-h%zpib)/float(h%nzp-1)
          do i=1,imx
            zcm(:) = flat(i,:,ixz)*1.e5
            do k=2,kmx-1
              pzps(k) = zcm(k+1)-zcm(k-1)
            enddo
            pzps(1) = 4.*zcm(2)-3.*zcm(1)-zcm(3)
            pzps(kmx) = 3.*zcm(kmx)-4.*zcm(kmx-1)+zcm(kmx-2)
!
! 3/21/06 btf: make pzps generic for 0.5 and 0.25 vertical resolutions:
            pzps(:) = pzps(:)/(2.*dlev)
            f(ix)%data(i,lat,:) = flat(i,:,ixw)*pzps(:)
          enddo
          f(ix)%data(imx,lat,:) = f(ix)%data(1,lat,:)     ! force periodic point
          f(ix)%data(:,lat,:) = f(ix)%data(:,lat,:) * .01 ! cm/s to m/s
!
! Terrestrial RHO = o2+o1+n2 (where o2,o1,n2 are in desired units):
!
        elseif (f(ix)%fname8 == 'RHO     ') then
          if (.not.h%ismtgcm.and..not.h%isjtgcm.and.
     |        trim(h%version)/='VTGCM') then
            f(ix)%data(:,lat,:) = fo2(:,:)+fo1(:,:)+fn2(:,:)
            select case (iden)
              case(0)
                f(ix)%units = 'MMR'
              case(1)
                f(ix)%units = 'CM3'
              case(2)
                f(ix)%units = 'CM3-MR'
              case(3)
                f(ix)%units = 'GM/CM3'
              case(4)
                f(ix)%units = 'PPBV'
              case default
                write(6,"('>>> WARNING mkderived RHO: unknown iden=',
     |            i4)") iden
            end select
!
! For jtgcm, rho = h+he+h2, where h2=1-h-h3. fh,fe,fh2 have been 
! converted to desired units above according to user flag iden.
!
          elseif (h%isjtgcm) then
            f(ix)%data(:,lat,:) = fh(:,:)+fhe(:,:)+fh2(:,:)
!
! For [mv]tgcm, rho = co2+n2+co+o
! (internal mkco2 defines co2 in f(ix)%data, w/ density conversion
!  done only if 2nd arg is > 0)
!
          else ! mtgcm or vtgcm
! Store co2 in f(ix)%data, in desired units:
            call mkco2(f(ix)%data(:,lat,:),1)
! Store co in fwk in desired units:
            if (allocated(fwk)) deallocate(fwk)
            allocate(fwk(imx,kmx,1),stat=ier)
            if (ier /= 0) 
     +        call allocerr(ier,"mkderived allocating fwk for co")
            ixco= ixfindc(h%fnames,h%nflds,'CO      ')
            call denconv(flat(:,:,ixco),fwk(:,:,1),1,
     +        barm,pkt,28.,imx,kmx,iden,'MMR',' ',f(ix)%units)
! Now rho = co2+n2+co+o (fn2 and fo1 already in desired units)
            f(ix)%data(:,lat,:) = f(ix)%data(:,lat,:)+
     +        fn2(:,:)+fwk(:,:,1)+fo1(:,:)
            deallocate(fwk)
! rho in kg/km3 for [mv]tgcm if mass density requested:
            if (iden==3) then
              f(ix)%data(:,lat,:) = f(ix)%data(:,lat,:)*1.e12 
              f(ix)%units = 'KG/KM3          '
            endif
          endif
!
! Save 3d parameters needed for ion velocities 
! (the drifts themselves are calculated by mkdrifts in a later 
!  lat loop in getflds)
!
        elseif (f(ix)%fname8 == 'UI      ' .or.
     +          f(ix)%fname8 == 'VI      ' .or.
     +          f(ix)%fname8 == 'WI      ' .or.
     +          f(ix)%fname8 == 'UIVI    ') then
          if (isavephi==0.and.ionvel>0) then
            call savephi(flat,h,imx,kmx,nlat,lat,ionvel)
            isavephi = 1
          endif
!
! UNVN: store vector sum of un+vn in the UNVN field slot:
!
        elseif (f(ix)%fname8 == 'UNVN    ') then
          ixu = ixfindc(h%fnames,h%nflds,'UN      ')
          ixv = ixfindc(h%fnames,h%nflds,'VN      ')
          if (ixu==0 .or. ixv==0) then
            write(6,"('>>> WARNING mkderived: need UN and VN to ',
     +        'make UNVN')")
            cycle fields_loop
          endif
          f(ix)%data(:,lat,:) = vecsum(flat(:,:,ixu),flat(:,:,ixv),
     +      size(flat,1),size(flat,2))
!
! UNVN:
! Define f(ixu)%data and f(ixv)%data if necessary, to be used by mkmaps 
!   when drawing UN+VN vector arrows. (If user also requested UN and
!   VN separately, they were allocated by allocfdat and will be defined 
!   by mkunderived, but in any case UN and VN should have been added 
!   to f(:) by allocf for UNVN)
!
          ifu = ixfindc(f%fname8,nf,'UN      ')
          if (ifu==0) then
            write(6,"('>>> mkderived: need UN for UNVN')")
            stop 'UN'
          endif
          ifv = ixfindc(f%fname8,nf,'VN      ')
          if (ifv==0) then
            write(6,"('>>> mkderived: need VN for UNVN')")
            stop 'VN'
          endif
          if (.not.associated(f(ifu)%data)) then ! UN was not requested
            allocate(f(ifu)%data(h%nlon-1,h%nlat,f(ifu)%nlev),
     +        stat=ier)
            if (ier /= 0) 
     +        call allocerr(ier,"mkderived allocating f(ifu)%data") 
            idefu = 1
          endif
          if (.not.associated(f(ifv)%data)) then	! VN was not requested
            allocate(f(ifv)%data(h%nlon-1,h%nlat,f(ifv)%nlev),
     +        stat=ier)
            if (ier /= 0) 
     +        call allocerr(ier,"mkderived allocating f(ifv)%data") 
            idefv = 1
          endif
          if (idefu > 0) f(ifu)%data(:,lat,:) = flat(:,:,ixu)
          if (idefv > 0) f(ifv)%data(:,lat,:) = flat(:,:,ixv)
!
! FOF2 (height-independent):
!
        elseif (f(ix)%fname8 == 'FOF2    ') then
          ixne = ixfindc(h%fnames,h%nflds,'NE      ')
          if (ixne==0 .or. ixz==0) then
            write(6,"('WARNING mkderived: need NE and Z to ',
     +        'make FOF2')")
            cycle fields_loop
          endif
          do i=1,imx
            call fof2int(flat(i,:,ixne),flat(i,:,ixz),kmx,hmf2,
     +        f(ix)%data(i,lat,1),nmf2,0,i,lat)
          enddo
!
! HMF2 (height-independent):
!
        elseif (f(ix)%fname8 == 'HMF2    ') then
          ixne = ixfindc(h%fnames,h%nflds,'NE      ')
          if (ixne==0 .or. ixz==0) then
            write(6,"('WARNING mkderived: need NE and Z to ',
     +        'make HMF2')")
            cycle fields_loop
          endif
          do i=1,imx
            call fof2int(flat(i,:,ixne),flat(i,:,ixz),kmx,
     +        f(ix)%data(i,lat,1),fof2,nmf2,0,i,lat)
          enddo
!
! NMF2 (NE at hmf2 heights):
!
        elseif (f(ix)%fname8 == 'NMF2    ') then
          ixne = ixfindc(h%fnames,h%nflds,'NE      ')
          if (ixne==0 .or. ixz==0) then
            write(6,"('WARNING mkderived: need NE and Z to ',
     +        'make NMF2')")
            cycle fields_loop
          endif
          do i=1,imx
            call fof2int(flat(i,:,ixne),flat(i,:,ixz),kmx,
     +        hmf2,fof2,f(ix)%data(i,lat,1),0,i,lat)
          enddo
!--------------------- begin emissions fields -----------------------
!
! E6300 emission (dependencies are t,o2,o,te,o2+,ne in cm3)
! (if ie6300 > 0, solred is called by mke6300)
!     subroutine mke6300(tn,xo2,xo,xn2,ht,te,xne,xo2p,id,fout,ie6300,
!    +  iyd,ut,glat,glon,f107d,f107a)
!
        elseif (f(ix)%fname8 == 'E6300   ') then
          ixte = ixfindc(h%fnames,h%nflds,'TE      ')
          ixo2p= ixfindc(h%fnames,h%nflds,'O2+     ')
          if (ixo2p <= 0) then
            ixo2p = ixfindc(h%fnames,h%nflds,'O2P     ')
            if (ixo2p <= 0) then
              write(6,"('>>> WARNING: Need either O+ or O2+ for',
     |          ' E6300')")
            else
!             write(6,"('Note: using O2P (alias for O2+) to calculate',
!    |          ' E6300: ixo2p=',i3)") ixo2p
            endif
          endif
          ixne = ixfindc(h%fnames,h%nflds,'NE      ')
          do i=1,imx
            call mke6300(flat(i,:,ixt),fo2_cm3(i,:),
     +        fo1_cm3(i,:),fn2_cm3(i,:),flat(i,:,ixz),
     +        flat(i,:,ixte),flat(i,:,ixne),flat(i,:,ixo2p),kmx,
     +        f(ix)%data(i,lat,:),ie6300,h%iyd,h%ut,gcmlat(lat),
     +        gcmlon(i),h%f107d,h%f107a)
          enddo        
!
! E5577 greenline emission (dependencies are t,o2,o1,te,ne,o2p,o21d):
! (components are specified by ie5577(5))
!     subroutine mke5577(tn,xo2,xo,xn2,ht,te,xne,xo2p,xo21d,
!    +  id,fout,ie5577,iyd,ut,glat,glon,f107d,f107a)
!
        elseif (f(ix)%fname8 == 'E5577   ') then
          ixte = ixfindc(h%fnames,h%nflds,'TE      ')
          if (ixte <= 0) then
            if (ie5577(2) > 0) then
              write(6,"('>>> mkderived: need TE for E5577 when ',
     |          ' ie5577(2) > 0')")
              call shutdown('need TE')
            else
              ixte = 1 ! dummy value, will not be used in mke5577
            endif
          endif
          ixne = ixfindc(h%fnames,h%nflds,'NE      ')
          if (ixne <= 0) then
            write(6,"('>>> mkderived: need NE for E5577')")
            call shutdown('need NE')
          endif
          ixo2p= ixfindc(h%fnames,h%nflds,'O2+     ')
          if (ixo2p <= 0) then
            ixo2p = ixfindc(h%fnames,h%nflds,'O2P     ')
            if (ixo2p <= 0.and.ie5577(2) > 0) then
              write(6,"('>>> mkderived: Need O2+ for E5577',
     |          ' when ie5577(2) > 0')")
              call shutdown('need O2+')
            elseif (ixo2p <= 0) then
              ixo2p = 1 ! dummy value, will not be used in mke5577
            endif
          endif
          ixo21d= ixfindc(h%fnames,h%nflds,'O21D    ')
          if (ixo21d > 0) then
            fik(:,:) = flat(:,:,ixo21d) 
          else
!
! If o21d is not available for E5577, use 1.e-20:
            if (lat == 1)
     |      write(6,"('Note: O21D not available: using O21D=1.e-20', 
     |        ' to calculate E5577')")
            fik(:,:) = 1.e-20
          endif

!         write(6,"('Calling mke5577: f107d,a=',2e12.4)") h%f107d,
!    |      h%f107a
!         write(6,"('Calling mke5577: ixt=',i3,' ixz=',i3,' ixte=',
!    |      i3,' ixne=',i3,' ixo2p=',i3)") ixt,ixz,ixte,ixne,ixo2p

          do i=1,imx
            call mke5577(flat(i,:,ixt),fo2_cm3(i,:),
     +        fo1_cm3(i,:),fn2_cm3(i,:),flat(i,:,ixz),
     +        flat(i,:,ixte),flat(i,:,ixne),flat(i,:,ixo2p),
     +        fik(i,:),kmx,f(ix)%data(i,lat,:),
     +        ie5577,h%iyd,h%ut,gcmlat(lat),gcmlon(i),h%f107d,
     +        h%f107a)
          enddo
!
! EO200 O2 0-0 band emission (dependencies are t,o2,o1):
!     subroutine mkeo200(tn,xo2,xo,xn2,id,fout)
!
! fo1d_cm3 is o1d in cm3:
!
	elseif (f(ix)%fname8 == 'EO200   ') then
          if (ieo200(3) /= 0) then
            ixo1d= ixfindc(h%fnames,h%nflds,'O1D     ')
            call denconv(flat(:,:,ixo1d),fo1d_cm3(:,:),1,
     +      barm,pkt,32.,imx,kmx,-1,'MMR','CM3',dunit)
!           flat_o1d = flat(:,:,ixo1d)
          endif
          do i=1,imx
            slt = fslt(slt,h%ut,gcmlon(i),1)        ! curr local time
            sza = getsza(h%mtime(1),slt,lat,gcmlon(i))    ! curr zen ang
!            print *, ' MKDERIVED: slt, sza = ', slt, sza
            call mkeo200(flat(i,:,ixt),h,fo2_cm3(i,:),
     +        fo1_cm3(i,:),fn2_cm3(i,:),fo1d_cm3(i,:),kmx,sza,
     +        ieo200,f(ix)%data(i,lat,:))
          enddo
!
! EOH83 OH 8-3 band emission (dependencies are t,o2,o1):
!     subroutine mkeoh83(tn,xo2,xo,xn2,id,fout)	! without ohrad
!
        elseif (f(ix)%fname8 == 'EOH83   ') then
          do i=1,imx
            call mkeoh83(flat(i,:,ixt),fo2_cm3(i,:),
     +        fo1_cm3(i,:),fn2_cm3(i,:),kmx,f(ix)%data(i,lat,:))
          enddo
!
! ECO215u CO2 15u emission (dependencies are t,o1,co2):
!     subroutine mkeco215(tn,xo,xco2,id,fout)
!
        elseif (f(ix)%fname8 == 'ECO215u ') then
          if (allocated(fwk)) deallocate(fwk)
          allocate(fwk(imx,kmx,1),stat=ier)
          if (ier /= 0) 
     +      call allocerr(ier,"mkderived allocating fwk for co2")
          ixco2= ixfindc(h%fnames,h%nflds,'CO2     ')
          call denconv(flat(:,:,ixco2),fwk(:,:,1),1,
     +      barm,pkt,44.,imx,kmx,-1,'MMR','CM3',dunit)
          do i=1,imx
            call mkeco215(flat(i,:,ixt),fo1_cm3(i,:),
     +        fwk(i,:,1),kmx,f(ix)%data(i,lat,:))
          enddo
          deallocate(fwk)
!
! ENO53u NO 53u emission (dependencies are t,o1,no):
!     subroutine mkeno53(tn,xo,xno,id,fout)
!
        elseif (f(ix)%fname8 == 'ENO53u  ') then
          if (allocated(fwk)) deallocate(fwk)
          allocate(fwk(imx,kmx,1),stat=ier)
          if (ier /= 0) 
     +      call allocerr(ier,"mkderived allocating fwk for no")
          ixno= ixfindc(h%fnames,h%nflds,'NO      ')
          call denconv(flat(:,:,ixno),fwk(:,:,1),1,
     +      barm,pkt,30.,imx,kmx,-1,'MMR','CM3',dunit)
          do i=1,imx
            call mkeno53(flat(i,:,ixt),fo1_cm3(i,:),
     +        fwk(i,:,1),kmx,f(ix)%data(i,lat,:))
          enddo
          deallocate(fwk)
!
! Strickland's excitation rates:
!
        elseif (trim(f(ix)%type) == 'EXCITED-STATE') then
          if (isaveex == 0) then
            if (allocated(fwk)) deallocate(fwk)
            allocate(fwk(imx,kmx,1),stat=ier)
            if (ier /= 0) 
     +        call allocerr(ier,"mkderived allocating fwk for no")
            fwk(:,:,1) = fo2_cm3(:,:)+fo1_cm3(:,:)+fn2_cm3(:,:)
            fields_loop1: do i=1,nf
              if (.not.f(i)%derived.or..not.associated(f(i)%data).or.
     +          trim(f(i)%type)/='EXCITED-STATE') cycle fields_loop1
              call exrate(f(i)%fname8,f(i)%data(:,lat,:),fwk(:,:,1),
     +          fn2_cm3,fo1_cm3,flat(:,:,ixz),imx,kmx,gcmlat(lat),h,hdr)
            enddo fields_loop1 ! (excitation states)
            isaveex = 1
          endif
!
! OH fields:
!
        elseif (trim(f(ix)%type) == 'OH-VIB' .or.
     +          trim(f(ix)%type) == 'OH-BAND') then
          if (isaveoh == 0) then
!
! Getoh needs o2,o,n2,h,o3,ho2,oh in cm3 (o2,o,n2 are already available
! as fo2_cm3, fo1_cm3, and fn2_cm3 from above):
!
            if (allocated(fwk)) deallocate(fwk)
            allocate(fwk(imx,kmx,4),stat=ier)	! h,o3,ho2,oh in cm3
            if (ier /= 0) 
     +        call allocerr(ier,"mkderived allocating fwk for oh")
! Store H cm3 in fwk(:,:,1)
            ixh= ixfindc(h%fnames,h%nflds,'H       ')
            if (ixh==0) write(6,"('>>> WARNING mkderived: need ',
     +        'H for OH')")
            call denconv(flat(:,:,ixh),fwk(:,:,1),1,
     +        barm,pkt,1.,imx,kmx,-1,'MMR','CM3',dunit)
! Store O3 cm3 in fwk(:,:,2)
            ixo3= ixfindc(h%fnames,h%nflds,'O3      ')
            if (ixo3==0) write(6,"('>>> WARNING mkderived: need ',
     +        'O3 for OH')")
            call denconv(flat(:,:,ixo3),fwk(:,:,2),1,
     +        barm,pkt,48.,imx,kmx,-1,'MMR','CM3',dunit)
! Store OH cm3 in fwk(:,:,4)
            ixoh= ixfindc(h%fnames,h%nflds,'OH      ')
            if (ixoh==0) write(6,"('>>> WARNING mkderived: need ',
     +        'OH for OH-V/B')")
            call denconv(flat(:,:,ixoh),fwk(:,:,4),1,
     +        barm,pkt,17.,imx,kmx,-1,'MMR','CM3',dunit)
!
! Store HO2 cm3 in fwk(:,:,3)
! 11/97 If ho2 is unavailable, use mkho2 to calculate ho2 from
!       t,o2,o,h,o3,oh,no (see allocfdat in getflds.f)
!     subroutine mkho2(tn,fo2,fo1,fn2,fh,fo3,foh,fno,fho2_out,kmx)
!
            ixho2= ixfindc(h%fnames,h%nflds,'HO2     ')
            if (ixho2==0) then
              ixno= ixfindc(h%fnames,h%nflds,'NO      ')
              ixo3= ixfindc(h%fnames,h%nflds,'O3      ')
              if (ixno==0.or.ixo3==0) then
                write(6,"('>>> WARNING mkderived: need HO2, or ',
     +            'NO and O3 to calculate OH-V/B')")
                stop 'OH'
              endif
! mkho2 returns column of ho2 in fwk(i,:,3)
              do i=1,imx
                call mkho2(flat(i,:,ixt),fo2_cm3(i,:),fo1_cm3(i,:),
     +            fn2_cm3(i,:),flat(i,:,ixh),flat(i,:,ixo3),
     +            flat(i,:,ixoh),flat(i,:,ixno),fwk(i,:,3),kmx)
              enddo
            else
              call denconv(flat(:,:,ixho2),fwk(:,:,3),1,
     +          barm,pkt,33.,imx,kmx,-1,'MMR','CM3',dunit)
            endif
!
! If iohglb>0, then calculate oh fields at global grid, otherwise,
!   calculate oh fields only at reqlocs locations. Note iohglb has
!   been set to 1 if either of the following are true:
! 1. if any requested locations are zonal means, global means, or local time
! 2. if ipltmaps>0, ipltlon>0, or ipltlat>0
!
            if (iohglb > 0) then
              call getoh(f,nf,flat,fo2_cm3,fo1_cm3,fn2_cm3,
     +          fwk(:,:,1),fwk(:,:,2),fwk(:,:,3),fwk(:,:,4),
     +          h,imx,kmx,nlat,lat,0,0)
            else	! at reqlocs locations only (pass -lat to getoh):
              do i=1,mxloc
                if (reqlocs(1,i)/=spval.and.reqlocs(2,i)/=spval) then
                  if (gcmlat(lat)==reqlocs(1,i)) then
		    reqlocs_tmp = reqlocs(2,i)    ! switch to kind=4 if default
                    call getoh(f,nf,flat,fo2_cm3,fo1_cm3,fn2_cm3,
     +                fwk(:,:,1),fwk(:,:,2),fwk(:,:,3),fwk(:,:,4),
     +                h,imx,kmx,nlat,-lat,lat,
     +                ixfind(gcmlon,nlon,reqlocs_tmp,dlon))
                  endif                  
                endif
              enddo
            endif
            isaveoh = 1
            deallocate(fwk)
          endif	! isaveoh==0
!
! Sodium fields:
!   Accumulate number of requested Na sp in nna, then call getna after
!   fields loop (see below):
!
        elseif (f(ix)%fname8(1:2) == 'Na') then
          if (nna+1 > mxnaf) then
            write(6,"('>>> WARNING mkderived: too many Na species ',
     +        'requested: mxnaf=',i2)") mxnaf
            f(ix)%data = 0.
          endif
          nna = nna+1
!
! CO2 is derived only if history is mtgcm or vtgcm (co2=1-o1-co-n2):
!   (see setmtgcm and setvtgcm in getflds.f)
! For mtgcm, n2 is on the histories. For vtgcm, n2 is in vtgcm_n2(:,:,lat))
! (mkco2 is internal, see below)
!
        elseif (f(ix)%fname8 == 'CO2     ') then
          call mkco2(f(ix)%data(:,lat,:),1)
!
! O/CO2 ratio 
        elseif (f(ix)%fname8 == 'O/CO2   ') then
          if (h%ismtgcm.or.trim(h%version)=='VTGCM') then
!
! For mtgcm or vtgcm histories: 
!   Get co2 in mass mixing ratio (mmr), and store in fwk(:,:,1).
!   (see comments re co2 above).
!
            if (allocated(fwk)) deallocate(fwk)
            allocate(fwk(imx,kmx,1),stat=ier)
            if (ier /= 0) 
     +        call allocerr(ier,"mkderived allocating fwk for co2")
            call mkco2(fwk(:,:,1),1)
!
! Form the o/co2 ratio (fo1 already in desired units):
            f(ix)%data(:,lat,:) = fo1(:,:)/fwk(:,:,1)
            deallocate(fwk)
!
! For non-mtgcm history, convert co2, then take fo1/co2 since fo1 is
! already converted:  
          else	! non-[mv]tgcm history
            ixco2= ixfindc(h%fnames,h%nflds,'CO2     ')
            if (ixco2==0) then
              write(6,"('>>> mkderived: need CO2 for O/CO2')")
              stop 'O/CO2'
            endif
            f(ix)%data(:,lat,:) = flat(:,:,ixco2)
!
! bf 7/14/03: must call denconv with isetout=1 for IBM's, when
! first arg is a pointer.
!
            call denconv(f(ix)%data(:,lat,:),f(ix)%data(:,lat,:),1,
     +        barm,pkt,44.,imx,kmx,iden,'MMR',' ',f(ix)%units)
            f(ix)%data(:,lat,:) = fo1(:,:)/f(ix)%data(:,lat,:)
          endif
!
! O/N2 ratio: 
        elseif (f(ix)%fname8 == 'O/N2    ') then
          f(ix)%data(:,lat,:) = fo1(:,:)/fn2(:,:)
!
! N2/O ratio: 
        elseif (f(ix)%fname8 == 'N2/O    ') then
          f(ix)%data(:,lat,:) = fn2(:,:)/fo1(:,:)
!
! O/(O2+N2) ratio: 
        elseif (f(ix)%fname8 == 'O/O2+N2 ') then
          f(ix)%data(:,lat,:) = fo1(:,:)/(fo2(:,:)+fn2(:,:))
!
! O/O2 ratio: 
        elseif (f(ix)%fname8 == 'O/O2    ') then
          f(ix)%data(:,lat,:) = fo1(:,:)/fo2(:,:)
!
! Eliassen-palm fluxes:
        elseif (trim(f(ix)%type)=='EPFLUX') then
          if (isaveepflux==0) then
            call mkrhokg ! get rho in kg/m3 into fwk(:,:,1)
!
! Save necessary fields for ep fluxes:
            call save_epv(flat,fwk(:,:,1),h,imx,kmx,nlat,gcmlat(lat),
     +        lat)
            deallocate(fwk)
            isaveepflux = 1
            mkepv = .true.
          endif
!
! Qbary:
! Requested fields may be QBARY, QH (horiz component), and/or QV (vertical)
!   (QH and QV have QBARY as a dependency).
! Save_qbary saves zonal means tn, un, and pottn for later calculation
!   of the QBARY fields themselves by calc_qbary (called from 
!   getflds if mkqbary is true). See qbary.f.
!
        elseif (trim(f(ix)%type)=='QBARY') then
#if defined(AIX) || defined(SUN) || defined(LINUX)
          write(6,"('>>> QBARY not available on IBM or SUN systems.')")
#else
          if (isaveqbary==0) then
            call mkrhokg ! get rho in kg/m3 into fwk(:,:,1)
            call save_qbary(flat,fwk(:,:,1),h,imx,kmx,nlat,gcmlat(lat),
     +        lat)
            deallocate(fwk)
            mkqbary = .true.
            isaveqbary = 1
          endif
#endif
!
! Total hydrogen: 2*h2o+h+oh+ho2+4*ch4+2*h2 (vol mix ratio)
!
        elseif (f(ix)%fname8 == 'HTOT    ') then
          if (allocated(fwk)) deallocate(fwk)
          allocate(fwk(imx,kmx,6),stat=ier)
          if (ier /= 0) 
     +      call allocerr(ier,"mkderived allocating 6 fwk for htot")
!
! Get indices in flat for htot dependencies: 
          ixh2o = ixfindc(h%fnames,h%nflds,'H2O     ')
          ixh =   ixfindc(h%fnames,h%nflds,'H       ')
          ixoh =  ixfindc(h%fnames,h%nflds,'OH      ')
          ixch4 = ixfindc(h%fnames,h%nflds,'CH4     ')
          ixh2  = ixfindc(h%fnames,h%nflds,'H2      ')
          if (ixh2o<=0.or.ixh<=0.or.ixoh<=0.or.ixch4<=0.or.ixh2<=0) then
            write(6,"('mkderived WARNING: missing dependencies',
     +        ' for HTOT: ixh2o,h,oh,ho2,ch4,h2=',6i3)") 
     +        ixh2o,ixh,ixoh,ixho2,ixch4,ixh2
          endif
!
! If ho2 is not available, use routine mkho2, which requires NO and O3.
! If ho2 is not available and either NO or O3 are not available, htot 
! is zeroed out:
          ixho2 = ixfindc(h%fnames,h%nflds,'HO2     ')
          if (ixho2 <= 0 ) then
            ixno  = ixfindc(h%fnames,h%nflds,'NO      ')
            ixo3  = ixfindc(h%fnames,h%nflds,'O3      ')
            if (ixno <= 0.or.ixo3 <= 0) then
              write(6,"('>>> WARNING: HO2 or NO and O3 needed for ',
     +          'HTOT: ixho2,no,o3=',3i3)") ixho2,ixno,ixo3
              write(6,"('    HTOT will be zero.')")
              f(ix)%data(:,lat,:) = 0.
              cycle fields_loop
            endif
            write(6,"('NOTE: ho2 is needed for HTOT, but ',
     +        'is unavailable on current history. Using mkho2..')")
            do i=1,imx
              call mkho2(flat(i,:,ixt),fo2_cm3(i,:),fo1_cm3(i,:),
     +          fn2_cm3(i,:),flat(i,:,ixh),flat(i,:,ixo3),
     +          flat(i,:,ixoh),flat(i,:,ixno),fwk(i,:,4),kmx)
            enddo
          else 
            call denconv(flat(:,:,ixho2),fwk(:,:,4),1,
     +        barm,pkt,33.,imx,kmx,-1,'MMR','CM3',dunit)
          endif
!
! Components must be in cm3 regardless of iden:
          call denconv(flat(:,:,ixh2o),fwk(:,:,1),1,
     +      barm,pkt,18.,imx,kmx,-1,'MMR','CM3',dunit)
          call denconv(flat(:,:,ixh)  ,fwk(:,:,2),1,
     +      barm,pkt,1. ,imx,kmx,-1,'MMR','CM3',dunit)
          call denconv(flat(:,:,ixoh) ,fwk(:,:,3),1,
     +      barm,pkt,17.,imx,kmx,-1,'MMR','CM3',dunit)
          call denconv(flat(:,:,ixch4),fwk(:,:,5),1,
     +      barm,pkt,16.,imx,kmx,-1,'MMR','CM3',dunit)
          call denconv(flat(:,:,ixh2) ,fwk(:,:,6),1,
     +      barm,pkt,2. ,imx,kmx,-1,'MMR','CM3',dunit)
!
! Calculate total hydrogen (volume mixing ratio):
! Total hydrogen: (2*h2o+h+oh+ho2+4*ch4+2*h2)/(o2+o+n2)
          f(ix)%data(:,lat,:) = 
     +      (2.*fwk(:,:,1)+fwk(:,:,2)+fwk(:,:,3)+fwk(:,:,4)+
     +      4.*fwk(:,:,5)+2.*fwk(:,:,6))/
     +      (fo2_cm3(:,:)+fo1_cm3(:,:)+fn2_cm3(:,:))
          deallocate(fwk)
!
! Total carbon = CO+CO2 volume mixing ratio:
        elseif (f(ix)%fname8 == 'CTOT    ') then
          if (allocated(fwk)) deallocate(fwk)
          allocate(fwk(imx,kmx,3),stat=ier)
          if (ier /= 0) 
     +      call allocerr(ier,"mkderived allocating 2 fwk for ctot")
          ixco =  ixfindc(h%fnames,h%nflds,'CO      ')
          ixco2 = ixfindc(h%fnames,h%nflds,'CO2     ')
          ixch4 = ixfindc(h%fnames,h%nflds,'CH4     ')
          if (ixco <= 0 .or. ixco2 <= 0 .or. ixch4 <= 0) then
            write(6,"('WARNING: need CO, CO2, and CH4 for CTOT: ',
     +        ' ixco,ixco2,ixch4=',2i3,' (will zero CTOT)')") 
     |        ixco,ixco2,ixch4
            f(ix)%data(:,lat,:) = 0.
            cycle fields_loop 
          endif
          call denconv(flat(:,:,ixco ),fwk(:,:,1),1,
     +      barm,pkt,28.,imx,kmx,-1,'MMR','CM3',dunit)
          call denconv(flat(:,:,ixco2),fwk(:,:,2),1,
     +      barm,pkt,44.,imx,kmx,-1,'MMR','CM3',dunit)
          call denconv(flat(:,:,ixch4),fwk(:,:,2),1,
     +      barm,pkt,16.,imx,kmx,-1,'MMR','CM3',dunit)
          f(ix)%data(:,lat,:) = (fwk(:,:,1)+fwk(:,:,2)+fwk(:,:,3))/
     +      (fo2_cm3(:,:)+fo1_cm3(:,:)+fn2_cm3(:,:))
          deallocate(fwk)
!
! Pressure (mb): use p=nkT, where n = rho (cm3), k=boltz
! Multiply by 1.e-3 to convert from dynes/cm2 to mb.
!
        elseif (f(ix)%fname8 == 'PMB     ') then
          f(ix)%data(:,lat,:) = ((fo2_cm3(:,:)+fo1_cm3(:,:)+
     |      fn2_cm3(:,:))*boltz*flat(:,:,ixt))*1.e-3
!
! Freeze point:
        elseif (f(ix)%fname8 == 'TNFP    ') then
          ixh2o = ixfindc(h%fnames,h%nflds,'H2O     ')
          dzp = (h%zpt-h%zpb)/float(h%nzp-1)
          do k=1,kmx
            zp = h%zpb+(k-1)*dzp
            press = p0*exp(-zp)*1.e-3
            do i=1,imx

              f(ix)%data(i,lat,k) = flat(i,k,ixt)-6077.4/(28.548-
     |          alog(flat(i,k,ixh2o)*barm(i,k)/rmh2o)-alog(press))

!             write(6,"('tnfp: lat=',i2,' k=',i2,' i=',i2,
!    |          ' tn=',e12.4,' h2o=',e12.4,' barm=',e12.4,' tnfp=',
!    |          e12.4)") lat,k,i,flat(i,k,ixt),flat(i,k,ixh2o),
!    |          barm(i,k),f(ix)%data(i,lat,k)

            enddo
          enddo
!
! jtgcm H2 = 1-HE-H (if jtgcm, H2 was made derived by sub setjtgcm)
! (jtgcm fh2 was defined at begining of this routine)
        elseif (f(ix)%fname8 == 'H2      '.and.h%isjtgcm) then
          f(ix)%data(:,lat,:) = fh2
!
! jtgcm ratio he/h2:
        elseif (f(ix)%fname8 == 'HE/H2   ') then
          f(ix)%data(:,lat,:) = fhe(:,:)/fh2(:,:)
!
! jtgcm ratio h/h2:
        elseif (f(ix)%fname8 == 'H/H2    ') then
          f(ix)%data(:,lat,:) = fh(:,:)/fh2(:,:)
!
! Unknown derived field:
        else
          write(6,"('>>> mkderived: unrecognized derived field ',
     +      a)") f(ix)%fname8
        endif
      enddo fields_loop
!
! Define requested Na species, if any:
      if (nna > 0) call getna(f,nf,flat,h,imx,kmx,lat)
!
      return
      contains
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine mkrhokg	! internal to mkderived
!
! Make rho in kg/m3 and store in fwk(imx,kmx,1):
! (used for epfluxes and qbary)
!
      if (allocated(fwk)) deallocate(fwk)
      allocate(fwk(imx,kmx,3),stat=ier)
      if (ier /= 0) 
     +  call allocerr(ier,"mkrhokg allocating fwk for rho")
!
! Get o2,o,n2 in gm/cm3 (store in fwk(:,:,1-3)):
!
      call denconv(flat(:,:,ixo2),fwk(:,:,1),1,
     +  barm,pkt,32.,imx,kmx,-1,'MMR','GM/CM3',dunit)
      call denconv(flat(:,:,ixo1),fwk(:,:,2),1,
     +  barm,pkt,16.,imx,kmx,-1,'MMR','GM/CM3',dunit)
      if (h%ismtgcm) then
        fwk(:,:,3) = flat(:,:,ixn2)
      else
        fwk(:,:,3)=max(.00001,(1.-flat(:,:,ixo2)-flat(:,:,ixo1)))
      endif
      call denconv(fwk(:,:,3),fwk(:,:,3),0,barm,pkt,28.,
     +  imx,kmx,-1,'MMR','GM/CM3',dunit)
!
! Sum o2,o,n2 and convert to kg/m3:
!
      fwk(:,:,1) = (fwk(:,:,1)+fwk(:,:,2)+fwk(:,:,3))*1000.
      end subroutine mkrhokg
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine mkco2(co2,idenconv)	! internal to mkderived
!
! Args:
      real,intent(out) :: co2(imx,kmx)
      integer,intent(in) :: idenconv
!
! Local:
      real :: hco2(imx,kmx)       ! co2 scale heights on mars
      real,parameter :: r=8.314e7 ! gas constant 
      real,parameter :: g=371.1   ! grav constant on mars
!
! Return co2 -- this is called only for [mv]tgcm
!
      ixco= ixfindc(h%fnames,h%nflds,'CO      ')
      if (ixco==0) then
        write(6,"('>>> mkderived: need CO to get CO2 ',
     +    'for [mv]tgcm: ixco=',i3)") ixco
        stop 'need CO'
      endif
      ixn2= ixfindc(h%fnames,h%nflds,'N2      ')
      if (ixn2==0.and.trim(h%version)/='VTGCM') then
        write(6,"('>>> mkderived: need N2 to get CO2 ',
     +    'for mtgcm: ixn2=',i3)") ixn2
        stop 'mtgcm N2'
      endif
      if (trim(h%version)/='VTGCM') then ! mtgcm
!
! As per Bougher (4-5/99), use "diffusive equilibrium extrapolation"
! if co2 mmr goes <= 0. This uses mean and co2 scale heights.
! H=R*T/g*Mbar, where R = gas constant = 8.314e7 and g (mars) = 350 (371?)
!
!       co2(:,:) = max(.00001,
!    +    (1.-flat(:,:,ixo1)-flat(:,:,ixco)-flat(:,:,ixn2)))
        co2(:,:) = (1.-flat(:,:,ixo1)-flat(:,:,ixco)-flat(:,:,ixn2))
        if (any(co2 <= 0.)) then
          co2(:,1:2) = max(.00001, ! insure bottom 2 levels are > 0
     +      (1.-flat(:,1:2,ixo1)-flat(:,1:2,ixco)-flat(:,1:2,ixn2)))
          hco2(:,:) = (r*flat(:,:,ixt))/(g*44.)
          do k=3,kmx
            do i=1,imx
              if (co2(i,k) <= 0.) then
                co2(i,k-1) = co2(i,k-2)*
     |                       exp(-(0.5*barm(i,k-1)/hco2(i,k-1)))
                co2(i,k)   = co2(i,k-1)*
     |                       exp(-(0.5*barm(i,k)/hco2(i,k)))
              endif
            enddo
          enddo
        endif
      else ! vtgcm
        if (h%isnew) then
          co2(:,:) = max(.00001,
     +      (1.-flat(:,:,ixo1)-flat(:,:,ixco)-flat(:,:,ixn2)))
        else
          co2(:,:) = max(.00001,
     +      (1.-flat(:,:,ixo1)-flat(:,:,ixco)-vtgcm_n2(:,:,lat)))
        endif
      endif
      if (idenconv > 0)
     +  call denconv(co2,co2,0,barm,pkt,44.,imx,kmx,iden,
     +    'MMR',' ',f(ix)%units)
      end subroutine mkco2
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      end subroutine mkderived
      end module mk_derived
