Bugfix: using namelist species files would sometimes set incorrect below-cloud...

Bugfix: using namelist species files would sometimes set incorrect below-cloud scavenging parameters when using more than 1 species. Update: per-species scavenging statistics are written out at end of simulation.

src/FLEXPART.f90

@@ -413,10 +413,15 @@ program flexpart
   call timemanager
 
 ! NIK 16.02.2005 
+  do i=1,nspec
     write(*,*) '**********************************************'
-  write(*,*) 'Total number of occurences of below-cloud scavenging', tot_blc_count
-  write(*,*) 'Total number of occurences of in-cloud    scavenging', tot_inc_count
+    write(*,*) 'Scavenging statistics for species ', species(i), ':'
+    write(*,*) 'Total number of occurences of below-cloud scavenging', &
+         & tot_blc_count(i)
+    write(*,*) 'Total number of occurences of in-cloud    scavenging', &
+         & tot_inc_count(i)
     write(*,*) '**********************************************'
+  end do
   
   write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
        &XPART MODEL RUN!'

src/FLEXPART_MPI.f90

@@ -461,26 +461,29 @@ program flexpart
 
 ! NIK 16.02.2005 
   if (lroot) then
-    call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, 1, MPI_INTEGER8, MPI_SUM, id_root, &
+    call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
          & mp_comm_used, mp_ierr)
-    call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, 1, MPI_INTEGER8, MPI_SUM, id_root, &
+    call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
          & mp_comm_used, mp_ierr)
   else
     if (mp_partgroup_pid.ge.0) then ! Skip for readwind process 
-      call MPI_Reduce(tot_blc_count, 0, 1, MPI_INTEGER8, MPI_SUM, id_root, &
+      call MPI_Reduce(tot_blc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
            & mp_comm_used, mp_ierr)
-      call MPI_Reduce(tot_inc_count, 0, 1, MPI_INTEGER8, MPI_SUM, id_root, &
+      call MPI_Reduce(tot_inc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
            & mp_comm_used, mp_ierr)
     end if
   end if
 
   if (lroot) then
+    do i=1,nspec
       write(*,*) '**********************************************'
+      write(*,*) 'Scavenging statistics for species ', species(i), ':'
       write(*,*) 'Total number of occurences of below-cloud scavenging', &
-         & tot_blc_count
+           & tot_blc_count(i)
       write(*,*) 'Total number of occurences of in-cloud    scavenging', &
-         & tot_inc_count
+           & tot_inc_count(i)
       write(*,*) '**********************************************'
+    end do
 
     write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
          &XPART MODEL RUN!'

src/com_mod.f90

@@ -139,7 +139,8 @@ module com_mod
   
 
 !NIK 16.02.2015
-  integer(selected_int_kind(16)) :: tot_blc_count=0, tot_inc_count=0
+  integer(selected_int_kind(16)), dimension(maxspec) :: tot_blc_count=0, &
+       &tot_inc_count=0
 
 
   !*********************************************************************
@@ -575,7 +576,8 @@ module com_mod
   integer :: numxgridn,numygridn
   real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn
   !real outheight(maxzgrid),outheighthalf(maxzgrid)
-  logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC
+  logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,WETDEPSPEC(maxspec),&
+       & OHREA,ASSSPEC
 
   ! numxgrid,numygrid       number of grid points in x,y-direction
   ! numxgridn,numygridn     number of grid points in x,y-direction for nested output grid
@@ -592,6 +594,7 @@ module com_mod
   ! DRYDEP                  .true., if dry deposition is switched on
   ! DRYDEPSPEC              .true., if dry deposition is switched on for that species
   ! WETDEP                  .true., if wet deposition is switched on
+  ! WETDEPSPEC              .true., if wet deposition is switched on for that species
   ! OHREA                   .true., if OH reaction is switched on
   ! ASSSPEC                 .true., if there are two species asscoiated
   !                    (i.e. transfer of mass between these two occurs

src/readreleases.f90

@@ -275,6 +275,7 @@ subroutine readreleases
   OHREA=.false.
   do i=1,maxspec
     DRYDEPSPEC(i)=.false.
+    WETDEPSPEC(i)=.false.
   end do
 
   if (readerror.ne.0) then
@@ -368,6 +369,7 @@ subroutine readreleases
     if ((dquer(i).le.0..and.(weta_gas(i).gt.0. .or. wetb_gas(i).gt.0.)) .or. &
          &(dquer(i).gt.0. .and. (crain_aero(i) .gt. 0. .or. csnow_aero(i).gt.0.)))  then
       WETDEP=.true.
+      WETDEPSPEC(i)=.true.
       if (lroot) then
         write (*,*) '  Below-cloud scavenging: ON'
 !  write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
@@ -379,6 +381,7 @@ subroutine readreleases
 ! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
     if (dquer(i).gt.0..and.(ccn_aero(i).gt.0. .or. in_aero(i).gt.0.))  then
       WETDEP=.true.
+      WETDEPSPEC(i)=.true.
       if (lroot) then
         write (*,*) '  In-cloud scavenging: ON'
 !        write (*,*) 'In-cloud scavenging coefficients: ',&
@@ -398,7 +401,7 @@ subroutine readreleases
       DRYDEPSPEC(i)=.true.
     endif
 
-  end do
+  end do ! end loop over species
 
   if (WETDEP.or.DRYDEP) DEP=.true.
 

src/readspecies.f90

@@ -196,8 +196,8 @@ subroutine readspecies(id_spec,pos_spec)
     decay(pos_spec)=pdecay
     weta_gas(pos_spec)=pweta_gas
     wetb_gas(pos_spec)=pwetb_gas
-    crain_aero=pcrain_aero
-    csnow_aero=pcsnow_aero
+    crain_aero(pos_spec)=pcrain_aero
+    csnow_aero(pos_spec)=pcsnow_aero
     ccn_aero(pos_spec)=pccn_aero
     in_aero(pos_spec)=pin_aero
     reldiff(pos_spec)=preldiff

src/wetdepo.f90

@@ -89,7 +89,7 @@ subroutine wetdepo(itime,ltsample,loutnext)
   real, parameter :: bcls(6) = (/22.7, 0.0, 0.0, 1321.0, 381.0, 0.0/) !now (Kyro et al 2009)
   real :: frac_act, liq_frac, dquer_m
 
-  integer :: blc_count, inc_count
+  integer(selected_int_kind(16)), dimension(nspec) :: blc_count, inc_count
   real    :: Si_dummy, wetscav_dummy
   logical :: readclouds_this_nest
 
@@ -106,8 +106,8 @@ subroutine wetdepo(itime,ltsample,loutnext)
 ! Loop over all particles
 !************************
 
-  blc_count=0
-  inc_count=0
+  blc_count(:)=0
+  inc_count(:)=0
 
   do jpart=1,numpart
 
@@ -257,6 +257,9 @@ subroutine wetdepo(itime,ltsample,loutnext)
       wetdeposit(ks)=0. 
       wetscav=0.
 
+! Cycle loop if wet deposition for the species is off
+      if (WETDEPSPEC(ks).eqv..false.) cycle
+
       if (ngrid.gt.0) then
         act_temp=ttn(ix,jy,hz,n,ngrid)
       else
@@ -273,13 +276,13 @@ subroutine wetdepo(itime,ltsample,loutnext)
 !******************************************************************
         if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
           !        if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then 
-          blc_count=blc_count+1
+          blc_count(ks)=blc_count(ks)+1
           wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
 
 ! For aerosols: if positive below-cloud parameters (Crain/Csnow or B), and dquer>0
 !*********************************************************************************
         else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then
-          blc_count=blc_count+1
+          blc_count(ks)=blc_count(ks)+1
 
 !NIK 17.02.2015
 ! For the calculation here particle size needs to be in meter and not um as dquer is
@@ -319,7 +322,7 @@ subroutine wetdepo(itime,ltsample,loutnext)
 ! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are
 ! given in species file, or if gas and positive Henry's constant
         if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then 
-          inc_count=inc_count+1
+          inc_count(ks)=inc_count(ks)+1
 ! if negative coefficients (turned off) set to zero for use in equation
           if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
           if (in_aero(ks).lt.0.) in_aero(ks)=0.
@@ -431,7 +434,7 @@ subroutine wetdepo(itime,ltsample,loutnext)
   end do ! all particles
 
 ! count the total number of below-cloud and in-cloud occurences:
-  tot_blc_count=tot_blc_count+blc_count
-  tot_inc_count=tot_inc_count+inc_count
+  tot_blc_count(1:nspec)=tot_blc_count(1:nspec)+blc_count(1:nspec)
+  tot_inc_count(1:nspec)=tot_inc_count(1:nspec)+inc_count(1:nspec)
 
 end subroutine wetdepo