| 1 | program make_composition_gen |
|---|
| 2 | implicit none |
|---|
| 3 | |
|---|
| 4 | ! ------------------------------------------------------------- |
|---|
| 5 | ! Purpose: to create "composition.in" file readable by kspectrum |
|---|
| 6 | ! Authour: Robin Wordsworth (2010) |
|---|
| 7 | ! ------------------------------------------------------------- |
|---|
| 8 | |
|---|
| 9 | include 'max.inc' |
|---|
| 10 | include 'formats.inc' |
|---|
| 11 | integer strlen,nmolec,nlev_mf,mol,imf,ind,indf |
|---|
| 12 | integer i,j,k,im,m |
|---|
| 13 | integer nb,na,ne |
|---|
| 14 | double precision alt(1:Nmax) |
|---|
| 15 | double precision pres(1:Nmax) |
|---|
| 16 | double precision temp(1:Nmax) |
|---|
| 17 | |
|---|
| 18 | double precision P_mf(1:Nmax) |
|---|
| 19 | double precision x(1:Nmax,1:Nmol_max) |
|---|
| 20 | double precision xlev(1:Nmol_max) |
|---|
| 21 | double precision xconst(1:Nmol_max) |
|---|
| 22 | character*100 composition_file,model_file |
|---|
| 23 | character*100 planet_descriptor |
|---|
| 24 | character*200 string |
|---|
| 25 | character*20 s |
|---|
| 26 | character*10 molec_names(1:Nmol_max),molname |
|---|
| 27 | character*1 spch |
|---|
| 28 | |
|---|
| 29 | integer iP, iT, iQ |
|---|
| 30 | double precision P(1:Nmax), T(1:Nmax), Q(1:Nmax) |
|---|
| 31 | |
|---|
| 32 | spch=' ' |
|---|
| 33 | |
|---|
| 34 | print*,'Name of atmosphere / planet:' |
|---|
| 35 | read*,planet_descriptor |
|---|
| 36 | |
|---|
| 37 | composition_file='./composition.in' |
|---|
| 38 | |
|---|
| 39 | ! ! nmolec=2 ! total number of molecules |
|---|
| 40 | |
|---|
| 41 | ! load temperature, pressure, mixing ratio data |
|---|
| 42 | open(9,file='T.dat') |
|---|
| 43 | read(9,*) iT |
|---|
| 44 | do i=1,iT |
|---|
| 45 | read(9,*) T(i) |
|---|
| 46 | enddo |
|---|
| 47 | close(9) |
|---|
| 48 | |
|---|
| 49 | open(10,file='p.dat') |
|---|
| 50 | read(10,*) iP |
|---|
| 51 | do i=1,iP |
|---|
| 52 | read(10,*) P(i) |
|---|
| 53 | P(i)=10.**P(i)/1013.25 |
|---|
| 54 | enddo |
|---|
| 55 | close(10) |
|---|
| 56 | |
|---|
| 57 | open(11,file='Q.dat') |
|---|
| 58 | read(11,*) nmolec |
|---|
| 59 | |
|---|
| 60 | print*,'nmolec=',nmolec |
|---|
| 61 | |
|---|
| 62 | do mol=1,nmolec |
|---|
| 63 | read(11,*) molec_names(mol) |
|---|
| 64 | enddo |
|---|
| 65 | read(11,*) iQ |
|---|
| 66 | do i=1,iQ |
|---|
| 67 | read(11,*) Q(i) |
|---|
| 68 | enddo |
|---|
| 69 | close(11) |
|---|
| 70 | |
|---|
| 71 | m=iP*iT*iQ ! total number of layers |
|---|
| 72 | |
|---|
| 73 | print*,'Temperature layers: ',iT |
|---|
| 74 | print*,'Pressure layers: ',iP |
|---|
| 75 | print*,'Mixing ratio layers: ',iQ |
|---|
| 76 | print*,'Total: ',m |
|---|
| 77 | |
|---|
| 78 | do mol=1,nmolec-1 |
|---|
| 79 | print*,'Please enter vmr of ',molec_names(mol) |
|---|
| 80 | read(*,*) xconst(mol) |
|---|
| 81 | end do |
|---|
| 82 | |
|---|
| 83 | |
|---|
| 84 | do i=1,iQ |
|---|
| 85 | do j=1,iP |
|---|
| 86 | do k=1,iT |
|---|
| 87 | |
|---|
| 88 | im = (i-1)*iP*iT + (j-1)*iT + k |
|---|
| 89 | |
|---|
| 90 | alt(im)=0.D0 |
|---|
| 91 | pres(im)=P(j) |
|---|
| 92 | temp(im)=T(k) |
|---|
| 93 | |
|---|
| 94 | do mol=1,nmolec-1 |
|---|
| 95 | x(im,mol)=xconst(mol)*(1.0-Q(i)) |
|---|
| 96 | enddo |
|---|
| 97 | x(im,nmolec)=Q(i) |
|---|
| 98 | |
|---|
| 99 | enddo |
|---|
| 100 | enddo |
|---|
| 101 | enddo |
|---|
| 102 | |
|---|
| 103 | open(12,file=composition_file(1:strlen(composition_file))) |
|---|
| 104 | string='Atmospheric composition input ' |
|---|
| 105 | string=string(1:strlen(string))//' data file for planet:' |
|---|
| 106 | string=string(1:strlen(string))//' '// |
|---|
| 107 | |
|---|
| 108 | & planet_descriptor(1:strlen(planet_descriptor)) |
|---|
| 109 | write(12,'(a)') string(1:strlen(string)) |
|---|
| 110 | write(12,39) 'Number of atmospheric levels: ',m |
|---|
| 111 | write(12,34) 'Number of molecules: ',nmolec |
|---|
| 112 | write(12,*) |
|---|
| 113 | string=' z (km) / P (atm) / T (K) *' |
|---|
| 114 | |
|---|
| 115 | |
|---|
| 116 | do mol=1,nmolec |
|---|
| 117 | |
|---|
| 118 | molname=molec_names(mol) |
|---|
| 119 | |
|---|
| 120 | nb=5 |
|---|
| 121 | if (strlen(molname).eq.1) then |
|---|
| 122 | na=7 |
|---|
| 123 | else if (strlen(molname).eq.2) then |
|---|
| 124 | na=6 |
|---|
| 125 | else if (strlen(molname).eq.3) then |
|---|
| 126 | na=5 |
|---|
| 127 | endif |
|---|
| 128 | s='*' |
|---|
| 129 | do ne=1,nb |
|---|
| 130 | s=s(1:strlen(s)-1)//spch//'*' |
|---|
| 131 | enddo |
|---|
| 132 | s=s(1:strlen(s)-1)//'/*' |
|---|
| 133 | do ne=1,na |
|---|
| 134 | s=s(1:strlen(s)-1)//spch//'*' |
|---|
| 135 | enddo |
|---|
| 136 | |
|---|
| 137 | string=string(1:strlen(string)-1) |
|---|
| 138 | & //s(1:strlen(s)-1) |
|---|
| 139 | & //'x[' |
|---|
| 140 | & //molname(1:strlen(molname)) |
|---|
| 141 | & //']*' |
|---|
| 142 | enddo |
|---|
| 143 | write(12,*) string(1:strlen(string)-1) |
|---|
| 144 | |
|---|
| 145 | do im=1,m |
|---|
| 146 | |
|---|
| 147 | do mol=1,nmolec |
|---|
| 148 | xlev(mol)=x(im,mol) |
|---|
| 149 | enddo ! mol |
|---|
| 150 | |
|---|
| 151 | write(12,50) alt(im),pres(im),temp(im),(xlev(mol),mol=1,nmolec) |
|---|
| 152 | |
|---|
| 153 | enddo |
|---|
| 154 | |
|---|
| 155 | write(12,*) |
|---|
| 156 | close(12) |
|---|
| 157 | |
|---|
| 158 | write(*,*) 'Output file successfully generated:' |
|---|
| 159 | write(*,*) composition_file(1:strlen(composition_file)) |
|---|
| 160 | |
|---|
| 161 | end |
|---|
