else         do k = 1, pver            do i = 1, ncol               sulfmix(i,k) = 0.            end do         end do      endif      if ( indirect ) then ! Method of Martin et. al.         do k=pver,1,-1            do i = 1,ncol               locrhoair(i,k) = pmid(i,k)/( Rdryair*t(i,k) )               lwcwat(i,k) = ( qm1(i,k,ixcldw)*(1.-fice(i,k))/max(0.01_r8,cld(i,k)) )* &                             locrhoair(i,k)!                 NOTE: 0.001 converts kg/m3 -> g/cm3               so4mass(i,k) = sulfmix(i,k)*locrhoair(i,k)*0.001               Aso4(i,k) = so4mass(i,k)*Acoef               if (Aso4(i,k) <= 280.0) then                  Aso4(i,k) = max(36.0_r8,Aso4(i,k))                  Ntot(i,k) = -1.15E-3*Aso4(i,k)**2 + 0.963*Aso4(i,k)+5.30                  rekappa = 0.80               else                  Aso4(i,k) = min(1500.0_r8,Aso4(i,k))                  Ntot(i,k) = -2.10E-4*Aso4(i,k)**2 + 0.568*Aso4(i,k)-27.9                  rekappa = 0.67               end if               if (land(i)) then ! Account for local background aerosol;                  bgaer = Cland*exp(-(zm(i,k)/Hland))                  Ntot(i,k) = max(bgaer,Ntot(i,k))               else                  bgaer = Cmarn*exp(-(zm(i,k)/Hmarn))                  Ntot(i,k) = max(bgaer,Ntot(i,k))               end if!               if (k == pver) then                  Ntotb = Ntot(i,k)               else                  Ntotb = Ntot(i,k+1)               end if!               relmod(i,k) = (( (recoef*lwcwat(i,k))/(rekappa*Ntotb))**reexp)*10000.0               relmod(i,k) = max(4.0_r8,relmod(i,k))               relmod(i,k) = min(20.0_r8,relmod(i,k))               if (cld(i,k) >= 0.01) then                  cldfrq(i,k) = 1.0               else                  cldfrq(i,k) = 0.0               end if               wrel(i,k) = relmod(i,k)*cldfrq(i,k)               wlwc(i,k) = lwcwat(i,k)*cldfrq(i,k)            end do         end do      else         do k = 1, pver            do i = 1, ncol               relmod(i,k) = rel(i,k)            end do         end do      end if!! Specify aerosol mass mixing ratio!      call aermix(lchnk   ,ncol    ,pnm     ,sulfmix ,aermmr  ,rh      )      call t_startf('radcswmx')      call radcswmx(lchnk   ,ncol    ,                            &                    pnm     ,pbr     ,qm1     ,rh      ,o3mmr   , &                    aermmr  ,cld     ,clwp    ,rel     ,rei     , &                    fice    ,eccf    ,coszrs  ,scon    ,solin   , &                    asdir   ,asdif   ,aldir   ,aldif   ,nmxrgn  , &                    pmxrgn  ,qrs     ,fsnt    ,fsntc   ,fsntoa  , &                    fsntoac ,fsnirt  ,fsnrtc  ,fsnirtsq,fsns    , &                    fsnsc   ,fsdsc   ,fsds    ,sols    ,soll    , &                    solsd   ,solld   )      call t_stopf('radcswmx')      call outfld('AERMMR  ',aermmr, pcols,lchnk)      call outfld('REL     ',relmod ,pcols,lchnk)      if ( indirect ) then         call outfld('MSO4    ',so4mass,pcols,lchnk)         call outfld('LWC     ',lwcwat ,pcols,lchnk)         call outfld('CLDFRQ  ',cldfrq ,pcols,lchnk)         call outfld('WREL    ',wrel   ,pcols,lchnk)         call outfld('WLWC    ',wlwc   ,pcols,lchnk)      end if! -- tls ---------------------------------------------------------------2!! Convert units of shortwave fields needed by rest of model from CGS to MKS!      do i=1,ncol         solin(i) = solin(i)*1.e-3         fsds(i)  = fsds(i)*1.e-3         fsnirt(i)= fsnirt(i)*1.e-3         fsnrtc(i)= fsnrtc(i)*1.e-3         fsnirtsq(i)= fsnirtsq(i)*1.e-3         fsnt(i)  = fsnt(i) *1.e-3         fsns(i)  = fsns(i) *1.e-3         fsntc(i) = fsntc(i)*1.e-3         fsnsc(i) = fsnsc(i)*1.e-3         fsdsc(i) = fsdsc(i)*1.e-3         fsntoa(i)=fsntoa(i)*1.e-3         fsntoac(i)=fsntoac(i)*1.e-3      end do!! Dump shortwave radiation information to history tape buffer (diagnostics)!      ftem(:ncol,:pver) = qrs(:ncol,:pver)/cpair      call outfld('QRS     ',ftem  ,pcols,lchnk)      call outfld('SOLIN   ',solin ,pcols,lchnk)      call outfld('FSDS    ',fsds  ,pcols,lchnk)      call outfld('FSNIRTOA',fsnirt,pcols,lchnk)      call outfld('FSNRTOAC',fsnrtc,pcols,lchnk)      call outfld('FSNRTOAS',fsnirtsq,pcols,lchnk)      call outfld('FSNT    ',fsnt  ,pcols,lchnk)      call outfld('FSNS    ',fsns  ,pcols,lchnk)      call outfld('FSNTC   ',fsntc ,pcols,lchnk)      call outfld('FSNSC   ',fsnsc ,pcols,lchnk)      call outfld('FSDSC   ',fsdsc ,pcols,lchnk)      call outfld('FSNTOA  ',fsntoa,pcols,lchnk)      call outfld('FSNTOAC ',fsntoac,pcols,lchnk)      call outfld('SOLS    ',sols  ,pcols,lchnk)      call outfld('SOLL    ',soll  ,pcols,lchnk)      call outfld('SOLSD   ',solsd ,pcols,lchnk)      call outfld('SOLLD   ',solld ,pcols,lchnk)!   end if!! Longwave radiation computation!   if (dolw) then!! Convert upward longwave flux units to CGS!      do i=1,ncol         lwupcgs(i) = lwup(i)*1000.      end do!! Do longwave computation. If not implementing greenhouse gas code then! first specify trace gas mixing ratios. If greenhouse gas code then:!  o ixtrcg   => indx of advected n2o tracer!  o ixtrcg+1 => indx of advected ch4 tracer!  o ixtrcg+2 => indx of advected cfc11 tracer!  o ixtrcg+3 => indx of advected cfc12 tracer!      if (trace_gas) then         call cnst_get_ind('N2O'  , in2o)         call cnst_get_ind('CH4'  , ich4)         call cnst_get_ind('CFC11', if11)         call cnst_get_ind('CFC12', if12)         call t_startf("radclwmx")         call radclwmx(lchnk   ,ncol    ,                            &                       lwupcgs ,t       ,qm1(1,1,1)       ,o3vmr ,   &                       pbr     ,pnm     ,pmln    ,piln    ,          &                       qm1(1,1,in2o)    ,qm1(1,1,ich4)    ,          &                       qm1(1,1,if11)    ,qm1(1,1,if12)    ,          &                       cld     ,emis    ,pmxrgn  ,nmxrgn  ,qrl     , &                       flns    ,flnt    ,flnsc   ,flntc   ,flwds   , &                       flut    ,flutc   )         call t_stopf("radclwmx")      else         call trcmix(lchnk   ,ncol    , &                     pmid    ,n2o     ,ch4     ,                     &                     cfc11   ,cfc12   )         call t_startf("radclwmx")         call radclwmx(lchnk     ,ncol    ,                            &                       lwupcgs   ,t       ,qm1(1,1,1)       ,o3vmr ,   &                       pbr       ,pnm     ,pmln    ,piln    ,          &                       n2o       ,ch4     ,cfc11   ,cfc12   ,          &                       cld       ,emis    ,pmxrgn  ,nmxrgn  ,qrl     , &                       flns      ,flnt    ,flnsc   ,flntc   ,flwds   , &                       flut      ,flutc   )         call t_stopf("radclwmx")      endif!! Convert units of longwave fields needed by rest of model from CGS to MKS!      do i=1,ncol         flnt(i)  = flnt(i)*1.e-3         flut(i)  = flut(i)*1.e-3         flutc(i) = flutc(i)*1.e-3         flns(i)  = flns(i)*1.e-3         flntc(i) = flntc(i)*1.e-3         flnsc(i) = flnsc(i)*1.e-3         flwds(i) = flwds(i)*1.e-3         lwcf(i)=flutc(i) - flut(i)         swcf(i)=fsntoa(i) - fsntoac(i)      end do!! Dump longwave radiation information to history tape buffer (diagnostics)!      call outfld('QRL     ',qrl/cpair ,pcols,lchnk)      call outfld('FLNT    ',flnt  ,pcols,lchnk)      call outfld('FLUT    ',flut  ,pcols,lchnk)      call outfld('FLUTC   ',flutc ,pcols,lchnk)      call outfld('FLNTC   ',flntc ,pcols,lchnk)      call outfld('FLNS    ',flns  ,pcols,lchnk)      call outfld('FLNSC   ',flnsc ,pcols,lchnk)      call outfld('LWCF    ',lwcf  ,pcols,lchnk)      call outfld('SWCF    ',swcf  ,pcols,lchnk)!   end if!   returnend subroutine radctl