Changeset 462f74b in flexpart.git for src/releaseparticles.f90

Timestamp:

Jul 12, 2016, 11:02:42 AM (8 years ago)

Author:

Sabine Eckhardt <sabine@…>

Branches:

master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie

Children:

842074e

Parents:

f28aa0a

Message:

first version of the backward scavenging

File:

• src/releaseparticles.f90 (modified) (5 diffs)

src/releaseparticles.f90

@@ -57 +57 @@
   real :: xaux,yaux,zaux,rfraction
   real :: topo,rhoaux(2),r,t,rhoout,ddx,ddy,rddx,rddy,p1,p2,p3,p4
+  real :: rhosum(nspec)
   real :: dz1,dz2,dz,xtn,ytn,xlonav,timecorrect(maxspec),press,pressold
   real :: presspart,average_timecorrect
@@ -86 +87 @@
   minpart=1
   do i=1,numpoint
+    do k=1,nspec
+      rhosum(k)=0
+    end do
     if ((itime.ge.ireleasestart(i)).and. &! are we within release interval?
          (itime.le.ireleaseend(i))) then
@@ -176 +180 @@
   ! scaled. Adjust the mass per particle by the species-dependent time correction factor
   ! divided by the species-average one
+  ! for the scavenging calculation the mass needs to be multiplied with rho of the particle layer and
+  ! divided by the sum of rho of all particles.
   !*****************************************************************************
             do k=1,nspec
                xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
                     *timecorrect(k)/average_timecorrect
-  !            write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k
+              if (SCAVDEP) then ! if there is no scavenging in wetdepo it will be set to 0
+!              if ( henry(k).gt.0 .or. &
+!                   crain_aero(k).gt.0. .or. csnow_aero(k).gt.0. .or. &
+!                   ccn_aero(k).gt.0. .or. in_aero(k).gt.0. )  then
+                xscav_frac1(ipart,k)=-1.
+!               write(*,*) '190: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i
+!               xscav_frac1(ipart,k)=(-1.)/real(npart(i)) &
+!                    *timecorrect(k)/average_timecorrect
+!                  else
+!                     xscav_frac1(ipart,k)=0
+!                  endif   
+               endif
   ! Assign certain properties to particle
   !**************************************
@@ -372 +389 @@
               do k=1,nspec
                 xmass1(ipart,k)=xmass1(ipart,k)*rhoout
+                if (SCAVDEP) then
+                     xscav_frac1(ipart,k)=xscav_frac1(ipart,k)
+!mctest                     xscav_frac1(ipart,k)=xscav_frac1(ipart,k)*rhoout
+                     rhosum(k)=rhosum(k)+rhoout
+!               write(*,*) '391: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i
+                endif
               end do
             endif
@@ -379 +402 @@
             goto 34      ! Storage space has been found, stop searching
           endif
-        end do
+        end do  ! i=1:numpoint
         if (ipart.gt.maxpart) goto 996
 
 34      minpart=ipart+1
-      end do
+      end do ! ipart=minpart,maxpart
+      if (SCAVDEP) then
+         do ipart=minpart,maxpart
+            do k=1,nspec
+              if (xscav_frac1(ipart,k).lt.0) then
+!mctest                   xscav_frac1(ipart,k)=xscav_frac1(ipart,k)/rhosum(k)
+!                 write(*,*) '409: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i
+               endif   
+            end do
+         end do
       endif
+      endif ! j=1,numrel
   end do