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