chem/module_radm.F

   1     MODULE module_radm
   2   USE module_data_radm2
   3   USE module_data_sorgam
   4     integer numchem
   5     parameter (numchem=numchem_radm)
   6 ! ..
   7     CONTAINS
   8       subroutine radm_driver(id,curr_secs,dtstep,config_flags,         &
   9                gmt,julday,t_phy,moist,p8w,t8w,                         &
  10                p_phy,chem,rho_phy,dz8w,z,z_at_w,vdrog3,         &
  11                ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,  &
  12                ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho,    &
  13                ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho,         &
  14                ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,&
  15                ids,ide, jds,jde, kds,kde,                              &
  16                ims,ime, jms,jme, kms,kme,                              &
  17                its,ite, jts,jte, kts,kte                               )
  18   USE module_configure
  19   USE module_state_description
  20   USE module_model_constants
  21    IMPLICIT NONE
  22
  23    INTEGER,      INTENT(IN   ) :: id,julday,                           &
  24                                   ids,ide, jds,jde, kds,kde,           &
  25                                   ims,ime, jms,jme, kms,kme,           &
  26                                   its,ite, jts,jte, kts,kte
  27    REAL(KIND=8), INTENT(IN   ) :: curr_secs
  28    REAL,         INTENT(IN   ) :: dtstep,gmt
  29 !
  30 ! advected moisture variables
  31 !
  32    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ),            &
  33          INTENT(IN ) ::                                   moist
  34 !
  35 !  advected chemical tracers
  36 !
  37    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),             &
  38          INTENT(INOUT ) ::                                chem
  39
  40 !
  41 ! arrays that hold photolysis rates
  42 !
  43    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                       &
  44          INTENT(INOUT ) ::                                             &
  45            ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,      &
  46            ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho,        &
  47            ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho,             &
  48            ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob
  49 !
  50 ! on input from met model
  51 !
  52    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         ,           &
  53           INTENT(IN   ) ::                                             &
  54                                                       t_phy,           &
  55                                                       p_phy,           &
  56                                                       dz8w,            &
  57                                                       z    ,           &
  58                                               t8w,p8w,z_at_w ,         &
  59                                                     rho_phy
  60 !
  61 ! for interaction with aerosols (really is output)
  62 !
  63       real ,      INTENT(INOUT) ::                                     &
  64                       vdrog3(ims:ime,kms:kme-0,jms:jme,ldrog)
  65       TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
  66
  67 ! ..
  68 ! .. Local Scalars ..
  69       REAL ::  clwchem,  dt60, dtcmax, dtcmin
  70       INTEGER :: i,j,k,iprt, jce, jcs, n, nr, ipr,jpr,nvr
  71 ! ..
  72 ! .. Local Arrays ..
  73       REAL :: p(kts:kte), rh(kts:kte), rj(kts:kte,nreacj),       &
  74         t(kts:kte), vcinp(kts:kte,numchem),wlc(kts:kte)
  75       real :: vdrog1(kts:kte,ldrog)
  76
  77 !
  78       integer(kind=8) :: ixhour
  79       INTEGER :: iaerosol_sorgam
  80       real(kind=8) :: xhour, xtime, xtimin
  81       real :: xmin
  82       xtime=curr_secs/60._8
  83       ixhour=int(gmt+.01,8)+int(xtime/60._8,8)
  84       xhour=real(ixhour,8)
  85       xmin=60.*gmt+real(xtime-xhour*60._8,8)
  86
  87       ipr=-10
  88       jpr=-10
  89       nvr=5
  90 !
  91 ! following is for combination radm/sorgam only, p_nu0 must be defined
  92 ! in that case
  93 !
  94       iaerosol_sorgam=0
  95       if(p_nu0.gt.1)iaerosol_sorgam=1
  96
  97 !
  98       chem=max(chem,epsilc)
  99       do 100 j=jts,jte
 100       do 100 i=its,ite
 101       vcinp=epsilc
 102       vdrog1=0.
 103       iprt=0
 104 !     if(xtime/60.ge.2.)then
 105 !     if((i.eq.12.and.j.eq.17).or.       &
 106 !         (i.eq.12.and.j.eq.7).or.       &
 107 !         (i.eq.1.and.j.eq.17))iprt=2
 108 !     endif
 109
 110 ! reorder
 111 !
 112 !       if(iprt.eq.2)print *,'k,chem(i,k,j,p_sulf),vcinp(k,lsulf)'
 113       do k=kts,kte
 114         vcinp(k,lso2)  =  max(chem(i,k,j,p_so2),epsilc)
 115         vcinp(k,Lsulf) =  max(chem(i,k,j,p_sulf),epsilc)
 116         vcinp(k,Lno2)  =  max(chem(i,k,j,p_no2),epsilc)
 117         vcinp(k,Lno)   =  max(chem(i,k,j,p_no),1.e-6)
 118 !       vcinp(k,Lno)   =  max(chem(i,k,j,p_no),epsilc)
 119         vcinp(k,Lo3)   =  max(chem(i,k,j,p_o3),epsilc)
 120         vcinp(k,Lhno3) =  max(chem(i,k,j,p_hno3),epsilc)
 121         vcinp(k,Lh2o2) =  max(chem(i,k,j,p_h2o2),epsilc)
 122         vcinp(k,Lald)  =  max(chem(i,k,j,p_ald),epsilc)
 123         vcinp(k,Lhcho) =  max(chem(i,k,j,p_hcho),epsilc)
 124         vcinp(k,Lop1)  =  max(chem(i,k,j,p_op1),epsilc)
 125         vcinp(k,Lop2)  =  max(chem(i,k,j,p_op2),epsilc)
 126         vcinp(k,Lpaa)  =  max(chem(i,k,j,p_paa),epsilc)
 127         vcinp(k,Lora1)  =  max(chem(i,k,j,p_ora1),epsilc)
 128         vcinp(k,Lora2)  =  max(chem(i,k,j,p_ora2),epsilc)
 129         vcinp(k,Lnh3)  =  max(chem(i,k,j,p_nh3),epsilc)
 130         vcinp(k,Ln2o5)  =  max(chem(i,k,j,p_n2o5),epsilc)
 131         vcinp(k,Lno3)  =  max(chem(i,k,j,p_no3),epsilc)
 132         vcinp(k,Lpan)  =  max(chem(i,k,j,p_pan),epsilc)
 133         vcinp(k,Lhc3)  =  max(chem(i,k,j,p_hc3),epsilc)
 134         vcinp(k,Lhc5)  =  max(chem(i,k,j,p_hc5),epsilc)
 135         vcinp(k,Lhc8)  =  max(chem(i,k,j,p_hc8),epsilc)
 136         vcinp(k,Leth)  =  max(chem(i,k,j,p_eth),epsilc)
 137         vcinp(k,Lco)  =  max(chem(i,k,j,p_co),epsilc)
 138         vcinp(k,Lol2)  =  max(chem(i,k,j,p_ol2),epsilc)
 139         vcinp(k,Lolt)  =  max(chem(i,k,j,p_olt),epsilc)
 140         vcinp(k,Loli)  =  max(chem(i,k,j,p_oli),epsilc)
 141         vcinp(k,Ltol)  =  max(chem(i,k,j,p_tol),epsilc)
 142         vcinp(k,Lxyl)  =  max(chem(i,k,j,p_xyl),epsilc)
 143         vcinp(k,Laco3)  =  max(chem(i,k,j,p_aco3),epsilc)
 144         vcinp(k,Ltpan)  =  max(chem(i,k,j,p_tpan),epsilc)
 145         vcinp(k,Lhono)  =  max(chem(i,k,j,p_hono),epsilc)
 146         vcinp(k,Lhno4)  =  max(chem(i,k,j,p_hno4),epsilc)
 147         vcinp(k,Lket)  =  max(chem(i,k,j,p_ket),epsilc)
 148         vcinp(k,Lgly)  =  max(chem(i,k,j,p_gly),epsilc)
 149         vcinp(k,Lmgly)  =  max(chem(i,k,j,p_mgly),epsilc)
 150         vcinp(k,Ldcb)  =  max(chem(i,k,j,p_dcb),epsilc)
 151         vcinp(k,Lonit)  =  max(chem(i,k,j,p_onit),epsilc)
 152         vcinp(k,Lcsl)  =  max(chem(i,k,j,p_csl),epsilc)
 153         vcinp(k,Lxyl)  =  max(chem(i,k,j,p_xyl),epsilc)
 154         vcinp(k,Liso)  =  max(chem(i,k,j,p_iso),epsilc)
 155         vcinp(k,Lho)  =  max(chem(i,k,j,p_ho),epsilc)
 156         vcinp(k,Lho2)  =  max(chem(i,k,j,p_ho2),epsilc)
 157 !       if(iprt.eq.2)then
 158 !        print *,k,chem(i,k,j,p_sulf),vcinp(k,lsulf)
 159 !       endif
 160         enddo
 161 !--- now do chemistry, need some input here
 162
 163       do k=kts,kte
 164          t(k) = t_phy(i,k,j)
 165          p(k) = .001*p_phy(i,k,j)
 166          rh(k) = .95
 167          rh(k) = MIN( .95, moist(i,k,j,p_qv) / &
 168                (3.80*exp(17.27*(t_phy(i,k,j)-273.)/ &
 169                (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))))
 170          rh(k)=max(.1,rh(k))
 171 !        wlc(k) = moist(i,k,j,p_qc)
 172          wlc(k) = 0.
 173       END DO
 174       dt60 = dtstep/60.
 175       xtimin = max(0._8,xtime-real(dt60,8))
 176       dtcmin = min(.05_8,xtime-xtimin)
 177       dtcmin = max(dtcmin,0.5/60.)
 178       dtcmax = min(5.,dt60)
 179       dtcmax = min(real(dtcmax,8),xtime-xtimin)
 180 !
 181 ! radm here is called with a vertical stack
 182 !
 183       jcs = kts
 184       jce = kte
 185
 186
 187 !
 188 ! fill photolysis rates for use in radm module
 189 !
 190       do k=kts,kte
 191         rj(k,1) = ph_no2(i,k,j)
 192         rj(k,2) = ph_o31d(i,k,j)
 193         rj(k,3) = ph_o33p(i,k,j)
 194         rj(k,4) = ph_hno2(i,k,j)
 195         rj(k,5) = ph_hno3(i,k,j)
 196         rj(k,6) = ph_hno4(i,k,j)
 197         rj(k,7) = ph_no3o2(i,k,j)
 198         rj(k,8) = ph_no3o(i,k,j)
 199         rj(k,9) = ph_h2o2(i,k,j)
 200         rj(k,10) = ph_ch2om(i,k,j)
 201         rj(k,11) = ph_ch2or(i,k,j)
 202         rj(k,12) = ph_ch3cho(i,k,j)
 203         rj(k,13) = ph_ch3o2h(i,k,j)
 204         rj(k,14) = ph_ch3coch3(i,k,j)
 205         rj(k,15) = ph_ch3coo2h(i,k,j)
 206         rj(k,16) = ph_ch3coc2h5(i,k,j)
 207         rj(k,17) = ph_hcocho(i,k,j)
 208         rj(k,18) = ph_hcochob(i,k,j)
 209         rj(k,19) = ph_ch3cocho(i,k,j)
 210         rj(k,20) = ph_hcochest(i,k,j)
 211         rj(k,21) = ph_ch3ono2(i,k,j)
 212       END DO
 213 !     print *,'before radm, i,j = ',i,j
 214 !     iprt=0
 215 !     if((i.eq.87.and.j.eq.15).or.(i.eq.87.and.j.eq.4))then
 216 !       iprt=1
 217 !     endif
 218       CALL radm(rj,wlc,vcinp,t,p,rh,xtime,xtimin,kts,kte,            &
 219                 iprt,dt60,dtcmax,dtcmin,vdrog1,iaerosol_sorgam)
 220 !     print *,'after  radm, i,j = ',i,j
 221       do k=kts,kte
 222         chem(i,k,j,p_so2)     = max(vcinp(k,lso2),epsilc)
 223         chem(i,k,j,p_sulf)     = max(vcinp(k,Lsulf),epsilc)
 224         chem(i,k,j,p_no2)     = max(vcinp(k,Lno2),epsilc)
 225         chem(i,k,j,p_no)     = max(vcinp(k,Lno),1.e-6)
 226         chem(i,k,j,p_o3)     = max(vcinp(k,Lo3),epsilc)
 227         chem(i,k,j,p_hno3)     = max(vcinp(k,Lhno3),epsilc)
 228         chem(i,k,j,p_h2o2)     = max(vcinp(k,Lh2o2),epsilc)
 229         chem(i,k,j,p_ald)     = max(vcinp(k,Lald),epsilc)
 230         chem(i,k,j,p_hcho)     = max(vcinp(k,Lhcho),epsilc)
 231         chem(i,k,j,p_op1)     = max(vcinp(k,Lop1),epsilc)
 232         chem(i,k,j,p_op2)     = max(vcinp(k,Lop2),epsilc)
 233         chem(i,k,j,p_paa)     = max(vcinp(k,Lpaa),epsilc)
 234         chem(i,k,j,p_ora1)     = max(vcinp(k,Lora1),epsilc)
 235         chem(i,k,j,p_ora2)     = max(vcinp(k,Lora2),epsilc)
 236         chem(i,k,j,p_nh3)     = max(vcinp(k,Lnh3),epsilc)
 237         chem(i,k,j,p_n2o5)     = max(vcinp(k,Ln2o5),epsilc)
 238         chem(i,k,j,p_no3)     = max(vcinp(k,Lno3),epsilc)
 239         chem(i,k,j,p_pan)     = max(vcinp(k,Lpan),epsilc)
 240         chem(i,k,j,p_hc3)     = max(vcinp(k,Lhc3),epsilc)
 241         chem(i,k,j,p_hc5)     = max(vcinp(k,Lhc5),epsilc)
 242         chem(i,k,j,p_hc8)     = max(vcinp(k,Lhc8),epsilc)
 243         chem(i,k,j,p_eth)     = max(vcinp(k,Leth),epsilc)
 244         chem(i,k,j,p_co)     = max(vcinp(k,Lco),epsilc)
 245         chem(i,k,j,p_ol2)     = max(vcinp(k,Lol2),epsilc)
 246         chem(i,k,j,p_olt)     = max(vcinp(k,Lolt),epsilc)
 247         chem(i,k,j,p_oli)     = max(vcinp(k,Loli),epsilc)
 248         chem(i,k,j,p_tol)     = max(vcinp(k,Ltol),epsilc)
 249         chem(i,k,j,p_xyl)     = max(vcinp(k,Lxyl),epsilc)
 250         chem(i,k,j,p_aco3)     = max(vcinp(k,Laco3),epsilc)
 251         chem(i,k,j,p_tpan)     = max(vcinp(k,Ltpan),epsilc)
 252         chem(i,k,j,p_hono)     = max(vcinp(k,Lhono),epsilc)
 253         chem(i,k,j,p_hno4)     = max(vcinp(k,Lhno4),epsilc)
 254         chem(i,k,j,p_ket)     = max(vcinp(k,Lket),epsilc)
 255         chem(i,k,j,p_gly)     = max(vcinp(k,Lgly),epsilc)
 256         chem(i,k,j,p_mgly)     = max(vcinp(k,Lmgly),epsilc)
 257         chem(i,k,j,p_dcb)     = max(vcinp(k,Ldcb),epsilc)
 258         chem(i,k,j,p_onit)     = max(vcinp(k,Lonit),epsilc)
 259         chem(i,k,j,p_csl)     = max(vcinp(k,Lcsl),epsilc)
 260         chem(i,k,j,p_iso)     = max(vcinp(k,Liso),epsilc)
 261         chem(i,k,j,p_ho)     = max(vcinp(k,Lho),epsilc)
 262         chem(i,k,j,p_ho2)     = max(vcinp(k,Lho2),epsilc)
 263         if(p_nu0.gt.1)then
 264         VDROG3(i,k,j,PXYL ) = VDROG1(k,PXYL )
 265         VDROG3(i,k,j,PTOL ) = VDROG1(k,PTOL )
 266         VDROG3(i,k,j,PCSL1) = VDROG1(k,PCSL1)
 267         VDROG3(i,k,j,PCSL2) = VDROG1(k,PCSL2)
 268         VDROG3(i,k,j,PHC8 ) = VDROG1(k,PHC8 )
 269         VDROG3(i,k,j,POLI1) = VDROG1(k,POLI1)
 270         VDROG3(i,k,j,POLI2) = VDROG1(k,POLI2)
 271         VDROG3(i,k,j,POLI3) = VDROG1(k,POLI3)
 272         VDROG3(i,k,j,POLT1) = VDROG1(k,POLT1)
 273         VDROG3(i,k,j,POLT2) = VDROG1(k,POLT2)
 274         VDROG3(i,k,j,POLT3) = VDROG1(k,POLT3)
 275         endif
 276       END DO
 277 !       if(iprt.eq.2)then
 278 !       print *,'after radm, k,chem(i,k,j,p_sulf),vcinp(k,lsulf)'
 279 !       do k=kts,kte
 280 !        print *,k,chem(i,k,j,p_sulf),vcinp(k,lsulf)
 281 !       enddo
 282 !       endif
 283 100   continue
 284
 285
 286 END SUBROUTINE radm_driver
 287
 288
 289       SUBROUTINE radm(rjj,wlcc,vcinp,tinp,pinp,rhinp,tstart,timemx,   &
 290                  jcs,jce,iprt,dt60,dtcmax,dtcmin,vdrog,iaerosol_sorgam)
 291          implicit none
 292 ! .. Parameters ..
 293         REAL, PARAMETER :: c302 = 5417.4, c303 = 19.83
 294 ! ..
 295 ! .. Scalar Arguments ..
 296         REAL(KIND=8),INTENT(IN) :: timemx, tstart
 297         REAL,INTENT(IN) :: dt60, dtcmax, dtcmin
 298         INTEGER, INTENT(IN) :: iprt, jce, jcs
 299 !
 300 !
 301
 302
 303          integer, intent (in) :: iaerosol_sorgam
 304          REAL,INTENT(IN) :: rjj(jcs:jce,nreacj),                       &
 305               wlcc(jcs:jce), tinp(jcs:jce),pinp(jcs:jce),rhinp(jcs:jce)
 306 ! ..
 307         real,intent (INOUT) :: vdrog(jcs:jce,ldrog),vcinp(jcs:jce,lspec)
 308 ! ..
 309 ! .. Local Scalars ..
 310         REAL :: dtc, r, timenow, tsqrd, xk0, xk2, xk3
 311         INTEGER :: i, ir, irdum, j, k, kdum, l, nr
 312 ! ..
 313 ! .. Local Arrays ..
 314         REAL :: prdrog(jcs:jce,ldrog)
 315         REAL :: aquad(jcs:jce), bquad(jcs:jce),                         &
 316           crj(jcs:jce,nreacj), crk(jcs:jce,nreack),                     &
 317           dum(jcs:jce), dvc(jcs:jce,ldiag), dvca(jcs:jce,ldiag),        &
 318           dvcg(jcs:jce,ldiag), h2o(jcs:jce,1),                          &
 319           loss(jcs:jce,lpred), lossl(jcs:jce,lump),                     &
 320           p(jcs:jce,1), patmot1(jcs:jce),                               &
 321           patmot2(jcs:jce), patmot3(jcs:jce),                           &
 322           pot(jcs:jce), prod(jcs:jce,lpred),                            &
 323           prodl(jcs:jce,lump), rh(jcs:jce),                             &
 324           rj(jcs:jce,nreacj), rk(jcs:jce,nreack),                       &
 325           t(jcs:jce,1), tin(jcs:jce), to300(jcs:jce),                   &
 326           vc(jcs:jce,1,lspec), vca(jcs:jce,1,lspec),                    &
 327           vcg(jcs:jce,1,lspec), vcl(jcs:jce,lump), wlc(jcs:jce)
 328 ! ..
 329 ! .. Intrinsic Functions ..
 330         INTRINSIC amax1, amin1, exp, log10
 331 ! ..
 332         IF (iprt==1) PRINT *, 'in radm ', jcs, jce, vcinp(jcs:jce,3), &
 333           vcinp(jcs:jce,lho2)
 334         IF (iprt==1) PRINT *, 'in radm ',  lspec, lho2
 335         IF (iprt==2) PRINT *, 'in radm ',  lsulf,vcinp(jcs:jcs+5,lsulf)
 336       r = 0.0820578
 337       do nr=1,ldrog
 338          do j=jcs,jce
 339           VDROG(j,nr)=0.
 340          enddo
 341       enddo
 342       DO nr = 1, nreacj
 343           DO j = jcs, jce
 344             rj(j,nr) = rjj(j,nr)
 345           END DO
 346       END DO
 347       DO j = jcs, jce
 348           wlc(j) = wlcc(j)
 349           t(j,1) = tinp(j)
 350           p(j,1) = pinp(j)
 351           rh(j) = rhinp(j)
 352       END DO
 353       DO l = 1, lspec
 354           DO j = jcs, jce
 355             vca(j,1,l) = epsilc
 356             vcg(j,1,l) = amax1(epsilc,vcinp(j,l))
 357             vc(j,1,l) = amax1(epsilc,vcinp(j,l))
 358           END DO
 359       END DO
 360         IF (iprt==1) PRINT *, ' radm', lho2, vc(jcs:jce,1,3), vc(jcs:jce,1,7), &
 361           vc(jcs:jce,1,lho2)
 362         DO l = 1, lpred
 363           DO j = jcs, jce
 364             prod(j,l) = 0.
 365             loss(j,l) = epsilc
 366           END DO
 367         END DO
 368         DO l = 1, nreacj
 369           DO j = jcs, jce
 370             crj(j,l) = 0.
 371           END DO
 372         END DO
 373         DO l = 1, ldiag
 374           DO j = jcs, jce
 375             dvca(j,l) = epsilc
 376             dvcg(j,l) = epsilc
 377             dvc(j,l) = epsilc
 378           END DO
 379         END DO
 380         DO l = 1, nreack
 381           DO j = jcs, jce
 382             rk(j,l) = 0.
 383             crk(j,l) = epsilc
 384           END DO
 385         END DO
 386         DO l = 1, lump
 387           DO j = jcs, jce
 388             vcl(j,l) = 1.e-9
 389             lossl(j,l) = epsilc
 390             prodl(j,l) = 0.
 391           END DO
 392         END DO
 393
 394         dtc = dtcmin
 395
 396         DO j = jcs, jce
 397           h2o(j,1) = .611E6*rh(j)*exp(c303-c302/t(j,1))/p(j,1)
 398         END DO
 399
 400         k = 1
 401         i = 1
 402         kdum = k
 403         DO j = jcs, jce
 404           tin(j) = 1./t(j,1) !RADM2.0 I --> IMRCHEM
 405 !RADM2.0 I --> IMRCHEM
 406           pot(j) = p(j,1)*tin(j)/101.3
 407 !RADM2.0 I --> IMRCHEM
 408           to300(j) = t(j,1)/300.
 409           patmot1(j) = const(1)*pot(j)
 410           patmot2(j) = const(2)*pot(j)
 411           patmot3(j) = const(3)*pot(j)*pot(j)
 412         END DO
 413         DO ir = 1, nreack
 414           DO j = jcs, jce
 415             rk(j,ir) = thafac(ir)*exp(-eor(ir)*tin(j))*patmot2(j)
 416           END DO
 417         END DO
 418         DO j = jcs, jce
 419 !3RD ORDER
 420           rk(j,16) = rk(j,16)*patmot3(j)/patmot2(j)*1.E-20
 421 !1ST ORDER
 422           rk(j,54) = rk(j,54)/patmot2(j)*60.
 423 !1ST ORDER
 424           rk(j,56) = rk(j,56)/patmot2(j)*60.
 425         END DO
 426         DO ir = 1, ntroe
 427           irdum = itroe(ir)
 428           DO j = jcs, jce
 429             aquad(j) = xk0300(ir)*to300(j)**(-xntroe(ir))
 430             aquad(j) = aquad(j)*patmot1(j)
 431             bquad(j) = xkf300(ir)*to300(j)**(-xmtroe(ir))
 432             bquad(j) = aquad(j)/bquad(j)
 433           END DO
 434           DO j = jcs, jce
 435             rk(j,irdum) = aquad(j)/(1.+bquad(j))*0.6**(1./(1.+(log10(bquad(j)) &
 436               )**2))
 437           END DO
 438           IF (ir>2) THEN
 439             DO j = jcs, jce
 440               rk(j,irdum) = rk(j,irdum)*patmot2(j)
 441             END DO
 442           ELSE
 443             DO j = jcs, jce
 444 !changed RADM2.0 IMRCHEM
 445               rk(j,irdum) = rk(j,irdum)/(afac(ir)*exp(bfac(ir)/t(j,1)))*60.
 446             END DO
 447           END IF
 448 !END DO 90 LOOP
 449         END DO
 450         DO j = jcs, jce
 451           tsqrd = t(j,1)*t(j,1) !was Imrchem 3d I --> IMRCHEM
 452           rk(j,30) = rk(j,30)*tsqrd
 453           rk(j,31) = rk(j,31)*tsqrd
 454           rk(j,50) = rk(j,50)*tsqrd
 455         END DO
 456         DO j = jcs, jce
 457           rk(j,1) = patmot1(j)*6.E-34*to300(j)**(-2.3)*patmot2(j)
 458           rk(j,12) = (2.2E-13*exp(620.*tin(j))+1.9E-33*patmot1(j)*exp(980.*tin &
 459             (j)))*patmot2(j)
 460 !         IF (iprt==1 .AND. j==jce) THEN
 461 !           PRINT *, j, tin(j), patmot1(j), patmot2(j), &
 462 !             1.9E-33*patmot1(j)*exp(980.*tin(j))
 463 !           PRINT *, rk(j,12), 2.2E-13*exp(620.*tin(j)), const(3), p(j,1)
 464 !         END IF
 465           xk0 = 7.2E-15*exp(785.*tin(j))
 466           xk2 = 4.1E-16*exp(1440.*tin(j))
 467           xk3 = 1.9E-33*exp(725.*tin(j))*patmot1(j)
 468           rk(j,25) = (xk0+xk3/(1.+xk3/xk2))*patmot2(j)
 469           rk(j,29) = (1.5E-13*(1.+2.439E-20*patmot1(j)))*patmot2(j)
 470           rk(j,13) = (3.08E-34*exp(2820.*tin(j))+2.66E-34*patmot1(j)*1.E-20* &
 471             exp(3180.*tin(j)))*patmot3(j)
 472         END DO
 473         DO j = jcs, jce
 474           dum(j) = amin1(rh(j),1.)
 475           dum(j) = amax1(dum(j),0.)
 476 !            RK(J,137) =  1./(600.*EXP(-(DUM(J)/.28)**2.8)+5.)
 477 !          RK(J,137)= CVMGP(0.2,0.0,DUM(J) - .70) ! HETEROGENOUS N2O5
 478           rk(j,23) = 0.0
 479 ! HOMOGENEOUS  N2O5
 480           IF (dum(j)-.7>=0.) THEN
 481             rk(j,137) = 0.2
 482           ELSE
 483             rk(j,137) = 0.
 484           END IF
 485         END DO
 486 !**
 487         DO j = jcs, jce
 488           vcl(j,lnox) = vc(j,1,lno) + vc(j,1,lno2)
 489           vcl(j,lhox) = max(1.e-9,vc(j,1,lho) + vc(j,1,lho2))
 490           vcl(j,lpao3) = vc(j,1,lpan) + vc(j,1,laco3)
 491           vcl(j,ln2n3) = vc(j,1,lno3) + vc(j,1,ln2o5)
 492         END DO
 493 !**********************************************************************
 494 !                  C H E M I C A L       S O L V E R
 495 !**********************************************************************
 496         timenow = 0.
 497 10      CONTINUE
 498
 499 ! Chemical solver
 500         CALL predraten(jcs,jce,iprt,crj,crk,rj,rk,vc,dvc,vca, &
 501                                        wlc,dvca,p,h2o,dvcg,t,r)
 502
 503         CALL producn(jcs,jce,iprt,crj,crk,loss,prod,prodl,lossl, &
 504                                            prdrog,iaerosol_sorgam)
 505
 506         CALL setdtc(jcs,jce,dtc,dtcmax,dtcmin,dt60,prod,loss,vc,timenow)
 507
 508         CALL integ1n(jcs,jce,iprt,dtc,vc,loss,prod,vcl,lossl,prodl,    &
 509                          rk,dvc,h2o,rj,vdrog,prdrog,iaerosol_sorgam)
 510
 511         timenow = timenow + dtc
 512         IF (iprt==2) PRINT *, 'end radm',  timenow,vc(jcs:jce,1,lsulf)
 513         IF ((timenow+0.001)<dt60) GO TO 10
 514 !********************END K - LEVEL LOOP
 515         IF (iprt==1) PRINT *, 'end radm',  vc(jcs:jce,1,3), vc(jcs:jce,1,7)
 516         DO l = 1, lspec
 517           DO j = jcs, jce
 518             vcinp(j,l) = amax1(epsilc,vc(j,1,l))
 519           END DO
 520         END DO
 521
 522
 523
 524         RETURN
 525 END SUBROUTINE radm
 526
 527       SUBROUTINE integ1n(jcs,jce,iprt,dtc,vc,loss,prod,vcl,lossl, &
 528               prodl,rk,dvc,h2o,rj,vdrog,prdrog,iaerosol_sorgam)
 529        implicit none
 530 !***********************************************************************
 531 !* DESCRIPTION:
 532 !* This is the chemical integration routine.  The routine employs the
 533 !* FOLLOWING APPROACH TO THE SOLUTION OF THE DIFFERENTIAL
 534 !* equations by the following algorithm:
 535 !*
 536 !*              dt * Pn + Yn (1 - dt/2 * (Ln/Yn))
 537 !*      Yn+1=   ---------------------------------
 538 !*                     (1 + dt/2 * (Ln/Yn))
 539 !*
 540 !***********************************************************************
 541 !**    Revision history:
 542 !**    NO.   DATE     WHO      WHAT
 543 !**    __    ____     ___      _______________________________________
 544 !**    03    89279    MCB      INTEGRATED VECTORIZATION CONTROL FROM
 545 !**                              JEFF BROOKS OF CRAY RESEARCH.  THIS
 546 !**                              INCLUDED PUSHING J LOOPS INSIDE THE
 547 !**                              IF BLOCKS AND INSIDE REST OF CODE FOR
 548 !**                              VECTORIZATION.  (THIS CLAIMS TO GIVE
 549 !**                              APPROXIMATELY A FACTOR OF 2.3 INCREASE
 550 !**                              IN EXECUTION TIME ALONE!)
 551 !**    02    89153    MCB      REMOVED CALL TO PREDRT2 - EMBEDDED CODE
 552 !**    01    89149    MCB       INTEGRATIED FINAL STAND ALONE VERSION
 553 !**
 554 !**    Input files: None
 555 !**
 556 !**    Output files:None
 557 !**
 558 !**    Called by: CHEM
 559 !**
 560 !**    CALLS THE FOLLOWING SUBROUTINES: NONE
 561 !**********************************************************************
 562 !*****THE FOLLOWING TEMPORARY ARRAYS WERE CREATED BY JEFF BROOKS FOR
 563 !*****VECTORIZABLE CODE.
 564 ! .. Scalar Arguments ..
 565         REAL, INTENT(IN) ::  dtc
 566         INTEGER, intent(in) :: iprt, jce, jcs
 567         integer, intent (in) :: iaerosol_sorgam
 568         real,intent(in)::prdrog(jcs:jce,ldrog)
 569         real,intent(inout)::vdrog(jcs:jce,ldrog)
 570
 571 ! ..
 572 ! .. Array Arguments ..
 573         REAL :: dvc(jcs:jce,ldiag), h2o(jcs:jce,1),                      &
 574           loss(jcs:jce,lpred), lossl(jcs:jce,lump), prod(jcs:jce,lpred), &
 575           prodl(jcs:jce,lump), rj(jcs:jce,nreacj), rk(jcs:jce,nreack),   &
 576           vc(jcs:jce,1,lspec), vcl(jcs:jce,lump)
 577 ! ..
 578 ! .. Local Scalars ..
 579         REAL :: eqno,alow
 580         INTEGER :: inumho, inumhox, iter, j, l, niter
 581 ! ..
 582 ! .. Local Arrays ..
 583         REAL :: crji(nreacj), crki(nreack), h2o2sv(jcs:jce),             &
 584                 hoxsv(jcs:jce),jb(jcs:jce), jbb(jcs:jce), jc(jcs:jce),   &
 585                 jloss(jcs:jce), jprod(jcs:jce), tauinv(jcs:jce),         &
 586                 taulinv(jcs:jce)
 587 ! ..
 588 ! .. Intrinsic Functions ..
 589         INTRINSIC exp, max
 590 ! ..
 591         alow=1.e-17
 592         DO j = jcs, jce
 593           jprod(j) = 0.
 594           jb(j) = 0.
 595           jc(j) = 0.
 596           jloss(j) = 0.
 597           jbb(j) = 0.
 598         END DO
 599 !     K      = 1
 600 !     if(iprt.eq.1)print *,'in integ1'
 601 !     if(iprt.eq.2)print *,'in integ1',lsulf,intgrt(lsulf),qdtc(lsulf)
 602         inumho = 3
 603         inumhox = 3
 604         DO l = 1, lpred
 605           IF (intgrt(l)==1) THEN
 606
 607             DO j = jcs, jce
 608 !              if(iprt.eq.1.and.(l.eq.3.or.l.eq.7))then
 609 !               print *,'i1',l,j,vc(j,1,l),loss(j,l),cmin(l)
 610 !              endif
 611               tauinv(j) = (loss(j,l)/vc(j,1,l))*0.5
 612               vc(j,1,l) = (vc(j,1,l)*(1.-dtc*tauinv(j))+prod(j,l)*dtc)/ &
 613                 (1+dtc*tauinv(j))
 614               vc(j,1,l) = max(cmin(l),vc(j,1,l))
 615 !              if(iprt.eq.1.and.(l.eq.3.or.l.eq.7))then
 616 !               print *,'i1',l,j,vc(j,1,l),dtc
 617 !              endif
 618 !              if(iprt.eq.2.and.(l.eq.lsulf))then
 619 !               print *,'i1',l,j,vc(j,1,l),dtc,loss(j,l),prod(j,l)
 620 !              endif
 621             END DO
 622           END IF
 623         END DO
 624
 625 !... do H2O2 exponentially and store old conc in H2O2SV for iteration
 626
 627         DO j = jcs, jce
 628           h2o2sv(j) = vc(j,1,lh2o2)
 629           tauinv(j) = loss(j,lh2o2)/vc(j,1,lh2o2)
 630           vc(j,1,lh2o2) = prod(j,lh2o2)/tauinv(j) + (h2o2sv(j)-prod(j,lh2o2)/ &
 631             tauinv(j))*exp(-dtc*tauinv(j))
 632           vc(j,1,lh2o2) = max(cmin(lh2o2),vc(j,1,lh2o2))
 633         END DO
 634
 635
 636 !*****LUMPED SPECIES INTEGRATION:
 637         DO l = 1, lump
 638           IF (l==lhox) THEN
 639             DO j = jcs, jce
 640               hoxsv(j) = vcl(j,lhox)
 641               taulinv(j) = lossl(j,l)/vcl(j,l)
 642
 643               vcl(j,l) = prodl(j,l)/taulinv(j) + (hoxsv(j)-prodl(j,l)/taulinv( &
 644                 j))*exp(-dtc*taulinv(j))
 645 !             vcl(j,l)=max(cmin(l),vcl(j,l))
 646               vcl(j,l)=max(1.e-9,vcl(j,l))
 647             END DO
 648           ELSE
 649             DO j = jcs, jce
 650               taulinv(j) = (lossl(j,l)/vcl(j,l))*0.5
 651               vcl(j,l) = (vcl(j,l)*(1.-dtc*taulinv(j))+prodl(j,l)*dtc)/ &
 652                 (1.+dtc*taulinv(j))
 653 !             vcl(j,l)=max(cmin(l),vcl(j,l))
 654               vcl(j,l)=max(1.e-9,vcl(j,l))
 655             END DO
 656           END IF
 657         END DO
 658
 659         DO j = jcs, jce
 660
 661 !*******[N2O5]:
 662           vc(j,1,ln2o5) = vcl(j,ln2n3) - vc(j,1,lno3)
 663           vc(j,1,ln2o5) = max(cmin(ln2o5),vc(j,1,ln2o5))
 664
 665 !*******[NO]:
 666 !       VC(J,1,LNO)  = VCL(J,LNOX) - VC(J,1,LNO2)
 667 !       VC(J,1,LNO)  = MAX(CMIN(LNO),VC(J,1,LNO))
 668 ! 8/31/00 Calculate steady state NO from NO2,O3 ...etc., renormalize NO,
 669 !23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
 670           jprod(j) = rj(j,1)*vc(j,1,lno2) + rj(j,4)*vc(j,1,lhono) + &
 671             rj(j,8)*vc(j,1,lno3)
 672           jloss(j) = rk(j,6)*vc(j,1,lo3) + rk(j,9)*vc(j,1,lho2) + &
 673             rk(j,15)*vc(j,1,lho) + rk(j,16)*vc(j,1,lno) + &
 674             rk(j,18)*vc(j,1,lno3) + rk(j,57)*dvc(j,lmo2) + &
 675             rk(j,58)*dvc(j,lhc3p) + rk(j,60)*dvc(j,lhc5p) + &
 676             rk(j,62)*dvc(j,lhc8p) + rk(j,64)*dvc(j,lol2p) + &
 677             rk(j,65)*dvc(j,loltp) + rk(j,66)*dvc(j,lolip) + &
 678             rk(j,67)*dvc(j,ltco3) + rk(j,68)*dvc(j,ltco3) + &
 679             rk(j,69)*dvc(j,ltolp) + rk(j,70)*dvc(j,lxylp) + &
 680             rk(j,71)*dvc(j,lethp) + rk(j,72)*dvc(j,lketp) + &
 681             rk(j,73)*dvc(j,loln) + rk(j,131)*dvc(j,lxo2)
 682           eqno = max(cmin(lno),jprod(j)/max(epsilc,jloss(j)))
 683           EQNO=MAX(1.e-6,JPROD(J)/MAX(EPSILC,JLOSS(J)))
 684 !*******[NO]:
 685           vc(j,1,lno) = max(1.e-6,eqno*vcl(j,lnox)/(eqno+vc(j,1,lno2)))
 686           vc(j,1,lno2) = vc(j,1,lno2)*vcl(j,lnox)/(eqno+vc(j,1,lno2))
 687
 688 !*******[PAN]:
 689           vc(j,1,lpan) = vcl(j,lpao3) - vc(j,1,laco3)
 690           vc(j,1,lpan) = max(cmin(lpan),vc(j,1,lpan))
 691 !         if(iprt.ge. 1.and.j.eq.jcs)then
 692 !            print *,'i1,no',j,vc(j,1,lno),vc(j,1,lno2),EQNO,JPROD(J),JLOSS(J),rj(j,1), &
 693 !                 vc(j,1,lno2),rj(j,8),vc(j,1,lno3),rj(j,4),vc(j,1,lhono)
 694 !         endif
 695         END DO
 696
 697         DO j = jcs, jce
 698           jprod(j) = rj(j,4)*vc(j,1,lhono) + rj(j,5)*vc(j,1,lhno3) + &
 699             2.*rj(j,9)*vc(j,1,lh2o2) + rj(j,13)*vc(j,1,lop1) + &
 700             rj(j,14)*vc(j,1,lop2) + rj(j,15)*vc(j,1,lpaa) + &
 701             2.E0*rk(j,5)*dvc(j,lo1d)*h2o(j,1) + .1E0*(rk(j,85)*vc(j,1,lolt)*vc &
 702             (j,1,lo3)+rk(j,87)*vc(j,1,liso)*vc(j,1,lo3)) + &
 703             .14E0*rk(j,86)*vc(j,1,loli)*vc(j,1,lo3)
 704
 705           jloss(j) = rk(j,7)*vc(j,1,lo3) + rk(j,14)*vc(j,1,lh2o2) + &
 706             rk(j,15)*vc(j,1,lno) + rk(j,24)*vc(j,1,lno2) + &
 707             rk(j,25)*vc(j,1,lhno3) + rk(j,26)*vc(j,1,lhno4) + &
 708             rk(j,28)*vc(j,1,lso2) + rk(j,29)*vc(j,1,lco) + rk(j,30)*ch4 + &
 709             rk(j,31)*vc(j,1,leth) + rk(j,32)*vc(j,1,lhc3) + &
 710             rk(j,33)*vc(j,1,lhc5) + rk(j,34)*vc(j,1,lhc8) + &
 711             rk(j,35)*vc(j,1,lol2) + rk(j,36)*vc(j,1,lolt) + &
 712             rk(j,37)*vc(j,1,loli) + rk(j,38)*vc(j,1,ltol) + &
 713             rk(j,39)*vc(j,1,lxyl) + rk(j,40)*vc(j,1,lcsl) + &
 714             rk(j,41)*vc(j,1,lhcho) + rk(j,42)*vc(j,1,lald) + &
 715             rk(j,43)*vc(j,1,lket) + rk(j,44)*vc(j,1,lgly) + &
 716             rk(j,45)*vc(j,1,lmgly) + rk(j,46)*vc(j,1,ldcb) + &
 717             .5E0*rk(j,47)*vc(j,1,lop1) + .5E0*rk(j,48)*vc(j,1,lop2) + &
 718             rk(j,49)*vc(j,1,lpaa) + rk(j,50)*vc(j,1,lpan) + &
 719             rk(j,51)*vc(j,1,lonit) + rk(j,52)*vc(j,1,liso)
 720           jbb(j) = rk(j,8)*vc(j,1,lo3) + rk(j,9)*vc(j,1,lno) + jloss(j)
 721         END DO
 722 !       if(iprt.eq.1)print *,'jloss '
 723
 724         DO niter = 1, inumhox
 725           DO j = jcs, jce
 726             crji(9) = rj(j,9)*vc(j,1,lh2o2)
 727             crki(10) = rk(j,10)*vc(j,1,lho2)*vc(j,1,lno2)
 728             crki(12) = rk(j,12)*vc(j,1,lho2)*vc(j,1,lho2)
 729             crki(13) = rk(j,13)*vc(j,1,lho2)*vc(j,1,lho2)*h2o(j,1)
 730             crki(14) = rk(j,14)*vc(j,1,lh2o2)*vc(j,1,lho)
 731             crki(15) = rk(j,15)*vc(j,1,lno)*vc(j,1,lho)
 732             crki(20) = rk(j,20)*vc(j,1,lno3)*vc(j,1,lho2)
 733             crki(24) = rk(j,24)*vc(j,1,lho)*vc(j,1,lno2)
 734             crki(25) = rk(j,25)*vc(j,1,lho)*vc(j,1,lhno3)
 735             crki(26) = rk(j,26)*vc(j,1,lho)*vc(j,1,lhno4)
 736             crki(27) = rk(j,27)*vc(j,1,lho)*vc(j,1,lho2)
 737             crki(30) = rk(j,30)*ch4*vc(j,1,lho)
 738             crki(31) = rk(j,31)*vc(j,1,leth)*vc(j,1,lho)
 739             crki(32) = rk(j,32)*vc(j,1,lhc3)*vc(j,1,lho)
 740             crki(33) = rk(j,33)*vc(j,1,lhc5)*vc(j,1,lho)
 741             crki(34) = rk(j,34)*vc(j,1,lhc8)*vc(j,1,lho)
 742             crki(35) = rk(j,35)*vc(j,1,lol2)*vc(j,1,lho)
 743             crki(36) = rk(j,36)*vc(j,1,lolt)*vc(j,1,lho)
 744             crki(37) = rk(j,37)*vc(j,1,loli)*vc(j,1,lho)
 745             crki(38) = rk(j,38)*vc(j,1,ltol)*vc(j,1,lho)
 746             crki(39) = rk(j,39)*vc(j,1,lxyl)*vc(j,1,lho)
 747             crki(40) = rk(j,40)*vc(j,1,lcsl)*vc(j,1,lho)
 748             crki(42) = rk(j,42)*vc(j,1,lald)*vc(j,1,lho)
 749             crki(43) = rk(j,43)*vc(j,1,lket)*vc(j,1,lho)
 750             crki(45) = rk(j,45)*vc(j,1,lmgly)*vc(j,1,lho)
 751             crki(46) = rk(j,46)*vc(j,1,ldcb)*vc(j,1,lho)
 752             crki(47) = rk(j,47)*vc(j,1,lop1)*vc(j,1,lho)
 753             crki(48) = rk(j,48)*vc(j,1,lop2)*vc(j,1,lho)
 754             crki(49) = rk(j,49)*vc(j,1,lpaa)*vc(j,1,lho)
 755             crki(50) = rk(j,50)*vc(j,1,lpan)*vc(j,1,lho)
 756             crki(51) = rk(j,51)*vc(j,1,lonit)*vc(j,1,lho)
 757             crki(52) = rk(j,52)*vc(j,1,liso)*vc(j,1,lho)
 758             crki(88) = rk(j,88)*vc(j,1,lho2)*dvc(j,lmo2)
 759             crki(89) = rk(j,89)*vc(j,1,lho2)*dvc(j,lethp)
 760             crki(90) = rk(j,90)*vc(j,1,lho2)*dvc(j,lhc3p)
 761             crki(91) = rk(j,91)*vc(j,1,lho2)*dvc(j,lhc5p)
 762             crki(92) = rk(j,92)*vc(j,1,lho2)*dvc(j,lhc8p)
 763             crki(93) = rk(j,93)*vc(j,1,lho2)*dvc(j,lol2p)
 764             crki(94) = rk(j,94)*vc(j,1,lho2)*dvc(j,loltp)
 765             crki(95) = rk(j,95)*vc(j,1,lho2)*dvc(j,lolip)
 766             crki(96) = rk(j,96)*vc(j,1,lho2)*dvc(j,lketp)
 767             crki(97) = rk(j,97)*vc(j,1,lho2)*vc(j,1,laco3)
 768             crki(98) = rk(j,98)*vc(j,1,lho2)*dvc(j,ltolp)
 769             crki(99) = rk(j,99)*vc(j,1,lho2)*dvc(j,lxylp)
 770             crki(100) = rk(j,100)*vc(j,1,lho2)*dvc(j,ltco3)
 771             crki(101) = rk(j,101)*vc(j,1,lho2)*dvc(j,loln)
 772             crki(127) = rk(j,127)*dvc(j,lxo2)*vc(j,1,lho2)
 773             crki(133) = rk(j,133)*dvc(j,lxno2)*vc(j,1,lho2)
 774
 775             lossl(j,lhox) = crki(10) + 2.*crki(12) + 2.*crki(13) + crki(15) + &
 776               crki(20) + crki(24) + crki(25) + crki(26) + 2.*crki(27) + &
 777               crki(30) + crki(31) + .83*crki(32) + crki(33) + crki(34) + &
 778               crki(35) + crki(36) + crki(37) + .75*crki(38) + .83*crki(39) + &
 779               1.8*crki(40) + crki(42) + crki(43) + crki(45) + crki(46) + &
 780               .5*crki(47) + .5*crki(48) + crki(49) + crki(50) + crki(51) + &
 781               crki(52) + crki(88) + crki(89) + crki(90) + crki(91) + &
 782               crki(92) + crki(93) + crki(94) + crki(95) + crki(96) + &
 783               crki(97) + crki(98) + crki(99) + crki(100) + crki(101) + &
 784               crki(127) + crki(133)
 785             lossl(j,lhox) = max(alow,lossl(j,lhox))
 786             prod(j,lh2o2) = crki(12) + crki(13)
 787             loss(j,lh2o2) = crji(9) + crki(14)
 788           END DO
 789
 790 !       if(iprt.eq.1)print *,'i1, hox,niter ',niter
 791           DO j = jcs, jce
 792 !********[HOX]:
 793             taulinv(j) = (lossl(j,lhox)/hoxsv(j))
 794             vcl(j,lhox) = prodl(j,lhox)/taulinv(j) + &
 795               (hoxsv(j)-prodl(j,lhox)/taulinv(j))*exp(-dtc*taulinv(j))
 796 !           vcl(j,lhox) = max(cmin(lhox),vcl(j,lhox))
 797               vcl(j,lhox)=max(1.e-9,vcl(j,lhox))
 798 !           if(iprt.eq.2.and.j.eq.jcs)print *,lossl(j,lhox),hoxsv(j),prodl(j,lhox),cmin(lhox),vcl(j,lhox)
 799             jc(j) = jprod(j) + (rk(j,9)*vc(j,1,lno)+rk(j,8)*vc(j,1,lo3))*vcl(j &
 800               ,lhox)
 801             jb(j) = jbb(j) + (rk(j,27)*vcl(j,lhox))
 802
 803 !********[H2O2]:
 804             tauinv(j) = loss(j,lh2o2)/h2o2sv(j)
 805             vc(j,1,lh2o2) = prod(j,lh2o2)/tauinv(j) + &
 806               (h2o2sv(j)-prod(j,lh2o2)/tauinv(j))*exp(-dtc*tauinv(j))
 807             vc(j,1,lh2o2) = max(cmin(lh2o2),vc(j,1,lh2o2))
 808
 809           END DO
 810
 811 !********[HO]:
 812 !       if(iprt.eq.1)print *,'i1, ho and ho2 '
 813           DO iter = 1, inumho
 814 !         if(iprt.eq.2)print *,'i1, ho and ho2 iter = ',iter
 815             DO j = jcs, jce
 816               vc(j,1,lho) = vc(j,1,lho) - (rk(j,27)*(vc(j,1, &
 817                 lho)**2)-jb(j)*vc(j,1,lho)+jc(j))/(2*rk(j,27)*vc(j,1,lho)-jb(j &
 818                 ))
 819             vc(j,1,lho) = max(cmin(lho),vc(j,1,lho))
 820 !           if(iprt.eq.2.and.j.eq.jcs)print *,vc(j,1,lho),cmin(lho),&
 821 !                                         rk(j,27),jb(j),jc(j)
 822             END DO
 823           END DO
 824
 825 !********[HO2]:
 826           DO j = jcs, jce
 827             vc(j,1,lho2) = vcl(j,lhox) - vc(j,1,lho)
 828             vc(j,1,lho2) = max(cmin(lho2),vc(j,1,lho2))
 829 !           if(iprt.eq.2.and.j.eq.jcs)print *,vc(j,1,lho2),vcl(j,lhox),vc(j,1,lho),cmin(lho2)
 830           END DO
 831
 832         END DO
 833 !
 834 !       if(iprt.eq.1)print *,'i1, aerosols '
 835         if(iaerosol_sorgam==1)then
 836 !       if(iprt.eq.2)print *,'i1, aerosols '
 837         DO L = 1, LDROG
 838            DO J = JCS, JCE
 839               VDROG( J, L ) = VDROG( J , L) + PRDROG( J, L ) * DTC
 840               VDROG( J, L ) = MAX( 0., VDROG( J, L ) )
 841            ENDDO
 842         ENDDO
 843         endif
 844 !
 845         RETURN
 846       END SUBROUTINE integ1n
 847       SUBROUTINE predraten(jcs,jce,iprt,crj,crk,rj,rk,vc,dvc,vca, &
 848                                            wlc,dvca,p,h2o,dvcg,t,r)
 849        implicit none
 850 ! .. Scalar Arguments ..
 851         REAL, INTENT(IN) :: r
 852         INTEGER, INTENT(IN) :: iprt, jce, jcs
 853
 854 ! ..
 855 ! .. Array Arguments ..
 856         REAL :: crj(jcs:jce,nreacj),  crk(jcs:jce,nreack),         &
 857           dvc(jcs:jce,ldiag), dvca(jcs:jce,ldiag),                 &
 858           dvcg(jcs:jce,ldiag), h2o(jcs:jce,1), p(jcs:jce,1),       &
 859           rj(jcs:jce,nreacj), rk(jcs:jce,nreack), t(jcs:jce,1),    &
 860           vc(jcs:jce,1,lspec), vca(jcs:jce,1,lspec), wlc(jcs:jce)
 861 ! ..
 862 ! .. Local Scalars ..
 863         INTEGER :: i, j, k, l
 864 ! ..
 865         k = 1
 866         i = 1
 867         DO j = jcs, jce
 868
 869           crj(j,1) = rj(j,1)*vc(j,1,lno2)
 870           crj(j,2) = rj(j,2)*vc(j,1,lo3)
 871           crj(j,3) = rj(j,3)*vc(j,1,lo3)
 872           crk(j,15) = rk(j,15)*vc(j,1,lno)*vc(j,1,lho)
 873           crj(j,4) = rj(j,4)*vc(j,1,lhono)
 874           crj(j,5) = rj(j,5)*vc(j,1,lhno3)
 875           crk(j,10) = rk(j,10)*vc(j,1,lho2)*vc(j,1,lno2)
 876           crj(j,6) = rj(j,6)*vc(j,1,lhno4)
 877           crj(j,7) = rj(j,7)*vc(j,1,lno3)
 878           crj(j,8) = rj(j,8)*vc(j,1,lno3)
 879           crj(j,9) = rj(j,9)*vc(j,1,lh2o2)
 880           crj(j,10) = rj(j,10)*vc(j,1,lhcho)
 881           crj(j,11) = rj(j,11)*vc(j,1,lhcho)
 882           crj(j,12) = rj(j,12)*vc(j,1,lald)
 883           crj(j,13) = rj(j,13)*vc(j,1,lop1)
 884           crj(j,14) = rj(j,14)*vc(j,1,lop2)
 885           crj(j,15) = rj(j,15)*vc(j,1,lpaa)
 886           crj(j,16) = rj(j,16)*vc(j,1,lket)
 887           crj(j,17) = rj(j,17)*vc(j,1,lgly)
 888           crj(j,18) = rj(j,18)*vc(j,k,lgly)
 889           crj(j,19) = rj(j,19)*vc(j,k,lmgly)
 890           crj(j,20) = rj(j,20)*vc(j,k,ldcb)
 891           crj(j,21) = rj(j,21)*vc(j,k,lonit)
 892         END DO
 893 !*
 894 !* end of photolysis equilimaxium species
 895         DO j = jcs, jce
 896           dvc(j,lo1d) = crj(j,2)/(rk(j,3)*n2+rk(j,4)*o2+rk(j,5)*h2o(j,1))
 897         END DO
 898         DO j = jcs, jce
 899           dvc(j,lo3p) = (crj(j,1)+crj(j,3)+crj(j,8)+rk(j,3)*dvc(j,lo1d)*n2+rk( &
 900             j,4)*dvc(j,lo1d)*o2)/(rk(j,2)*vc(j,1,lno2)+rk(j,1)*o2)
 901         END DO
 902         DO j = jcs, jce
 903           crk(j,1) = rk(j,1)*dvc(j,lo3p)*o2
 904           crk(j,2) = rk(j,2)*dvc(j,lo3p)*vc(j,1,lno2)
 905           crk(j,3) = rk(j,3)*dvc(j,lo1d)*n2
 906           crk(j,4) = rk(j,4)*dvc(j,lo1d)*o2
 907           crk(j,5) = rk(j,5)*dvc(j,lo1d)*h2o(j,1)
 908           crk(j,6) = rk(j,6)*vc(j,k,lo3)*vc(j,1,lno)
 909           crk(j,7) = rk(j,7)*vc(j,k,lo3)*vc(j,1,lho)
 910           crk(j,8) = rk(j,8)*vc(j,k,lo3)*vc(j,1,lho2)
 911           crk(j,9) = rk(j,9)*vc(j,k,lho2)*vc(j,1,lno)
 912           crk(j,11) = rk(j,11)*vc(j,1,lhno4)
 913           crk(j,12) = rk(j,12)*vc(j,1,lho2)*vc(j,1,lho2)
 914           crk(j,13) = rk(j,13)*vc(j,1,lho2)*vc(j,1,lho2)*h2o(j,1)
 915 !     if(iprt.eq.1.and.j.eq.jce)then
 916 !      print *,'in predate ',CRK(J, 12),CRK(J, 13),RK(J, 12),RK(J,13)
 917 !      print *,VC(J,1, LHO2),H2O(J,1)
 918 !     endif
 919           crk(j,14) = rk(j,14)*vc(j,1,lh2o2)*vc(j,1,lho)
 920           crk(j,16) = rk(j,16)*vc(j,1,lno)*vc(j,1,lno)*o2
 921           crk(j,17) = rk(j,17)*vc(j,1,lo3)*vc(j,1,lno2)
 922           crk(j,22) = rk(j,22)*vc(j,1,ln2o5)
 923           crk(j,23) = rk(j,23)*vc(j,1,ln2o5)*h2o(j,1)
 924           crk(j,24) = rk(j,24)*vc(j,1,lho)*vc(j,1,lno2)
 925           crk(j,25) = rk(j,25)*vc(j,1,lho)*vc(j,1,lhno3)
 926           crk(j,26) = rk(j,26)*vc(j,1,lho)*vc(j,1,lhno4)
 927           crk(j,27) = rk(j,27)*vc(j,1,lho)*vc(j,1,lho2)
 928           crk(j,28) = rk(j,28)*vc(j,1,lho)*vc(j,1,lso2)
 929 !     if(iprt.eq.2.and.j.eq.jcs)then
 930 !      print *,'in predate ',crk(j,28),rk(j,28),vc(j,1,lho),vc(j,1,lso2)
 931 !     endif
 932           crk(j,29) = rk(j,29)*vc(j,1,lco)*vc(j,1,lho)
 933           crk(j,30) = rk(j,30)*ch4*vc(j,1,lho)
 934           crk(j,31) = rk(j,31)*vc(j,1,leth)*vc(j,1,lho)
 935           crk(j,32) = rk(j,32)*vc(j,1,lhc3)*vc(j,1,lho)
 936           crk(j,33) = rk(j,33)*vc(j,1,lhc5)*vc(j,1,lho)
 937           crk(j,34) = rk(j,34)*vc(j,1,lhc8)*vc(j,1,lho)
 938           crk(j,35) = rk(j,35)*vc(j,1,lol2)*vc(j,1,lho)
 939           crk(j,36) = rk(j,36)*vc(j,1,lolt)*vc(j,1,lho)
 940           crk(j,37) = rk(j,37)*vc(j,1,loli)*vc(j,1,lho)
 941           crk(j,38) = rk(j,38)*vc(j,1,ltol)*vc(j,1,lho)
 942           crk(j,39) = rk(j,39)*vc(j,1,lxyl)*vc(j,1,lho)
 943           crk(j,40) = rk(j,40)*vc(j,1,lcsl)*vc(j,1,lho)
 944           crk(j,41) = rk(j,41)*vc(j,1,lhcho)*vc(j,1,lho)
 945           crk(j,42) = rk(j,42)*vc(j,1,lald)*vc(j,1,lho)
 946           crk(j,43) = rk(j,43)*vc(j,1,lket)*vc(j,1,lho)
 947           crk(j,44) = rk(j,44)*vc(j,1,lgly)*vc(j,1,lho)
 948           crk(j,45) = rk(j,45)*vc(j,1,lmgly)*vc(j,1,lho)
 949           crk(j,46) = rk(j,46)*vc(j,1,ldcb)*vc(j,1,lho)
 950           crk(j,47) = rk(j,47)*vc(j,1,lop1)*vc(j,1,lho)
 951           crk(j,48) = rk(j,48)*vc(j,1,lop2)*vc(j,1,lho)
 952           crk(j,49) = rk(j,49)*vc(j,1,lpaa)*vc(j,1,lho)
 953           crk(j,50) = rk(j,50)*vc(j,1,lpan)*vc(j,1,lho)
 954           crk(j,18) = rk(j,18)*vc(j,1,lno3)*vc(j,1,lno)
 955           crk(j,19) = rk(j,19)*vc(j,1,lno3)*vc(j,1,lno2)
 956           crk(j,20) = rk(j,20)*vc(j,1,lno3)*vc(j,1,lho2)
 957           crk(j,21) = rk(j,21)*vc(j,1,lno3)*vc(j,1,lno2)
 958           crk(j,51) = rk(j,51)*vc(j,1,lonit)*vc(j,1,lho)
 959           crk(j,52) = rk(j,52)*vc(j,1,liso)*vc(j,1,lho)
 960           crk(j,53) = rk(j,53)*vc(j,1,laco3)*vc(j,1,lno2)
 961           crk(j,54) = rk(j,54)*vc(j,1,lpan)
 962           crk(j,56) = rk(j,56)*vc(j,1,ltpan)
 963           crk(j,78) = rk(j,78)*vc(j,1,ldcb)*vc(j,1,lno3)
 964         END DO
 965         DO j = jcs, jce
 966           dvc(j,ltco3) = (crj(j,20)+crk(j,56)+0.9*crk(j,40)+crk(j,46)+crk(j,78 &
 967             ))/(rk(j,55)*vc(j,k,lno2)+rk(j,68)*vc(j,k,lno)+rk(j,100)*vc(j,k, &
 968             lho2)+rk(j,126)*vc(j,k,laco3)+rk(j,114)*dvc(j,lmo2))
 969         END DO
 970         DO j = jcs, jce
 971           crk(j,55) = rk(j,55)*dvc(j,ltco3)*vc(j,k,lno2)
 972
 973           dvc(j,lhc3p) = (.83*crk(j,32)+0.5*crk(j,48)+crk(j,51))/ &
 974             (rk(j,58)*vc(j,k,lno)+rk(j,90)*vc(j,k,lho2)+ &
 975             rk(j,116)*vc(j,k,laco3)+rk(j,104)*dvc(j,lmo2))
 976           crk(j,58) = rk(j,58)*dvc(j,lhc3p)*vc(j,k,lno)
 977         END DO
 978         DO j = jcs, jce
 979           dvc(j,lhc5p) = crk(j,33)/(rk(j,60)*vc(j,k,lno)+rk(j,91)*vc(j,k,lho2) &
 980             +rk(j,117)*vc(j,k,laco3)+rk(j,105)*dvc(j,lmo2))
 981           crk(j,60) = rk(j,60)*dvc(j,lhc5p)*vc(j,k,lno)
 982           dvc(j,lhc8p) = crk(j,34)/(rk(j,62)*vc(j,k,lno)+rk(j,92)*vc(j,k,lho2) &
 983             +rk(j,118)*vc(j,k,laco3)+rk(j,106)*dvc(j,lmo2))
 984           crk(j,62) = rk(j,62)*dvc(j,lhc8p)*vc(j,k,lno)
 985           dvc(j,lol2p) = crk(j,35)/(rk(j,64)*vc(j,k,lno)+rk(j,93)*vc(j,k,lho2) &
 986             +rk(j,119)*vc(j,k,laco3)+rk(j,107)*dvc(j,lmo2))
 987           crk(j,64) = rk(j,64)*dvc(j,lol2p)*vc(j,k,lno)
 988           dvc(j,loltp) = (crk(j,36)+crk(j,52))/(rk(j,65)*vc(j,k,lno)+rk(j,94)* &
 989             vc(j,k,lho2)+rk(j,120)*vc(j,k,laco3)+rk(j,108)*dvc(j,lmo2))
 990           crk(j,65) = rk(j,65)*dvc(j,loltp)*vc(j,k,lno)
 991           dvc(j,lolip) = crk(j,37)/(rk(j,66)*vc(j,k,lno)+rk(j,95)*vc(j,k,lho2) &
 992             +rk(j,121)*vc(j,k,laco3)+rk(j,109)*dvc(j,lmo2))
 993           crk(j,66) = rk(j,66)*dvc(j,lolip)*vc(j,k,lno)
 994           crk(j,67) = rk(j,67)*vc(j,k,laco3)*vc(j,k,lno)
 995           crk(j,68) = rk(j,68)*dvc(j,ltco3)*vc(j,k,lno)
 996           dvc(j,ltolp) = 0.75*crk(j,38)/(rk(j,69)*vc(j,k,lno)+rk(j,98)*vc(j,k, &
 997             lho2)+rk(j,124)*vc(j,k,laco3)+rk(j,112)*dvc(j,lmo2))
 998           crk(j,69) = rk(j,69)*dvc(j,ltolp)*vc(j,k,lno)
 999           dvc(j,lxylp) = 0.83*crk(j,39)/(rk(j,70)*vc(j,k,lno)+rk(j,99)*vc(j,k, &
1000             lho2)+rk(j,125)*vc(j,k,laco3)+rk(j,113)*dvc(j,lmo2))
1001
1002           crk(j,70) = rk(j,70)*dvc(j,lxylp)*vc(j,k,lno)
1003           dvc(j,lethp) = (crj(j,16)+crk(j,31))/(rk(j,71)*vc(j,k,lno)+rk(j,89)* &
1004             vc(j,k,lho2)+rk(j,115)*vc(j,k,laco3)+rk(j,103)*dvc(j,lmo2))
1005
1006           crk(j,71) = rk(j,71)*dvc(j,lethp)*vc(j,k,lno)
1007           dvc(j,lketp) = crk(j,43)/(rk(j,72)*vc(j,k,lno)+rk(j,96)*vc(j,k,lho2) &
1008             +rk(j,122)*vc(j,k,laco3)+rk(j,110)*dvc(j,lmo2))
1009
1010         END DO
1011         DO j = jcs, jce
1012           crk(j,72) = rk(j,72)*dvc(j,lketp)*vc(j,k,lno)
1013           crk(j,80) = rk(j,80)*vc(j,k,lol2)*vc(j,k,lno3)
1014           crk(j,81) = rk(j,81)*vc(j,k,lolt)*vc(j,k,lno3)
1015           crk(j,82) = rk(j,82)*vc(j,k,loli)*vc(j,k,lno3)
1016           crk(j,83) = rk(j,83)*vc(j,k,liso)*vc(j,k,lno3)
1017           dvc(j,loln) = (crk(j,80)+crk(j,81)+crk(j,82)+crk(j,83))/ &
1018             (rk(j,73)*vc(j,k,lno)+rk(j,101)*vc(j,k,lho2)+ &
1019             rk(j,138)*dvc(j,lmo2)+rk(j,139)*vc(j,k,laco3)+ &
1020             rk(j,140)*dvc(j,loln))
1021 !         IF (iprt==1 .AND. j==jce-1) THEN
1022 !           PRINT *, dvc(j,loln), crk(j,80), crk(j,81), crk(j,82), crk(j,83), &
1023 !             rk(j,73), vc(j,1,lno), rk(j,101), vc(j,1,lho2), rk(j,138), &
1024 !             dvc(j,lmo2), rk(j,139), vc(j,1,laco3), rk(j,140), dvc(j,loln)
1025 !           PRINT *, dvc(j,loln), crk(j,80), crk(j,81), crk(j,82), crk(j,83), &
1026 !             rk(j,73), vc(j,k,lno), rk(j,101), vc(j,k,lho2), rk(j,138), &
1027 !             dvc(j,lmo2), rk(j,139), vc(j,k,laco3), rk(j,140), dvc(j,loln)
1028 !         END IF
1029
1030         END DO
1031         DO j = jcs, jce
1032           crk(j,73) = rk(j,73)*dvc(j,loln)*vc(j,k,lno)
1033           crk(j,74) = rk(j,74)*vc(j,k,lhcho)*vc(j,k,lno3)
1034           crk(j,75) = rk(j,75)*vc(j,k,lald)*vc(j,k,lno3)
1035           crk(j,76) = rk(j,76)*vc(j,k,lgly)*vc(j,k,lno3)
1036           crk(j,77) = rk(j,77)*vc(j,k,lmgly)*vc(j,k,lno3)
1037           crk(j,79) = rk(j,79)*vc(j,k,lcsl)*vc(j,k,lno3)
1038           crk(j,84) = rk(j,84)*vc(j,k,lol2)*vc(j,k,lo3)
1039           crk(j,85) = rk(j,85)*vc(j,k,lolt)*vc(j,k,lo3)
1040           crk(j,86) = rk(j,86)*vc(j,k,loli)*vc(j,k,lo3)
1041           crk(j,87) = rk(j,87)*vc(j,k,liso)*vc(j,k,lo3)
1042           crk(j,89) = rk(j,89)*vc(j,k,lho2)*dvc(j,lethp)
1043           crk(j,90) = rk(j,90)*vc(j,k,lho2)*dvc(j,lhc3p)
1044           crk(j,91) = rk(j,91)*vc(j,k,lho2)*dvc(j,lhc5p)
1045           crk(j,92) = rk(j,92)*vc(j,k,lho2)*dvc(j,lhc8p)
1046           crk(j,93) = rk(j,93)*vc(j,k,lho2)*dvc(j,lol2p)
1047           crk(j,94) = rk(j,94)*vc(j,k,lho2)*dvc(j,loltp)
1048           crk(j,95) = rk(j,95)*vc(j,k,lho2)*dvc(j,lolip)
1049           crk(j,96) = rk(j,96)*vc(j,k,lho2)*dvc(j,lketp)
1050           crk(j,97) = rk(j,97)*vc(j,k,lho2)*vc(j,k,laco3)
1051           crk(j,98) = rk(j,98)*vc(j,k,lho2)*dvc(j,ltolp)
1052           crk(j,99) = rk(j,99)*vc(j,k,lho2)*dvc(j,lxylp)
1053           crk(j,100) = rk(j,100)*vc(j,k,lho2)*dvc(j,ltco3)
1054           crk(j,101) = rk(j,101)*vc(j,k,lho2)*dvc(j,loln)
1055           crk(j,115) = rk(j,115)*dvc(j,lethp)*vc(j,k,laco3)
1056           crk(j,116) = rk(j,116)*dvc(j,lhc3p)*vc(j,k,laco3)
1057           crk(j,117) = rk(j,117)*dvc(j,lhc5p)*vc(j,k,laco3)
1058           crk(j,118) = rk(j,118)*dvc(j,lhc8p)*vc(j,k,laco3)
1059           crk(j,119) = rk(j,119)*dvc(j,lol2p)*vc(j,k,laco3)
1060           crk(j,120) = rk(j,120)*dvc(j,loltp)*vc(j,k,laco3)
1061           crk(j,121) = rk(j,121)*dvc(j,lolip)*vc(j,k,laco3)
1062           crk(j,122) = rk(j,122)*dvc(j,lketp)*vc(j,k,laco3)
1063           crk(j,123) = rk(j,123)*vc(j,k,laco3)*vc(j,k,laco3)
1064           crk(j,124) = rk(j,124)*vc(j,k,laco3)*dvc(j,ltolp)
1065           crk(j,125) = rk(j,125)*vc(j,k,laco3)*dvc(j,lxylp)
1066           crk(j,126) = rk(j,126)*vc(j,k,laco3)*dvc(j,ltco3)
1067         END DO
1068         DO j = jcs, jce
1069           dvc(j,lxo2) = (0.25*crk(j,33)+0.75*crk(j,34)+0.9*crk(j,40)+crk(j,50) &
1070             +2.0*crk(j,68)+2.0*crk(j,126)+crk(j,114))/ &
1071             (rk(j,131)*vc(j,k,lno)+rk(j,127)*vc(j,k,lho2)+ &
1072             rk(j,128)*dvc(j,lmo2)+rk(j,129)*vc(j,k,laco3))
1073
1074           crk(j,127) = rk(j,127)*dvc(j,lxo2)*vc(j,k,lho2)
1075           crk(j,129) = rk(j,129)*dvc(j,lxo2)*vc(j,k,laco3)
1076           crk(j,130) = 0.
1077           crk(j,131) = rk(j,131)*dvc(j,lxo2)*vc(j,k,lno)
1078           dvc(j,lxno2) = crk(j,79)/(rk(j,132)*vc(j,k,lno2)+rk(j,133)*vc(j,k, &
1079             lho2)+rk(j,134)*dvc(j,lmo2)+rk(j,135)*vc(j,k,laco3))
1080         END DO
1081         DO j = jcs, jce
1082           crk(j,132) = rk(j,132)*dvc(j,lxno2)*vc(j,k,lno2)
1083           crk(j,133) = rk(j,133)*dvc(j,lxno2)*vc(j,k,lho2)
1084           crk(j,135) = rk(j,135)*dvc(j,lxno2)*vc(j,k,laco3)
1085           crk(j,136) = 0.
1086           crk(j,137) = rk(j,137)*vc(j,k,ln2o5)
1087           crk(j,138) = rk(j,138)*dvc(j,lmo2)*dvc(j,loln)
1088           crk(j,139) = rk(j,139)*vc(j,k,laco3)*dvc(j,loln)
1089           crk(j,140) = rk(j,140)*dvc(j,loln)*dvc(j,loln)
1090 !         IF (iprt==1 .AND. j==jce-1) THEN
1091 !           PRINT *, 'crk(j,140),RK(J, 140),DVC(J,LOLN)'
1092 !           PRINT *, crk(j,140), rk(j,140), dvc(j,loln)
1093 !         END IF
1094
1095           dvc(j,lmo2) = (crj(j,12)+crj(j,15)+crk(j,30)+.5*crk(j,47)+crk(j,67)+ &
1096             .22*crk(j,85)+.31*crk(j,86)+.22*crk(j,87)+.5*(crk(j,111)+crk(j, &
1097             115)+crk(j,116)+crk(j,117)+crk(j,118)+crk(j,119)+crk(j,120)+crk(j, &
1098             121)+crk(j,122))+2.0*crk(j,123)+crk(j,124)+crk(j,125)+crk(j,126)+ &
1099             crk(j,129)+crk(j,135)+.5*crk(j,139))/(rk(j,57)*vc(j,k,lno)+rk(j,88 &
1100             )*vc(j,k,lho2)+rk(j,102)*dvc(j,lmo2)+rk(j,103)*dvc(j,lethp)+ &
1101             rk(j,104)*dvc(j,lhc3p)+rk(j,105)*dvc(j,lhc5p)+ &
1102             rk(j,106)*dvc(j,lhc8p)+rk(j,107)*dvc(j,lol2p)+ &
1103             rk(j,108)*dvc(j,loltp)+rk(j,109)*dvc(j,lolip)+ &
1104             rk(j,110)*dvc(j,lketp)+rk(j,111)*vc(j,k,laco3)+ &
1105             rk(j,112)*dvc(j,ltolp)+rk(j,113)*dvc(j,lxylp)+ &
1106             rk(j,114)*dvc(j,ltco3)+rk(j,128)*dvc(j,lxo2)+ &
1107             rk(j,134)*dvc(j,lxno2)+rk(j,138)*dvc(j,loln))
1108         END DO
1109         DO j = jcs, jce
1110           crk(j,57) = rk(j,57)*dvc(j,lmo2)*vc(j,k,lno)
1111           crk(j,88) = rk(j,88)*vc(j,k,lho2)*dvc(j,lmo2)
1112           crk(j,102) = rk(j,102)*dvc(j,lmo2)*dvc(j,lmo2)
1113           crk(j,103) = rk(j,103)*dvc(j,lmo2)*dvc(j,lethp)
1114           crk(j,104) = rk(j,104)*dvc(j,lmo2)*dvc(j,lhc3p)
1115           crk(j,105) = rk(j,105)*dvc(j,lmo2)*dvc(j,lhc5p)
1116           crk(j,106) = rk(j,106)*dvc(j,lmo2)*dvc(j,lhc8p)
1117           crk(j,107) = rk(j,107)*dvc(j,lmo2)*dvc(j,lol2p)
1118           crk(j,108) = rk(j,108)*dvc(j,lmo2)*dvc(j,loltp)
1119           crk(j,109) = rk(j,109)*dvc(j,lmo2)*dvc(j,lolip)
1120           crk(j,110) = rk(j,110)*dvc(j,lmo2)*dvc(j,lketp)
1121           crk(j,111) = rk(j,111)*dvc(j,lmo2)*vc(j,k,laco3)
1122           crk(j,112) = rk(j,112)*dvc(j,lmo2)*dvc(j,ltolp)
1123           crk(j,113) = rk(j,113)*dvc(j,lmo2)*dvc(j,lxylp)
1124           crk(j,114) = rk(j,114)*dvc(j,lmo2)*dvc(j,ltco3)
1125           crk(j,128) = rk(j,128)*dvc(j,lxo2)*dvc(j,lmo2)
1126           crk(j,134) = rk(j,134)*dvc(j,lxno2)*dvc(j,lmo2)
1127         END DO
1128 !       if(iprt.eq.2)then
1129 !         print *,'list of crks ',vc(jcs,1,lho2),vc(jcs,1,lho)
1130 !         do l=1,140
1131 !           print *,'crk at l ',l,crk(1,l)
1132 !         enddo
1133 !         print *,'list of crjs '
1134 !         do l=1,21
1135 !           print *,'crj at l ',l,crj(1,l)
1136 !         enddo
1137 !       endif
1138
1139         RETURN
1140       END SUBROUTINE predraten
1141       SUBROUTINE producn(jcs,jce,iprt,crj,crk,loss,prod,prodl,lossl, &
1142                                                prdrog,iaerosol_sorgam)
1143 !***********************************************************************
1144 !**    DESCRIPTION:
1145 !**     Subroutine PRODUC derives the differential equations from the
1146 !**     reaction rates computed in the PREDRATE subroutine.
1147 !***********************************************************************
1148 ! .. Scalar Arguments ..
1149         implicit none
1150         INTEGER,INTENT(IN)   :: iprt, jce, jcs
1151         real , intent (out)  :: prdrog(jcs:jce,ldrog)
1152         integer, intent (in) :: iaerosol_sorgam
1153
1154 ! ..
1155 ! .. Array Arguments ..
1156         REAL :: crj(jcs:jce,nreacj), crk(jcs:jce,nreack),            &
1157           loss(jcs:jce,lpred), lossl(jcs:jce,lump),                  &
1158           prod(jcs:jce,lpred), prodl(jcs:jce,lump)
1159 ! ..
1160 ! .. Local Scalars ..
1161         REAL    :: alow
1162         INTEGER :: j, l
1163 ! ..
1164 ! .. Intrinsic Functions ..
1165         INTRINSIC max
1166 ! ..
1167         alow = 1.e-17
1168         PRDROG = 0.
1169         DO l = 1, lump
1170           DO j = jcs, jce
1171             prodl(j,l) = 0.0
1172           END DO
1173         END DO
1174         DO l = 1, lpred
1175           DO j = jcs, jce
1176             prod(j,l) = 0.0
1177           END DO
1178         END DO
1179
1180         DO j = jcs, jce
1181           prodl(j,lpao3) = crj(j,16) + crj(j,19) + .02*crj(j,20) + crk(j,42) + &
1182             crk(j,45) + crk(j,49) + .05*crk(j,68) + crk(j,75) + crk(j,77) + &
1183             .03*crk(j,114)
1184         END DO
1185
1186         DO j = jcs, jce
1187           lossl(j,lpao3) = crk(j,67) + crk(j,97) + crk(j,111) + crk(j,115) + &
1188             crk(j,116) + crk(j,117) + crk(j,118) + crk(j,119) + crk(j,120) + &
1189             crk(j,121) + crk(j,122) + 2.*crk(j,123) + crk(j,124) + &
1190             crk(j,125) + crk(j,129) + crk(j,135) + crk(j,50) + &
1191             .95*crk(j,126) + crk(j,139)
1192           lossl(j,lpao3) = max(alow,lossl(j,lpao3))
1193         END DO
1194
1195         DO j = jcs, jce
1196           loss(j,lo3) = crk(j,2) + crk(j,5) + crk(j,6) + crk(j,7) + crk(j,8) + &
1197             crk(j,17) + crk(j,84) + crk(j,85) + crk(j,86) + crk(j,87)
1198         END DO
1199
1200         DO j = jcs, jce
1201           prod(j,lo3) = crj(j,1) + crj(j,8)
1202 !         PROD(J,LO3) = CRK(J,1)
1203         END DO
1204
1205         DO j = jcs, jce
1206           prod(j,lno2) = crj(j,5) + crj(j,6) + crj(j,7) + crj(j,21) + &
1207             crk(j,6) + crk(j,9) + crk(j,11) + 2.*crk(j,16) + 2.*crk(j,18) + &
1208             crk(j,19) + crk(j,22) + crk(j,26) + crk(j,51) + crk(j,54) + &
1209             crk(j,56) + crk(j,57) + .964*crk(j,58) + .92*crk(j,60) + &
1210             .76*crk(j,62) + crk(j,64) + crk(j,65) + crk(j,66) + crk(j,67) + &
1211             crk(j,68) + crk(j,69) + crk(j,70) + crk(j,71) + crk(j,72) + &
1212             crk(j,131) + crk(j,138) + crk(j,139) + 2.0*crk(j,140)
1213 !          IF (iprt==1 .AND. j==jce) THEN
1214 !!       print *,'CRK(J,58),CRK(J,60),CRK(J,62),CRK(J,64),CRK(J,66),'
1215 !!    1         ,'CRK(J,71),CRK(J,138),CRK(J,139),CRK(J,140)'
1216 !            PRINT *, 'CRK(J,58),CRK(J,60),CRK(J,62),CRK(J,64),..'
1217 !            PRINT *, crk(j,58), crk(j,60), crk(j,62), crk(j,64), crk(j,65), &
1218 !              crk(j,66), crk(j,138), crk(j,139), crk(j,140)
1219 !          END IF
1220         END DO
1221
1222         DO j = jcs, jce
1223           loss(j,lno2) = crj(j,1) + crk(j,2) + crk(j,10) + crk(j,17) + &
1224             crk(j,19) + crk(j,21) + crk(j,24) + crk(j,53) + crk(j,55) + &
1225             crk(j,132)
1226 !        if(iprt.eq.1.and.j.eq.jce)then
1227 !           print *,LNO2,LOSS(J,LNO2),CRJ(J,  9),CRK(J, 14)
1228 !        endif
1229         END DO
1230
1231         DO j = jcs, jce
1232           prodl(j,lnox) = crj(j,5) + crj(j,6) + crj(j,8) + crj(j,21) + &
1233             crk(j,11) + crk(j,18) + crk(j,22) + crk(j,26) + crk(j,54) + &
1234             crk(j,56) + crk(j,73) + crj(j,4) + crj(j,7) + crk(j,19) + &
1235             crk(j,51)
1236         END DO
1237
1238         DO j = jcs, jce
1239           lossl(j,lnox) = crk(j,10) + crk(j,17) + crk(j,21) + crk(j,24) + &
1240             crk(j,53) + crk(j,55) + crk(j,132) + crk(j,15) + .036*crk(j,58) + &
1241             .08*crk(j,60) + .24*crk(j,62)
1242         END DO
1243
1244 !....... NTOTAL=NOX+HNO3+NO3+2*N2O5+HONO+HNO4+PAN+TPAN+ONIT
1245 !....... ---->           TAKE THE COMMENTED STUFF
1246 !........ UNCOMMENTED IS: NTOTAL=NOX+HNO3+NO3+2*N2O5+HONO+HNO4
1247 !-------> DON'T FORGET TO CHANGE INTEG AND CHEM
1248
1249
1250 !      LOSSL(J,LNTOTAL)= CRK(J,80)+CRK(J,81)+
1251 !     1                  CRK(J,82)+CRK(J,83)
1252
1253 !      LOSSL(J,LNTOTAL)= CRK(J,53)+CRK(J,55)+
1254 !     1                  CRK(J,132)+
1255 !     1                  CRK(J,80)+CRK(J,81)+
1256 !     1                  CRK(J,82)+CRK(J,83)
1257 !      PRODL(J,LNTOTAL)=
1258 !     1                 CRK(J,20)+
1259 !     1                 CRK(J,50)+CRK(J,54)+CRK(J,56)+CRK(J,73)+
1260 !     1                 CRK(J,51)
1261
1262 !      PRODL(J,LNTOTAL)= CRK(J,101) + CRK(J,73)
1263
1264         DO j = jcs, jce
1265           prodl(j,ln2n3) = crk(j,17) + crk(j,25) + crk(j,50)
1266         END DO
1267
1268         DO j = jcs, jce
1269           lossl(j,ln2n3) = crk(j,23) + crj(j,7) + crj(j,8) + crk(j,18) + &
1270             crk(j,19) + crk(j,20) + crk(j,74) + crk(j,75) + crk(j,76) + &
1271             crk(j,77) + crk(j,78) + crk(j,79) + crk(j,80) + crk(j,81) + &
1272             crk(j,82) + crk(j,83) + crk(j,137)
1273         END DO
1274
1275         DO j = jcs, jce
1276           loss(j,lpan) = crk(j,50) + crk(j,54)
1277         END DO
1278
1279         DO j = jcs, jce
1280           prod(j,lpan) = crk(j,53)
1281         END DO
1282
1283         DO j = jcs, jce
1284           loss(j,lhno3) = crj(j,5) + crk(j,25)
1285         END DO
1286
1287         DO j = jcs, jce
1288           prod(j,lhno3) = crk(j,20) + 2.D0*crk(j,23) + crk(j,24) + crk(j,74) + &
1289             crk(j,75) + crk(j,76) + crk(j,77) + crk(j,78) + crk(j,79) + &
1290             2.*crk(j,137)
1291         END DO
1292
1293         DO j = jcs, jce
1294           loss(j,lh2o2) = max(alow,crj(j,9) + crk(j,14) )
1295 !        if(iprt.eq.1.and.j.eq.jce)then
1296 !           print *,LH2O2,LOSS(J,LH2O2),CRJ(J,  9),CRK(J, 14)
1297 !        endif
1298         END DO
1299
1300         DO j = jcs, jce
1301           prod(j,lh2o2) = crk(j,12) + crk(j,13)
1302 !        if(iprt.eq.1.and.j.eq.jce)then
1303 !           print *,LH2O2,prod(J,LH2O2),CRK(J, 12),CRK(J, 13)
1304 !        endif
1305         END DO
1306
1307         DO j = jcs, jce
1308           loss(j,lhcho) = crj(j,10) + crj(j,11) + crk(j,41) + crk(j,74)
1309         END DO
1310
1311         DO j = jcs, jce
1312           prod(j,lhcho) = crj(j,13) + .13*crj(j,17) + .45*crj(j,18) + &
1313             .009*crk(j,32) + .5*crk(j,47) + crk(j,50) + crk(j,57) + &
1314             .09*crk(j,58) + .04*crk(j,62) + 1.6*crk(j,64) + crk(j,65) + &
1315             .28*crk(j,66) + crk(j,73) + crk(j,84) + .53*crk(j,85) + &
1316             .18*crk(j,86) + .53*crk(j,87) + 1.5*crk(j,102) + .75*crk(j,103) + &
1317             .75*crk(j,104) + .77*crk(j,105) + .80*crk(j,106) + &
1318             1.55*crk(j,107) + 1.25*crk(j,108) + .89*crk(j,109) + &
1319             .75*crk(j,110) + crk(j,111) + crk(j,112) + crk(j,113) + &
1320             .5*crk(j,114) + .8*crk(j,119) + .5*crk(j,120) + .14*crk(j,121) + &
1321             crk(j,128) + crk(j,134) + 1.75*crk(j,138) + crk(j,139) + &
1322             2.0*crk(j,140)
1323         END DO
1324
1325         DO j = jcs, jce
1326           prod(j,lhono) = crk(j,15)
1327         END DO
1328
1329         DO j = jcs, jce
1330           loss(j,lhono) = crj(j,4)
1331         END DO
1332
1333         DO j = jcs, jce
1334           prod(j,lhno4) = crk(j,10)
1335         END DO
1336
1337         DO j = jcs, jce
1338           loss(j,lhno4) = crj(j,6) + crk(j,11) + crk(j,26)
1339         END DO
1340
1341         DO j = jcs, jce
1342           prod(j,ln2o5) = crk(j,21)
1343         END DO
1344
1345         DO j = jcs, jce
1346           loss(j,ln2o5) = crk(j,22) + crk(j,23) + crk(j,137)
1347         END DO
1348
1349         DO j = jcs, jce
1350           prod(j,lno3) = crk(j,17) + crk(j,22) + crk(j,25) + crk(j,50)
1351         END DO
1352
1353         DO j = jcs, jce
1354           loss(j,lno3) = crj(j,7) + crj(j,8) + crk(j,18) + crk(j,19) + &
1355             crk(j,20) + crk(j,21) + crk(j,74) + crk(j,75) + crk(j,76) + &
1356             crk(j,77) + crk(j,78) + crk(j,79) + crk(j,80) + crk(j,81) + &
1357             crk(j,82) + crk(j,83)
1358         END DO
1359
1360         DO j = jcs, jce
1361           loss(j,lco) = crk(j,29)
1362         END DO
1363
1364         DO j = jcs, jce
1365           prod(j,lco) = crj(j,10) + crj(j,11) + crj(j,12) + 1.87*crj(j,17) + &
1366             1.55*crj(j,18) + crj(j,19) + crk(j,41) + 2.*crk(j,44) + &
1367             crk(j,45) + .95*crk(j,68) + crk(j,74) + 2.*crk(j,76) + crk(j,77) + &
1368             .42*crk(j,84) + .33*crk(j,85) + .23*crk(j,86) + .33*crk(j,87) + &
1369             .475*crk(j,114) + .95*crk(j,126)
1370         END DO
1371
1372         DO j = jcs, jce
1373           loss(j,lald) = crj(j,12) + crk(j,42) + crk(j,75)
1374         END DO
1375
1376         DO j = jcs, jce
1377           prod(j,lald) = crj(j,14) + .075*crk(j,32) + .2*crj(j,21) + &
1378             .5*crk(j,48) + .75*crk(j,58) + .38*crk(j,60) + .35*crk(j,62) + &
1379             .2*crk(j,64) + crk(j,65) + 1.45*crk(j,66) + crk(j,73) + &
1380             crk(j,71) + .5*crk(j,85) + .72*crk(j,86) + .5*crk(j,87) + &
1381             .75*crk(j,103) + .15*crk(j,104) + .41*crk(j,105) + &
1382             .46*crk(j,106) + .35*crk(j,107) + .75*crk(j,108) + &
1383             .725*crk(j,109) + crk(j,115) + .2*crk(j,116) + .14*crk(j,117) + &
1384             .1*crk(j,118) + .6*crk(j,119) + crk(j,120) + .725*crk(j,121) + &
1385             crk(j,138) + crk(j,139) + 2.0*crk(j,140)
1386         END DO
1387
1388         DO j = jcs, jce
1389           loss(j,lop1) = crj(j,13) + crk(j,47)
1390         END DO
1391
1392         DO j = jcs, jce
1393           prod(j,lop1) = crk(j,88)
1394         END DO
1395
1396         DO j = jcs, jce
1397           loss(j,lop2) = crj(j,14) + crk(j,48)
1398         END DO
1399
1400         DO j = jcs, jce
1401           prod(j,lop2) = crk(j,89) + crk(j,90) + crk(j,91) + crk(j,92) + &
1402             crk(j,93) + crk(j,94) + crk(j,95) + crk(j,96) + crk(j,98) + &
1403             crk(j,99) + crk(j,100) + crk(j,127) + crk(j,133)
1404         END DO
1405
1406         DO j = jcs, jce
1407           loss(j,lpaa) = crj(j,15) + crk(j,49)
1408         END DO
1409
1410         DO j = jcs, jce
1411           prod(j,lpaa) = crk(j,97)
1412         END DO
1413
1414         DO j = jcs, jce
1415           loss(j,lket) = crj(j,16) + crk(j,43)
1416         END DO
1417
1418         DO j = jcs, jce
1419           prod(j,lket) = .8*crj(j,21) + .025*crk(j,32) + .25*crk(j,58) + &
1420             .69*crk(j,60) + 1.06*crk(j,62) + .10*crk(j,66) + .10*crk(j,86) + &
1421             .6*crk(j,104) + .75*crk(j,105) + 1.39*crk(j,106) + &
1422             .55*crk(j,109) + .8*crk(j,116) + .86*crk(j,117) + .9*crk(j,118) + &
1423             .55*crk(j,121)
1424         END DO
1425
1426         DO j = jcs, jce
1427           loss(j,lgly) = crj(j,17) + crj(j,18) + crk(j,44) + crk(j,76)
1428         END DO
1429
1430         DO j = jcs, jce
1431           prod(j,lgly) = .89*crk(j,68) + .16*crk(j,69) + .16*crk(j,112) + &
1432             .44*crk(j,114) + .2*crk(j,124) + .89*crk(j,126)
1433         END DO
1434
1435         DO j = jcs, jce
1436           loss(j,lmgly) = crj(j,19) + crk(j,45) + crk(j,77)
1437         END DO
1438
1439         DO j = jcs, jce
1440           prod(j,lmgly) = .11*crk(j,68) + .17*crk(j,69) + .450*crk(j,70) + &
1441             crk(j,72) + .75*crk(j,110) + .17*crk(j,112) + .45*crk(j,113) + &
1442             .05*crk(j,114) + crk(j,122) + .8*crk(j,124) + crk(j,125) + &
1443             .11*crk(j,126)
1444         END DO
1445
1446         DO j = jcs, jce
1447           loss(j,ldcb) = crj(j,20) + crk(j,46) + crk(j,78)
1448         END DO
1449
1450         DO j = jcs, jce
1451           loss(j,ldcb) = max(alow,loss(j,ldcb))
1452         END DO
1453
1454         DO j = jcs, jce
1455           prod(j,ldcb) = .70*crk(j,69) + .806*crk(j,70) + .7*crk(j,112) + &
1456             .806*crk(j,113) + crk(j,124) + crk(j,125)
1457         END DO
1458
1459         DO j = jcs, jce
1460           loss(j,lonit) = crj(j,21) + crk(j,51)
1461         END DO
1462
1463         DO j = jcs, jce
1464           prod(j,lonit) = .036*crk(j,58) + .08*crk(j,60) + .24*crk(j,62) + &
1465             crk(j,101) + crk(j,132)
1466         END DO
1467
1468         DO j = jcs, jce
1469           loss(j,lso2) = crk(j,28)
1470         END DO
1471
1472         DO j = jcs, jce
1473           loss(j,lsulf) = 0.
1474         END DO
1475
1476         DO j = jcs, jce
1477           prod(j,lsulf) = crk(j,28)
1478 !         if(iprt==2)print *,' j,prod = ',j,prod(j,lsulf)
1479         END DO
1480
1481         DO j = jcs, jce
1482           loss(j,leth) = crk(j,31)
1483         END DO
1484
1485         DO j = jcs, jce
1486           loss(j,lhc3) = crk(j,32)
1487         END DO
1488
1489         DO j = jcs, jce
1490           loss(j,lhc5) = crk(j,33)
1491         END DO
1492
1493         DO j = jcs, jce
1494           loss(j,lhc8) = crk(j,34)
1495         END DO
1496
1497         DO j = jcs, jce
1498           loss(j,lol2) = crk(j,35) + crk(j,80) + crk(j,84)
1499         END DO
1500
1501         DO j = jcs, jce
1502           loss(j,lolt) = crk(j,36) + crk(j,81) + crk(j,85)
1503         END DO
1504
1505         DO j = jcs, jce
1506           loss(j,loli) = crk(j,37) + crk(j,82) + crk(j,86)
1507         END DO
1508
1509         DO j = jcs, jce
1510           loss(j,ltol) = crk(j,38)
1511         END DO
1512
1513         DO j = jcs, jce
1514           loss(j,lcsl) = crk(j,40) + .5*crk(j,79)
1515         END DO
1516
1517         DO j = jcs, jce
1518           prod(j,lcsl) = .25*crk(j,38) + .17*crk(j,39)
1519         END DO
1520
1521         DO j = jcs, jce
1522           loss(j,lxyl) = crk(j,39)
1523         END DO
1524
1525         DO j = jcs, jce
1526           loss(j,laco3) = crk(j,53) + crk(j,67) + crk(j,97) + crk(j,111) + &
1527             crk(j,115) + crk(j,116) + crk(j,117) + crk(j,118) + crk(j,119) + &
1528             crk(j,120) + crk(j,121) + crk(j,122) + 2.*crk(j,123) + &
1529             crk(j,124) + crk(j,125) + .95*crk(j,126) + crk(j,129) + &
1530             crk(j,135) + crk(j,139)
1531         END DO
1532
1533         DO j = jcs, jce
1534           prod(j,laco3) = crj(j,16) + crj(j,19) + .02*crj(j,20) + crk(j,42) + &
1535             crk(j,45) + crk(j,49) + crk(j,54) + .05*crk(j,68) + crk(j,75) + &
1536             crk(j,77) + .03*crk(j,114)
1537         END DO
1538
1539         DO j = jcs, jce
1540           loss(j,liso) = crk(j,52) + crk(j,83) + crk(j,87)
1541         END DO
1542
1543         DO j = jcs, jce
1544           loss(j,ltpan) = crk(j,56)
1545         END DO
1546
1547         DO j = jcs, jce
1548           prod(j,ltpan) = crk(j,55)
1549         END DO
1550
1551         DO j = jcs, jce
1552           loss(j,lora1) = 1.E-27
1553         END DO
1554
1555         DO j = jcs, jce
1556           prod(j,lora1) = .4*crk(j,84) + .06*crk(j,86) + .2*crk(j,85) + &
1557             .2*crk(j,87)
1558         END DO
1559
1560         DO j = jcs, jce
1561           loss(j,lora2) = 1.E-27
1562         END DO
1563
1564         DO j = jcs, jce
1565           prod(j,lora2) = .2*crk(j,85) + .29*crk(j,86) + .2*crk(j,87) + &
1566             .5*crk(j,111) + .5*crk(j,114) + .5*crk(j,115) + .5*crk(j,116) + &
1567             .5*crk(j,117) + .5*crk(j,118) + .5*crk(j,119) + .5*crk(j,120) + &
1568             .5*crk(j,121) + .5*crk(j,122) + .5*crk(j,139)
1569         END DO
1570
1571         DO j = jcs, jce
1572           lossl(j,lhox) = crk(j,15) + crk(j,24) + crk(j,25) + crk(j,26) + &
1573             crk(j,27) + crk(j,30) + crk(j,31) + .83*crk(j,32) + crk(j,33) + &
1574             crk(j,34) + crk(j,35) + crk(j,36) + crk(j,37) + .75*crk(j,38) + &
1575             .83*crk(j,39) + 1.8*crk(j,40) + crk(j,42) + crk(j,43) + &
1576             crk(j,45) + crk(j,46) + crk(j,49) + crk(j,50) + crk(j,51) + &
1577             crk(j,52) + crk(j,10) + 2.*crk(j,12) + 2.*crk(j,13) + crk(j,20) + &
1578             crk(j,27) + crk(j,88) + crk(j,89) + crk(j,90) + crk(j,91) + &
1579             .5*crk(j,47) + .5*crk(j,48) + crk(j,92) + crk(j,93) + crk(j,94) + &
1580             crk(j,95) + crk(j,96) + crk(j,97) + crk(j,98) + crk(j,99) + &
1581             crk(j,100) + crk(j,101) + crk(j,127) + crk(j,133)
1582           lossl(j,lhox) = max(alow,lossl(j,lhox))
1583         END DO
1584
1585         DO j = jcs, jce
1586           prodl(j,lhox) = crj(j,4) + crj(j,5) + crj(j,6) + 2.*crj(j,9) + &
1587             crj(j,13) + crj(j,14) + crj(j,15) + 2.*crk(j,5) + 2.*crj(j,11) + &
1588             crj(j,12) + crj(j,13) + crj(j,14) + .8*crj(j,18) + crj(j,19) + &
1589             .98*crj(j,20) + crj(j,21) + crk(j,11) + crk(j,57) + &
1590             .964*crk(j,58) + .92*crk(j,60) + .76*crk(j,62) + crk(j,64) + &
1591             crk(j,65) + crk(j,66) + .92*crk(j,68) + crk(j,69) + crk(j,70) + &
1592             crk(j,71) + crk(j,72) + crk(j,74) + crk(j,76) + .12*crk(j,84) + &
1593             .33*crk(j,85) + .40*crk(j,86) + .33*crk(j,87) + crk(j,102) + &
1594             crk(j,103) + crk(j,104) + crk(j,105) + crk(j,106) + crk(j,107) + &
1595             crk(j,108) + crk(j,109) + crk(j,110) + .5*crk(j,111) + &
1596             2.0*crk(j,112) + 2.*crk(j,113) + .46*crk(j,114) + .5*crk(j,115) + &
1597             .5*crk(j,116) + .5*crk(j,117) + .5*crk(j,118) + .5*crk(j,119) + &
1598             .5*crk(j,120) + .5*crk(j,121) + .5*crk(j,122) + crk(j,124) + &
1599             crk(j,125) + .92*crk(j,126) + crk(j,128) + crk(j,134) + &
1600             .5*crk(j,138)
1601         END DO
1602
1603 !      DO 850 L=1 ,LPRED
1604 !         DO 850 J=JCS,JCE
1605 !            PROD(J,L)= PROD(J,L) + PRODS(J,L)
1606 !850   CONTINUE
1607
1608 !      DO 900 J=JCS,JCE
1609 !        PRODL(J,LNOX)   = PRODL(J,LNOX)+ PRODS(J,LNO) + PRODS(J,LNO2)
1610 !        PRODL(J,LNTOTAL)= PRODL(J,LNTOTAL)+PRODS(J,LNO)+PRODS(J,LNO2)
1611 !900   CONTINUE
1612       DO J = JCS, JCE
1613          PRDROG(J,PXYL)  =        CRK(J, 39)
1614          PRDROG(J,PTOL)  =        CRK(J, 38)
1615          PRDROG(J,PCSL1) =        CRK(J, 40)
1616          PRDROG(J,PCSL2) = 0.50 * CRK(J, 79)
1617          PRDROG(J,PHC8)  =        CRK(J, 34)
1618          PRDROG(J,POLI1) =        CRK(J, 37)
1619          PRDROG(J,POLI2) =        CRK(J, 82)
1620          PRDROG(J,POLI3) =        CRK(J, 86)
1621          PRDROG(J,POLT1) =        CRK(J, 36)
1622          PRDROG(J,POLT2) =        CRK(J, 81)
1623          PRDROG(J,POLT3) =        CRK(J, 85)
1624 !
1625 ! next lines for radm only, RACM would be different
1626 !
1627          PRDROG(J,PAPI1) =        0.
1628          PRDROG(J,PAPI2) =        0.
1629          PRDROG(J,PAPI3) =        0.
1630          PRDROG(J,PLIM1) =        0.
1631          PRDROG(J,PLIM2) =        0.
1632          PRDROG(J,PLIM3) =        0.
1633       ENDDO
1634
1635         RETURN
1636       END SUBROUTINE producn
1637       SUBROUTINE setdtc(jcs,jce,dtc,dtcmax,dtcmin,dt60,prod,loss,vc,   &
1638                                                              timenow )
1639         implicit none
1640         REAL, PARAMETER :: huge=1.e10
1641 ! .. Scalar Arguments ..
1642         REAL,  intent(in)    :: dt60, dtcmax, dtcmin, timenow
1643         INTEGER, intent(in)  :: jce, jcs
1644         REAL,  intent(in)    :: loss(jcs:jce,lpred),                    &
1645                                 prod(jcs:jce,lpred), vc(jcs:jce,1,lspec)
1646         real,  intent(inout) :: dtc
1647
1648 ! ..
1649 ! ..
1650 ! .. Local Scalars ..
1651         INTEGER :: j, k, l
1652 ! ..
1653 ! .. Local Arrays ..
1654         REAL :: dtlsp(lspec), dum(jcs:jce)
1655 ! ..
1656 ! .. Intrinsic Functions ..
1657         INTRINSIC abs, max, min
1658 ! ..
1659 ! ..
1660         k = 1
1661
1662         DO l = 1, lspec
1663           dtlsp(l) = huge
1664         END DO
1665         DO l = 1, lpred
1666           IF (qdtc(l)==1) THEN
1667             DO j = jcs, jce
1668               dum(j) = prod(j,l) - loss(j,l)
1669 !             dum(j) = max(abs(dum(j)),epsilc)
1670               dum(j) = max(abs(dum(j)),1.e-30)
1671               dum(j) = .02*vc(j,1,l)/dum(j)
1672 !              DUM(J) = CVMGP(DUM(J),HUGE,VC(J,1,L)-epsilc*100.)
1673               IF (vc(j,1,l)-1.e-10>=0.) THEN
1674                 dum(j) = dum(j)
1675               ELSE
1676                 dum(j) = huge
1677               END IF
1678
1679             END DO
1680             DO j = jcs, jce
1681               dtlsp(l) = min(dtlsp(l),dum(j))
1682             END DO
1683           END IF
1684         END DO
1685 !       IF (dtc<=dtcmax*.9) THEN
1686 !         dtc = dtc*1.1
1687 !       ELSE
1688 !         dtc = dtcmax
1689 !       END IF
1690           dtc = dtcmax
1691         DO l = 1, lpred
1692           IF (qdtc(l)==1) THEN
1693             IF (dtlsp(l)<dtc) THEN
1694               dtc = dtlsp(l)
1695             END IF
1696           END IF
1697         END DO
1698         IF (dtc<dtcmin) THEN
1699           dtc = dtcmin
1700         END IF
1701         IF ((timenow+dtc)>dt60) dtc = dt60 - timenow
1702         RETURN
1703       END SUBROUTINE setdtc
1704       SUBROUTINE chemin
1705        implicit none
1706 ! .. Scalar Arguments ..
1707         RETURN
1708       END SUBROUTINE chemin
1709
1710     END MODULE module_radm