c kspectrum (http://www.meso-star.com/en_Products.html) - This file is part of kspectrum c Copyright (C) 2008-2015 - Méso-Star - Vincent Eymet c c This file must be used under the terms of the CeCILL license. c This source file is licensed as described in the file COPYING, which c you should have received as part of this distribution. The terms c are also available at c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt c subroutine read_composition(nis,Hnmol,Hmol_index,Hmol_niso, & Hiso_index,Hmol_name,Hiso_abundance,Hiso_Qref, & Hiso_g,Hiso_mass, & planet_descriptor,m,alt,pres,temp, & nmol,mol_index,mol_niso,iso_index, & mol_name,abundances,iso_Qref,iso_g,iso_mass,x) implicit none c c Purpose: to read the atmospheric composition, pressure and temperature profiles c c Input: c + nis: option for selecting isotopes when not specified c + Hnmol: number of molecules used in the HITRAN04.par database c + Hmol_index(mol) is the index of the molecule (index 'mol') used in the HITRAN04.par database c + Hmol_niso(mol) is the number of isotopes for HITRAN molecule index 'mol' c + Hiso_index(mol,iso) is the index (used in the "molparam.in" file) of the isotope index 'iso' for molecule index 'mol' c + Hmol_name(mol) is a character string that describes HITRAN molecule 'mol' c + Hiso_abundance(mol,iso) is the abundance of isotope 'iso' for each HITRAN molecule 'mol' c + Hiso_Qref(mol,iso) are reference (@ T=Tref) partition functions c + Hiso_g(mol,iso) are state independent degeneracy factors c + Hiso_mass(mol,iso) are molar masses [g/mol] c c Output: c + planet_descriptor: character string that identifies the atmosphere c + m: number of atmospheric levels c + alt: altitude at each level (m) c + pres: pressure at each level (atm) c + temp: temperature at each level (K) c + nmol: number of molecules radiatively active c + mol_index: indexes (using HITRAN convention) of the nmol molecules c + mol_niso: number of isotopes for each molecule c + iso_index: indexes (using HITRAN convention) of isotopes c + mol_name: name of each molecule c + abundances(i,mol,iso) is the abundance of isotope 'iso' for molecule 'mol', at level 'i' c + iso_Qref(mol,iso) are reference (@ T=Tref) partition functions c + iso_g(mol,iso) are state independent degeneracy factors c + iso_mass(mol,iso) are molar masses [g/mol] c + x: molar fraction (mixratio) of each molecule, at each atmospheric level c include 'max.inc' include 'formats.inc' include 'descriptors.inc' integer nis integer Hnmol,nsmol,nsmol_identified,sm,im integer mol,iso,nmol,ios,m,i,c1,cstart,csf integer hdof,rnmol,mf,hm,tm,itm integer amolf,amol_ind,Hamol_ind,iso_ind,imolf integer isof,sm1,add integer nsm(1:Nmol_max) integer ind_mol_identified(1:Nmol_max,1:4) integer Hmol_index(1:Nmol_max) integer Hmol_niso(1:Nmol_max) integer spec(1:Nmol_max),mt(1:Nmol_max) integer Hiso_index(1:Nmol_max,1:Niso_max) integer mol_index(1:Nmol_max) integer mol_niso(1:Nmol_max) integer iso_index(1:Nmol_max,1:Niso_max) integer shift double precision alt(1:Nmax) double precision pres(1:Nmax) double precision temp(1:Nmax) double precision Hiso_abundance(1:Nmol_max,1:Niso_max) double precision Hiso_Qref(1:Nmol_max,1:Niso_max) integer Hiso_g(1:Nmol_max,1:Niso_max) double precision Hiso_mass(1:Nmol_max,1:Niso_max) double precision abundances(1:Nmax,1:Nmol_max,1:Niso_max) double precision iso_Qref(1:Nmol_max,1:Niso_max) integer iso_g(1:Nmol_max,1:Niso_max) double precision iso_mass(1:Nmol_max,1:Niso_max) double precision x(1:Nmax,1:Nmol_max) double precision smol_mass(1:Nmol_max) double precision sum_x,sum_others,sum_Hisoab double precision sum_ab character*(Nchar_mx) input_file character*(Nchar_mx) string character*(Nchar_mx) file_composition,planet_descriptor character*(Nchar_mx) file_info character*(Nchar_mx) strtmp1 character*10 smol_name(1:Nmol_max,1:2) character*10 Hmol_name(1:Nmol_max) character*10 mol_name(1:Nmol_max) character*10 cmol,hmol,smol,smol1,amol,imol character*10 tmol character*1 char logical mpf_ex c t* variables: read from "molparam.in" integer tNlev integer tNmol integer tmol_index(1:Nmol_max) integer tmol_niso(1:Nmol_max) character*10 tmol_name(1:Nmol_max) integer tiso_index(1:Nmol_max,1:Niso_max) double precision tiso_abundance(1:Nmax,1:Nmol_max,1:Niso_max) double precision tiso_Qref(1:Nmol_max,1:Niso_max) integer tiso_g(1:Nmol_max,1:Niso_max) double precision tiso_mass(1:Nmol_max,1:Niso_max) c label integer strlen character*(Nchar_mx) label label='subroutine read_composition' input_file='./data/molparam.in' inquire(file=input_file(1:strlen(input_file)),exist=mpf_ex) if (mpf_ex) then call read_isotopic_composition(input_file, & tNlev,tNmol,tmol_name,tmol_index,tmol_niso,tiso_index, & tiso_abundance,tiso_Qref,tiso_g,tiso_mass) else tNmol=0 endif file_composition='./data/composition.in' open(12,file=file_composition(1:strlen(file_composition)), & status='old',iostat=ios) if (ios.ne.0) then ! file not found call error(label) write(*,*) 'Data file could not be found:' write(*,*) file_composition(1:strlen(file_composition)) stop else write(*,*) 'Reading atmospheric composition from file:' write(*,*) file_composition(1:strlen(file_composition)) endif read(12,'(a)') string csf=0 do c1=1,strlen(string)-1 char=string(c1:c1+1) if (char.eq.':') then csf=1 cstart=c1+2 goto 123 endif enddo 123 continue if (csf.eq.0) then call error(label) write(*,*) 'Planet descritor could not be found in file:' write(*,*) file_composition(1:strlen(file_composition)) stop endif planet_descriptor=string(cstart:strlen(string)) read(12,*) read(12,*) m read(12,*) read(12,*) nmol read(12,*) read(12,'(a)') string call identify_molecules(string,nmol,mol_name) do i=1,m read(12,50) alt(i),pres(i),temp(i),(x(i,mol),mol=1,nmol) enddo close(12) if ((mpf_ex).and.(tNlev.ne.m)) then call error(label) write(*,*) 'Number of levels in composition file:',m write(*,*) 'Number of levels in isotopic composition file:' & ,tNlev write(*,*) 'should be the same' stop endif c identify corresponding molecules in the list of HITRAN molecules hdof=0 rnmol=0 nsmol_identified=0 do mol=1,nmol nsm(mol)=0 enddo do mol=1,nmol mf=0 cmol=mol_name(mol) c look for molecule "mol" in the list of HITRAN molecules do hm=1,Hnmol hmol=Hmol_name(hm) if (cmol(1:strlen(cmol)).eq.hmol(1:strlen(hmol))) then mf=1 mol_index(mol)=Hmol_index(hm) c spec(mol)=0 ! molecule 'mol' is not a special molecule goto 124 endif enddo ! hm 124 continue c If identification of molecule "mol" failed if (mf.eq.0) then call error(label) write(*,*) 'molecule: "',cmol(1:strlen(cmol)), & '" could not be identified' write(*,*) 'within the list of HITRAN molecules' stop endif enddo ! mol c Set number of isotopes mol_niso and indexes of isotopes iso_index do mol=1,nmol c each molecule index "mol" must be identified within the "molparam.in" file mf=0 cmol=mol_name(mol) do tm=1,tNmol tmol=tmol_name(tm) if (cmol(1:strlen(cmol)).eq.tmol(1:strlen(tmol))) then mf=1 itm=tm goto 125 endif enddo ! hm 125 continue if (mf.eq.1) then ! If molecule "mol" was identified in "molparam.in" mol_niso(mol)=0 do iso=1,tmol_niso(itm) mol_niso(mol)=mol_niso(mol)+1 iso_index(mol,mol_niso(mol))=iso do i=1,m abundances(i,mol,mol_niso(mol))= & tiso_abundance(i,itm,iso) enddo ! i iso_Qref(mol,mol_niso(mol))=tiso_Qref(itm,iso) iso_g(mol,mol_niso(mol))=tiso_g(itm,iso) iso_mass(mol,mol_niso(mol))=tiso_mass(itm,iso) enddo ! iso else if (mf.eq.0) then ! If identification of molecule "mol" failed if (nis.eq.1) then ! use only main isotopologue mol_niso(mol)=1 iso_index(mol,1)=1 do i=1,m abundances(i,mol,1)=1.0D+0 enddo ! i iso_Qref(mol,1)=Hiso_Qref(mol_index(mol), & iso_index(mol,1)) iso_g(mol,1)=Hiso_g(mol_index(mol), & iso_index(mol,1)) iso_mass(mol,1)=Hiso_mass(mol_index(mol), & iso_index(mol,1)) else ! use all known isotopologues mol_niso(mol)=Hmol_niso(mol_index(mol)) do iso=1,mol_niso(mol) iso_index(mol,iso)=Hiso_index(mol_index(mol) & ,iso) do i=1,m abundances(i,mol,iso)= & Hiso_abundance(mol_index(mol), & iso_index(mol,iso)) enddo ! i iso_Qref(mol,iso)=Hiso_Qref(mol_index(mol), & iso_index(mol,iso)) iso_g(mol,iso)=Hiso_g(mol_index(mol), & iso_index(mol,iso)) iso_mass(mol,iso)=Hiso_mass(mol_index(mol), & iso_index(mol,iso)) enddo ! iso endif ! nis else ! mf!=0 and mf!= 1 call error(label) write(*,*) 'mf=',mf write(*,*) 'should be equal to 0 or 1' stop endif ! mf=0 enddo ! mol c Debug file_info='./results/composition_info.txt' open(21,file=file_info(1:strlen(file_info))) write(21,*) 'Number of molecules identified:',nmol do mol=1,nmol write(21,*) 'molecule',mol,' /',nmol & ,'; ',mol_name(mol) & ,' (',mol_index(mol),')' write(21,*) 'Number of isotopes:',mol_niso(mol) do iso=1,mol_niso(mol) write(21,*) 'isotope',iso,' abundances(',mol,',', & iso,')=',abundances(1,mol,iso), & ' iso_index=',iso_index(mol,iso), & ' iso_mass=',iso_mass(mol,iso) enddo if (mol_niso(mol).gt.1) then sum_x=0.0D+0 do iso=1,mol_niso(mol) sum_x=sum_x+ & abundances(1,mol,iso) enddo write(21,*) 'sum of abundances=',sum_x endif write(21,*) enddo close(21) write(*,*) 'Composition information file was generated:' write(*,*) file_info(1:strlen(file_info)) c Debug return end subroutine find_hm(mol,Hnmol,Hmol_index,mol_index,hm) implicit none include 'max.inc' integer mol,Hnmol integer Hmol_index(1:Nmol_max) integer mol_index(1:Nmol_max) integer hm integer molec,hmf hmf=0 do molec=1,Hnmol if (mol_index(mol).eq.Hmol_index(molec)) then hm=molec hmf=1 goto 123 endif enddo ! hm 123 continue if (hmf.eq.0) then write(*,*) 'Error from routine find_hm:' write(*,*) 'hm could not be found for:' write(*,*) 'mol=',mol write(*,*) 'mol_index(',mol,')=',mol_index(mol) write(*,*) 'Hnmol=',Hnmol do molec=1,Hnmol write(*,*) 'Hmol_index(',molec,')=',Hmol_index(molec) enddo stop endif return end subroutine find_hiso(mol,iso,hm,Hmol_niso,Hiso_index,iso_index, & hiso) implicit none include 'max.inc' integer mol,iso,hm integer Hmol_niso(1:Nmol_max) integer Hiso_index(1:Nmol_max,1:Niso_max) integer iso_index(1:Nmol_max,1:Niso_max) integer hiso integer isotop,hisof hisof=0 do isotop=1,Hmol_niso(hm) if (iso_index(mol,iso).eq.Hiso_index(hm,isotop)) then hiso=isotop hisof=1 goto 123 endif enddo ! hiso 123 continue if (hisof.eq.0) then write(*,*) 'Error from routine find_hiso:' write(*,*) 'hiso could not be found for:' write(*,*) 'mol=',mol write(*,*) 'iso=',iso write(*,*) 'iso_index(',mol,',',iso,')=',iso_index(mol,iso) write(*,*) 'Hmol_niso(',hm,')=',Hmol_niso(hm) do isotop=1,Hmol_niso(hm) write(*,*) 'Hiso_index(',hm,',',isotop,')=', & Hiso_index(hm,isotop) enddo stop endif return end