Commit f75967de authored by Espen Sollum's avatar Espen Sollum
Browse files

Fixed a bug with number of occurances of wet scavenging written to stdout

parent d6a09777
......@@ -435,22 +435,22 @@ program flexpart
call mpif_finalize
stop
endif
call timemanager
! NIK 16.02.2005
if (lroot) then
call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, 1, mp_pp, MPI_SUM, id_root, &
call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, 1, MPI_INTEGER, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, 1, mp_pp, MPI_SUM, id_root, &
call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, 1, MPI_INTEGER, 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, tot_blc_count, 1, mp_pp, MPI_SUM, id_root, &
call MPI_Reduce(tot_blc_count, tot_blc_count, 1, MPI_INTEGER, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
call MPI_Reduce(tot_inc_count, tot_inc_count, 1, mp_pp, MPI_SUM, id_root, &
call MPI_Reduce(tot_inc_count, tot_inc_count, 1, MPI_INTEGER, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
end if
end if
......
......@@ -131,9 +131,12 @@ module com_mod
!ZHG SEP 2015 wheather or not to read clouds from GRIB
logical :: readclouds
!ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed)
logical :: sumclouds
!NIK 16.02.2015
integer :: tot_blc_count, tot_inc_count
integer :: tot_blc_count=0, tot_inc_count=0
......@@ -745,7 +748,7 @@ module com_mod
! files in cases where differences with MPI version is minor (eso)
!*****************************************************************
integer :: mpi_mode=0 ! .gt. 0 if running MPI version
logical :: lroot=.true. ! true if serial version, or if MPI and root process
logical :: lroot=.true. ! true if serial version, or if MPI .and. root process
contains
subroutine com_mod_allocate(nmpart)
......
......@@ -137,9 +137,9 @@ subroutine timemanager
!ZHG 2015
!CGZ-lifetime: set lifetime to 0
checklifetime(:,:)=0
species_lifetime(:,:)=0
print*, 'Initialized lifetime'
! checklifetime(:,:)=0
! species_lifetime(:,:)=0
! print*, 'Initialized lifetime'
!CGZ-lifetime: set lifetime to 0
......@@ -421,9 +421,9 @@ subroutine timemanager
!CGZ-lifetime: output species lifetime
!ZHG
write(*,*) 'Overview species lifetime in days', &
real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
write(*,*) 'all info:',species_lifetime
! write(*,*) 'Overview species lifetime in days', &
! real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
! write(*,*) 'all info:',species_lifetime
!ZHG
!CGZ-lifetime: output species lifetime
......@@ -593,12 +593,12 @@ subroutine timemanager
xmass1(j,ks)/xmass(npoint(j),ks))
!ZHG 2015
!CGZ-lifetime: Check mass fraction left/save lifetime
if(real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.0.01.and.checklifetime(j,ks).eq.0.)then
! if(real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.0.01.and.checklifetime(j,ks).eq.0.)then
!Mass below 1% of initial >register lifetime
checklifetime(j,ks)=abs(itra1(j)-itramem(j))
species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
species_lifetime(ks,2)= species_lifetime(ks,2)+1
endif
! checklifetime(j,ks)=abs(itra1(j)-itramem(j))
! species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
! species_lifetime(ks,2)= species_lifetime(ks,2)+1
! endif
!CGZ-lifetime: Check mass fraction left/save lifetime
!ZHG 2015
else
......
......@@ -150,9 +150,9 @@ subroutine timemanager
end if ! (lroot)
!CGZ-lifetime: set lifetime to 0
checklifetime(:,:)=0
species_lifetime(:,:)=0
print*, 'Initialized lifetime'
! checklifetime(:,:)=0
! species_lifetime(:,:)=0
! print*, 'Initialized lifetime'
!CGZ-lifetime: set lifetime to 0
......@@ -563,7 +563,7 @@ subroutine timemanager
endif
!CGZ-lifetime: output species lifetime
! if (lroot) then
if (lroot) then
! write(*,*) 'Overview species lifetime in days', &
! real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
! write(*,*) 'all info:',species_lifetime
......@@ -572,7 +572,7 @@ subroutine timemanager
! if (verbosity.gt.0) then
! write (*,*) 'timemanager> starting simulation'
! end if
! end if
end if
45 format(i13,' SECONDS SIMULATED: ',i13, ' PARTICLES: Uncertainty: ',3f7.3)
46 format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
......@@ -756,13 +756,13 @@ subroutine timemanager
xmass1(j,ks)/xmass(npoint(j),ks))
!CGZ-lifetime: Check mass fraction left/save lifetime
if(lroot.and.real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.0.01.and.checklifetime(j,ks).eq.0.)then
! if(lroot.and.real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.0.01.and.checklifetime(j,ks).eq.0.)then
!Mass below 1% of initial >register lifetime
checklifetime(j,ks)=abs(itra1(j)-itramem(j))
! checklifetime(j,ks)=abs(itra1(j)-itramem(j))
species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
species_lifetime(ks,2)= species_lifetime(ks,2)+1
endif
! species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
! species_lifetime(ks,2)= species_lifetime(ks,2)+1
! endif
!CGZ-lifetime: Check mass fraction left/save lifetime
endif
......
......@@ -567,7 +567,7 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!***********************************************************************************
if (readclouds) then !HG METHOD
if (readclouds) then !HG METHOD
! The method is loops all grids vertically and constructs the 3D matrix for clouds
! Cloud top and cloud bottom gid cells are assigned as well as the total column
! cloud water. For precipitating grids, the type and whether it is in or below
......@@ -576,104 +576,105 @@ if (readclouds) then !HG METHOD
! to scavenging. Also clouds that are not precipitating are defined which may be
! to include future cloud processing by non-precipitating-clouds.
!***********************************************************************************
if (readclouds) write(*,*) 'using cloud water from ECMWF'
clw(:,:,:,n)=0
icloud_stats(:,:,:,n)=0
clouds(:,:,:,n)=0
do jy=0,nymin1
do ix=0,nxmin1
lsp=lsprec(ix,jy,1,n)
convp=convprec(ix,jy,1,n)
prec=lsp+convp
tot_cloud_h=0
! Find clouds in the vertical
do kz=1, nz-1 !go from top to bottom
if (clwc(ix,jy,kz,n).gt.0) then
! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n) ! Column cloud water [m3/m3]
icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz)) ! Cloud BOT height stats [m]
write(*,*) 'using cloud water from ECMWF'
clw(:,:,:,n)=0
icloud_stats(:,:,:,n)=0
clouds(:,:,:,n)=0
do jy=0,nymin1
do ix=0,nxmin1
lsp=lsprec(ix,jy,1,n)
convp=convprec(ix,jy,1,n)
prec=lsp+convp
tot_cloud_h=0
! Find clouds in the vertical
do kz=1, nz-1 !go from top to bottom
if (clwc(ix,jy,kz,n).gt.0) then
! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n) ! Column cloud water [m3/m3]
icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz)) ! Cloud BOT height stats [m]
!ZHG 2015 extra for testing
clh(ix,jy,kz,n)=height(kz+1)-height(kz)
icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz)) ! Cloud TOP height [m]
clh(ix,jy,kz,n)=height(kz+1)-height(kz)
icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz)) ! Cloud TOP height [m]
!ZHG
endif
end do
endif
end do
! If Precipitation. Define removal type in the vertical
if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
! If Precipitation. Define removal type in the vertical
if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
do kz=nz,1,-1 !go Bottom up!
if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
clouds(ix,jy,kz,n)=1 ! is a cloud
if (lsp.ge.convp) then
clouds(ix,jy,kz,n)=3 ! lsp in-cloud
else
clouds(ix,jy,kz,n)=2 ! convp in-cloud
endif ! convective or large scale
elseif((clw(ix,jy,kz,n).le.0) .and. (icloud_stats(ix,jy,3,n).ge.height(kz)) ) then ! is below cloud
if (lsp.ge.convp) then
clouds(ix,jy,kz,n)=5 ! lsp dominated washout
else
clouds(ix,jy,kz,n)=4 ! convp dominated washout
endif ! convective or large scale
endif
do kz=nz,1,-1 !go Bottom up!
if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
clouds(ix,jy,kz,n)=1 ! is a cloud
if (lsp.ge.convp) then
clouds(ix,jy,kz,n)=3 ! lsp in-cloud
else
clouds(ix,jy,kz,n)=2 ! convp in-cloud
endif ! convective or large scale
elseif((clw(ix,jy,kz,n).le.0) .and. (icloud_stats(ix,jy,3,n).ge.height(kz)) ) then ! is below cloud
if (lsp.ge.convp) then
clouds(ix,jy,kz,n)=5 ! lsp dominated washout
else
clouds(ix,jy,kz,n)=4 ! convp dominated washout
endif ! convective or large scale
endif
if (height(kz).ge. 19000) then ! set a max height for removal
clouds(ix,jy,kz,n)=0
endif !clw>0
end do !nz
endif ! precipitation
end do
end do
if (height(kz).ge. 19000) then ! set a max height for removal
clouds(ix,jy,kz,n)=0
endif !clw>0
end do !nz
endif ! precipitation
end do
end do
!**************************************************************************
else ! use old definitions
else ! use old definitions
!**************************************************************************
! create a cloud and rainout/washout field, clouds occur where rh>80%
! total cloudheight is stored at level 0
if (.not.readclouds) write(*,*) 'using cloud water from Parameterization'
do jy=0,nymin1
do ix=0,nxmin1
! if (.not.readclouds) write(*,*) 'using cloud water from Parameterization'
write(*,*) 'using cloud water from Parameterization'
do jy=0,nymin1
do ix=0,nxmin1
! OLD METHOD
rain_cloud_above(ix,jy)=0
lsp=lsprec(ix,jy,1,n)
convp=convprec(ix,jy,1,n)
cloudsh(ix,jy,n)=0
do kz_inv=1,nz-1
kz=nz-kz_inv+1
pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
clouds(ix,jy,kz,n)=0
if (rh.gt.0.8) then ! in cloud
if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
rain_cloud_above(ix,jy)=1
cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
height(kz)-height(kz-1)
if (lsp.ge.convp) then
clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
else
clouds(ix,jy,kz,n)=2 ! convp dominated rainout
rain_cloud_above(ix,jy)=0
lsp=lsprec(ix,jy,1,n)
convp=convprec(ix,jy,1,n)
cloudsh(ix,jy,n)=0
do kz_inv=1,nz-1
kz=nz-kz_inv+1
pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
clouds(ix,jy,kz,n)=0
if (rh.gt.0.8) then ! in cloud
if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
rain_cloud_above(ix,jy)=1
cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
height(kz)-height(kz-1)
if (lsp.ge.convp) then
clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
else
clouds(ix,jy,kz,n)=2 ! convp dominated rainout
endif
else ! no precipitation
clouds(ix,jy,kz,n)=1 ! cloud
endif
else ! no precipitation
clouds(ix,jy,kz,n)=1 ! cloud
endif
else ! no cloud
if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
if (lsp.ge.convp) then
clouds(ix,jy,kz,n)=5 ! lsp dominated washout
else
clouds(ix,jy,kz,n)=4 ! convp dominated washout
else ! no cloud
if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
if (lsp.ge.convp) then
clouds(ix,jy,kz,n)=5 ! lsp dominated washout
else
clouds(ix,jy,kz,n)=4 ! convp dominated washout
endif
endif
endif
endif
end do
end do
!END OLD METHOD
end do
end do
endif !readclouds
end do
endif !readclouds
!********* TEST ***************
! WRITE OUT SOME TEST VARIABLES
......@@ -837,5 +838,5 @@ endif !readclouds
! sum(rho(:,:,:,n)*rho(:,:,:,n)),sum(drhodz(:,:,:,n)*drhodz(:,:,:,n)),&
! sum(ww(:,:,:,n)*ww(:,:,:,n)), &
! sum(clouds(:,:,:,n)), sum(cloudsh(:,:,n)),sum(idx),sum(pinmconv)
end subroutine verttransform
end subroutine verttransform
......@@ -72,7 +72,7 @@ subroutine wetdepo(itime,ltsample,loutnext)
integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
integer :: ks, kp
integer :: n1,n2, icbot,ictop, indcloud !TEST
! integer :: n1,n2, icbot,ictop, indcloud !TEST
real :: S_i, act_temp, cl, cle ! in cloud scavenging
real :: clouds_h ! cloud height for the specific grid point
real :: xtn,ytn,lsp,convp,cc,grfraction(3),prec(3),wetscav,totprec
......@@ -105,7 +105,6 @@ subroutine wetdepo(itime,ltsample,loutnext)
! Loop over all particles
!************************
blc_count=0
inc_count=0
......@@ -184,7 +183,7 @@ subroutine wetdepo(itime,ltsample,loutnext)
n=memind(2)
if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
n=memind(1)
n=memind(1)
if (ngrid.eq.0) then
clouds_v=clouds(ix,jy,hz,n)
......@@ -198,7 +197,7 @@ subroutine wetdepo(itime,ltsample,loutnext)
! if there is no precipitation or the particle is above the clouds no
! scavenging is done
if (clouds_v.le.1) goto 20
if (clouds_v.le.1) goto 20
! 1) Parameterization of the the area fraction of the grid cell where the
! precipitation occurs: the absolute limit is the total cloud cover, but
......@@ -310,6 +309,7 @@ subroutine wetdepo(itime,ltsample,loutnext)
if (clouds_v.lt.4) then ! In-cloud
! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are given in species file
if (weta_in(ks).gt.0. .or. wetb_in(ks).gt.0.) then
inc_count=inc_count+1
! if negative coefficients (turned off) set to zero for use in equation
if (weta_in(ks).lt.0.) weta_in(ks)=0.
if (wetb_in(ks).lt.0.) wetb_in(ks)=0.
......@@ -354,9 +354,9 @@ subroutine wetdepo(itime,ltsample,loutnext)
! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
!OLD
if (readclouds) then
wetscav=S_i*(prec(1)/3.6E6)
wetscav=S_i*(prec(1)/3.6E6)
else
wetscav=S_i*(prec(1)/3.6E6)/clouds_h
wetscav=S_i*(prec(1)/3.6E6)/clouds_h
endif
!ZHG 2015 TEST
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment