| 1 | c************************************************************************** |
|---|
| 2 | c |
|---|
| 3 | subroutine nltecool(ngrid,nlayer,nq,pplay,pt,pq,dtnlte) |
|---|
| 4 | c |
|---|
| 5 | c This code was designed as a delivery for the "Martian Environment Models" |
|---|
| 6 | c project ( ESA contract 11369/95/nl/jg CCN2 ) |
|---|
| 7 | c Computes non-LTE heating rates from CO2 emission at 15 um |
|---|
| 8 | c in the Martian upper atmosphere. |
|---|
| 9 | c Uses a simplified model consisting of two excited levels with two |
|---|
| 10 | c emission bands, one of them stronger than the other, which correspond |
|---|
| 11 | c to the behaviours of the 626 fundamental band and the isotopic fund.bands. |
|---|
| 12 | c It uses a cool-to-space approximation with tabulated escape functions. |
|---|
| 13 | c These escape functions have been precomputed for the strong and weak bands, |
|---|
| 14 | c and are given as a function of pressure in separate files. |
|---|
| 15 | c The output values are the heating rates (actually, cooling, since they |
|---|
| 16 | c are always negative) for the two bands, i.e., the total cooling is the |
|---|
| 17 | c sum of them. |
|---|
| 18 | c Miguel A. Lopez Valverde |
|---|
| 19 | c Instituto de Astrofisica de Andalucia (CSIC), Granada, Spain |
|---|
| 20 | c |
|---|
| 21 | c Version 1b. See description above. 22-March-2000. |
|---|
| 22 | c Adapted as a subroutine for use in GCM -- PLR/SRL 6/2000 |
|---|
| 23 | c Version 1c. Inclusion of VMR in the tabulation of escape functions. |
|---|
| 24 | c Table contains now only 1 input file -- Miguel 11/Jul/2000 |
|---|
| 25 | c Version 1d data contained in original input file "nlte_escape.dat" |
|---|
| 26 | c now stored in include file "nltedata.h" Y.Wanherdrick 09/2000 |
|---|
| 27 | |
|---|
| 28 | c jul 2011 fgg Modified to allow variable O |
|---|
| 29 | c |
|---|
| 30 | c*************************************************************************** |
|---|
| 31 | |
|---|
| 32 | implicit none |
|---|
| 33 | |
|---|
| 34 | #include "nltedata.h" ! (Equivalent to the reading of the "nlte_escape.dat" file) |
|---|
| 35 | #include "dimensions.h" |
|---|
| 36 | #include "dimphys.h" |
|---|
| 37 | #include "chimiedata.h" |
|---|
| 38 | #include "conc.h" !Added to have "dynamic composition" in the scheme |
|---|
| 39 | #include "tracer.h" !" |
|---|
| 40 | #include "callkeys.h" |
|---|
| 41 | |
|---|
| 42 | c Input and output variables |
|---|
| 43 | c |
|---|
| 44 | integer ngrid ! no. of horiz. gridpoints |
|---|
| 45 | integer nlayer ! no. of atmospheric layers |
|---|
| 46 | integer nq ! no. of tracers |
|---|
| 47 | real pplay(ngrid,nlayer) ! input pressure grid |
|---|
| 48 | real pt(ngrid,nlayer) ! input temperatures |
|---|
| 49 | real pq(ngrid,nlayer,nq) ! input mmrs |
|---|
| 50 | real dtnlte(ngrid,nlayer) ! output temp. tendencies |
|---|
| 51 | |
|---|
| 52 | c |
|---|
| 53 | c Standard atmosphere variables |
|---|
| 54 | c |
|---|
| 55 | real nt ! number density [cm-3] |
|---|
| 56 | real co2(nlayer) ! " of CO2 |
|---|
| 57 | real o3p(nlayer) ! " of atomic oxygen |
|---|
| 58 | real n2co(nlayer) ! " of N2 + CO |
|---|
| 59 | real pyy(nlayer) ! auxiliary pressure grid |
|---|
| 60 | |
|---|
| 61 | c |
|---|
| 62 | c Vectors and indexes for the tabulation of escape functions and VMR |
|---|
| 63 | c |
|---|
| 64 | c np ! # data points in tabulation |
|---|
| 65 | c pnb(np) ! Pressure in tabulation |
|---|
| 66 | c ef1(np) ! Esc.funct.#1, tabulated |
|---|
| 67 | c ef2(np) ! Esc.funct.#2, tabulated |
|---|
| 68 | c co2vmr(np) ! CO2 VMR tabulated |
|---|
| 69 | c o3pvmr(np) ! CO2 VMR tabulated |
|---|
| 70 | c n2covmr(np) ! N2+CO VMR tabulated |
|---|
| 71 | real escf1(nlayer) ! Esc.funct.#1, interpolated |
|---|
| 72 | real escf2(nlayer) ! Esc.funct.#2, interpolated |
|---|
| 73 | |
|---|
| 74 | |
|---|
| 75 | c |
|---|
| 76 | c Local Constants |
|---|
| 77 | c |
|---|
| 78 | real nu1, nu2 ! freq. of energy levels |
|---|
| 79 | real imr1, imr2 ! isotopic abundances |
|---|
| 80 | real hplanck, gamma, vlight ! physical constants |
|---|
| 81 | real ee |
|---|
| 82 | real rfvt ! collisional rate |
|---|
| 83 | real rfvto3p ! " |
|---|
| 84 | real rfvv ! " |
|---|
| 85 | |
|---|
| 86 | c |
|---|
| 87 | c Local variables for the main loop |
|---|
| 88 | c |
|---|
| 89 | real n1, n2, co2t ! ground populations |
|---|
| 90 | real l1, p1, p12 ! prod & losses |
|---|
| 91 | real l2, p2, p21 |
|---|
| 92 | real tt ! dummy variable |
|---|
| 93 | real c1, c2 ! molecular constants |
|---|
| 94 | real ae1, ae2 ! einstein spontaneous emission |
|---|
| 95 | real a1, a2, a12, a21 |
|---|
| 96 | real pl1, pl2 |
|---|
| 97 | real el1, el2 |
|---|
| 98 | real hr1, hr2 ! heating rate due to each band |
|---|
| 99 | real hr(nlayer) ! total heating rate |
|---|
| 100 | |
|---|
| 101 | c |
|---|
| 102 | c Indexes |
|---|
| 103 | c |
|---|
| 104 | integer i |
|---|
| 105 | integer j,ii |
|---|
| 106 | |
|---|
| 107 | c |
|---|
| 108 | c Rate coefficients |
|---|
| 109 | c |
|---|
| 110 | real k19xca, k19xcb |
|---|
| 111 | real k19cap1, k19cap2 |
|---|
| 112 | real k19cbp1, k19cbp2 |
|---|
| 113 | real d19c, d19cp1, d19cp2 |
|---|
| 114 | real k20xc, k20cp1, k20cp2 |
|---|
| 115 | real k21xc, k21cp2 |
|---|
| 116 | |
|---|
| 117 | logical firstcall |
|---|
| 118 | data firstcall/.true./ |
|---|
| 119 | save firstcall,ef1,ef2,co2vmr,n2covmr,o3pvmr,pnb |
|---|
| 120 | |
|---|
| 121 | c |
|---|
| 122 | c Data |
|---|
| 123 | c |
|---|
| 124 | data nu1, nu2, hplanck, gamma, vlight, ee/ |
|---|
| 125 | 1 667.38, 662.3734, 6.6261e-27, 1.191e-5, 3.e10, 1.438769/ |
|---|
| 126 | |
|---|
| 127 | c************************************************************************* |
|---|
| 128 | c PROGRAM STARTS |
|---|
| 129 | c************************************************************************* |
|---|
| 130 | |
|---|
| 131 | imr1 = 0.987 |
|---|
| 132 | imr2 = 0.00408 + 0.0112 |
|---|
| 133 | rfvt = 0.1 |
|---|
| 134 | rfvto3p = 1.0 |
|---|
| 135 | rfvv = 0.1 |
|---|
| 136 | |
|---|
| 137 | if(firstcall) then |
|---|
| 138 | |
|---|
| 139 | do i=1,np |
|---|
| 140 | pnb(i)=1.0e-4*exp(pnb(i)) ! p into Pa |
|---|
| 141 | end do |
|---|
| 142 | |
|---|
| 143 | firstcall = .false. |
|---|
| 144 | |
|---|
| 145 | endif |
|---|
| 146 | |
|---|
| 147 | c |
|---|
| 148 | c MAIN LOOP, for each gridpoint and altitude: |
|---|
| 149 | c |
|---|
| 150 | do j=1,ngrid ! loop over grid points |
|---|
| 151 | c |
|---|
| 152 | c set up local pressure grid |
|---|
| 153 | c |
|---|
| 154 | do ii=1,nlayer |
|---|
| 155 | pyy(ii)=pplay(j,ii) |
|---|
| 156 | enddo |
|---|
| 157 | ! |
|---|
| 158 | ! Interpolate escape functions and VMR to the desired grid |
|---|
| 159 | ! |
|---|
| 160 | call interp1(escf2,pyy,nlayer,ef2,pnb,np) |
|---|
| 161 | call interp1(escf1,pyy,nlayer,ef1,pnb,np) |
|---|
| 162 | if(nltemodel.eq.0) then |
|---|
| 163 | call interp3(co2,o3p,n2co,pyy,nlayer, |
|---|
| 164 | & co2vmr,o3pvmr,n2covmr,pnb,np) |
|---|
| 165 | endif |
|---|
| 166 | |
|---|
| 167 | do i=1,nlayer ! loop over layers |
|---|
| 168 | C |
|---|
| 169 | C test if p lies outside range (p > 3.5 Pa) |
|---|
| 170 | C changed to 1 Pa since transition will always be higher than this |
|---|
| 171 | C |
|---|
| 172 | if(pyy(i) .gt. 1.0 .or. pyy(i) .lt. 4.0e-6) then |
|---|
| 173 | hr(i)=0.0 |
|---|
| 174 | dtnlte(j,i)=0.0 |
|---|
| 175 | else |
|---|
| 176 | c |
|---|
| 177 | c if(pt(j,i).lt.1.0)print*,pt(j,i) |
|---|
| 178 | nt = pyy(i)/(1.381e-17*pt(j,i)) ! nt in cm-3 |
|---|
| 179 | !Dynamic composition |
|---|
| 180 | if(nltemodel.eq.1) then |
|---|
| 181 | co2(i)=pq(j,i,igcm_co2)*mmean(j,i)/mmol(igcm_co2) |
|---|
| 182 | o3p(i)=pq(j,i,igcm_o)*mmean(j,i)/mmol(igcm_o) |
|---|
| 183 | n2co(i)=pq(j,i,igcm_co)*mmean(j,i)/mmol(igcm_co) + |
|---|
| 184 | $ pq(j,i,igcm_n2)*mmean(j,i)/mmol(igcm_n2) |
|---|
| 185 | endif |
|---|
| 186 | |
|---|
| 187 | !Mixing ratio to density |
|---|
| 188 | co2(i)=co2(i)*nt ! CO2 density in cm-3 |
|---|
| 189 | o3p(i)=o3p(i)*nt ! O3p density in cm-3 |
|---|
| 190 | n2co(i)=n2co(i)*nt ! N2+CO in cm-3 |
|---|
| 191 | c molecular populations |
|---|
| 192 | n1 = co2(i) * imr1 |
|---|
| 193 | n2 = co2(i) * imr2 |
|---|
| 194 | co2t = n1 + n2 |
|---|
| 195 | |
|---|
| 196 | c intermediate collisional rates |
|---|
| 197 | tt = pt(j,i)*pt(j,i) |
|---|
| 198 | |
|---|
| 199 | if (pt(j,i).le.175.0) then |
|---|
| 200 | k19xca = 3.3e-15 |
|---|
| 201 | k19xcb = 7.6e-16 |
|---|
| 202 | else |
|---|
| 203 | k19xca = 4.2e-12 * exp( -2988.0/pt(j,i) + 303930.0/tt) |
|---|
| 204 | k19xcb = 2.1e-12 * exp( -2659.0/pt(j,i) + 223052.0/tt) |
|---|
| 205 | endif |
|---|
| 206 | k19xca = k19xca * rfvt |
|---|
| 207 | k19xcb = k19xcb * rfvt |
|---|
| 208 | k19cap1 = k19xca * 2.0 * exp( -ee*nu1/pt(j,i) ) |
|---|
| 209 | k19cap2 = k19xca * 2.0 * exp( -ee*nu2/pt(j,i) ) |
|---|
| 210 | k19cbp1 = k19xcb * 2.0 * exp( -ee*nu1/pt(j,i) ) |
|---|
| 211 | k19cbp2 = k19xcb * 2.0 * exp( -ee*nu2/pt(j,i) ) |
|---|
| 212 | d19c = k19xca*co2t + k19xcb*n2co(i) |
|---|
| 213 | d19cp1 = k19cap1*co2t + k19cbp1*n2co(i) |
|---|
| 214 | d19cp2 = k19cap2*co2t + k19cbp2*n2co(i) |
|---|
| 215 | ! |
|---|
| 216 | k20xc = 3.e-12 * rfvto3p |
|---|
| 217 | k20cp1 = k20xc * 2.0 * exp( -ee/pt(j,i) * nu1 ) |
|---|
| 218 | k20cp2 = k20xc * 2.0 * exp( -ee/pt(j,i) * nu2 ) |
|---|
| 219 | ! |
|---|
| 220 | k21xc = 2.49e-11 * 0.5 * rfvv |
|---|
| 221 | k21cp2 = k21xc * exp( - ee/pt(j,i) * (nu2-nu1) ) |
|---|
| 222 | ! |
|---|
| 223 | l1 = d19c + k20xc*o3p(i) + k21cp2*n2 |
|---|
| 224 | p1 = ( d19cp1 + k20cp1*o3p(i) ) * n1 |
|---|
| 225 | p12 = k21xc*n1 |
|---|
| 226 | ! |
|---|
| 227 | l2 = d19c + k20xc*o3p(i) + k21xc*n1 |
|---|
| 228 | p2 = ( d19cp2 + k20cp2*o3p(i) ) * n2 |
|---|
| 229 | p21 = k21cp2*n2 |
|---|
| 230 | |
|---|
| 231 | c radiative rates |
|---|
| 232 | ae1 = 1.3546 * 1.66 / 4.0 * escf1(i) |
|---|
| 233 | ae2 = ( 1.3452 + 1.1878 ) * 1.66 / 4.0 * escf2(i) |
|---|
| 234 | l1 = l1 + ae1 |
|---|
| 235 | l2 = l2 + ae2 |
|---|
| 236 | |
|---|
| 237 | c solving the system |
|---|
| 238 | c1 = gamma*nu1**3. * 0.5 |
|---|
| 239 | c2 = gamma*nu2**3. * 0.5 |
|---|
| 240 | a1 = c1 * p1 / (n1*l1) |
|---|
| 241 | a2 = c2 * p2 / (n2*l2) |
|---|
| 242 | a12 = (nu1/nu2)**3. * n2/n1 * p12/l1 |
|---|
| 243 | a21 = (nu2/nu1)**3. * n1/n2 * p21/l2 |
|---|
| 244 | el2 = (a2 + a21 * a1 ) / ( 1.0 - a21 * a12 ) |
|---|
| 245 | el1 = a1 + a12 * el2 |
|---|
| 246 | pl1 = el1 * n1 / c1 |
|---|
| 247 | pl2 = el2 * n2 / c2 |
|---|
| 248 | |
|---|
| 249 | c heating rate |
|---|
| 250 | hr1 = - hplanck*vlight * nu1 * ae1 * pl1 |
|---|
| 251 | hr2 = - hplanck*vlight * nu2 * ae2 * pl2 |
|---|
| 252 | hr(i) = hr1 + hr2 |
|---|
| 253 | dtnlte(j,i)=0.1*hr(i)*pt(j,i)/(4.4*pyy(i)) ! dtnlte in K s-1 |
|---|
| 254 | c write(7,25)pxx(i),hr1,hr2,hr(i),qt |
|---|
| 255 | c 25 format(' ',1p5e12.4) |
|---|
| 256 | |
|---|
| 257 | endif |
|---|
| 258 | |
|---|
| 259 | enddo ! end loop over layers |
|---|
| 260 | enddo ! end loop over grid points |
|---|
| 261 | c close(7) |
|---|
| 262 | c |
|---|
| 263 | return |
|---|
| 264 | end |
|---|
| 265 | |
|---|
| 266 | c*********************************************************************** |
|---|
| 267 | |
|---|
| 268 | subroutine interp1(escout,p,nlayer,escin,pin,nl) |
|---|
| 269 | C |
|---|
| 270 | C subroutine to perform linear interpolation in pressure from 1D profile |
|---|
| 271 | C escin(nl) sampled on pressure grid pin(nl) to profile |
|---|
| 272 | C escout(nlayer) on pressure grid p(nlayer). |
|---|
| 273 | C |
|---|
| 274 | real escout(nlayer),p(nlayer) |
|---|
| 275 | real escin(nl),pin(nl),wm,wp |
|---|
| 276 | integer nl,nlayer,n1,n,nm,np |
|---|
| 277 | do n1=1,nlayer |
|---|
| 278 | if(p(n1) .gt. 3.5 .or. p(n1) .lt. 4.0e-6) then |
|---|
| 279 | escout(n1) = 0.0 |
|---|
| 280 | else |
|---|
| 281 | do n = 1,nl-1 |
|---|
| 282 | if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then |
|---|
| 283 | nm=n |
|---|
| 284 | np=n+1 |
|---|
| 285 | wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np)) |
|---|
| 286 | wp=1.0 - wm |
|---|
| 287 | endif |
|---|
| 288 | enddo |
|---|
| 289 | escout(n1) = escin(nm)*wm + escin(np)*wp |
|---|
| 290 | endif |
|---|
| 291 | enddo |
|---|
| 292 | return |
|---|
| 293 | end |
|---|
| 294 | |
|---|
| 295 | c*********************************************************************** |
|---|
| 296 | |
|---|
| 297 | subroutine interp3(esco1,esco2,esco3,p,nlayer, |
|---|
| 298 | 1 esci1,esci2,esci3,pin,nl) |
|---|
| 299 | C |
|---|
| 300 | C subroutine to perform 3 simultaneous linear interpolations in pressure from |
|---|
| 301 | C 1D profiles esci1-3(nl) sampled on pressure grid pin(nl) to 1D profiles |
|---|
| 302 | C esco1-3(nlayer) on pressure grid p(ngrid,nlayer). |
|---|
| 303 | C |
|---|
| 304 | real esco1(nlayer),esco2(nlayer),esco3(nlayer),p(nlayer) |
|---|
| 305 | real esci1(nl), esci2(nl), esci3(nl), pin(nl),wm,wp |
|---|
| 306 | integer nl,nlayer,n1,n,nm,np |
|---|
| 307 | do n1=1,nlayer |
|---|
| 308 | if (p(n1).gt. 3.5 .or. p(n1) .lt. 4.0e-6) then |
|---|
| 309 | esco1(n1)=0.0 |
|---|
| 310 | esco2(n1)=0.0 |
|---|
| 311 | esco3(n1)=0.0 |
|---|
| 312 | else |
|---|
| 313 | do n = 1,nl-1 |
|---|
| 314 | if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then |
|---|
| 315 | nm=n |
|---|
| 316 | np=n+1 |
|---|
| 317 | wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np)) |
|---|
| 318 | wp=1.0 - wm |
|---|
| 319 | endif |
|---|
| 320 | enddo |
|---|
| 321 | esco1(n1) = esci1(nm)*wm + esci1(np)*wp |
|---|
| 322 | esco2(n1) = esci2(nm)*wm + esci2(np)*wp |
|---|
| 323 | esco3(n1) = esci3(nm)*wm + esci3(np)*wp |
|---|
| 324 | endif |
|---|
| 325 | enddo |
|---|
| 326 | return |
|---|
| 327 | end |
|---|
