[3012] | 1 | MODULE nltecool_mod |
---|
| 2 | |
---|
| 3 | IMPLICIT NONE |
---|
[3016] | 4 | |
---|
| 5 | ! Escape functions and VMRs from tabulated values. |
---|
| 6 | ! Origin: nlte_escape.dat file form Miguel L. |
---|
| 7 | ! (Y. Wanherdrick, 09/2000) |
---|
| 8 | integer,parameter :: np=68 ! # data points in tabulation |
---|
[3012] | 9 | |
---|
[3016] | 10 | real, save :: pnb(np) = [ ! Pressure (exponent) in tabulation |
---|
| 11 | & 12.0000, 11.0000, 10.8000, |
---|
| 12 | & 10.6000, 10.40000, 10.20000, |
---|
| 13 | & 10.00000, 9.80000, 9.60000, |
---|
| 14 | & 9.40000, 9.20000, 9.00000, |
---|
| 15 | & 8.80000, 8.60000, 8.40000, |
---|
| 16 | & 8.20000, 8.00000, 7.80000, |
---|
| 17 | & 7.60000, 7.40000, 7.20000, |
---|
| 18 | & 7.00000, 6.80000, 6.60000, |
---|
| 19 | & 6.40000, 6.20000, 6.00000, |
---|
| 20 | & 5.80000, 5.60000, 5.40000, |
---|
| 21 | & 5.20000, 5.00000, 4.80000, |
---|
| 22 | & 4.60000, 4.40000, 4.20000, |
---|
| 23 | & 4.00000, 3.80000, 3.60000, |
---|
| 24 | & 3.40000, 3.20000, 3.00000, |
---|
| 25 | & 2.80000, 2.60000, 2.40000, |
---|
| 26 | & 2.20000, 2.00000, 1.80000, |
---|
| 27 | & 1.60000, 1.40000, 1.20000, |
---|
| 28 | & 1.00000, 0.800000, 0.599999, |
---|
| 29 | & 0.400000, 0.200000, 0., |
---|
| 30 | & -0.200000, -0.400001, -0.600000, |
---|
| 31 | & -0.800000, -1.00000, -1.20000, |
---|
| 32 | & -1.40000, -1.60000, -1.80000, |
---|
| 33 | & -2.00000, -3.00000 ] ! values modified in nltecool |
---|
| 34 | !$OMP THREADPRIVATE(pnb) |
---|
| 35 | |
---|
| 36 | real,parameter :: ef1(np) = [ ! Esc.funct.#1, tabulated |
---|
| 37 | & 4.58112E-04, 4.58112E-04, 4.58112E-04, |
---|
| 38 | & 4.58112E-04, 4.58112E-04, 4.58112E-04, |
---|
| 39 | & 4.58112E-04, 4.58112E-04, 4.58112E-04, |
---|
| 40 | & 4.58112E-04, 4.58112E-04, 4.58112E-04, |
---|
| 41 | & 4.58112E-04, 4.58112E-04, 4.58112E-04, |
---|
| 42 | & 4.58112E-04, 4.58112E-04, 4.58112E-04, |
---|
| 43 | & 4.58112E-04, 4.58112E-04, 4.58112E-04, |
---|
| 44 | & 4.58112E-04, 4.61707E-04, 4.76886E-04, |
---|
| 45 | & 4.95638E-04, 5.20935E-04, 5.55511E-04, |
---|
| 46 | & 6.01219E-04, 6.63734E-04, 7.50691E-04, |
---|
| 47 | & 8.63474E-04, 1.00900E-03, 1.19642E-03, |
---|
| 48 | & 1.42690E-03, 1.71398E-03, 2.06663E-03, |
---|
| 49 | & 2.48974E-03, 3.01578E-03, 3.64350E-03, |
---|
| 50 | & 4.40323E-03, 5.32066E-03, 6.40456E-03, |
---|
| 51 | & 7.72069E-03, 9.25684E-03, 1.10905E-02, |
---|
| 52 | & 1.32374E-02, 1.57643E-02, 1.87388E-02, |
---|
| 53 | & 2.22072E-02, 2.63099E-02, 3.10614E-02, |
---|
| 54 | & 3.66948E-02, 4.32373E-02, 5.15022E-02, |
---|
| 55 | & 6.21455E-02, 7.77212E-02, 9.92027E-02, |
---|
| 56 | & 0.131155, 0.179470, 0.258913, |
---|
| 57 | & 0.380549, 0.530450, 0.643180, |
---|
| 58 | & 0.741061, 0.826336, 0.922787, |
---|
| 59 | & 0.997203, 1.00000 ] |
---|
| 60 | |
---|
| 61 | real,parameter :: ef2(np) = [ ! Esc.funct.#2, tabulated |
---|
| 62 | & 1.98992E-03, 1.98992E-03, 1.98992E-03, |
---|
| 63 | & 1.98992E-03, 1.98992E-03, 1.98992E-03, |
---|
| 64 | & 1.98992E-03, 1.98992E-03, 1.98992E-03, |
---|
| 65 | & 1.98992E-03, 1.98992E-03, 2.01376E-03, |
---|
| 66 | & 2.09450E-03, 2.22993E-03, 2.42056E-03, |
---|
| 67 | & 2.68018E-03, 3.04398E-03, 3.43896E-03, |
---|
| 68 | & 3.80282E-03, 4.20622E-03, 4.76121E-03, |
---|
| 69 | & 8.01698E-03, 1.19947E-02, 1.69149E-02, |
---|
| 70 | & 2.24497E-02, 2.85244E-02, 3.54813E-02, |
---|
| 71 | & 4.39264E-02, 5.46248E-02, 6.75367E-02, |
---|
| 72 | & 8.29931E-02, 1.01717E-01, 0.123422, |
---|
| 73 | & 0.148468, 0.177096, 0.208816, |
---|
| 74 | & 0.244003, 0.282013, 0.322559, |
---|
| 75 | & 0.365542, 0.410518, 0.457384, |
---|
| 76 | & 0.505358, 0.553627, 0.600472, |
---|
| 77 | & 0.644807, 0.687185, 0.727429, |
---|
| 78 | & 0.764734, 0.798562, 0.828699, |
---|
| 79 | & 0.854797, 0.877717, 0.897874, |
---|
| 80 | & 0.915258, 0.929904, 0.942381, |
---|
| 81 | & 0.952906, 0.962173, 0.970191, |
---|
| 82 | & 0.976437, 0.981501, 0.985406, |
---|
| 83 | & 0.988560, 0.991111, 0.993653, |
---|
| 84 | & 0.995561, 1.00000 ] |
---|
| 85 | |
---|
| 86 | real,parameter :: co2vmr(np) = [ ! CO2 VMR tabulated |
---|
| 87 | & 0.950000, 0.950000, 0.950000, |
---|
| 88 | & 0.950000, 0.950000, 0.950000, |
---|
| 89 | & 0.950000, 0.950000, 0.950000, |
---|
| 90 | & 0.950000, 0.950000, 0.950000, |
---|
| 91 | & 0.950000, 0.950000, 0.950000, |
---|
| 92 | & 0.950000, 0.950000, 0.950000, |
---|
| 93 | & 0.950000, 0.950000, 0.950000, |
---|
| 94 | & 0.950000, 0.950000, 0.950000, |
---|
| 95 | & 0.950000, 0.950000, 0.950000, |
---|
| 96 | & 0.950000, 0.950000, 0.950000, |
---|
| 97 | & 0.950000, 0.950000, 0.950000, |
---|
| 98 | & 0.950000, 0.950000, 0.950000, |
---|
| 99 | & 0.950000, 0.950000, 0.950000, |
---|
| 100 | & 0.950000, 0.950000, 0.950000, |
---|
| 101 | & 0.950000, 0.950000, 0.950000, |
---|
| 102 | & 0.950000, 0.950000, 0.950000, |
---|
| 103 | & 0.950000, 0.950000, 0.950000, |
---|
| 104 | & 0.950000, 0.950000, 0.950000, |
---|
| 105 | & 0.950000, 0.950000, 0.950000, |
---|
| 106 | & 0.950000, 0.950000, 0.950000, |
---|
| 107 | & 0.949619, 0.947694, 0.945830, |
---|
| 108 | & 0.944016, 0.940557, 0.937068, |
---|
| 109 | & 0.932366, 0.893661 ] |
---|
| 110 | |
---|
| 111 | real,parameter :: o3pvmr(np) = [ ! O3p VMR tabulated |
---|
| 112 | & 5.06756E-08, 9.16539E-07, 1.68217E-06, |
---|
| 113 | & 3.00843E-06, 5.03151E-06, 8.07489E-06, |
---|
| 114 | & 1.23137E-05, 1.79029E-05, 2.45308E-05, |
---|
| 115 | & 3.27431E-05, 4.26692E-05, 5.44396E-05, |
---|
| 116 | & 6.78865E-05, 8.33147E-05, 1.00148E-04, |
---|
| 117 | & 1.18846E-04, 1.39681E-04, 1.64909E-04, |
---|
| 118 | & 1.93617E-04, 2.25161E-04, 2.60834E-04, |
---|
| 119 | & 3.01501E-04, 3.44953E-04, 3.91011E-04, |
---|
| 120 | & 4.40377E-04, 4.90820E-04, 5.43200E-04, |
---|
| 121 | & 5.95335E-04, 6.45420E-04, 6.93166E-04, |
---|
| 122 | & 7.43729E-04, 7.93710E-04, 8.44394E-04, |
---|
| 123 | & 8.94318E-04, 9.44732E-04, 9.94964E-04, |
---|
| 124 | & 1.04901E-03, 1.10008E-03, 1.16302E-03, |
---|
| 125 | & 1.22989E-03, 1.30026E-03, 1.37131E-03, |
---|
| 126 | & 1.45556E-03, 1.55186E-03, 1.66328E-03, |
---|
| 127 | & 1.77802E-03, 1.91546E-03, 2.07503E-03, |
---|
| 128 | & 2.24903E-03, 2.47117E-03, 2.71728E-03, |
---|
| 129 | & 2.99739E-03, 3.33582E-03, 3.73507E-03, |
---|
| 130 | & 4.20819E-03, 4.76887E-03, 5.42558E-03, |
---|
| 131 | & 6.20815E-03, 7.14473E-03, 8.28545E-03, |
---|
| 132 | & 9.51779E-03, 1.08140E-02, 1.22359E-02, |
---|
| 133 | & 1.36870E-02, 1.51495E-02, 1.67196E-02, |
---|
| 134 | & 1.85485E-02, 3.36252E-02 ] |
---|
| 135 | |
---|
| 136 | real, parameter :: n2covmr(np) = [ ! N2+CO VMR tabulated |
---|
| 137 | & 2.71412E-02, 2.71464E-02, 2.71490E-02, |
---|
| 138 | & 2.71523E-02, 2.71558E-02, 2.71617E-02, |
---|
| 139 | & 2.71672E-02, 2.71749E-02, 2.71837E-02, |
---|
| 140 | & 2.71943E-02, 2.72058E-02, 2.72189E-02, |
---|
| 141 | & 2.72326E-02, 2.72483E-02, 2.72661E-02, |
---|
| 142 | & 2.72848E-02, 2.73054E-02, 2.73279E-02, |
---|
| 143 | & 2.73514E-02, 2.73775E-02, 2.74048E-02, |
---|
| 144 | & 2.74345E-02, 2.74672E-02, 2.75021E-02, |
---|
| 145 | & 2.75404E-02, 2.75826E-02, 2.76340E-02, |
---|
| 146 | & 2.77013E-02, 2.78220E-02, 2.79707E-02, |
---|
| 147 | & 2.81759E-02, 2.84339E-02, 2.87587E-02, |
---|
| 148 | & 2.91600E-02, 2.96561E-02, 3.02558E-02, |
---|
| 149 | & 3.09922E-02, 3.18062E-02, 3.27010E-02, |
---|
| 150 | & 3.35635E-02, 3.44388E-02, 3.53327E-02, |
---|
| 151 | & 3.62143E-02, 3.70941E-02, 3.79315E-02, |
---|
| 152 | & 3.87626E-02, 3.95524E-02, 4.03071E-02, |
---|
| 153 | & 4.10071E-02, 4.16229E-02, 4.21231E-02, |
---|
| 154 | & 4.25167E-02, 4.27964E-02, 4.29773E-02, |
---|
| 155 | & 4.30488E-02, 4.29638E-02, 4.28049E-02, |
---|
| 156 | & 4.26788E-02, 4.26822E-02, 4.29426E-02, |
---|
| 157 | & 4.34634E-02, 4.42559E-02, 4.53038E-02, |
---|
| 158 | & 4.65879E-02, 4.80262E-02, 4.96303E-02, |
---|
| 159 | & 5.14885E-02, 6.91651E-02 ] |
---|
| 160 | |
---|
[3012] | 161 | CONTAINS |
---|
[38] | 162 | c************************************************************************** |
---|
| 163 | c |
---|
[414] | 164 | subroutine nltecool(ngrid,nlayer,nq,pplay,pt,pq,dtnlte) |
---|
[38] | 165 | c |
---|
| 166 | c This code was designed as a delivery for the "Martian Environment Models" |
---|
| 167 | c project ( ESA contract 11369/95/nl/jg CCN2 ) |
---|
| 168 | c Computes non-LTE heating rates from CO2 emission at 15 um |
---|
| 169 | c in the Martian upper atmosphere. |
---|
| 170 | c Uses a simplified model consisting of two excited levels with two |
---|
| 171 | c emission bands, one of them stronger than the other, which correspond |
---|
| 172 | c to the behaviours of the 626 fundamental band and the isotopic fund.bands. |
---|
| 173 | c It uses a cool-to-space approximation with tabulated escape functions. |
---|
| 174 | c These escape functions have been precomputed for the strong and weak bands, |
---|
| 175 | c and are given as a function of pressure in separate files. |
---|
| 176 | c The output values are the heating rates (actually, cooling, since they |
---|
| 177 | c are always negative) for the two bands, i.e., the total cooling is the |
---|
| 178 | c sum of them. |
---|
| 179 | c Miguel A. Lopez Valverde |
---|
| 180 | c Instituto de Astrofisica de Andalucia (CSIC), Granada, Spain |
---|
| 181 | c |
---|
| 182 | c Version 1b. See description above. 22-March-2000. |
---|
| 183 | c Adapted as a subroutine for use in GCM -- PLR/SRL 6/2000 |
---|
| 184 | c Version 1c. Inclusion of VMR in the tabulation of escape functions. |
---|
| 185 | c Table contains now only 1 input file -- Miguel 11/Jul/2000 |
---|
| 186 | c Version 1d data contained in original input file "nlte_escape.dat" |
---|
| 187 | c now stored in include file "nltedata.h" Y.Wanherdrick 09/2000 |
---|
[414] | 188 | |
---|
| 189 | c jul 2011 fgg Modified to allow variable O |
---|
[38] | 190 | c |
---|
| 191 | c*************************************************************************** |
---|
| 192 | |
---|
[1036] | 193 | use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_n2, mmol |
---|
[1047] | 194 | use conc_mod, only: mmean |
---|
[38] | 195 | implicit none |
---|
| 196 | |
---|
[3012] | 197 | include "callkeys.h" |
---|
[414] | 198 | |
---|
[38] | 199 | c Input and output variables |
---|
| 200 | c |
---|
[3012] | 201 | integer,intent(in) :: ngrid ! no. of horiz. gridpoints |
---|
| 202 | integer,intent(in) :: nlayer ! no. of atmospheric layers |
---|
| 203 | integer,intent(in) :: nq ! no. of tracers |
---|
| 204 | real,intent(in) :: pplay(ngrid,nlayer) ! input pressure grid (Pa) |
---|
| 205 | real,intent(in) :: pt(ngrid,nlayer) ! input temperatures (K) |
---|
| 206 | real,intent(in) :: pq(ngrid,nlayer,nq) ! input mmrs (kg/kg_air) |
---|
| 207 | real,intent(out) :: dtnlte(ngrid,nlayer) ! output temp. tendencies (K/s) |
---|
[38] | 208 | |
---|
| 209 | c |
---|
| 210 | c Standard atmosphere variables |
---|
| 211 | c |
---|
| 212 | real nt ! number density [cm-3] |
---|
| 213 | real co2(nlayer) ! " of CO2 |
---|
| 214 | real o3p(nlayer) ! " of atomic oxygen |
---|
| 215 | real n2co(nlayer) ! " of N2 + CO |
---|
| 216 | real pyy(nlayer) ! auxiliary pressure grid |
---|
| 217 | |
---|
| 218 | c |
---|
| 219 | c Vectors and indexes for the tabulation of escape functions and VMR |
---|
| 220 | c |
---|
| 221 | c np ! # data points in tabulation |
---|
| 222 | c pnb(np) ! Pressure in tabulation |
---|
| 223 | c ef1(np) ! Esc.funct.#1, tabulated |
---|
| 224 | c ef2(np) ! Esc.funct.#2, tabulated |
---|
| 225 | c co2vmr(np) ! CO2 VMR tabulated |
---|
| 226 | c o3pvmr(np) ! CO2 VMR tabulated |
---|
| 227 | c n2covmr(np) ! N2+CO VMR tabulated |
---|
| 228 | real escf1(nlayer) ! Esc.funct.#1, interpolated |
---|
| 229 | real escf2(nlayer) ! Esc.funct.#2, interpolated |
---|
| 230 | |
---|
| 231 | |
---|
| 232 | c |
---|
| 233 | c Local Constants |
---|
| 234 | c |
---|
| 235 | real nu1, nu2 ! freq. of energy levels |
---|
| 236 | real imr1, imr2 ! isotopic abundances |
---|
| 237 | real hplanck, gamma, vlight ! physical constants |
---|
| 238 | real ee |
---|
| 239 | real rfvt ! collisional rate |
---|
| 240 | real rfvto3p ! " |
---|
| 241 | real rfvv ! " |
---|
| 242 | |
---|
| 243 | c |
---|
| 244 | c Local variables for the main loop |
---|
| 245 | c |
---|
| 246 | real n1, n2, co2t ! ground populations |
---|
| 247 | real l1, p1, p12 ! prod & losses |
---|
| 248 | real l2, p2, p21 |
---|
| 249 | real tt ! dummy variable |
---|
| 250 | real c1, c2 ! molecular constants |
---|
| 251 | real ae1, ae2 ! einstein spontaneous emission |
---|
| 252 | real a1, a2, a12, a21 |
---|
| 253 | real pl1, pl2 |
---|
| 254 | real el1, el2 |
---|
| 255 | real hr1, hr2 ! heating rate due to each band |
---|
| 256 | real hr(nlayer) ! total heating rate |
---|
| 257 | |
---|
| 258 | c |
---|
| 259 | c Indexes |
---|
| 260 | c |
---|
| 261 | integer i |
---|
| 262 | integer j,ii |
---|
| 263 | |
---|
| 264 | c |
---|
| 265 | c Rate coefficients |
---|
| 266 | c |
---|
| 267 | real k19xca, k19xcb |
---|
| 268 | real k19cap1, k19cap2 |
---|
| 269 | real k19cbp1, k19cbp2 |
---|
| 270 | real d19c, d19cp1, d19cp2 |
---|
| 271 | real k20xc, k20cp1, k20cp2 |
---|
| 272 | real k21xc, k21cp2 |
---|
| 273 | |
---|
[3016] | 274 | logical,save :: firstcall=.true. |
---|
| 275 | !$OMP THREADPRIVATE(firstcall) |
---|
[38] | 276 | |
---|
| 277 | c |
---|
| 278 | c Data |
---|
| 279 | c |
---|
| 280 | data nu1, nu2, hplanck, gamma, vlight, ee/ |
---|
| 281 | 1 667.38, 662.3734, 6.6261e-27, 1.191e-5, 3.e10, 1.438769/ |
---|
| 282 | |
---|
| 283 | c************************************************************************* |
---|
| 284 | c PROGRAM STARTS |
---|
| 285 | c************************************************************************* |
---|
| 286 | |
---|
| 287 | imr1 = 0.987 |
---|
| 288 | imr2 = 0.00408 + 0.0112 |
---|
| 289 | rfvt = 0.1 |
---|
| 290 | rfvto3p = 1.0 |
---|
| 291 | rfvv = 0.1 |
---|
| 292 | |
---|
| 293 | if(firstcall) then |
---|
| 294 | |
---|
| 295 | do i=1,np |
---|
| 296 | pnb(i)=1.0e-4*exp(pnb(i)) ! p into Pa |
---|
| 297 | end do |
---|
| 298 | |
---|
| 299 | firstcall = .false. |
---|
| 300 | |
---|
| 301 | endif |
---|
| 302 | c |
---|
| 303 | c MAIN LOOP, for each gridpoint and altitude: |
---|
| 304 | c |
---|
| 305 | do j=1,ngrid ! loop over grid points |
---|
| 306 | c |
---|
| 307 | c set up local pressure grid |
---|
| 308 | c |
---|
| 309 | do ii=1,nlayer |
---|
| 310 | pyy(ii)=pplay(j,ii) |
---|
| 311 | enddo |
---|
| 312 | ! |
---|
| 313 | ! Interpolate escape functions and VMR to the desired grid |
---|
| 314 | ! |
---|
| 315 | call interp1(escf2,pyy,nlayer,ef2,pnb,np) |
---|
| 316 | call interp1(escf1,pyy,nlayer,ef1,pnb,np) |
---|
[414] | 317 | if(nltemodel.eq.0) then |
---|
| 318 | call interp3(co2,o3p,n2co,pyy,nlayer, |
---|
| 319 | & co2vmr,o3pvmr,n2covmr,pnb,np) |
---|
| 320 | endif |
---|
[38] | 321 | |
---|
| 322 | do i=1,nlayer ! loop over layers |
---|
| 323 | C |
---|
| 324 | C test if p lies outside range (p > 3.5 Pa) |
---|
| 325 | C changed to 1 Pa since transition will always be higher than this |
---|
| 326 | C |
---|
| 327 | if(pyy(i) .gt. 1.0 .or. pyy(i) .lt. 4.0e-6) then |
---|
| 328 | hr(i)=0.0 |
---|
| 329 | dtnlte(j,i)=0.0 |
---|
| 330 | else |
---|
| 331 | c |
---|
| 332 | c if(pt(j,i).lt.1.0)print*,pt(j,i) |
---|
| 333 | nt = pyy(i)/(1.381e-17*pt(j,i)) ! nt in cm-3 |
---|
[414] | 334 | !Dynamic composition |
---|
| 335 | if(nltemodel.eq.1) then |
---|
[472] | 336 | co2(i)=pq(j,i,igcm_co2)*mmean(j,i)/mmol(igcm_co2) |
---|
| 337 | o3p(i)=pq(j,i,igcm_o)*mmean(j,i)/mmol(igcm_o) |
---|
| 338 | n2co(i)=pq(j,i,igcm_co)*mmean(j,i)/mmol(igcm_co) + |
---|
| 339 | $ pq(j,i,igcm_n2)*mmean(j,i)/mmol(igcm_n2) |
---|
[414] | 340 | endif |
---|
| 341 | |
---|
| 342 | !Mixing ratio to density |
---|
[38] | 343 | co2(i)=co2(i)*nt ! CO2 density in cm-3 |
---|
| 344 | o3p(i)=o3p(i)*nt ! O3p density in cm-3 |
---|
| 345 | n2co(i)=n2co(i)*nt ! N2+CO in cm-3 |
---|
| 346 | c molecular populations |
---|
| 347 | n1 = co2(i) * imr1 |
---|
| 348 | n2 = co2(i) * imr2 |
---|
| 349 | co2t = n1 + n2 |
---|
| 350 | |
---|
| 351 | c intermediate collisional rates |
---|
| 352 | tt = pt(j,i)*pt(j,i) |
---|
| 353 | |
---|
| 354 | if (pt(j,i).le.175.0) then |
---|
| 355 | k19xca = 3.3e-15 |
---|
| 356 | k19xcb = 7.6e-16 |
---|
| 357 | else |
---|
| 358 | k19xca = 4.2e-12 * exp( -2988.0/pt(j,i) + 303930.0/tt) |
---|
| 359 | k19xcb = 2.1e-12 * exp( -2659.0/pt(j,i) + 223052.0/tt) |
---|
| 360 | endif |
---|
| 361 | k19xca = k19xca * rfvt |
---|
| 362 | k19xcb = k19xcb * rfvt |
---|
| 363 | k19cap1 = k19xca * 2.0 * exp( -ee*nu1/pt(j,i) ) |
---|
| 364 | k19cap2 = k19xca * 2.0 * exp( -ee*nu2/pt(j,i) ) |
---|
| 365 | k19cbp1 = k19xcb * 2.0 * exp( -ee*nu1/pt(j,i) ) |
---|
| 366 | k19cbp2 = k19xcb * 2.0 * exp( -ee*nu2/pt(j,i) ) |
---|
| 367 | d19c = k19xca*co2t + k19xcb*n2co(i) |
---|
| 368 | d19cp1 = k19cap1*co2t + k19cbp1*n2co(i) |
---|
| 369 | d19cp2 = k19cap2*co2t + k19cbp2*n2co(i) |
---|
| 370 | ! |
---|
| 371 | k20xc = 3.e-12 * rfvto3p |
---|
| 372 | k20cp1 = k20xc * 2.0 * exp( -ee/pt(j,i) * nu1 ) |
---|
| 373 | k20cp2 = k20xc * 2.0 * exp( -ee/pt(j,i) * nu2 ) |
---|
| 374 | ! |
---|
| 375 | k21xc = 2.49e-11 * 0.5 * rfvv |
---|
| 376 | k21cp2 = k21xc * exp( - ee/pt(j,i) * (nu2-nu1) ) |
---|
| 377 | ! |
---|
| 378 | l1 = d19c + k20xc*o3p(i) + k21cp2*n2 |
---|
| 379 | p1 = ( d19cp1 + k20cp1*o3p(i) ) * n1 |
---|
| 380 | p12 = k21xc*n1 |
---|
| 381 | ! |
---|
| 382 | l2 = d19c + k20xc*o3p(i) + k21xc*n1 |
---|
| 383 | p2 = ( d19cp2 + k20cp2*o3p(i) ) * n2 |
---|
| 384 | p21 = k21cp2*n2 |
---|
| 385 | |
---|
| 386 | c radiative rates |
---|
| 387 | ae1 = 1.3546 * 1.66 / 4.0 * escf1(i) |
---|
| 388 | ae2 = ( 1.3452 + 1.1878 ) * 1.66 / 4.0 * escf2(i) |
---|
| 389 | l1 = l1 + ae1 |
---|
| 390 | l2 = l2 + ae2 |
---|
| 391 | |
---|
| 392 | c solving the system |
---|
| 393 | c1 = gamma*nu1**3. * 0.5 |
---|
| 394 | c2 = gamma*nu2**3. * 0.5 |
---|
| 395 | a1 = c1 * p1 / (n1*l1) |
---|
| 396 | a2 = c2 * p2 / (n2*l2) |
---|
| 397 | a12 = (nu1/nu2)**3. * n2/n1 * p12/l1 |
---|
| 398 | a21 = (nu2/nu1)**3. * n1/n2 * p21/l2 |
---|
| 399 | el2 = (a2 + a21 * a1 ) / ( 1.0 - a21 * a12 ) |
---|
| 400 | el1 = a1 + a12 * el2 |
---|
| 401 | pl1 = el1 * n1 / c1 |
---|
| 402 | pl2 = el2 * n2 / c2 |
---|
| 403 | |
---|
| 404 | c heating rate |
---|
| 405 | hr1 = - hplanck*vlight * nu1 * ae1 * pl1 |
---|
| 406 | hr2 = - hplanck*vlight * nu2 * ae2 * pl2 |
---|
| 407 | hr(i) = hr1 + hr2 |
---|
| 408 | dtnlte(j,i)=0.1*hr(i)*pt(j,i)/(4.4*pyy(i)) ! dtnlte in K s-1 |
---|
| 409 | c write(7,25)pxx(i),hr1,hr2,hr(i),qt |
---|
| 410 | c 25 format(' ',1p5e12.4) |
---|
| 411 | |
---|
| 412 | endif |
---|
| 413 | |
---|
| 414 | enddo ! end loop over layers |
---|
| 415 | enddo ! end loop over grid points |
---|
| 416 | c close(7) |
---|
| 417 | c |
---|
| 418 | |
---|
[3012] | 419 | end subroutine nltecool |
---|
| 420 | |
---|
[38] | 421 | c*********************************************************************** |
---|
| 422 | |
---|
| 423 | subroutine interp1(escout,p,nlayer,escin,pin,nl) |
---|
| 424 | C |
---|
| 425 | C subroutine to perform linear interpolation in pressure from 1D profile |
---|
| 426 | C escin(nl) sampled on pressure grid pin(nl) to profile |
---|
| 427 | C escout(nlayer) on pressure grid p(nlayer). |
---|
| 428 | C |
---|
| 429 | real escout(nlayer),p(nlayer) |
---|
| 430 | real escin(nl),pin(nl),wm,wp |
---|
| 431 | integer nl,nlayer,n1,n,nm,np |
---|
| 432 | do n1=1,nlayer |
---|
| 433 | if(p(n1) .gt. 3.5 .or. p(n1) .lt. 4.0e-6) then |
---|
| 434 | escout(n1) = 0.0 |
---|
| 435 | else |
---|
| 436 | do n = 1,nl-1 |
---|
| 437 | if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then |
---|
| 438 | nm=n |
---|
| 439 | np=n+1 |
---|
| 440 | wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np)) |
---|
| 441 | wp=1.0 - wm |
---|
| 442 | endif |
---|
| 443 | enddo |
---|
| 444 | escout(n1) = escin(nm)*wm + escin(np)*wp |
---|
| 445 | endif |
---|
| 446 | enddo |
---|
| 447 | |
---|
[3012] | 448 | end subroutine interp1 |
---|
| 449 | |
---|
[38] | 450 | c*********************************************************************** |
---|
| 451 | |
---|
| 452 | subroutine interp3(esco1,esco2,esco3,p,nlayer, |
---|
| 453 | 1 esci1,esci2,esci3,pin,nl) |
---|
| 454 | C |
---|
| 455 | C subroutine to perform 3 simultaneous linear interpolations in pressure from |
---|
| 456 | C 1D profiles esci1-3(nl) sampled on pressure grid pin(nl) to 1D profiles |
---|
| 457 | C esco1-3(nlayer) on pressure grid p(ngrid,nlayer). |
---|
| 458 | C |
---|
| 459 | real esco1(nlayer),esco2(nlayer),esco3(nlayer),p(nlayer) |
---|
| 460 | real esci1(nl), esci2(nl), esci3(nl), pin(nl),wm,wp |
---|
| 461 | integer nl,nlayer,n1,n,nm,np |
---|
| 462 | do n1=1,nlayer |
---|
| 463 | if (p(n1).gt. 3.5 .or. p(n1) .lt. 4.0e-6) then |
---|
| 464 | esco1(n1)=0.0 |
---|
| 465 | esco2(n1)=0.0 |
---|
| 466 | esco3(n1)=0.0 |
---|
| 467 | else |
---|
| 468 | do n = 1,nl-1 |
---|
| 469 | if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then |
---|
| 470 | nm=n |
---|
| 471 | np=n+1 |
---|
| 472 | wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np)) |
---|
| 473 | wp=1.0 - wm |
---|
| 474 | endif |
---|
| 475 | enddo |
---|
| 476 | esco1(n1) = esci1(nm)*wm + esci1(np)*wp |
---|
| 477 | esco2(n1) = esci2(nm)*wm + esci2(np)*wp |
---|
| 478 | esco3(n1) = esci3(nm)*wm + esci3(np)*wp |
---|
| 479 | endif |
---|
| 480 | enddo |
---|
[3012] | 481 | |
---|
| 482 | end subroutine interp3 |
---|
| 483 | |
---|
| 484 | END MODULE nltecool_mod |
---|