c****************************************************************
      SUBROUTINE REMSTR(aa,bb,out)
      CHARACTER (LEN=80)      ::    aa,out
      CHARACTER (LEN=1)       ::    bb              !only single chars! 

      out = ''
      i = 1
      j = 0
      DO WHILE (i.LE.LEN(aa))
         IF (aa(i:i).NE.bb) THEN
            j = j + 1
            out(j:j) = aa(i:i)
         END IF
         i = i + 1
      END DO

      END SUBROUTINE
c****************************************************************
c
      PROGRAM ATMBIN
c     
c     convert binary - ascii, and MARCS to SIU binary
c     
c      INCLUDE 'atmos_definitions.inc'
      PARAMETER(NDMAX=80, NATOM=92, NELMAX=92, NMOLMAX=14)
      PARAMETER(NR_ION=3, MAXBREAK=2, MAX_ETA_CONT=2)
c
c define parameters for MARCS models
c
      PARAMETER(ndep_m = 80, nion_m = 48, nel_m = 16, nlev_m=21)
      PARAMETER(nmol_m = 0, nbreak_m = 2, nel_new_m = 63)

      REAL*8 boltz
      INTEGER atomnr(nel_m), ionst(nel_m), break(nbreak_m)
      CHARACTER filein*80,fileout*80, foutbin*80
      CHARACTER filein1*80,fileout1*80

      CHARACTER*60 out
      CHARACTER*1  bb    
      DATA boltz/1.380658d-16/
      DATA atomnr /1,2,6,7,8,10,11,12,13,14,16,19,20,24,26,28/ ! not used
      DATA ionst /2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3/             ! not used
      DATA break /15, 50/               ! break points in tau scale for MARCS!
c
      INTEGER*4 n_unit, n_eli
      INTEGER*4 ndepth                                   
      INTEGER*4 nmol                               
      INTEGER*4 nel                                
      INTEGER*4 nbreak                        
c
      STRUCTURE /atm_head2_type/
        CHARACTER  date*30, user*20              
        REAL*4 teff, logg, x, y, z, xi, alpha
        INTEGER*2     status            ! Atm Status       
        INTEGER*2     eps_nr            ! Nr of atm element abundances
        REAL*4        eps(NATOM)        ! atm element abundances
        REAL*4        eps_dev(NATOM)    ! temp. deviation of element abundances
        INTEGER*2     ref_nr            ! Nr of Ref.-character (atm)
        CHARACTER*1000 ref              ! atm references   
        INTEGER*2      cmt_nr           ! Nr of atm comment character
        CHARACTER*1000 cmt              ! atm comment

        INTEGER*4   nscale(MAXBREAK)    ! breakpoints of tau-scale
        INTEGER*4   idx_el(NELMAX)      ! atomic-numbers of used elements
        INTEGER*4   nion  (NELMAX)      ! Ionisation steps taken into account
        INTEGER*4   eta_cont_nr         ! Amount of continuum opacities
        CHARACTER*2 cel (NELMAX)        ! Name of elements concerned
        CHARACTER*4 cmol(NMOLMAX)       ! name of molecules concerned
      END STRUCTURE
c
      STRUCTURE /atm_type/
        RECORD /atm_head2_type/ h
         INTEGER*4   ndepth              ! Nr of depthpoints
         INTEGER*4   nmol                ! Nr of different molecules
         INTEGER*4   nelem               ! Nr of different molecule formingelements
         INTEGER*4   nel                 ! Nr of different elements
         INTEGER*4   nbreak              ! Nr of breakpoints in tau-scale
         INTEGER*4   nscale(MAXBREAK)    ! breakpoints of tau-scale
         INTEGER*4   idx_el(NELMAX)      ! atomic-numbers of used elements
         INTEGER*4   nion  (NELMAX)      ! Ionisation steps taken into account
        REAL*4  tau(NDMAX)              ! tau(reference) depth scale
        REAL*4  logtau(NDMAX)           ! log(tau(reference))
        REAL*4  t  (NDMAX)              ! Temperature structure
        REAL*4  pe (NDMAX)              ! electron pressure per depth
        REAL*4  pg (NDMAX)              ! gas pressure          
        REAL*4  pp (NDMAX,NR_ION)       ! partial pressures of each ion class
        REAL*4  pk (NDMAX)              ! stratif. of Pnuc      
        REAL*4  rho(NDMAX)              ! mass density stratif. 
        REAL*4  pmol(NDMAX,NMOLMAX)     ! molecular partial pressures / partitionfunction
        REAL*4  umol(NDMAX,NMOLMAX)     ! molecular partial pressures / partitionfunction
        REAL*4  pion(NDMAX,NR_ION,NATOM)! atomic partial pressures    / partitionfunction
        REAL*4  uion(NDMAX,NR_ION,NATOM) ! atomic partial pressures    / partitionfunction
        REAL*4  kappa_ref(NDMAX)        ! reference kappa (i.e. kappa(5000A))
        REAL*4  etasig  (NDMAX,MAX_ETA_CONT)! scattering coefficient/kapref
        REAL*4  eta_cont(NDMAX,MAX_ETA_CONT)! continuum opacities
        REAL*8  eta_cont_w(MAX_ETA_CONT) ! Wavelength points of continuum opacities
      END STRUCTURE
      RECORD/atm_type/atm_int(1)
c
c DEFINITIONS FOR MARCS
c

      REAL*4  ABUN(nel_new_m), ab_m(nel_new_m)  ! Abundance of 16 elements on 12-hydrogen scale
      INTEGER*2 ATNR(nel_m)           ! Atomic number of the elements
      INTEGER*2 INST(nel_m)           ! Number of ionisation stages considered for each element
      CHARACTER IONS(nion_m)          ! String array containing the name of all ions
      REAL*4  KREF(ndep_m)            ! Kappa @500nm 
      REAL*4  ABNLTE(nlev_m, ndep_m)  ! Non-LTE abundance for each depth point and atomic level
      INTEGER*4 NDEP_MARCS                   ! n depth points (already defined in the parameter)
      REAL*4  ABLTE(nlev_m, ndep_m)          ! LTE abundance for each depth point and atomic level      
      REAL*8  PARTEL(nion_m, ndep_m)  ! Partition function for each depth point and ion
      REAL*4  PELE(ndep_m)                   ! Electron pressure
      REAL*8  PGAS(ndep_m)                   ! Gas pressure 
      REAL*8  PRESEL(nion_m, ndep_m)  ! Partial pressure for each depth point and ion
      REAL*4  RHOG(ndep_m)                   ! Density
      REAL*4  TAU5(ndep_m)            ! Tau @500nm
      REAL*4  TEMP(ndep_m)            ! Temperature  

      REAL*4  pgas_m(ndep_m)            ! p gas MARCS in single precision
      REAL*4  part_m(ndep_m,NR_ION,nel_m)  ! partition functions transposed
      REAL*4  pres_m(ndep_m,NR_ION,nel_m)  ! partial pressures transposed

      INTEGER*4 at_idx(nel_m)           ! Atomic number of the elements
      INTEGER*4 at_nion(nel_m)          ! Number of ionisation stages considered for each element

      OPEN (4,file='model_lst.inp',status='old')
      OPEN (6,file='model.lst',status='UNKNOWN')
      READ (4,711) filein1, fileout1
      READ (4,711)
      READ (4,711) filein, fileout
      READ (4,711) foutbin
      PRINT 711,   filein, fileout
      PRINT 711,   foutbin
 

  711 FORMAT(A80)
8112  FORMAT(5x,f6.0,7x,f5.2,2f8.3) 
11    FORMAT(f7.1,f8.1,2f8.3) 
5888  FORMAT(i3,E13.6,F9.1,2E13.6,f13.5,E13.6)  
712   FORMAT(10I4)
713   FORMAT(7F13.6,/,7F13.6)                  !adapt format to the number of molecules!
714   FORMAT(I4,14F13.6)  
715   FORMAT(1X,A5,3X,A6,8X,A1,8X,A10,3(8X,A2),12X,A3)

      i = 1
      filein = TRIM(ADJUSTR(filein))

c      IF (filein(LEN(filein)-2:LEN(filein)).EQ.'dat')  THEN     !assume all MAFAGS bin files are dat

c have to go through this step now to put the header to the MODEL, which
c is now just copy-paste from analogous MAFAGS model --> should be updated later     

        OPEN (2,file = filein1,status='old',
     *    FORM = 'UNFORMATTED')
       OPEN (27, FILE = fileout1,STATUS = 'UNKNOWN', FORM = 'FORMATTED')
        READ(2)
     *     ndepth,nel,nmol,nbreak,
     *    (atm_int(i).idx_el(n)    ,n=1,nel   ),
     *    (atm_int(i).nion(n)      ,n=1,nel   ),
     *    (atm_int(i).nscale(n)    ,n=1,nbreak),
     *    (atm_int(i).tau(n)       ,n=1,ndepth),
     *    (atm_int(i).t(n)         ,n=1,ndepth),
     *    (atm_int(i).kappa_ref(n) ,n=1,ndepth),
     *    (atm_int(i).pe(n)        ,n=1,ndepth),
     *    (atm_int(i).pg(n)        ,n=1,ndepth),
     *    (atm_int(i).pk(n)        ,n=1,ndepth),
     *    (atm_int(i).rho(n)       ,n=1,ndepth),
     *   ((atm_int(i).pmol(n,k)  ,n=1,ndepth),k=1,nmol),
     *   ((atm_int(i).umol(n,k)  ,n=1,ndepth),k=1,nmol),
     *  (((atm_int(i).pion(n,m,k),n=1,ndepth),m=1,
     *                               atm_int(i).nion(k)),k=1,nel),
     *  (((atm_int(i).uion(n,m,k),n=1,ndepth),m=1,
     *                               atm_int(i).nion(k)),k=1,nel)

        atm_int(i).ndepth      = ndepth 
        atm_int(i).nel         = nel
        atm_int(i).nmol        = nmol
        atm_int(i).nbreak      = nbreak

        READ(2)  atm_int(i).h.date,atm_int(i).h.user,
     *                atm_int(i).h.ref_nr,atm_int(i).h.teff,
     *                atm_int(i).h.logg,atm_int(i).h.x,atm_int(i).h.y,
     *                atm_int(i).h.z,atm_int(i).h.alpha,
     *                atm_int(i).h.eps_nr,atm_int(i).h.status,
     *                atm_int(i).h.cmt_nr,
     *                atm_int(i).h.ref(:atm_int(i).h.ref_nr),
     *               (atm_int(i).h.eps(k),k=1,atm_int(i).h.eps_nr),
     *                atm_int(i).h.cmt(:atm_int(i).h.cmt_nr)

        WRITE(27,"('N depth points',I4)") ndepth
        WRITE(27,"('N of elements',I4)") nel
        WRITE(27,"('N of molecules',I4)") nmol
        WRITE(27,"('N of break points in tau scale',I4)") nbreak
         WRITE(27,"('Atomic numbers of elements')")
         WRITE(27,712) (atm_int(i).idx_el(n), n=1, nel)
         WRITE(27,"('Number of ionization ',
     *                     'stages for each element')")
         WRITE(27,712) (atm_int(i).nion(n), n=1, nel)
         WRITE(27,"('Logarithmic element ',
     *            'abundances, A = log(N_el/N_H) + 12')")
         WRITE(27,"(10(F5.2,2X))") (atm_int(i).h.eps(k),
     *                           k=1,atm_int(i).h.eps_nr)
         WRITE(27,"('Break poins in tau-scale', 10I4)")
     *             (atm_int(i).nscale(n), n=1, nbreak)
         WRITE(27,715) 'Depth','lgTau5','T','Kappa_ref',
     *                 'Pe','Pg','Pk','Rho'
        DO n=1, ndepth
c
c  TEMPERATURE IS NOT IN LOG SCALE IN THE ASCII FILE (SIU uses LOG!)
c
         WRITE(27,714) n, atm_int(i).tau(n),
     *                 10**(atm_int(i).t(n)), atm_int(i).kappa_ref(n),
     *                 atm_int(i).pe(n), atm_int(i).pg(n),
     *                 atm_int(i).pk(n), atm_int(i).rho(n)
        ENDDO
c
         WRITE(27,"(/,'Partial pressures ',
     *             'for ',
     *              I2,' molecules')") nmol
         WRITE(27,"('H2,CH,NH,OH,C2,CN,CO,N2,NO,O2,MgH,SiH,CaH,TiO')")
        DO n=1, ndepth
        WRITE(27,714) n, (atm_int(i).pmol(n,k),k =1, nmol)
        ENDDO
          WRITE(27,"('Partition functions for molecules')")
        DO n=1, ndepth
        WRITE(27,714) n, (atm_int(i).umol(n,k),k =1, nmol)
        ENDDO
        n_eli = nel*3-3
        WRITE(27,"(/,'Partial pressures ',
     *             'and partition functions for ',
     *           I2,' elements (up to ',I2,' ionization stages)')")
     *             nel, atm_int(i).nion(1)
        DO n=1, ndepth
           WRITE(27,"(I4,300F13.6)")
     *                 n,((atm_int(i).pion(n,m,k),
     *                 m=1,atm_int(i).nion(k)), k=1,nel)
        ENDDO

        DO n=1, ndepth
           WRITE(27,"(I4,300F13.6)") n,((atm_int(i).uion(n,m,k),
     *                 m=1,atm_int(i).nion(k)), k=1,nel)
        ENDDO

        WRITE(27,*)  atm_int(i).h.date,'x'
        WRITE(27,*)  atm_int(i).h.user,'x'
        WRITE(27,*)   atm_int(i).h.ref_nr,'x'
          WRITE(27,*) atm_int(i).h.teff,'x'
          WRITE(27,*) atm_int(i).h.logg,'x'
          WRITE(27,*) atm_int(i).h.x,'x'
          WRITE(27,*) atm_int(i).h.y,'x'
          WRITE(27,*) atm_int(i).h.z,'x'
          WRITE(27,*) atm_int(i).h.alpha,'x'
          WRITE(27,*) atm_int(i).h.eps_nr,'x'
          WRITE(27,*) atm_int(i).h.status,'x'
         WRITE(27,*)  atm_int(i).h.cmt_nr,'x'
          WRITE(27,*) atm_int(i).h.ref(:atm_int(i).h.ref_nr),'x'
       
        CLOSE(27)

c       GOTO 991
c
c ------------------------------ MARCS MODELS ------------------
c
                                          !adapted for KL output from MARCS
c DO NOT USE i subscript - at the moment is used to designate atm_int(i)
c
         OPEN(3, file = filein, form = "unformatted")  
         OPEN(8, FILE=fileout,STATUS ='UNKNOWN',FORM='FORMATTED')

c   ndep,tau5,temp,pele,rhog,kref,pgas,abun,inst,atnr,partel,presel,n,nstar 

         READ(3) 
     *     NDEP_MARCS,              ! n depths
     *    (TAU5(n), n=1,ndep_m),
     *    (TEMP(n), n=1,ndep_m),
     *    (PELE(n), n=1,ndep_m),
     *    (RHOG(n), n=1,ndep_m),
     *    (KREF(n), n=1,ndep_m),
     *    (PGAS(n), n=1,ndep_m), 
     *    (ABUN(n), n=1,nel_new_m),    ! abundances for each el
     *    (INST(n), n=1,nel_m),    ! n ion stages for each el
     *    (ATNR(n), n=1,nel_m),    ! atomic numbers 
     *    ((PARTEL(m,n), m=1, nion_m), n=1,ndep_m),
     *    ((PRESEL(m,n), m=1, nion_m), n=1,ndep_m)

         m = 0
         DO i = 1, nel_m
          DO j = 1, INST(i)
          m = m + 1            ! nion_m = 48 (16 elements, each 3 ion stages)
            DO k = 1, ndep_m
              part_m(k,j,i) = PARTEL(m,k)
              pres_m(k,j,i) = PRESEL(m,k)
            ENDDO
          ENDDO
         ENDDO

c add relative eps_dev for elements not used in calculations
c only used in previous MULTI models with abundances for 16 elements specified in the ATNR array
c        m = 1
c         DO i = 1, nel_new_m ! 63 elements = assume array index equal to Z
c            IF (i.eq.ATNR(m)) THEN
c            ab_m(i) =  ABUN(m)
c            m = m + 1        ! m runs from 1 to nel_m (=16)
c            ELSE 
c            ab_m(i)= 0.00
c            ENDIF
c         ENDDO

C         (((PARTEL(n,m,k),  n=1,ndep_m), m=1, INST(k)), k=1,nel_m),
C         (((PRESEL(n,m,k),  n=1,ndep_m), m=1, INST(k)), k=1,nel_m)

        WRITE(8,"('N depth points',I4)") ndep_m
        WRITE(8,"('N of elements',I4)") nel_m
        WRITE(8,"('N of molecules',I4)") nmol_m             !set to 0 - unknown
        WRITE(8,"('N of break points in tau scale',I4)") nbreak_m     !set to 0 - unknown
         WRITE(8,"('Atomic numbers of elements')")
         WRITE(8,712) (ATNR(n), n=1, nel_m)
         WRITE(8,"('Number of ionization ',
     *                     'stages for each element')")
         WRITE(8,712) (INST(n), n=1, nel_m)
         WRITE(8,"('Logarithmic element ',
     *            'abundances, A = log(N_el/N_H) + 12')")
         WRITE(8,"(10(F5.2,2X))") (ABUN(n), n=1,nel_new_m)
         WRITE(8,"('Break poins in tau-scale', 10I4)")
     *             (break(n), n=1, nbreak_m)
         WRITE(8,715) 'Depth','lgTau5','T','Kappa_ref',
     *                 'Pe','Pg','Pk','Rho'
          
c
c  TEMPERATURE IS NOT IN LOG SCALE IN THE ASCII FILE (SIU uses LOG!)
c
        DO n = 1, ndep_m
         WRITE(8,714) n, TAU5(n), 10**(TEMP(n)), KREF(n),
     *                 PELE(n), PGAS(n),
     *                 PGAS(n), RHOG(n)
        ENDDO
         
        WRITE(8,"(/,'Partial pressures ',
     *             'and partition functions for ',
     *           I2,' elements (up to ',I2,' ionization stages)')")
     *             nel_m, ionst(2)                 !for He, for H only 2 ion stages!

        DO n = 1, ndep_m
           WRITE(8,"(I4,300F13.6)") 
     *                       n, ((pres_m(n,m,k),
     *                       m=1, INST(k)), k=1,nel_m)
        ENDDO

        DO n = 1, ndep_m
           WRITE(8,"(I4,300F13.6)") 
     *                       n, ((part_m(n,m,k), 
     *                       m=1, INST(k)), k=1,nel_m)
        ENDDO

c -------------------------------------------------
c
        i = 1 ! assign the number to atm_int(i).xxx

        atm_int(i).h.date = 'December 2010'
        atm_int(i).h.user = 'created from MARCS'
        atm_int(i).h.eps_nr = nel_m
      
        DO j = 1, nel_m                ! change from double to single precision
           at_idx(j) = ATNR(j)
           at_nion(j) = INST(j)
        END DO      
        
        atm_int(i).h.eps_nr = nel_new_m
        DO j = 1, atm_int(i).h.eps_nr
           atm_int(i).h.eps(j) = ABUN(j)
        ENDDO

        DO j = 1, ndep_m                ! change from double to single precision
        pgas_m(j) = PGAS(j)
        END DO


        WRITE(8,*)  'HEADER FROM MAFAGS-OS SUN'
        WRITE(8,*)  atm_int(i).h.date,'x'
        WRITE(8,*)  atm_int(i).h.user,'x'
        WRITE(8,*)   atm_int(i).h.ref_nr,'x'
          WRITE(8,*) atm_int(i).h.teff,'x'
          WRITE(8,*) atm_int(i).h.logg,'x'
          WRITE(8,*) atm_int(i).h.x,'x'
          WRITE(8,*) atm_int(i).h.y,'x'
          WRITE(8,*) atm_int(i).h.z,'x'
          WRITE(8,*) atm_int(i).h.alpha,'x'
          WRITE(8,*) atm_int(i).h.eps_nr,'x'
          WRITE(8,*) atm_int(i).h.status,'x'
         WRITE(8,*)  atm_int(i).h.cmt_nr,'x'
          WRITE(8,*) atm_int(i).h.ref(:atm_int(i).h.ref_nr),'x'
        CLOSE(8)

        OPEN(9, FILE=foutbin,STATUS ='UNKNOWN',FORM='UNFORMATTED')

        WRITE(9)
     *     ndep_m,nel_m,nmol_m,nbreak_m,
     *    (at_idx(n),  n=1,nel_m),
     *    (at_nion(n), n=1,nel_m),
     *    (break(n), n=1,nbreak_m),  ! user-defined 
     *    (TAU5(n),  n=1,ndep_m),
     *    (TEMP(n),  n=1,ndep_m),
     *    (KREF(n),  n=1,ndep_m),
     *    (PELE(n),  n=1,ndep_m),
     *    (pgas_m(n),  n=1,ndep_m),
     *    (pgas_m(n),  n=1,ndep_m),
     *    (RHOG(n),  n=1,ndep_m),
     *   (((pres_m(n,m,k),n=1,ndep_m), m=1, INST(k)),k=1,nel_m),     !nion - parameter - header!
     *   (((part_m(n,m,k),n=1,ndep_m), m=1, INST(k)),k=1,nel_m)
 
        WRITE(9)  atm_int(i).h.date,atm_int(i).h.user,
     *                atm_int(i).h.ref_nr,atm_int(i).h.teff,
     *                atm_int(i).h.logg,atm_int(i).h.x,atm_int(i).h.y,
     *                atm_int(i).h.z,atm_int(i).h.alpha,
     *                atm_int(i).h.eps_nr,atm_int(i).h.status,
     *                atm_int(i).h.cmt_nr,
     *                atm_int(i).h.ref(:atm_int(i).h.ref_nr),
     *               (atm_int(i).h.eps(k),k=1,atm_int(i).h.eps_nr),
     *                atm_int(i).h.cmt(:atm_int(i).h.cmt_nr)


        CLOSE(9) !binary model atmosphere!
 991    CONTINUE
                                     
        CLOSE(2)
        CLOSE(3)

        STOP
        END
