| [1047] | 1 |       SUBROUTINE concentrations(ngrid,nlayer,nq, | 
|---|
 | 2 |      &                          pplay,pt,pdt,pq,pdq,ptimestep) | 
|---|
| [38] | 3 |                                               | 
|---|
| [1036] | 4 |       use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d, | 
|---|
 | 5 |      &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, | 
|---|
 | 6 |      &                      igcm_oh, igcm_ho2, igcm_n2, igcm_ar, | 
|---|
 | 7 |      &                      igcm_h2o_vap, igcm_n, igcm_no, igcm_no2, | 
|---|
 | 8 |      &                      igcm_n2d, igcm_co2plus, igcm_oplus, | 
|---|
 | 9 |      &                      igcm_o2plus, igcm_coplus, igcm_cplus, | 
|---|
 | 10 |      &                      igcm_nplus, igcm_noplus, igcm_n2plus, | 
|---|
| [1660] | 11 |      &                      igcm_hplus, igcm_hco2plus, mmol, | 
|---|
| [2158] | 12 |      &                      igcm_he, igcm_elec | 
|---|
| [1047] | 13 |       use conc_mod, only: mmean, Akknew, rnew, cpnew | 
|---|
| [370] | 14 |       implicit none | 
|---|
 | 15 |  | 
|---|
| [618] | 16 | !======================================================================= | 
|---|
 | 17 | ! CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R | 
|---|
 | 18 | ! | 
|---|
| [1047] | 19 | ! mmean(ngrid,nlayer)   amu | 
|---|
 | 20 | ! cpnew(ngrid,nlayer)   J/kg/K | 
|---|
 | 21 | ! rnew(ngrid,nlayer)    J/kg/K | 
|---|
 | 22 | ! akknew(ngrid,nlayer)  coefficient of thermal concduction | 
|---|
| [618] | 23 | ! | 
|---|
 | 24 | ! version: April 2012 - Franck Lefevre | 
|---|
 | 25 | !======================================================================= | 
|---|
| [370] | 26 |  | 
|---|
| [618] | 27 | !     declarations | 
|---|
| [370] | 28 |   | 
|---|
| [38] | 29 | #include "callkeys.h" | 
|---|
 | 30 |  | 
|---|
| [618] | 31 | !     input/output | 
|---|
| [38] | 32 |  | 
|---|
| [1047] | 33 |       integer,intent(in) :: ngrid ! number of atmospheric columns | 
|---|
 | 34 |       integer,intent(in) :: nlayer ! number of atmospheric layers | 
|---|
| [1036] | 35 |       integer,intent(in) :: nq ! number of tracers | 
|---|
| [1047] | 36 |       real,intent(in) :: pplay(ngrid,nlayer) | 
|---|
 | 37 |       real,intent(in) :: pt(ngrid,nlayer) | 
|---|
 | 38 |       real,intent(in) :: pdt(ngrid,nlayer) | 
|---|
 | 39 |       real,intent(in) :: pq(ngrid,nlayer,nq) | 
|---|
 | 40 |       real,intent(in) :: pdq(ngrid,nlayer,nq) | 
|---|
| [1036] | 41 |       real,intent(in) :: ptimestep | 
|---|
| [38] | 42 |  | 
|---|
| [618] | 43 | !     local variables | 
|---|
| [38] | 44 |  | 
|---|
| [618] | 45 |       integer       :: i, l, ig, iq | 
|---|
| [1036] | 46 |       integer, save :: nbq | 
|---|
 | 47 |       integer,allocatable,save :: niq(:) | 
|---|
 | 48 |       real          :: ni(nq), ntot | 
|---|
| [1047] | 49 |       real          :: zq(ngrid, nlayer, nq) | 
|---|
 | 50 |       real          :: zt(ngrid, nlayer) | 
|---|
| [1036] | 51 |       real,allocatable,save    :: aki(:) | 
|---|
 | 52 |       real,allocatable,save    :: cpi(:) | 
|---|
| [370] | 53 |  | 
|---|
 | 54 |       logical, save :: firstcall = .true. | 
|---|
 | 55 |  | 
|---|
| [38] | 56 |       if (firstcall) then | 
|---|
 | 57 |  | 
|---|
| [1036] | 58 |          ! allocate local saved arrays: | 
|---|
 | 59 |          allocate(aki(nq)) | 
|---|
 | 60 |          allocate(cpi(nq)) | 
|---|
 | 61 |          allocate(niq(nq)) | 
|---|
| [618] | 62 | !        find index of chemical tracers to use | 
|---|
 | 63 | !        initialize thermal conductivity and specific heat coefficients | 
|---|
 | 64 | !        !? values are estimated | 
|---|
| [38] | 65 |  | 
|---|
| [618] | 66 |          nbq = 0 ! to count number of tracers used in this subroutine | 
|---|
| [370] | 67 |  | 
|---|
| [618] | 68 |          if (igcm_co2 /= 0) then | 
|---|
 | 69 |             nbq = nbq + 1 | 
|---|
 | 70 |             niq(nbq) = igcm_co2 | 
|---|
 | 71 |             aki(nbq) = 3.072e-4 | 
|---|
 | 72 |             cpi(nbq) = 0.834e3 | 
|---|
 | 73 |          end if | 
|---|
 | 74 |          if (igcm_co /= 0) then | 
|---|
 | 75 |             nbq = nbq + 1 | 
|---|
 | 76 |             niq(nbq) = igcm_co | 
|---|
 | 77 |             aki(nbq) = 4.87e-4 | 
|---|
 | 78 |             cpi(nbq) = 1.034e3 | 
|---|
 | 79 |          end if | 
|---|
 | 80 |          if (igcm_o /= 0) then | 
|---|
 | 81 |             nbq = nbq + 1 | 
|---|
 | 82 |             niq(nbq) = igcm_o | 
|---|
 | 83 |             aki(nbq) = 7.59e-4 | 
|---|
 | 84 |             cpi(nbq) = 1.3e3 | 
|---|
 | 85 |          end if | 
|---|
 | 86 |          if (igcm_o1d /= 0) then | 
|---|
 | 87 |             nbq = nbq + 1 | 
|---|
 | 88 |             niq(nbq) = igcm_o1d | 
|---|
 | 89 |             aki(nbq) = 7.59e-4  !? | 
|---|
 | 90 |             cpi(nbq) = 1.3e3    !? | 
|---|
 | 91 |          end if | 
|---|
 | 92 |          if (igcm_o2 /= 0) then | 
|---|
 | 93 |             nbq = nbq + 1 | 
|---|
 | 94 |             niq(nbq) = igcm_o2 | 
|---|
 | 95 |             aki(nbq) = 5.68e-4 | 
|---|
 | 96 |             cpi(nbq) = 0.9194e3 | 
|---|
 | 97 |          end if | 
|---|
 | 98 |          if (igcm_o3 /= 0) then | 
|---|
 | 99 |             nbq = nbq + 1 | 
|---|
 | 100 |             niq(nbq) = igcm_o3 | 
|---|
 | 101 |             aki(nbq) = 3.00e-4  !? | 
|---|
 | 102 |             cpi(nbq) = 0.800e3  !? | 
|---|
 | 103 |          end if | 
|---|
 | 104 |          if (igcm_h /= 0) then | 
|---|
 | 105 |             nbq = nbq + 1 | 
|---|
 | 106 |             niq(nbq) = igcm_h | 
|---|
 | 107 |             aki(nbq) = 0.0 | 
|---|
 | 108 |             cpi(nbq) = 20.780e3 | 
|---|
 | 109 |          end if | 
|---|
 | 110 |          if (igcm_h2 /= 0) then | 
|---|
 | 111 |             nbq = nbq + 1 | 
|---|
 | 112 |             niq(nbq) = igcm_h2 | 
|---|
 | 113 |             aki(nbq) = 36.314e-4 | 
|---|
 | 114 |             cpi(nbq) = 14.266e3 | 
|---|
 | 115 |          end if | 
|---|
 | 116 |          if (igcm_oh /= 0) then | 
|---|
 | 117 |             nbq = nbq + 1 | 
|---|
 | 118 |             niq(nbq) = igcm_oh | 
|---|
 | 119 |             aki(nbq)  = 7.00e-4 !? | 
|---|
 | 120 |             cpi(nbq)  = 1.045e3 | 
|---|
 | 121 |          end if | 
|---|
 | 122 |          if (igcm_ho2 /= 0) then | 
|---|
 | 123 |             nbq = nbq + 1 | 
|---|
 | 124 |             niq(nbq) = igcm_ho2 | 
|---|
 | 125 |             aki(nbq) = 0.0 | 
|---|
 | 126 |             cpi(nbq) = 1.065e3  !? | 
|---|
 | 127 |          end if | 
|---|
 | 128 |          if (igcm_n2 /= 0) then | 
|---|
 | 129 |             nbq = nbq + 1 | 
|---|
 | 130 |             niq(nbq) = igcm_n2 | 
|---|
 | 131 |             aki(nbq) = 5.6e-4 | 
|---|
 | 132 |             cpi(nbq) = 1.034e3 | 
|---|
 | 133 |          end if | 
|---|
 | 134 |          if (igcm_ar /= 0) then | 
|---|
 | 135 |             nbq = nbq + 1 | 
|---|
 | 136 |             niq(nbq) = igcm_ar | 
|---|
| [1660] | 137 |             aki(nbq) = 3.4e-4 | 
|---|
 | 138 |             cpi(nbq) = 5.2e2 | 
|---|
 | 139 |             ! aki/cpi from Leslie A. Guildner, | 
|---|
 | 140 |             ! JOURNAL OF RESEARCH of the National Bureau of Standards-  | 
|---|
 | 141 |             ! A. Physics and Chemistry | 
|---|
 | 142 |             ! Vol. 79A, No.2, March-April 1975 | 
|---|
| [618] | 143 |          end if | 
|---|
 | 144 |          if (igcm_h2o_vap /= 0) then | 
|---|
 | 145 |             nbq = nbq + 1 | 
|---|
 | 146 |             niq(nbq) = igcm_h2o_vap | 
|---|
 | 147 |             aki(nbq) = 0.0 | 
|---|
 | 148 |             cpi(nbq) = 1.870e3 | 
|---|
 | 149 |          end if | 
|---|
| [635] | 150 |          if (igcm_n /= 0) then | 
|---|
 | 151 |             nbq = nbq + 1 | 
|---|
 | 152 |             niq(nbq) = igcm_n | 
|---|
 | 153 |             aki(nbq) = 0.0 | 
|---|
 | 154 |             cpi(nbq) = 0.0 | 
|---|
 | 155 |          endif | 
|---|
 | 156 |          if(igcm_no /= 0) then | 
|---|
 | 157 |             nbq = nbq + 1 | 
|---|
 | 158 |             niq(nbq) = igcm_no | 
|---|
 | 159 |             aki(nbq) = 0.0 | 
|---|
 | 160 |             cpi(nbq) = 0.0 | 
|---|
 | 161 |          endif | 
|---|
 | 162 |          if(igcm_no2 /= 0) then | 
|---|
 | 163 |             nbq = nbq + 1 | 
|---|
 | 164 |             niq(nbq) = igcm_no2 | 
|---|
 | 165 |             aki(nbq) = 0.0 | 
|---|
 | 166 |             cpi(nbq) = 0.0 | 
|---|
 | 167 |          endif | 
|---|
 | 168 |          if(igcm_n2d /= 0) then | 
|---|
 | 169 |             nbq = nbq + 1 | 
|---|
 | 170 |             niq(nbq) = igcm_n2d | 
|---|
 | 171 |             aki(nbq) = 0.0 | 
|---|
 | 172 |             cpi(nbq) = 0.0 | 
|---|
 | 173 |          endif | 
|---|
| [1660] | 174 |          if(igcm_he /= 0) then | 
|---|
 | 175 |             nbq = nbq + 1 | 
|---|
 | 176 |             niq(nbq) = igcm_he | 
|---|
 | 177 |             aki(nbq) = 30.e-4 | 
|---|
 | 178 |             cpi(nbq) = 5.2e3 | 
|---|
 | 179 |             ! aki/cpi from Leslie A. Guildner, | 
|---|
 | 180 |             ! JOURNAL OF RESEARCH of the National Bureau of Standards-  | 
|---|
 | 181 |             ! A. Physics and Chemistry | 
|---|
 | 182 |             ! Vol. 79A, No.2, March-April 1975 | 
|---|
 | 183 |          endif | 
|---|
| [635] | 184 |          if(igcm_co2plus /= 0) then | 
|---|
 | 185 |             nbq = nbq + 1 | 
|---|
 | 186 |             niq(nbq) = igcm_co2plus | 
|---|
 | 187 |             aki(nbq) = 0.0 | 
|---|
 | 188 |             cpi(nbq) = 0.0 | 
|---|
 | 189 |          endif | 
|---|
 | 190 |          if(igcm_oplus /= 0) then | 
|---|
 | 191 |             nbq = nbq + 1 | 
|---|
 | 192 |             niq(nbq) = igcm_oplus | 
|---|
 | 193 |             aki(nbq) = 0.0 | 
|---|
 | 194 |             cpi(nbq) = 0.0 | 
|---|
 | 195 |          endif | 
|---|
 | 196 |          if(igcm_o2plus /= 0) then | 
|---|
 | 197 |             nbq = nbq + 1 | 
|---|
 | 198 |             niq(nbq) = igcm_o2plus | 
|---|
 | 199 |             aki(nbq) = 0.0 | 
|---|
 | 200 |             cpi(nbq) = 0.0 | 
|---|
 | 201 |          endif | 
|---|
 | 202 |          if(igcm_coplus /= 0) then | 
|---|
 | 203 |             nbq = nbq + 1 | 
|---|
 | 204 |             niq(nbq) = igcm_coplus | 
|---|
 | 205 |             aki(nbq) = 0.0 | 
|---|
 | 206 |             cpi(nbq) = 0.0 | 
|---|
 | 207 |          endif | 
|---|
 | 208 |          if(igcm_cplus /= 0) then | 
|---|
 | 209 |             nbq = nbq + 1 | 
|---|
 | 210 |             niq(nbq) = igcm_cplus | 
|---|
 | 211 |             aki(nbq) = 0.0 | 
|---|
 | 212 |             cpi(nbq) = 0.0 | 
|---|
 | 213 |          endif | 
|---|
 | 214 |          if(igcm_nplus /= 0) then | 
|---|
 | 215 |             nbq = nbq + 1 | 
|---|
 | 216 |             niq(nbq) = igcm_nplus | 
|---|
 | 217 |             aki(nbq) = 0.0 | 
|---|
 | 218 |             cpi(nbq) = 0.0 | 
|---|
 | 219 |          endif | 
|---|
 | 220 |          if(igcm_noplus /= 0) then | 
|---|
 | 221 |             nbq = nbq + 1 | 
|---|
 | 222 |             niq(nbq) = igcm_noplus | 
|---|
 | 223 |             aki(nbq) = 0.0 | 
|---|
 | 224 |             cpi(nbq) = 0.0 | 
|---|
 | 225 |          endif | 
|---|
 | 226 |          if(igcm_n2plus /= 0) then | 
|---|
 | 227 |             nbq = nbq + 1 | 
|---|
 | 228 |             niq(nbq) = igcm_n2plus | 
|---|
 | 229 |             aki(nbq) = 0.0 | 
|---|
 | 230 |             cpi(nbq) = 0.0 | 
|---|
 | 231 |          endif | 
|---|
 | 232 |          if(igcm_hplus /= 0) then | 
|---|
 | 233 |             nbq = nbq + 1 | 
|---|
 | 234 |             niq(nbq) = igcm_hplus | 
|---|
 | 235 |             aki(nbq) = 0.0 | 
|---|
 | 236 |             cpi(nbq) = 0.0 | 
|---|
 | 237 |          endif | 
|---|
 | 238 |          if(igcm_hco2plus /= 0) then | 
|---|
 | 239 |             nbq = nbq + 1 | 
|---|
 | 240 |             niq(nbq) = igcm_hco2plus | 
|---|
 | 241 |             aki(nbq) = 0.0 | 
|---|
 | 242 |             cpi(nbq) = 0.0 | 
|---|
 | 243 |          endif | 
|---|
| [2158] | 244 |          if(igcm_elec /= 0) then | 
|---|
 | 245 |             nbq = nbq + 1 | 
|---|
 | 246 |             niq(nbq) = igcm_elec | 
|---|
 | 247 |             aki(nbq) = 0.0 | 
|---|
 | 248 |             cpi(nbq) = 0.0 | 
|---|
 | 249 |          endif | 
|---|
| [635] | 250 |           | 
|---|
| [1036] | 251 |          ! tell the world about it: | 
|---|
 | 252 |          write(*,*) "concentrations: firstcall, nbq=",nbq | 
|---|
 | 253 |          write(*,*) "  niq(1:nbq)=",niq(1:nbq) | 
|---|
 | 254 |          write(*,*) "  aki(1:nbq)=",aki(1:nbq) | 
|---|
 | 255 |          write(*,*) "  cpi(1:nbq)=",cpi(1:nbq) | 
|---|
| [38] | 256 |  | 
|---|
| [618] | 257 |          firstcall = .false. | 
|---|
| [38] | 258 |  | 
|---|
| [618] | 259 |       end if ! if (firstcall) | 
|---|
| [38] | 260 |  | 
|---|
| [618] | 261 | !     update temperature | 
|---|
| [38] | 262 |  | 
|---|
| [1047] | 263 |       do l = 1,nlayer | 
|---|
 | 264 |          do ig = 1,ngrid | 
|---|
| [370] | 265 |             zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep | 
|---|
 | 266 |          end do | 
|---|
 | 267 |       end do | 
|---|
| [618] | 268 |  | 
|---|
 | 269 | !     update tracers | 
|---|
 | 270 |  | 
|---|
| [1047] | 271 |       do l = 1,nlayer | 
|---|
 | 272 |          do ig = 1,ngrid | 
|---|
| [618] | 273 |             do i = 1,nbq | 
|---|
 | 274 |                iq = niq(i)  | 
|---|
 | 275 |                zq(ig,l,iq) = max(1.e-30, pq(ig,l,iq) | 
|---|
 | 276 |      $                                 + pdq(ig,l,iq)*ptimestep) | 
|---|
| [370] | 277 |             end do | 
|---|
 | 278 |          end do | 
|---|
 | 279 |       end do | 
|---|
| [618] | 280 |  | 
|---|
 | 281 | !     mmean : mean molecular mass | 
|---|
 | 282 | !     rnew  : specific gas constant | 
|---|
 | 283 |  | 
|---|
 | 284 |       mmean(:,:)  = 0. | 
|---|
 | 285 |  | 
|---|
| [1047] | 286 |       do l = 1,nlayer | 
|---|
 | 287 |          do ig = 1,ngrid | 
|---|
| [618] | 288 |             do i = 1,nbq | 
|---|
 | 289 |                iq = niq(i)  | 
|---|
 | 290 |                mmean(ig,l) = mmean(ig,l) + zq(ig,l,iq)/mmol(iq) | 
|---|
| [370] | 291 |             end do | 
|---|
 | 292 |             mmean(ig,l) = 1./mmean(ig,l) | 
|---|
 | 293 |             rnew(ig,l) = 8.314/mmean(ig,l)*1.e3     ! J/kg/K             | 
|---|
 | 294 |          end do | 
|---|
 | 295 |       end do | 
|---|
| [618] | 296 |  | 
|---|
 | 297 | !     cpnew  : specicic heat | 
|---|
 | 298 | !     akknew : thermal conductivity cofficient | 
|---|
 | 299 |        | 
|---|
 | 300 |       cpnew(:,:)  = 0. | 
|---|
 | 301 |       akknew(:,:) = 0. | 
|---|
 | 302 |  | 
|---|
| [1047] | 303 |       do l = 1,nlayer | 
|---|
 | 304 |          do ig = 1,ngrid | 
|---|
| [370] | 305 |             ntot = pplay(ig,l)/(1.381e-23*zt(ig,l))*1.e-6  ! in #/cm3 | 
|---|
| [618] | 306 |             do i = 1,nbq | 
|---|
 | 307 |                iq = niq(i)  | 
|---|
 | 308 |                ni(iq) = ntot*zq(ig,l,iq)*mmean(ig,l)/mmol(iq) | 
|---|
| [1035] | 309 |                cpnew(ig,l) = cpnew(ig,l) + ni(iq)*cpi(i) | 
|---|
 | 310 |                akknew(ig,l) = akknew(ig,l) + ni(iq)*aki(i) | 
|---|
| [370] | 311 |             end do  | 
|---|
 | 312 |             cpnew(ig,l) = cpnew(ig,l)/ntot | 
|---|
 | 313 |             akknew(ig,l) = akknew(ig,l)/ntot | 
|---|
 | 314 |          end do | 
|---|
| [618] | 315 | !        print*, l, mmean(1,l), cpnew(1,l), rnew(1,l) | 
|---|
| [370] | 316 |       end do | 
|---|
| [38] | 317 |  | 
|---|
 | 318 |       return | 
|---|
| [635] | 319 |       end  | 
|---|