!
      module den_convert
      use hist,only: history
      use get_hist,only: vtgcm_n2
      implicit none
!
! Density species conversion.
! dunits(ndunits) is list of valid density unit strings:
! MMR    = mass mixing ratios
! CM3    = number densities
! CM3-MR = number density mixing ratios
! GM/CM3 = mass density 
! Order of dunits corresonds to iden=0,1,2,3, for
!  conversion from history:
!   iden = 0 -> MMR to MMR (no change from history)
!   iden = 1 -> MMR to CM3
!   iden = 2 -> MMR to CM3-MR
!   iden = 3 -> MMR to GM/CM3
!
      real,parameter :: boltz=1.3805e-16
      integer,parameter :: ndunits=4, len_units=16
      character(len=len_units),save :: dunits(ndunits) = 
     +  (/'MMR     ','CM3     ','CM3-MR  ','GM/CM3  '/)
!
! (barm is different for mtgcm, vtgcm, jtgcm, see below)
      real,allocatable,save ::
     +        barm(:,:),	! 1./(o2/32.+o1/16.+n2/28.)
     +        pkt(:,:)		! p0*e(-z)/kT
      contains
!-------------------------------------------------------------------
      subroutine mkdenparms(flat,h,imx,kmx,lat)
      use proc,only: p0
!
! Define barm and pkt at current latitude for use in density 
!   conversions (if vtgcm, lat is only 1->18 and kmx=venus_kmx=35):
! 2/2/98: check for mtgcm (mars) histories.
! 8/15/02: check for jtgcm (jupiter) histories.
! 
! Args:
      integer,intent(in) :: imx,kmx,lat
      type(history),intent(in) :: h
      real,intent(in) :: flat(imx,kmx,h%nflds)
!
! Locals:
      integer :: ixo1,ixo2,ixt,k,ier,ixco,ixn2,ixh,ixhe
      real,save :: dz
      real :: fmin,fmax
!
! Externals:
      integer,external :: ixfindc
!
! Make necessary allocations:
      if (.not.allocated(barm)) then
        allocate(barm(imx,kmx),stat=ier)
        if (ier /= 0) call allocerr(ier,"mkdenparms allocate barm")
      endif
      if (.not.allocated(pkt)) then
        allocate(pkt(imx,kmx),stat=ier)
        if (ier /= 0) call allocerr(ier,"mkdenparms allocate pkt")
      endif
      dz = (h%zpt-h%zpb)/float(kmx-1)
      ixo1 = ixfindc(h%fnames,h%nflds,'O1      ')
      ixo2 = ixfindc(h%fnames,h%nflds,'O2      ')
      ixt  = ixfindc(h%fnames,h%nflds,'TN      ')
      if (.not.h%ismtgcm.and..not.h%isjtgcm.and..not.h%isvtgcm) then
        do k=1,kmx
          barm(:,k)=1./(flat(:,k,ixo2)/32.+flat(:,k,ixo1)/16.+
     +      max(.00001,(1.-flat(:,k,ixo2)-flat(:,k,ixo1)))/28.)
          pkt(:,k)=p0*exp(-(h%zpb+(k-1)*dz))/(boltz*flat(:,k,ixt))
        enddo
!
! For jtgcm, barm = 1./(he/4 + h/1 + h2/2.2)
      elseif (h%isjtgcm) then
        ixh  = ixfindc(h%fnames,h%nflds,'H       ')
        if (ixh==0) write(6,"('>>> WARNING mkdenparms: need H to ',
     |      'calculate barm for jtgcm history.')") 
        ixhe = ixfindc(h%fnames,h%nflds,'HE      ')
        if (ixhe==0) write(6,"('>>> WARNING mkdenparms: need HE to ',
     |      'calculate barm for jtgcm history.')") 
        do k=1,kmx
          barm(:,k) = 1./(flat(:,k,ixhe)/4. + flat(:,k,ixh)/1. +
     |      (1.-flat(:,k,ixhe)-flat(:,k,ixh))/2.2)
          pkt(:,k)=p0*exp(-(h%zpb+(k-1)*dz))/(boltz*flat(:,k,ixt))
        enddo 
      else	! mtgcm or vtgcm history
!
! For [mv]tgcm histories, barm = 1/(o/16+(co+n2)/28.+co2/44.) where
!   co2=1-o-co-n2 (mmr) (n2 is on mtgcm histories, and for vtgcm
!   histories, n2 is obtained from an external dataset, currently
!   defined in vtgcm_n2())
! vtgcm_n2(imxp3,venus_,kmx,venus_nlat) is read from external data 
! 
        ixco = ixfindc(h%fnames,h%nflds,'CO      ')
        if (ixco==0) then
          write(6,"('>>> WARNING mkdenparms: need CO to ',
     +      'calculate barm for [mv]tgcm history.')") 
        endif
        if (h%ismtgcm) then ! is mtgcm
          ixn2 = ixfindc(h%fnames,h%nflds,'N2      ')
          if (ixn2==0) then
            write(6,"('>>> WARNING mkdenparms: need N2 to ',
     +        'calculate barm for mtgcm history.')") 
          endif
          do k=1,kmx
            barm(:,k) = 1./(flat(:,k,ixo1)/16.+(flat(:,k,ixco)+
     +        flat(:,k,ixn2))/28.+max(.00001,(1.-flat(:,k,ixo1)-
     +        flat(:,k,ixco)-flat(:,k,ixn2)))/44.)
            pkt(:,k)=p0*exp(-(h%zpb+(k-1)*dz))/(boltz*flat(:,k,ixt))
          enddo
!
! If old vtgcm history, vtgcm_n2 was read from external data file
! by getvtgcm (gethist.F). If new vtgcm history, N2 was read from
! source history, so is in flat(...ixn2).
!
        else
          ixn2 = ixfindc(h%fnames,h%nflds,'N2      ')
          if (ixn2==0) then ! old history, use vtgcm_n2
            write(6,"('>>> WARNING mkdenparms: ixn2=',i3,' will use',
     +        'vtgcm_n2 to calculate barm for vtgcm history.')") ixn2
            do k=1,kmx
              barm(:,k) = 1./(flat(:,k,ixo1)/16.+(flat(:,k,ixco)+
     |          vtgcm_n2(1:imx,k,lat))/28.+max(.00001,
     |          (1.-flat(:,k,ixo1)-flat(:,k,ixco)-vtgcm_n2(:,k,lat)))/
     |           44.)
              pkt(:,k)=p0*exp(-(h%zpb+(k-1)*dz))/(boltz*flat(:,k,ixt))
            enddo ! k=1,kmx
          else             ! new history, use flat(..ixn2)
            do k=1,kmx
              barm(:,k) = 1./(flat(:,k,ixo1)/16.+(flat(:,k,ixco)+
     |          flat(:,k,ixn2))/28.+max(.00001,
     |          (1.-flat(:,k,ixo1)-flat(:,k,ixco)-flat(:,k,ixn2)))/
     |           44.)
              pkt(:,k)=p0*exp(-(h%zpb+(k-1)*dz))/(boltz*flat(:,k,ixt))
            enddo ! k=1,kmx
          endif ! ixn2
        endif
      endif
      end subroutine mkdenparms
!-------------------------------------------------------------------
      subroutine denconvert(f,nflds,flat,h,imx,kmx,lat,iden)
      use fields,only: field
!
! Do unit conversion of density fields in f(nflds) according to iden.
! (this routine called only if iden > 0)
! (Called from getflds)
!
! Args:
      integer,intent(in) :: nflds	! number of fields in f
      type(field),intent(inout) :: f(nflds)	! fields to be defined
      type(history) :: h		! the current history
      integer,intent(in) :: imx,kmx	! number of lons and zp levels
      real,intent(inout) :: flat(imx,kmx,h%nflds) ! the current lat slice
      integer,intent(in) :: lat		! current lat index
      integer,intent(in) :: iden	! density conversion flag
!
! Locals:
      integer :: ixo3,ixo2,ixo1,ixno,ixno2,ixh,ixho2,ixoh
      integer :: ix,k,id1,id2
      real :: fmin,fmax
      real,allocatable :: ftmp(:,:)
!
! Externals:
      integer,external :: ixfindc
!
! Use f(nflds) fields to do density conversion:
!
      do ix=1,nflds
        if (f(ix)%type .eq. 'DENSITY         '.and.
     +      .not. f(ix)%derived .and.
     +      f(ix)%wt > 0. .and.
     +      associated(f(ix)%data)) then
!
! 4/03 bf:
! For temp array. Note calls to denconv must have isetout=1 for AIX
! (Setting isetout=0 (i.e., overwriting input array) worked on SGI, 
!  but input array was unchanged under AIX). 
!
          id1 = size(f(ix)%data(:,lat,1))
          id2 = size(f(ix)%data(1,lat,:))
          allocate(ftmp(id1,id2))
!
          select case (f(ix)%fname8)
! OX (O+O3):
          case ('OX      ') 
            ixo3  = ixfindc(h%fnames,h%nflds,'O3      ')
            ixo1  = ixfindc(h%fnames,h%nflds,'O1      ')
            if (ixo3 == 0 .or. ixo1 == 0) then
              if (lat == 1)
     +          write(6,"('>>> WARNING: need O3 and O1 to do OX ',
     +            'density conversion -- OX will be in MMR.')")
              f(ix)%units = 'MMR'
            else
              f(ix)%data(:,lat,:) = f(ix)%data(:,lat,:) * f(ix)%wt *
     +           (flat(:,:,ixo3)/48.+flat(:,:,ixo1)/16.)/
     +           (flat(:,:,ixo3)+flat(:,:,ixo1))
              call denconv(f(ix)%data(:,lat,:),ftmp,
     +          1,barm,pkt,f(ix)%wt,imx,kmx,iden,'MMR',' ',f(ix)%units)
              f(ix)%data(:,lat,:) = ftmp(:,:)

            endif
! NOZ (NO+NO2):
          case ('NOZ     ') 
            ixno  = ixfindc(h%fnames,h%nflds,'NO      ')
            ixno2 = ixfindc(h%fnames,h%nflds,'NO2     ')
            if (ixno == 0 .or.ixno2 == 0) then
              if (lat == 1)
     +          write(6,"('>>> WARNING: need NO, NO2 to do NOZ ',
     +            'density conversion -- NOZ will be in MMR.')")
              f(ix)%units = 'MMR'
            else
              f(ix)%data(:,lat,:) = f(ix)%data(:,lat,:) * f(ix)%wt *
     +           (flat(:,:,ixno)/30.+flat(:,:,ixno2)/46.)/
     +           (flat(:,:,ixno)+flat(:,:,ixno2))
              call denconv(f(ix)%data(:,lat,:),ftmp,
     +          1,barm,pkt,f(ix)%wt,imx,kmx,iden,'MMR',' ',f(ix)%units)
              f(ix)%data(:,lat,:) = ftmp(:,:)

            endif
! HOX (OH+HO2+H):
          case ('HOX     ') 
            ixoh  = ixfindc(h%fnames,h%nflds,'OH      ')
            ixho2 = ixfindc(h%fnames,h%nflds,'HO2     ')
            ixh   = ixfindc(h%fnames,h%nflds,'H       ')
            if (ixoh == 0 .or.ixho2 == 0 .or. ixh == 0) then
              if (lat == 1)
     +          write(6,"('>>> WARNING: need OH, HO2, and H to do HOX ',
     +            'density conversion -- HOX will be in MMR.')")
              f(ix)%units = 'MMR'
            else
              f(ix)%data(:,lat,:) = f(ix)%data(:,lat,:) * f(ix)%wt *
     +           (flat(:,:,ixoh)/17.+flat(:,:,ixho2)/33.+flat(:,:,ixh))/ 
     +           (flat(:,:,ixoh)+flat(:,:,ixho2)+flat(:,:,ixh))
              call denconv(f(ix)%data(:,lat,:),ftmp,
     +          1,barm,pkt,f(ix)%wt,imx,kmx,iden,'MMR',' ',f(ix)%units)
              f(ix)%data(:,lat,:) = ftmp(:,:)
            endif
!
! HTOT and CTOT remain volume mixing ratios, as defined in mkderived:
          case ('HTOT    ')
          case ('CTOT    ')
! Normal (non-family) density species:
          case default
            call denconv(f(ix)%data(:,lat,:),ftmp,
     +        1,barm,pkt,f(ix)%wt,imx,kmx,iden,'MMR',' ',f(ix)%units)
            f(ix)%data(:,lat,:) = ftmp(:,:)
          end select
          if (allocated(ftmp)) deallocate(ftmp)
        endif
      enddo
      end subroutine denconvert
!
!-------------------------------------------------------------------
!
      subroutine denconv(fin,fout,isetfout,barm,pkt,w,
     +  id1,id2,iden,units_frm,units_to,units_out)
!
! Given field fin(id1,id2) in units units_frm, convert to
!   units_to, placing converted field in fout if isetfout, or
!   in fin otherwise.
! On input:
!   fin(id1,id2)	input field
!   fout(id1,id2) 	output field (if isetfout>0)
!   isetfout (int) 	fout flag
!
! bf 7/14/03: Calling this routine with isetfout==0 when fin is a
! pointer (e.g., f(ix)%data) does NOT work on the IBM. Must call with
! isetfout==1 in this case (if fin is a local array, isetfout==0 does work)
!
!   barm(id1,id2) 	1./(o2/32+o1/16+n2/28) (made from mmr)
!   pkt(id1,id2)	p0*e(-z)/kT   
!   w			molecular weight of input field
!   id1,id2		dimensions of input (and output) arrays
!   iden		if 0 <= iden <= 3, is density conversion flag
!                       (if iden < 0, then units_frm and units_to are used)
!   units_frm,units_to	if iden < 0, is units from and to for conversion, 
!                         otherwise these are ignored
!                         (must be one of dunits)
!   units_out		set to units_to if no error, otherwise blanks
! On output:
!   If isetfout, then converted field is in fout, otherwise is in fin
!     (i.e., fin is changed if isetfout==0)
!   units_out has been set to units_to if no error, blank otherwise
!   All other args unchanged.
!
! Args:
      integer,intent(in) :: id1,id2,iden,isetfout
      real,dimension(id1,id2),intent(inout) :: fin
      real,dimension(id1,id2),intent(in) :: barm,pkt
      real,dimension(id1,id2),intent(out) :: fout
      real,intent(in) :: w
      character(len=*),intent(in) :: units_frm,units_to
      character(len=*),intent(out) :: units_out
!
! Local:
      real,dimension(id1,id2) :: f	     ! auto-array
      character(len=len_units) :: ufrm,uto   ! local unit strings
      character(len=len_units*2+4) :: ufrmto ! e.g. "MMR_to_CM3"
      integer :: len_frm,len_to,len_frmto    ! trimmed char lengths
      integer :: ix
      real :: fmin,fmax
!
! Externals:
      integer,external :: ixfindc
!
      len_frm = len_trim(units_frm)
      len_to = len_trim(units_to)
!
! Use iden flag if set, then units_frm should be MMR:
!   iden = 0 -> leave species as on history (most are mass mixing ratio)
!   iden = 1 -> convert species to number densities (cm3)
!   iden = 2 -> convert species to number density mixing ratios
!   iden = 3 -> convert species to mass density (gm/cm3)
!
      if (iden >= 0) then
        if (len_frm > 0 .and. units_frm(1:len_frm) /= 'MMR') then
          write(6,"('>>> WARNING denconv: iden=',i2,' so expected',
     +      ' units_frm to be MMR. Instead units_frm=',a)")
     +      iden,units_frm
        endif
        ufrm = 'MMR'
        select case (iden)	! could use uto=dunits(iden+1)
        case(0)			! MMR to MMR (no change)
          uto = 'MMR'
        case(1)			! MMR to CM3
          uto = 'CM3'
        case(2)			! MMR to CM3-MR
          uto = 'CM3-MR'
        case(3)			! MMR to GM/CM3
          uto = 'GM/CM3'
        case default
          write(6,"('>>> WARNING denconv: unknown iden value=',i4)")
     +      iden
        end select
        len_to = len_trim(uto)
        len_frm = len_trim(ufrm)
      else			! iden < 0 -> use units_frm,to
!
! Maybe put this validation stuff in internal routine.
        if (len_frm > len_units) then
          write(6,"('>>> WARNING denconv: units_frm too long',
     +      ' (must be <= ',i2,')')") len_units
          len_frm = len_units
          return
        endif
        if (len_to > len_units) then
          write(6,"('>>> WARNING denconv: units_to too long',
     +      ' (must be <= ',i2,')')") len_units
          len_to = len_units
          return
        endif
        ix = ixfindc(dunits,ndunits,units_frm)
        if (ix == 0) then
          write(6,"('>>> WARNING denconv: unrecognized units_frm: ',
     +      a)") units_frm 
          write(6,"('  valid density units are:',/(4a16))") dunits
          return
        endif   
        ix = ixfindc(dunits,ndunits,units_to)
        if (ix == 0) then
          write(6,"('>>> WARNING denconv: unrecognized units_to: ',
     +      a)") units_to 
          write(6,"('  valid units are:',/(4a16))") dunits
          return
        endif   
        ufrm = units_frm(1:len_frm)
        uto  = units_to(1:len_to)
      endif	! iden >= 0
      if (ufrm(1:len_frm) == uto(1:len_to)) then
        f(:,:) = fin(:,:)
        goto 100
      else
        write(ufrmto,"(a,'_to_',a)") ufrm(1:len_frm),uto(1:len_to)
      endif
      len_frmto = len_trim(ufrmto)
!
! Validate w:
      if (w < 1.e-20 .or. w > 80.) then
        write(6,"('>>> WARNING denconv: bad w (mol weight)=',e12.4,
     +    ' no conversion done')") w
        return
      endif
!
! Now have units from and to in ufrmto. Do the conversion:
!
      if (len_frmto > 0) then
        select case (ufrmto(1:len_frmto))
          case ('MMR_to_CM3') 
            f = fin * pkt * barm / w 
          case ('MMR_to_CM3-MR') 
            f = fin * barm / w
          case ('MMR_to_GM/CM3') 
            f = fin * pkt * barm * 1.66e-24
          case ('CM3_to_MMR') 
            f = fin * w / (pkt*barm)
          case default
            write(6,"('>>> WARNING denconv: unrecognized density',
     +        ' conversion: ',a,' (iden=',i3,')')") 
     +        ufrmto(1:len_frmto),iden
        end select
      endif ! len_frmto > 0
!
! Return converted field according to isetfout:
!
 100  continue 
      if (isetfout > 0) then
        fout = f
      else
        fin = f
      endif
      write(units_out,"(a)") uto(1:len_to)
      return
      end subroutine denconv
      end module den_convert
