program make_composition implicit none c c Purpose: to create input data file "composition_venus.in" c include 'max.inc' include 'formats.inc' integer m,i,j,nlev_mf,mol,imf,ind,iso logical indf integer band double precision alt(1:Nmax) double precision pres(1:Nmax) double precision adia(1:Nmax) double precision temp(1:Nmax) double precision Tsol double precision Tesp double precision dens_mix(1:Nmax) double precision grav(1:Nmax) double precision zeta(1:Nmax) double precision capa(1:Nmax) double precision P_mf(1:Nmax) double precision x(1:Nmax,1:Nmol_max) double precision nu_lo(1:Nbmx) double precision nu_hi(1:Nbmx) double precision nu_inf,nu_tr,nu_sup double precision dnu_lw,dnu_sw integer Nband,Nband_lw,Nband_sw character*(Nchar_mx) composition_file,model_file,nu_file character*(Nchar_mx) input_file character*(Nchar_mx) iso_comp_file character*(Nchar_mx) planet_descriptor c (P,T) levels double precision DeltaT integer iprange,itrange integer NPrange,NTrange(1:Nrange_mx) double precision Prange(1:Nrange_mx) double precision Trange_limits(1:Nrange_mx,1:2) double precision Trange(1:Nrange_mx,1:Nrange_mx) integer Nlev,ilev double precision Plev(1:Nmax) double precision Tlev(1:Nmax) double precision zlev(1:Nmax) double precision xlev(1:Nmax,1:Nmol_max) double precision x_HDO(1:Nmax) double precision abundance_HDO(1:Nmax) logical hdo_found,h2o_found integer ihdo,ih2o logical imol_found integer imol logical iso161_found,iso162_found integer i161,i162 integer conv_err c HITRAN data integer HNmol integer Hmol_index(1:Nmol_max) integer Hmol_niso(1:Nmol_max) character*10 Hmol_name(1:Nmol_max) integer Hiso_index(1:Nmol_max,1:Niso_max) 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) c abundances of isotopologues integer Nmol integer mol_index(1:Nmol_max) integer mol_niso(1:Nmol_max) character*10 mol_name(1:Nmol_max) integer iso_index(1:Nmol_max,1:Niso_max) double precision iso_abundance(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) c label integer strlen character*(Nchar_mx) label label='program make_composition' planet_descriptor='VENUS 140km' c get molecules names and indexes from HITRAN input_file='./molparam_hitran2008.txt' call read_molparam(input_file, & HNmol,Hmol_index,Hmol_niso,Hmol_name, & Hiso_index,Hiso_abundance,Hiso_Qref,Hiso_g,Hiso_mass) c Prange: range of pressure levels for spectra generation NPrange=14 if (NPrange.gt.Nrange_mx) then call error(label) write(*,*) 'NPrange=',NPrange write(*,*) 'while Nrange_mx=',Nrange_mx stop endif Prange(1)=1.0D+7 ! Pa do iprange=2,NPrange Prange(iprange)=Prange(iprange-1)/10.0D+0 ! Pa enddo ! iprange c Trange: range of temperature levels for spectra generation DeltaT=50.0D+0 ! K c Prange(1)=10^7 Pa Trange_limits(1,1)=450.0D+0 ! K Trange_limits(1,2)=750.0D+0 ! K c Prange(2)=10^6 Pa Trange_limits(2,1)=300.0D+0 ! K Trange_limits(2,2)=750.0D+0 ! K c Prange(3)=10^5 Pa Trange_limits(3,1)=200.0D+0 ! K Trange_limits(3,2)=550.0D+0 ! K c Prange(4)=10^4 Pa Trange_limits(4,1)=150.0D+0 ! K Trange_limits(4,2)=400.0D+0 ! K c Prange(5)=10^3 Pa Trange_limits(5,1)=150.0D+0 ! K Trange_limits(5,2)=300.0D+0 ! K c Prange(6)=10^2 Pa Trange_limits(6,1)=100.0D+0 ! K Trange_limits(6,2)=300.0D+0 ! K c Prange(7)=10^1 Pa Trange_limits(7,1)=100.0D+0 ! K Trange_limits(7,2)=300.0D+0 ! K c Prange(8)=10^0 Pa Trange_limits(8,1)=100.0D+0 ! K Trange_limits(8,2)=300.0D+0 ! K c Prange(9)=10^-1 Pa Trange_limits(9,1)=100.0D+0 ! K Trange_limits(9,2)=300.0D+0 ! K c Prange(10)=10^-2 Pa Trange_limits(10,1)=100.0D+0 ! K Trange_limits(10,2)=300.0D+0 ! K c Prange(11)=10^-3 Pa Trange_limits(11,1)=100.0D+0 ! K Trange_limits(11,2)=300.0D+0 ! K c Prange(12)=10^-4 Pa Trange_limits(12,1)=100.0D+0 ! K Trange_limits(12,2)=300.0D+0 ! K c Prange(13)=10^-5 Pa Trange_limits(13,1)=100.0D+0 ! K Trange_limits(13,2)=300.0D+0 ! K c Prange(14)=10^-6 Pa Trange_limits(14,1)=100.0D+0 ! K Trange_limits(14,2)=300.0D+0 ! K do iprange=1,NPrange NTrange(iprange)=int((Trange_limits(iprange,2) & -Trange_limits(iprange,1))/DeltaT+1) if (NTrange(iprange).gt.Nrange_mx) then call error(label) write(*,*) 'NTrange(',iprange,')=',NTrange(iprange) write(*,*) 'while Nrange_mx=',Nrange_mx stop endif do itrange=1,NTrange(iprange) Trange(iprange,itrange)=Trange_limits(iprange,1) & +(itrange-1)*DeltaT enddo ! itrange c Debug c write(*,*) 'Prange(',iprange,')=',Prange(iprange) c write(*,*) 'Tmin=',Trange_limits(iprange,1), c & ' Tmax=',Trange_limits(iprange,2), c & ' NTrange(',iprange,')=',NTrange(iprange) c write(*,*) (Trange(iprange,itrange), c & itrange=1,NTrange(iprange)) c Debug enddo ! iprange Nlev=0 do iprange=1,NPrange Nlev=Nlev+NTrange(iprange) enddo ! iprange if (Nlev.gt.Nmax) then call error(label) write(*,*) 'Nlev=',Nlev write(*,*) 'while Nmax=',Nmax stop endif ilev=0 do iprange=1,NPrange do itrange=1,NTrange(iprange) ilev=ilev+1 Plev(ilev)=Prange(iprange) ! Pa Tlev(ilev)=Trange(iprange,itrange) ! K zlev(ilev)=0.0D+0 ! m enddo ! itrange enddo ! iprange c Reading physical model of the atmosphere call read_model(m,alt,pres,temp,Tsol,Tesp, & dens_mix,adia,grav,zeta,capa) call read_mixratios(Nmol,nlev_mf,P_mf,x,mol_name) c 2013-10-09 c Identify HDO molecule index hdo_found=.false. do mol=1,Nmol if (mol_name(mol).eq.'HDO') then hdo_found=.true. ihdo=mol goto 112 endif enddo ! mol 112 continue if (.not.hdo_found) then call error(label) write(*,*) 'HDO could not be identified among:' do mol=1,Nmol write(*,*) mol,mol_name(mol)(1:strlen(mol_name(mol))) enddo ! mol endif c Identify H2O molecule index h2o_found=.false. do mol=1,Nmol if (mol_name(mol).eq.'H2O') then h2o_found=.true. ih2o=mol goto 113 endif enddo ! mol 113 continue if (.not.h2o_found) then call error(label) write(*,*) 'H2O could not be identified among:' do mol=1,Nmol write(*,*) mol,mol_name(mol)(1:strlen(mol_name(mol))) enddo ! mol endif c Interpolation of species concentrations for every (P,T) level do ilev=1,Nlev indf=.false. do imf=1,nlev_mf-1 if ((Plev(ilev).gt.P_mf(imf)).and. & (Plev(ilev).lt.P_mf(imf+1))) then indf=.true. ind=imf goto 321 endif enddo ! imf 321 continue c if (indf) then do mol=1,Nmol xlev(ilev,mol)=x(ind,mol) & +(x(ind+1,mol)-x(ind,mol)) & /(dlog(P_mf(ind+1))-dlog(P_mf(ind))) & *(dlog(Plev(ilev))-dlog(P_mf(ind))) enddo ! mol else ! indf=.false. if (Plev(ilev).gt.P_mf(nlev_mf)) then do mol=1,Nmol xlev(ilev,mol)=x(nlev_mf,mol) enddo endif if (Plev(ilev).lt.P_mf(1)) then do mol=1,Nmol xlev(ilev,mol)=x(1,mol) enddo ! mol endif ! Plev(ilev)