Aerosol interaction vn14.2

src/activation.F90

@@ -261,14 +261,12 @@ subroutine activate(dt, cloud_mass, cloud_number, w, rho, dnumber, dmac, T, p,
 
         activated_arg = sum(dnccn_all(:))
 
+        dactive = zero_real_wp
+        dmass_d = zero_real_wp
+
         if (Smax_act > smax_act_min .and. .not. l_warm) then
 
-          if (iopt_inuc < 4) then
-            ! For lower-order ice nucleation options, need to initialise
-            ! dactive and dmass_d
-            dactive = zero_real_wp
-            dmass_d = zero_real_wp
-          else
+          if (iopt_inuc > 3) then
             ! For higher-order ice nucleation schemes, dactive and dmass_d
             ! can be based on dustphys
             dactive = 0.01*dustphys%N(1)

src/aerosol_routines.F90

@@ -14,6 +14,9 @@ module aerosol_routines
   use lognormal_funcs, only: MNtoRm ! DPG - added this for MNtoRm since was
                                     ! causing circular conflicts as wanted
                                     ! to use it in which_mode_to_use.F90
+  USE umprintmgr,  ONLY: umPrint, umMessage
+
+
 
   implicit none
 
@@ -209,8 +212,8 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, &
 
     character(len=*), parameter :: RoutineName='EXAMINE_AEROSOL'
 
-    real, parameter :: aero_mact_mean_max=2.3e-15 !kg 3xsigma+mean(0.3micron) density 1777
-    real, parameter :: dust_mact_mean_max=3.1e-13 !kg 3xsigma+mean(1micron) density 1777
+    real, parameter :: aero_mact_mean_max=1e30!3e-13 !kg -make bigger for coarse soluble aerosol
+    real, parameter :: dust_mact_mean_max=1e30!3e-13 !kg 3xsigma+mean(1micron) density 1777
     ! When l_bypass_whichmode is True (or no larger modes exist to move aerosol to)
     ! aerosol mean sizes in hydrometeors can become so large that they lead to instability
     ! during hydrometeor sedimentation. Limiting the maximum mean mass stops these edge
@@ -315,7 +318,18 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, &
             aeroact(k)%mact=mac
             aeroact(k)%rcrit=0.0
             aeroact(k)%mact_mean=aeroact(k)%mact /(aeroact(k)%nact + epsilon(ntot))
-            aeroact(k)%mact_mean=min(aero_mact_mean_max, aeroact(k)%mact_mean)
+            if (aeroact(k)%mact_mean .gt. aero_mact_mean_max)then
+!                               WRITE(umMessage, *) 'PRFx1 ',aeroact(k)%mact_mean,mac,ntot,ntot*aero_mact_mean_max/mac
+!                               CALL umPrint( umMessage) 
+              aeroact(k)%mact_mean=min(aero_mact_mean_max, aeroact(k)%mact_mean)
+!              aeroact(k)%mact=aero_mact_mean_max*aeroact(k)%nact
+!              mac=aeroact(k)%mact
+
+              aeroact(k)%nact=aeroact(k)%mact/aero_mact_mean_max
+              ntot=aeroact(k)%nact
+
+
+            endif
             ! Get mean radius of distribution
             rm_arc=MNtoRm(aeroact(k)%mact,aeroact(k)%nact,density,sigma_arc)
 
@@ -340,15 +354,15 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, &
             aeroact(k)%rcrit2=rcrit2
             aeroact(k)%mact2=mact2
             aeroact(k)%mact2_mean=aeroact(k)%mact2/(aeroact(k)%nact2 + epsilon(ntot))
-            aeroact(k)%mact2_mean=min(aero_mact_mean_max, aeroact(k)%mact2_mean)
+!            aeroact(k)%mact2_mean=min(aero_mact_mean_max, aeroact(k)%mact2_mean)
 
           end if
 
           aeroact(k)%nact1=max(0.0_wp, aeroact(k)%nact-aeroact(k)%nact2)
           aeroact(k)%mact1=max(0.0_wp, aeroact(k)%mact-aeroact(k)%mact2)
           aeroact(k)%rcrit1=0.0
           aeroact(k)%mact1_mean=aeroact(k)%mact1/(aeroact(k)%nact1+epsilon(mar))
-          aeroact(k)%mact1_mean=min(aero_mact_mean_max, aeroact(k)%mact1_mean)
+!          aeroact(k)%mact1_mean=min(aero_mact_mean_max, aeroact(k)%mact1_mean)
 
           if (cloud_number > epsilon(1.0_wp)) then
             aeroact(k)%nratio1=max(0.0,min(1.0,aeroact(k)%nact1/cloud_number))
@@ -457,14 +471,27 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, &
             dustact(k)%mact=mad
             dustact(k)%rcrit=0.0
             dustact(k)%mact_mean=dustact(k)%mact/(dustact(k)%nact+epsilon(mad))
-            dustact(k)%mact_mean=min(dust_mact_mean_max, dustact(k)%mact_mean)
+!            dustact(k)%mact_mean=min(dust_mact_mean_max, dustact(k)%mact_mean)
+            if (dustact(k)%mact_mean .gt. dust_mact_mean_max)then
+!                               WRITE(umMessage, *) 'PRFx2 ',dustact(k)%mact_mean,mad,nitot,nitot*dust_mact_mean_max/mad
+!                               CALL umPrint( umMessage) 
+              dustact(k)%mact_mean=min(dust_mact_mean_max, dustact(k)%mact_mean)
+!              dustact(k)%mact=dust_mact_mean_max*dustact(k)%nact
+!              mad=dustact(k)%mact
+
+              dustact(k)%nact=dustact(k)%mact/dust_mact_mean_max
+              nitot=dustact(k)%nact
+
+
+            endif
 
             if (snow_number > epsilon(1.0_wp)) then
               dustact(k)%nact2=nitot*ratio_s
               dustact(k)%rcrit2=0.0
               dustact(k)%mact2=mad *ratio_s
               dustact(k)%mact2_mean=dustact(k)%mact2/(dustact(k)%nact2+epsilon(nhtot))
-              dustact(k)%mact2_mean=min(dust_mact_mean_max, dustact(k)%mact2_mean)
+!              dustact(k)%mact2_mean=min(dust_mact_mean_max, dustact(k)%mact2_mean)
+
               dustact(k)%nratio2=max(0.0,min(1.0,dustact(k)%nact2/(snow_number+epsilon(nhtot)) ))
             end if
 
@@ -473,7 +500,7 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, &
               dustact(k)%rcrit3=0.0
               dustact(k)%mact3=mad *ratio_g
               dustact(k)%mact3_mean=dustact(k)%mact3/(dustact(k)%nact3+epsilon(nhtot))
-              dustact(k)%mact3_mean=min(dust_mact_mean_max, dustact(k)%mact3_mean)
+!              dustact(k)%mact3_mean=min(dust_mact_mean_max, dustact(k)%mact3_mean)
               dustact(k)%nratio3=max(0.0,min(1.0,dustact(k)%nact3/(graupel_number+epsilon(nhtot)) ))
             end if
 
@@ -482,7 +509,8 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, &
                dustact(k)%mact1=max(0.0_wp, dustact(k)%mact-dustact(k)%mact2-dustact(k)%mact3)
                dustact(k)%rcrit1=0.0
                dustact(k)%mact1_mean=dustact(k)%mact1/(dustact(k)%nact1+epsilon(mar))
-               dustact(k)%mact1_mean=min(dust_mact_mean_max, dustact(k)%mact1_mean)
+!               dustact(k)%mact1_mean=min(dust_mact_mean_max, dustact(k)%mact1_mean)
+
                dustact(k)%nratio1=max(0.0,min(1.0,dustact(k)%nact1/(ice_number+epsilon(nhtot)) ))
             end if
 
@@ -559,14 +587,25 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, &
               aeroice(k)%mact=maai
               aeroice(k)%rcrit=0.0
               aeroice(k)%mact_mean=aeroice(k)%mact/(aeroice(k)%nact+epsilon(nhtot))
-              aeroice(k)%mact_mean=min(aero_mact_mean_max, aeroice(k)%mact_mean)
+!              aeroice(k)%mact_mean=min(aero_mact_mean_max, aeroice(k)%mact_mean)
+              if (aeroice(k)%mact_mean .gt. aero_mact_mean_max)then
+!                               WRITE(umMessage, *) 'PRFx3 ',aeroice(k)%mact_mean,maai,nitot,nitot*aero_mact_mean_max/maai
+!                               CALL umPrint( umMessage) 
+                aeroice(k)%mact_mean=min(aero_mact_mean_max, aeroice(k)%mact_mean)
+!                aeroice(k)%mact=aero_mact_mean_max*aeroice(k)%nact
+!                maai=aeroice(k)%mact
+              aeroice(k)%nact=aeroice(k)%mact/aero_mact_mean_max
+              nitot=aeroice(k)%nact
+
+              endif
+
 
               if (ratio_s > epsilon(1.0_wp)) then
                 aeroice(k)%nact2=nitot*ratio_s
                 aeroice(k)%rcrit2=0.0
                 aeroice(k)%mact2=maai*ratio_s
                 aeroice(k)%mact2_mean=aeroice(k)%mact2/(aeroice(k)%nact2+epsilon(nhtot))
-                aeroice(k)%mact2_mean=min(aero_mact_mean_max, aeroice(k)%mact2_mean)
+!                aeroice(k)%mact2_mean=min(aero_mact_mean_max, aeroice(k)%mact2_mean)
                 aeroice(k)%nratio2=max(0.0, min(1.0,aeroice(k)%nact2/(snow_number+epsilon(nhtot)) ))
               end if
 
@@ -647,6 +686,16 @@ subroutine examine_aerosol(aerofields, qfields, aerophys, aerochem, aeroact, &
               dustliq(k)%rcrit=0.0
               dustliq(k)%mact_mean=dustliq(k)%mact/(dustliq(k)%nact+epsilon(nhtot))
               dustliq(k)%mact_mean=min(dust_mact_mean_max, dustliq(k)%mact_mean)
+              if (dustliq(k)%mact_mean .gt. dust_mact_mean_max)then
+!                               WRITE(umMessage, *) 'PRFx4 ',dustliq(k)%mact_mean,madl,ntot,ntot*aero_mact_mean_max/madl
+!                               CALL umPrint( umMessage) 
+                dustliq(k)%mact_mean=min(dust_mact_mean_max, dustliq(k)%mact_mean)
+!                dustliq(k)%mact=dust_mact_mean_max*dustliq(k)%nact
+!                madl=dustliq(k)%mact
+                dustliq(k)%nact=dustliq(k)%mact/dust_mact_mean_max
+                ntot=dustliq(k)%nact
+
+              endif
 
               if (ratio_r > epsilon(1.0_wp)) then
                 dustliq(k)%nact2=ntot*ratio_r

src/mphys_tidy.F90

@@ -2,7 +2,8 @@ module mphys_tidy
   use variable_precision, only: wp
   use process_routines, only: process_rate,  process_name
   use aerosol_routines, only: aerosol_active
-  use thresholds, only: thresh_tidy, thresh_atidy
+  use thresholds, only: thresh_tidy, thresh_atidy,&
+      aeromass_small, aeronumber_small
   use passive_fields, only: exner, pressure
   use mphys_switches, only:                       &
        i_qv, i_ql, i_nl, i_qr, i_nr, i_m3r, i_th, &
@@ -14,6 +15,8 @@ module mphys_tidy
        i_an2, i_am2, i_am4, i_am5, l_warm,        &
        i_an6, i_am6, i_am7, i_am8, i_am9,         &
        i_an11, i_an12,                            &
+       i_an1,i_am1, i_an2, i_am2, i_an3, i_am3,   &
+       i_am6, i_an6, i_am10, i_an10,              &
        l_process, ntotalq, ntotala,               &
        i_qstart, i_nstart, i_m3start,             &
        l_separate_rain, l_tidy_conserve_E, l_tidy_conserve_q, &
@@ -208,12 +211,16 @@ subroutine qtidy(ixy_inner, dt, nz, qfields, procs, aerofields, aeroact, dustact
     !--------------------------------------------------------------------------
     IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle)
 
-    l_qsig(0)=.false.
-    l_qsneg(0)=.false.
+    l_qsig(:)=.false.
+    l_qsneg(:)=.false. !should this be l_qsneg(:)
 
     do k = 1, nz
     nr_reset=.false.
+    qr_reset=.false.
     m3r_reset=.false.
+    nl_reset=.false.
+    ql_reset=.false.
+    ni_reset=.false.
     qi_reset=.false.
     qs_reset=.false.
     ns_reset=.false.
@@ -865,12 +872,88 @@ subroutine tidy_ain(qfields, aerofields)
       end if
       if (i_an11 > 0) then
         if (((qfields(k, i_ql) + qfields(k,i_qr) + qfields(k,i_qi) + qfields(k,i_qs) + qfields(k,i_qg) &
-                                                                 <=0.0 .and. aerofields(k,i_an12)>0.0) &
+                                                                 <=0.0 .and. aerofields(k,i_an11)>0.0) &
                                                             .or. aerofields(k,i_an11) < 0.0)) then
           aerofields(k,i_an11)=0.0
         end if
       end if
 
+
+
+!        aerofields(:, i_am1)=a1(ks:ke,i,j)      AitkenSolMass         
+!        aerofields(:, i_an1)=a2(ks:ke,i,j)       AitkenSolNumber         
+!        aerofields(:, i_am2)=a3(ks:ke,i,j)      AccumSolMass         
+!        aerofields(:, i_an2)=a4(ks:ke,i,j)       AccumSolNumber
+!        aerofields(:, i_am3)=a5(ks:ke,i,j)      CoarseSolMass
+!        aerofields(:, i_an3)=a6(ks:ke,i,j)       CoarseSolNumber
+!        aerofields(:, i_am4)=a7(ks:ke,i,j)      ActSolLiq_casim
+!        aerofields(:, i_am5)=a8(ks:ke,i,j)      ActSolRain_casim              <-this is not currently available
+!        aerofields(:, i_am6)=a9(ks:ke,i,j)      CoarseDustMass
+!        aerofields(:, i_an6)=a10(ks:ke,i,j)     CoarseDustNumber
+
+       IF (i_am1 > 0) then
+         if (aerofields(k, i_am1) < aeromass_small) then
+           aerofields(k, i_am1)=0.0
+           aerofields(k, i_an1)=0.0
+         end if
+       END IF
+       IF (i_am2 > 0) then
+         if (aerofields(k, i_am2) < aeromass_small) then
+           aerofields(k, i_am2)=0.0
+           aerofields(k, i_an2)=0.0
+         end if
+       END IF
+       IF (i_am3 > 0) then
+         if (aerofields(k, i_am3) < aeromass_small) then
+           aerofields(k, i_am3)=0.0
+           aerofields(k, i_an3)=0.0
+         end if
+       END IF
+       IF (i_am6 > 0) then
+         if (aerofields(k, i_am6) < aeromass_small) then
+           aerofields(k, i_am6)=0.0
+           aerofields(k, i_an6)=0.0
+         end if
+       END IF
+       IF (i_an10 > 0) then
+         if (aerofields(k, i_am10) < aeromass_small) then
+           aerofields(k, i_am10)=0.0
+           aerofields(k, i_an10)=0.0
+         end if
+       END IF
+
+       IF (i_am1 > 0) then
+         if (aerofields(k, i_an1) < aeronumber_small) then
+           aerofields(k, i_am1)=0.0
+           aerofields(k, i_an1)=0.0
+         end if
+       END IF
+       IF (i_am2 > 0) then
+         if (aerofields(k, i_an2) < aeronumber_small) then
+           aerofields(k, i_am2)=0.0
+           aerofields(k, i_an2)=0.0
+         end if
+       END IF
+       IF (i_am3 > 0) then
+         if (aerofields(k, i_an3) < aeronumber_small) then
+           aerofields(k, i_am3)=0.0
+           aerofields(k, i_an3)=0.0
+         end if
+       END IF
+       IF (i_am6 > 0) then
+         if (aerofields(k, i_an6) < aeronumber_small) then
+           aerofields(k, i_am6)=0.0
+           aerofields(k, i_an6)=0.0
+         end if
+       END IF
+       IF (i_am10 > 0) then
+         if (aerofields(k, i_an10) < aeronumber_small) then
+           aerofields(k, i_am10)=0.0
+           aerofields(k, i_an10)=0.0
+         end if
+       END IF
+
+
     end do
 
     IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle)