| 1 | !========================================================================== |
|---|
| 2 | |
|---|
| 3 | subroutine photolysis(nlayer, & |
|---|
| 4 | lswitch, press, temp, sza, tauref, & |
|---|
| 5 | zmmean, dist_sol, rmco2, rmo3, v_phot) |
|---|
| 6 | |
|---|
| 7 | !========================================================================== |
|---|
| 8 | |
|---|
| 9 | use comcstfi_h |
|---|
| 10 | use photolysis_mod, only : nb_phot_max |
|---|
| 11 | |
|---|
| 12 | implicit none |
|---|
| 13 | |
|---|
| 14 | #include "chimiedata.h" |
|---|
| 15 | |
|---|
| 16 | !========================================================================== |
|---|
| 17 | ! input: |
|---|
| 18 | !========================================================================== |
|---|
| 19 | |
|---|
| 20 | integer, intent(in) :: nlayer ! number of atmospheric layers |
|---|
| 21 | integer :: lswitch ! interface level between chemistries |
|---|
| 22 | real :: press(nlayer) ! pressure (hPa) |
|---|
| 23 | real :: temp(nlayer) ! temperature (K) |
|---|
| 24 | real :: sza ! solar zenith angle (deg) |
|---|
| 25 | real :: tauref ! optical depth at 7 hpa |
|---|
| 26 | real :: zmmean(nlayer) ! mean molecular mass (g) |
|---|
| 27 | real :: dist_sol ! sun distance (AU) |
|---|
| 28 | real :: rmco2(nlayer) ! co2 mixing ratio |
|---|
| 29 | real :: rmo3(nlayer) ! ozone mixing ratio |
|---|
| 30 | |
|---|
| 31 | !========================================================================== |
|---|
| 32 | ! output: interpolated photodissociation rates (s-1) |
|---|
| 33 | !========================================================================== |
|---|
| 34 | |
|---|
| 35 | real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot |
|---|
| 36 | |
|---|
| 37 | !========================================================================== |
|---|
| 38 | ! local: |
|---|
| 39 | !========================================================================== |
|---|
| 40 | |
|---|
| 41 | integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp, & |
|---|
| 42 | iozo, isza, itau, it, l |
|---|
| 43 | |
|---|
| 44 | integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, & |
|---|
| 45 | j_o3_o, j_h2o, j_hdo, j_h2o2, j_ho2, j_no, j_no2, & |
|---|
| 46 | j_hno3, j_hno4 |
|---|
| 47 | |
|---|
| 48 | real :: col(nlayer) ! overhead air column (molecule cm-2) |
|---|
| 49 | real :: colo3(nlayer) ! overhead ozone column (molecule cm-2) |
|---|
| 50 | real :: poids(2,2,2,2,2) ! 5D interpolation weights |
|---|
| 51 | real :: tref ! temperature at 1.9 hPa in the gcm (K) |
|---|
| 52 | real :: table_temp(ntemp) ! temperatures at 1.9 hPa in jmars (K) |
|---|
| 53 | real :: cinf, csup, cicol, ciozo, cisza, citemp, citau |
|---|
| 54 | real :: colo3min, dp, coef |
|---|
| 55 | real :: ratio_o3(nlayer) |
|---|
| 56 | real :: tau |
|---|
| 57 | real :: j(nlayer,nd) |
|---|
| 58 | |
|---|
| 59 | !========================================================================== |
|---|
| 60 | ! day/night criterion |
|---|
| 61 | !========================================================================== |
|---|
| 62 | |
|---|
| 63 | if (sza <= 95.) then |
|---|
| 64 | |
|---|
| 65 | !========================================================================== |
|---|
| 66 | ! temperatures at 1.9 hPa in lookup table |
|---|
| 67 | !========================================================================== |
|---|
| 68 | |
|---|
| 69 | table_temp(1) = 226.2 |
|---|
| 70 | table_temp(2) = 206.2 |
|---|
| 71 | table_temp(3) = 186.2 |
|---|
| 72 | table_temp(4) = 169.8 |
|---|
| 73 | |
|---|
| 74 | !========================================================================== |
|---|
| 75 | ! interpolation in solar zenith angle |
|---|
| 76 | !========================================================================== |
|---|
| 77 | |
|---|
| 78 | indsza = nsza - 1 |
|---|
| 79 | do isza = 1,nsza |
|---|
| 80 | if (szatab(isza) >= sza) then |
|---|
| 81 | indsza = isza - 1 |
|---|
| 82 | indsza = max(indsza, 1) |
|---|
| 83 | exit |
|---|
| 84 | end if |
|---|
| 85 | end do |
|---|
| 86 | cisza = (sza - szatab(indsza)) & |
|---|
| 87 | /(szatab(indsza + 1) - szatab(indsza)) |
|---|
| 88 | |
|---|
| 89 | !========================================================================== |
|---|
| 90 | ! interpolation in dust (tau) |
|---|
| 91 | !========================================================================== |
|---|
| 92 | |
|---|
| 93 | tau = min(tauref, tautab(ntau)) |
|---|
| 94 | tau = max(tau, tautab(1)) |
|---|
| 95 | |
|---|
| 96 | indtau = ntau - 1 |
|---|
| 97 | do itau = 1,ntau |
|---|
| 98 | if (tautab(itau) >= tau) then |
|---|
| 99 | indtau = itau - 1 |
|---|
| 100 | indtau = max(indtau, 1) |
|---|
| 101 | exit |
|---|
| 102 | end if |
|---|
| 103 | end do |
|---|
| 104 | citau = (tau - tautab(indtau)) & |
|---|
| 105 | /(tautab(indtau + 1) - tautab(indtau)) |
|---|
| 106 | |
|---|
| 107 | !========================================================================== |
|---|
| 108 | ! co2 and ozone columns |
|---|
| 109 | !========================================================================== |
|---|
| 110 | |
|---|
| 111 | ! co2 column at model top (molecule.cm-2) |
|---|
| 112 | |
|---|
| 113 | col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100. & |
|---|
| 114 | /(zmmean(lswitch-1)*g) |
|---|
| 115 | |
|---|
| 116 | ! ozone column at model top |
|---|
| 117 | |
|---|
| 118 | colo3(lswitch-1) = 0. |
|---|
| 119 | |
|---|
| 120 | ! co2 and ozone columns for other levels (molecule.cm-2) |
|---|
| 121 | |
|---|
| 122 | do l = lswitch-2,1,-1 |
|---|
| 123 | dp = (press(l) - press(l+1))*100. |
|---|
| 124 | col(l) = col(l+1) + (rmco2(l+1) + rmco2(l))*0.5 & |
|---|
| 125 | *6.022e22*dp/(zmmean(l)*g) |
|---|
| 126 | col(l) = min(col(l), colairtab(1)) |
|---|
| 127 | colo3(l) = colo3(l+1) + (rmo3(l+1) + rmo3(l))*0.5 & |
|---|
| 128 | *6.022e22*dp/(zmmean(l)*g) |
|---|
| 129 | end do |
|---|
| 130 | |
|---|
| 131 | ! ratio ozone column/minimal theoretical column (0.1 micron-atm) |
|---|
| 132 | |
|---|
| 133 | ! ro3 = 7.171e-10 is the o3 mixing ratio for a uniform |
|---|
| 134 | ! profile giving a column 0.1 micron-atmosphere at |
|---|
| 135 | ! a surface pressure of 10 hpa. |
|---|
| 136 | |
|---|
| 137 | do l = 1,lswitch-1 |
|---|
| 138 | colo3min = col(l)*7.171e-10 |
|---|
| 139 | ratio_o3(l) = colo3(l)/colo3min |
|---|
| 140 | ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo)*10.) |
|---|
| 141 | ratio_o3(l) = max(ratio_o3(l), 1.) |
|---|
| 142 | end do |
|---|
| 143 | |
|---|
| 144 | !========================================================================== |
|---|
| 145 | ! temperature dependence |
|---|
| 146 | !========================================================================== |
|---|
| 147 | |
|---|
| 148 | ! 1) search for temperature at 1.9 hPa (tref): vertical interpolation |
|---|
| 149 | |
|---|
| 150 | tref = temp(1) |
|---|
| 151 | do l = (lswitch-1)-1,1,-1 |
|---|
| 152 | if (press(l) > 1.9) then |
|---|
| 153 | cinf = (press(l) - 1.9)/(press(l) - press(l+1)) |
|---|
| 154 | csup = 1. - cinf |
|---|
| 155 | tref = cinf*temp(l+1) + csup*temp(l) |
|---|
| 156 | exit |
|---|
| 157 | end if |
|---|
| 158 | end do |
|---|
| 159 | |
|---|
| 160 | ! 2) interpolation in temperature |
|---|
| 161 | |
|---|
| 162 | tref = min(tref,table_temp(1)) |
|---|
| 163 | tref = max(tref,table_temp(ntemp)) |
|---|
| 164 | |
|---|
| 165 | do it = 2, ntemp |
|---|
| 166 | if (table_temp(it) <= tref) then |
|---|
| 167 | citemp = (log(tref) - log(table_temp(it))) & |
|---|
| 168 | /(log(table_temp(it-1)) - log(table_temp(it))) |
|---|
| 169 | indtemp = it - 1 |
|---|
| 170 | exit |
|---|
| 171 | end if |
|---|
| 172 | end do |
|---|
| 173 | |
|---|
| 174 | !========================================================================== |
|---|
| 175 | ! loop over vertical levels |
|---|
| 176 | !========================================================================== |
|---|
| 177 | |
|---|
| 178 | do l = 1,lswitch-1 |
|---|
| 179 | |
|---|
| 180 | ! interpolation in air column |
|---|
| 181 | |
|---|
| 182 | indcol = nz - 1 |
|---|
| 183 | do icol = 1,nz |
|---|
| 184 | if (colairtab(icol) < col(l)) then |
|---|
| 185 | indcol = icol - 1 |
|---|
| 186 | exit |
|---|
| 187 | end if |
|---|
| 188 | end do |
|---|
| 189 | cicol = (log(col(l)) - log(colairtab(indcol + 1))) & |
|---|
| 190 | /(log(colairtab(indcol)) - log(colairtab(indcol + 1))) |
|---|
| 191 | |
|---|
| 192 | ! interpolation in ozone column |
|---|
| 193 | |
|---|
| 194 | indozo = nozo - 1 |
|---|
| 195 | do iozo = 1,nozo |
|---|
| 196 | if (table_ozo(iozo)*10. >= ratio_o3(l)) then |
|---|
| 197 | indozo = iozo - 1 |
|---|
| 198 | indozo = max(indozo, 1) |
|---|
| 199 | exit |
|---|
| 200 | end if |
|---|
| 201 | end do |
|---|
| 202 | ciozo = (ratio_o3(l) - table_ozo(indozo)*10.) & |
|---|
| 203 | /(table_ozo(indozo + 1)*10. - table_ozo(indozo)*10.) |
|---|
| 204 | |
|---|
| 205 | ! 4-dimensional interpolation weights |
|---|
| 206 | |
|---|
| 207 | ! poids(temp,sza,co2,o3,tau) |
|---|
| 208 | |
|---|
| 209 | poids(1,1,1,1,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau) |
|---|
| 210 | poids(1,1,1,2,1) = citemp*(1.-cisza)*cicol*ciozo*(1.-citau) |
|---|
| 211 | poids(1,1,2,1,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau) |
|---|
| 212 | poids(1,1,2,2,1) = citemp*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau) |
|---|
| 213 | poids(1,2,1,1,1) = citemp*cisza*cicol*(1.-ciozo)*(1.-citau) |
|---|
| 214 | poids(1,2,1,2,1) = citemp*cisza*cicol*ciozo*(1.-citau) |
|---|
| 215 | poids(1,2,2,1,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau) |
|---|
| 216 | poids(1,2,2,2,1) = citemp*cisza*(1.-cicol)*ciozo*(1.-citau) |
|---|
| 217 | poids(2,1,1,1,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau) |
|---|
| 218 | poids(2,1,1,2,1) = (1.-citemp)*(1.-cisza)*cicol*ciozo*(1.-citau) |
|---|
| 219 | poids(2,1,2,1,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau) |
|---|
| 220 | poids(2,1,2,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau) |
|---|
| 221 | poids(2,2,1,1,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)*(1.-citau) |
|---|
| 222 | poids(2,2,1,2,1) = (1.-citemp)*cisza*cicol*ciozo*(1.-citau) |
|---|
| 223 | poids(2,2,2,1,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau) |
|---|
| 224 | poids(2,2,2,2,1) = (1.-citemp)*cisza*(1.-cicol)*ciozo*(1.-citau) |
|---|
| 225 | ! |
|---|
| 226 | poids(1,1,1,1,2) = citemp*(1.-cisza)*cicol*(1.-ciozo)*citau |
|---|
| 227 | poids(1,1,1,2,2) = citemp*(1.-cisza)*cicol*ciozo*citau |
|---|
| 228 | poids(1,1,2,1,2) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau |
|---|
| 229 | poids(1,1,2,2,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo*citau |
|---|
| 230 | poids(1,2,1,1,2) = citemp*cisza*cicol*(1.-ciozo)*citau |
|---|
| 231 | poids(1,2,1,2,2) = citemp*cisza*cicol*ciozo*citau |
|---|
| 232 | poids(1,2,2,1,2) = citemp*cisza*(1.-cicol)*(1.-ciozo)*citau |
|---|
| 233 | poids(1,2,2,2,2) = citemp*cisza*(1.-cicol)*ciozo*citau |
|---|
| 234 | poids(2,1,1,1,2) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*citau |
|---|
| 235 | poids(2,1,1,2,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo*citau |
|---|
| 236 | poids(2,1,2,1,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau |
|---|
| 237 | poids(2,1,2,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*citau |
|---|
| 238 | poids(2,2,1,1,2) = (1.-citemp)*cisza*cicol*(1.-ciozo)*citau |
|---|
| 239 | poids(2,2,1,2,2) = (1.-citemp)*cisza*cicol*ciozo*citau |
|---|
| 240 | poids(2,2,2,1,2) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*citau |
|---|
| 241 | poids(2,2,2,2,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo*citau |
|---|
| 242 | |
|---|
| 243 | ! 4-dimensional interpolation in the lookup table |
|---|
| 244 | |
|---|
| 245 | do ij = 1,nd |
|---|
| 246 | j(l,ij) = & |
|---|
| 247 | poids(1,1,1,1,1)*jphot(indtemp,indsza,indcol,indozo,indtau,ij) & |
|---|
| 248 | + poids(1,1,1,2,1)*jphot(indtemp,indsza,indcol,indozo+1,indtau,ij) & |
|---|
| 249 | + poids(1,1,2,1,1)*jphot(indtemp,indsza,indcol+1,indozo,indtau,ij) & |
|---|
| 250 | + poids(1,1,2,2,1)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau,ij) & |
|---|
| 251 | + poids(1,2,1,1,1)*jphot(indtemp,indsza+1,indcol,indozo,indtau,ij) & |
|---|
| 252 | + poids(1,2,1,2,1)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau,ij) & |
|---|
| 253 | + poids(1,2,2,1,1)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau,ij) & |
|---|
| 254 | + poids(1,2,2,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,ij) & |
|---|
| 255 | + poids(2,1,1,1,1)*jphot(indtemp+1,indsza,indcol,indozo,indtau,ij) & |
|---|
| 256 | + poids(2,1,1,2,1)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau,ij) & |
|---|
| 257 | + poids(2,1,2,1,1)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau,ij) & |
|---|
| 258 | + poids(2,1,2,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,ij) & |
|---|
| 259 | + poids(2,2,1,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau,ij) & |
|---|
| 260 | + poids(2,2,1,2,1)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,ij) & |
|---|
| 261 | + poids(2,2,2,1,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,ij) & |
|---|
| 262 | + poids(2,2,2,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,ij) & |
|---|
| 263 | ! |
|---|
| 264 | + poids(1,1,1,1,2)*jphot(indtemp,indsza,indcol,indozo,indtau+1,ij) & |
|---|
| 265 | + poids(1,1,1,2,2)*jphot(indtemp,indsza,indcol,indozo+1,indtau+1,ij) & |
|---|
| 266 | + poids(1,1,2,1,2)*jphot(indtemp,indsza,indcol+1,indozo,indtau+1,ij) & |
|---|
| 267 | + poids(1,1,2,2,2)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,ij) & |
|---|
| 268 | + poids(1,2,1,1,2)*jphot(indtemp,indsza+1,indcol,indozo,indtau+1,ij) & |
|---|
| 269 | + poids(1,2,1,2,2)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,ij) & |
|---|
| 270 | + poids(1,2,2,1,2)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,ij) & |
|---|
| 271 | + poids(1,2,2,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,ij) & |
|---|
| 272 | + poids(2,1,1,1,2)*jphot(indtemp+1,indsza,indcol,indozo,indtau+1,ij) & |
|---|
| 273 | + poids(2,1,1,2,2)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,ij) & |
|---|
| 274 | + poids(2,1,2,1,2)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,ij) & |
|---|
| 275 | + poids(2,1,2,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,ij) & |
|---|
| 276 | + poids(2,2,1,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,ij) & |
|---|
| 277 | + poids(2,2,1,2,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,ij) & |
|---|
| 278 | + poids(2,2,2,1,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,ij) & |
|---|
| 279 | + poids(2,2,2,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,ij) |
|---|
| 280 | end do |
|---|
| 281 | |
|---|
| 282 | ! correction for sun distance |
|---|
| 283 | |
|---|
| 284 | do ij = 1,nd |
|---|
| 285 | j(l,ij) = j(l,ij)*(1.52/dist_sol)**2. |
|---|
| 286 | end do |
|---|
| 287 | |
|---|
| 288 | !========================================================================== |
|---|
| 289 | ! end of loop over vertical levels |
|---|
| 290 | !========================================================================== |
|---|
| 291 | |
|---|
| 292 | end do |
|---|
| 293 | |
|---|
| 294 | else |
|---|
| 295 | |
|---|
| 296 | !========================================================================== |
|---|
| 297 | ! night |
|---|
| 298 | !========================================================================== |
|---|
| 299 | |
|---|
| 300 | j(:,:) = 0. |
|---|
| 301 | |
|---|
| 302 | end if |
|---|
| 303 | |
|---|
| 304 | ! photodissociation rates numbering in the lookup table |
|---|
| 305 | |
|---|
| 306 | ! jmars.20140930 |
|---|
| 307 | |
|---|
| 308 | j_o2_o = 1 ! o2 + hv -> o + o |
|---|
| 309 | j_o2_o1d = 2 ! o2 + hv -> o + o(1d) |
|---|
| 310 | j_co2_o = 3 ! co2 + hv -> co + o |
|---|
| 311 | j_co2_o1d = 4 ! co2 + hv -> co + o(1d) |
|---|
| 312 | j_o3_o1d = 5 ! o3 + hv -> o2 + o(1d) |
|---|
| 313 | j_o3_o = 6 ! o3 + hv -> o2 + o |
|---|
| 314 | j_h2o = 7 ! h2o + hv -> h + oh |
|---|
| 315 | j_h2o2 = 8 ! h2o2 + hv -> oh + oh |
|---|
| 316 | j_ho2 = 9 ! ho2 + hv -> oh + o |
|---|
| 317 | j_no = 10 ! no + hv -> n + o |
|---|
| 318 | j_no2 = 11 ! no2 + hv -> no + o |
|---|
| 319 | j_hno3 = 12 ! hno3 + hv -> no2 + oh |
|---|
| 320 | j_hno4 = 13 ! hno4 + hv -> no2 + ho2 |
|---|
| 321 | |
|---|
| 322 | ! fill v_phot array |
|---|
| 323 | |
|---|
| 324 | v_phot(:,:) = 0. |
|---|
| 325 | |
|---|
| 326 | do l = 1,lswitch-1 |
|---|
| 327 | v_phot(l, 1) = j(l,j_o2_o) |
|---|
| 328 | v_phot(l, 2) = j(l,j_o2_o1d) |
|---|
| 329 | v_phot(l, 3) = j(l,j_co2_o) |
|---|
| 330 | v_phot(l, 4) = j(l,j_co2_o1d) |
|---|
| 331 | v_phot(l, 5) = j(l,j_o3_o1d) |
|---|
| 332 | v_phot(l, 6) = j(l,j_o3_o) |
|---|
| 333 | v_phot(l, 7) = j(l,j_h2o) |
|---|
| 334 | v_phot(l, 8) = j(l,j_h2o2) |
|---|
| 335 | v_phot(l, 9) = j(l,j_ho2) |
|---|
| 336 | v_phot(l,10) = 0. ! h2 missing in lookup table |
|---|
| 337 | v_phot(l,11) = j(l,j_no) |
|---|
| 338 | v_phot(l,12) = j(l,j_no2) |
|---|
| 339 | v_phot(l,13) = 0. ! n2 missing in lookup table |
|---|
| 340 | end do |
|---|
| 341 | |
|---|
| 342 | return |
|---|
| 343 | end |
|---|
| 344 | |
|---|
| 345 | !***************************************************************** |
|---|
