!
      module mk_auxflds
      use hist,only: history 
      implicit none
!
      contains
!-------------------------------------------------------------------
      subroutine wrdat_bf(lu,h,xylocations,prof_all,aux_all,nx,maxy,
     |  ny,maxz,nz,yy,ut,proclab)
      use proc
      use input
      implicit none
!
! Write ascii data file in Bougher format
!
! On input: zprange used
!   lu = logical unit for connection to file (open statement by this routine)
!   h = history structure with volume info for output
!   prof_all = collect profile data for output
!   aux_all = auxillary fields derived for Bougher output
!   nx,ny,nz = dimensions of profiles (locs,zps,flds)
!   maxy,maxz = 2nd and 3rd dimensions for prof_all array
!   yy = array of vertical coords
!   ut = model time (for computing local times as output
!   proclab = used to create file name (will be "proclab".dat)
!
! 9/14/04 btf: Using xylocations instead of xylocs (passed from mkxyloc).
!   Also fixed print of date.
!
! Args:
      integer,intent(in) :: lu,maxy,maxz,nx,ny,nz
      type(history),intent(in) :: h
      real,intent(in) :: ut,yy(ny),prof_all(nx,maxy,maxz) ! (loc,zpht,fld)
      real,intent(in) :: aux_all(mxloc,100,10)
      character(len=*),intent(in) :: proclab
      real,intent(in) :: xylocations(2,mxloc)
!
      logical,external :: isslt
      real,external :: fslt
      integer,external :: ixfindc
!
! Locals:
      character(len=8) :: date
      character(len=80) :: fname
      integer :: iloc,j,jj,k,nlocs,ix0
      real :: lon,slt
      real :: wrflds(mxloc,100,mxfproc)
!     real,pointer :: dum=>NULL()
      real :: dum=0.
!
! Get current date
      call date_and_time(date) ! ccyymmdd format
      write(6,"('wrdat_bf: date=',a)") date

      if (len_trim(bf_format)==0) then
        write(6,"(/,'>>> wrdat_bf: bf_format is empty.')")
        write(6,"('You must provide format bf_format in ',
     |    'namelist inputs in order to write data file.')")
        stop 'need bf_format in namelist input'
      endif
!
! Open file if necessary:
      fname = ' '
      write(fname,"(a,'.dat')") proclab
!     open(lu,file=fname,status='NEW')
      open(lu,file=fname,status='UNKNOWN')
      write(6,"('wrascii: opened file ',a,' with lu ',i3)")
     +  trim(fname),lu
!
! Build array for output
!
      jj = 0
      field_loop: do j=1,mxfproc
        if(aux_fields(j)=='') cycle field_loop
	if(aux_fields(j)(1:3)=='AUX') then
	  jj = jj + 1
	  read(aux_fields(j)(4:5),"(i2)") ix0
	  write(6,*) ' write field ',jj,' == ',
     +  	   aux_fields(j), j,ix0
	  do iloc=1,nx
	    wrflds(iloc,:,jj) = aux_all(iloc,:,ix0)
	  enddo
        else
          ix0 = ixfindc(cfields,nz,aux_fields(j))
	  if(ix0 > 0) then 
	    jj = jj + 1
	    write(6,*) ' write field ',jj,' == ',
     +  	     aux_fields(j), j,ix0
	    do iloc=1,nx
	      wrflds(iloc,:,jj) = prof_all(iloc,:,ix0)
	    enddo
	  else
	    write(6,*) ' field not found in aux_all or prof_all ',
     +  	     aux_fields(j)
	  endif
	endif
      enddo field_loop
!
      nlocs = 0
      do iloc=1,nx
        if (xylocations(1,iloc)/=spval.and.xylocations(2,iloc)/=spval) 
     +    nlocs = nlocs+1
      enddo

!
! Output data in Bougher-format
!
!   Print header info:
!
      write(lu,1001)
      write(lu,1002) date(5:6),date(7:8),date(1:4),
     +               h%version,trim(h%mssvol),h%mtime,
     +               trim(mgcm_fspath)
      write(lu,1001)
      write(lu,1003) (aux_fields(j),j=1,jj)

!1007 format(//,'ALTITUDE SCALE: '//,(6f10.2))
      write(lu,"('ALTITUDE SCALE: ',/,(6e10.2))") xyloc_htscale
!     write(lu,1007) xyloc_htscale

      write(lu,1004)
      do iloc=1,nlocs
        if (isslt(xylocs(2,iloc),slt)) then	! is local time
          write(6,"('mkauxflds slt: iloc=',i3,' xylocations=',2f8.2)")
     |      iloc,xylocations(:,iloc)
          write(lu,"(2f8.2,' (lat,slt)')") xylocations(:,iloc)
        else                                    ! is longitude
          write(6,"('mkauxflds lon: iloc=',i3,' xylocations=',2f8.2)")
     |      iloc,xylocations(:,iloc)
          write(lu,"(2f8.2,' (lat,lon)')") xylocations(:,iloc)
        endif

      enddo ! iloc
      write(lu,1001)
!
!   Print data:
!
      loc_loop: do iloc=1,nx
        if (xylocs(1,iloc)==spval.or.xylocs(2,iloc)==spval) 
     +    cycle loc_loop
        if (isslt(xylocs(2,iloc),slt)) then	! is local time
!          xyprof%glon = fslt(slt,h%ut,dum,3)
!          ixlon = ixfind(gcmlon,nlon,xyprof%glon,dlon)
!          xyprof%glon = gcmlon(ixlon)
        else				! is regular longitude
!          xyprof%glon = xylocations(2,iloc)
!          ixlon = ixfind(gcmlon,nlon,xyprof%glon,dlon)
           lon = xylocs(2,iloc)
           slt = fslt(dum,ut,lon,1)
        endif
	write(lu,1000) xylocations(1,iloc), slt
        ht_loop: do k=1,ny
          write(lu,bf_format) (wrflds(iloc,k,j),j=1,jj)
        enddo ht_loop
      enddo loc_loop
!
 1000 format(//20x,'LAT=',f10.2,2x,'LT=',f10.2//)
 1001 format(/,'**********************************************',/)
 1002 format(/,'DATE: ',a2,'/',a2,'/',a4,//,
     +       a10,a,'  (DAY,HR,MIN= ',3i4,')',//,3x,a)
 1003 format('FIELDS:',//,(5a10,/))
 1004 format(//,'LOCATIONS:',//)
!1005 format(f6.2,a10,/)
!1006 format(f6.2,f6.2,/)
!1007 format(//,'ALTITUDE SCALE: '//,(6f10.2))
!
      return
      end subroutine wrdat_bf
c-------------------------------------------------------------------
      subroutine mkauxdev(prof_all,ny,nf,nflds,f,aux_all)
c     Model Atmospheres auxillary quantities incl. (AUX1-7 FLDS)
c     Aerobraking auxillary quantities incl. DSHT  (AUX8   FLD)
      use proc
      use fields,only: field,mxnaf
      use input
      implicit none
!
! Args:
      real,intent(in) :: prof_all(mxloc,100,nf)
      real,intent(out) :: aux_all(mxloc,100,10)
      integer,intent(in) :: ny,nf,nflds
      type(field),intent(in) :: f(nf)
!
! Locals:
      integer :: iloc,ix,k,nfd,i,j,ifu,ifv,ier,nna
      integer :: idefu=0,idefv=0
      integer ::
     +  ixo2,ixo1,ixt,ixz,ixu,ixv,ixw,ixne,ixte,ixo2p,ixo21d,ixco2,
     +  ixn2,ixno,ixh,ixo3,ixho2,ixoh,ixco,ixh2o,ixch4,ixh2
      integer :: iepvy,iepvz,iepvdiv
      real,allocatable :: fwk(:,:,:)
      character(len=16) :: dunit
      real :: fmin,fmax,dht
      real :: reqlocs_tmp
!
! Externals:
      integer,external :: ixfindc
!
! get indices for required fields
!
       ixt   = ixfindc(cfields,nf,'TN      ')
       ixo1  = ixfindc(cfields,nf,'O1      ')
       ixco  = ixfindc(cfields,nf,'CO      ')
       ixn2  = ixfindc(cfields,nf,'N2       ')
       ixo2  = ixfindc(cfields,nf,'O2       ')
       ixco2 = ixfindc(cfields,nf,'CO2     ')
       ixz   = ixfindc(cfields,nf,'Z       ')

      if (ixt==0.or.ixo1==0.or.ixco==0.or.ixn2==0.or.ixo2==0.or.
     |  ixco2==0.or.ixz==0) then
        write(6,"(/,'>>> mkauxdev: one or more of the following fields',
     |    ' is not available in namelist cfields: ',
     |    'TN,O1,CO,N2,O2,CO2,Z')")
        write(6,"('These fields must be available to make ',
     |    'aux fields.')")
        write(6,"('ixt,ixo1,ixco,ixn2,ixo2,ixco2,ixz=',7i4)")
     |    ixt,ixo1,ixco,ixn2,ixo2,ixco2,ixz
        stop 'fields missing for mkauxdev'
      endif
!
! Auxillary Bougher fields (model atmospheres plus AB versions)
! Requires Altitude be in km units
!
      loc_loop: do iloc=1,mxloc
        if (xylocs(1,iloc)==spval.or.xylocs(2,iloc)==spval)
     +                                          cycle loc_loop
!     Full height loop --------------
	ht_loop: do k=1,ny
! Subscript #3 of the array PROF_ALL has value 0 which is less than the lower bound of 1
          aux_all(iloc,k,1)= prof_all(iloc,k,ixo1)             ! XTOT
     +  		   + prof_all(iloc,k,ixco)
     +  		   + prof_all(iloc,k,ixn2)
     +  	           + prof_all(iloc,k,ixco2)
          aux_all(iloc,k,2)=(prof_all(iloc,k,ixo1)*16.         ! XMASS
     +  		   + prof_all(iloc,k,ixco)*28.
     +  		   + prof_all(iloc,k,ixn2)*28.
     +  	           + prof_all(iloc,k,ixco2)*44.)/
     +                       aux_all(iloc,k,1)
          aux_all(iloc,k,3)= aux_all(iloc,k,1) * 1.38E-16      ! PR(microbar)
     +                     * prof_all(iloc,k,ixt)
          aux_all(iloc,k,4)= alog(1.2e-03/aux_all(iloc,k,3))   ! ALOG(PR)
          aux_all(iloc,k,5)= aux_all(iloc,k,3) *               ! RHO (KG/KM3)
     +                       aux_all(iloc,k,2) /
     +                       (8.314e+07*prof_all(iloc,k,ixt))
     +                       * 1.0E+12
          aux_all(iloc,k,6)= 371.1 * (3388.25/                 ! GZ (cm/s2)
     +   		     (3388.25+prof_all(iloc,k,ixz)))**2.0
          aux_all(iloc,k,7)= 8.314E+07*prof_all(iloc,k,ixt)/   ! SHT (KM)
     +  		     (aux_all(iloc,k,2)*aux_all(iloc,k,6))
     +  		     * 1.0E-05
          aux_all(iloc,k,8)= 0.0                               ! DSHT(km)
        enddo ht_loop
!     Partial height loop --------------
!     DSHT : CENTER DIFFERENCING MOST ACCURATE / VARIABLE KM GRID
       do k=2,ny-1
	  DHT = (prof_all(iloc,k+1,ixz) - prof_all(iloc,k-1,ixz)) ! DHT(km)
          aux_all(iloc,k,8)= -1./log(aux_all(iloc,k+1,5)/
     +          aux_all(iloc,k-1,5))*DHT                       ! DSHT (km)
       enddo
      enddo loc_loop
!
      return
      end subroutine mkauxdev
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      end module mk_auxflds
