Skip to content
GitLab
Commit d6a09777 authored by Espen Sollum's avatar Espen Sollum
Browse files

Updates to Henrik's wet depo scheme

parent 88d8c3d8
Hide whitespace changes
Inline Side-by-side
src/FLEXPART_MPI.f90
View file @ d6a09777
......@@ -164,10 +164,11 @@ program flexpart
if (ldirect.le.0) then
write(*,FMT='(80("#"))')
write(*,*) '#### FLEXPART_MPI> ERROR: ', &
& 'MPI version not (yet) working with backward runs. ',&
& 'Use the serial version instead.'
& 'MPI version not (yet) working with backward runs. '
write(*,*) '#### Use the serial version instead.'
write(*,FMT='(80("#"))')
stop
! call mpif_finalize
! stop
end if
......
src/com_mod.f90
View file @ d6a09777
......@@ -129,7 +129,7 @@ module com_mod
! gdomainfill .T., if domain-filling is global, .F. if not
!hg
!ZHG SEP 2015 wheather or not to read clouds from GRIB
logical :: readclouds
!NIK 16.02.2015
......@@ -346,9 +346,11 @@ module com_mod
real :: ww(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: tt(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: qv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
!hg adding cloud water
real :: clwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
!ZHG adding cloud water
real :: clwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem) !liquid [kg/kg]
real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem) !ice [kg/kg]
real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem) !combined [m3/m3]
real :: pv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: rho(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: drhodz(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
......@@ -361,8 +363,10 @@ module com_mod
!scavenging NIK, PS
integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
integer :: cloudsh(0:nxmax-1,0:nymax-1,numwfmem)
! integer :: icloudbot(0:nxmax-1,0:nymax-1,numwfmem)
! integer :: icloudthck(0:nxmax-1,0:nymax-1,numwfmem)
!ZHG Sep 2015
real :: icloud_stats(0:nxmax-1,0:nymax-1,5,numwfmem)
real :: clh(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
! uu,vv,ww [m/2] wind components in x,y and z direction
......@@ -404,8 +408,8 @@ module com_mod
real :: hmix(0:nxmax-1,0:nymax-1,1,numwfmem)
real :: tropopause(0:nxmax-1,0:nymax-1,1,numwfmem)
real :: oli(0:nxmax-1,0:nymax-1,1,numwfmem)
! real :: diffk(0:nxmax-1,0:nymax-1,1,numwfmem) this is not in use?
! real :: diffk(0:nxmax-1,0:nymax-1,1,numwfmem) ESO: this is not in use?
! logical :: beneath_cloud=.true.
! ps surface pressure
! sd snow depth
! msl mean sea level pressure
......@@ -666,6 +670,10 @@ module com_mod
integer(kind=2), allocatable, dimension(:) :: cbt
!CGZ-lifetime
real, allocatable, dimension(:,:) ::checklifetime, species_lifetime
!CGZ-lifetime
! numpart actual number of particles in memory
! itra1 (maxpart) [s] temporal positions of the particles
! npoint(maxpart) indicates the release point of each particle
......@@ -758,7 +766,9 @@ module com_mod
allocate(itra1(nmpart),npoint(nmpart),nclass(nmpart),&
& idt(nmpart),itramem(nmpart),itrasplit(nmpart),&
& xtra1(nmpart),ytra1(nmpart),ztra1(nmpart),&
& xmass1(nmpart, maxspec))
& xmass1(nmpart, maxspec),&
& checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
allocate(uap(nmpart),ucp(nmpart),uzp(nmpart),us(nmpart),&
& vs(nmpart),ws(nmpart),cbt(nmpart))
......
src/gridcheck.f90
View file @ d6a09777
......@@ -194,12 +194,12 @@ subroutine gridcheck
isec1(6)=132 ! indicatorOfParameter
elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
isec1(6)=133 ! indicatorOfParameter
!hg
!ZHG FOR CLOUDS FROM GRIB
elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
isec1(6)=246 ! indicatorOfParameter
elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
isec1(6)=247 ! indicatorOfParameter
!hg end
!ZHG end
elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
isec1(6)=134 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
......
src/mpi_mod.f90
View file @ d6a09777
......@@ -88,7 +88,7 @@ module mpi_mod
! MPI tags/requests for send/receive operation
integer :: tm1
integer, parameter :: nvar_async=27 !29 :DBG:
integer, parameter :: nvar_async=26 !27 !29 :DBG:
!integer, dimension(:), allocatable :: tags
integer, dimension(:), allocatable :: reqs
......@@ -118,7 +118,7 @@ module mpi_mod
logical, parameter :: mp_dev_mode = .false.
logical, parameter :: mp_dbg_out = .false.
logical, parameter :: mp_time_barrier=.true.
logical, parameter :: mp_measure_time=.true.
logical, parameter :: mp_measure_time=.false.
logical, parameter :: mp_exact_numpart=.true.
! for measuring CPU/Wall time
......@@ -206,7 +206,7 @@ contains
write(*,*) 'amount of memory is being allocated.'
write(*,FMT='(80("#"))')
! Force "syncronized" version if all processes will call getfields
else if (.not.lmp_sync.and.mp_np.lt.read_grp_min) then
else if ((.not.lmp_sync.and.mp_np.lt.read_grp_min).or.(mp_np.eq.1)) then
if (lroot) then
write(*,FMT='(80("#"))')
write(*,*) '#### mpi_mod::mpif_init> WARNING: ', &
......@@ -962,10 +962,13 @@ contains
! cloud water/ice:
if (readclouds) then
call MPI_Bcast(clwc(:,:,:,li:ui),d3s1,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
if (mp_ierr /= 0) goto 600
call MPI_Bcast(ciwc(:,:,:,li:ui),d3s1,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2_size1*5,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
if (mp_ierr /= 0) goto 600
! call MPI_Bcast(clwc(:,:,:,li:ui),d3s1,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
! if (mp_ierr /= 0) goto 600
! call MPI_Bcast(ciwc(:,:,:,li:ui),d3s1,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
! if (mp_ierr /= 0) goto 600
end if
! 2D fields
......@@ -1188,7 +1191,6 @@ contains
! mind -- index where to place new fields
!
! TODO
! Transfer cloud water/ice
!
!*******************************************************************************
use com_mod
......@@ -1333,18 +1335,20 @@ contains
! Send cloud water if it exists. Increment counter always (as on receiving end)
if (readclouds) then
i=i+1
call MPI_Isend(clwc(:,:,:,mind),d3s1,mp_pp,dest,tm1,&
call MPI_Isend(icloud_stats(:,:,:,mind),d2s1*5,mp_pp,dest,tm1,&
&MPI_COMM_WORLD,reqs(i),mp_ierr)
if (mp_ierr /= 0) goto 600
i=i+1
call MPI_Isend(ciwc(:,:,:,mind),d3s1,mp_pp,dest,tm1,&
&MPI_COMM_WORLD,reqs(i),mp_ierr)
if (mp_ierr /= 0) goto 600
! else
! i=i+2
end if
! call MPI_Isend(clwc(:,:,:,mind),d3s1,mp_pp,dest,tm1,&
! &MPI_COMM_WORLD,reqs(i),mp_ierr)
! if (mp_ierr /= 0) goto 600
! i=i+1
! call MPI_Isend(ciwc(:,:,:,mind),d3s1,mp_pp,dest,tm1,&
! &MPI_COMM_WORLD,reqs(i),mp_ierr)
! if (mp_ierr /= 0) goto 600
end if
end do
if (mp_measure_time) call mpif_mtime('commtime',1)
......@@ -1370,7 +1374,6 @@ contains
! memstat -- input, used to resolve windfield index being received
!
! TODO
! Transfer cloud water/ice
!
!*******************************************************************************
use com_mod
......@@ -1390,12 +1393,6 @@ contains
! integer :: d3s1,d3s2,d2s1,d2s2
!*******************************************************************************
! :TODO: don't need these
! d3s1=d3_size1
! d3s2=d3_size2
! d2s1=d2_size1
! d2s2=d2_size2
! At the time this immediate receive is posted, memstat is the state of
! windfield indices at the previous/current time. From this, the future
! state is deduced.
......@@ -1544,13 +1541,19 @@ contains
! For now assume that data at all steps either have or do not have water
if (readclouds) then
j=j+1
call MPI_Irecv(clwc(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
&MPI_COMM_WORLD,reqs(j),mp_ierr)
if (mp_ierr /= 0) goto 600
j=j+1
call MPI_Irecv(ciwc(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
&MPI_COMM_WORLD,reqs(j),mp_ierr)
if (mp_ierr /= 0) goto 600
call MPI_Irecv(icloud_stats(:,:,:,mind),d2s1*5,mp_pp,id_read,MPI_ANY_TAG,&
&MPI_COMM_WORLD,reqs(j),mp_ierr)
if (mp_ierr /= 0) goto 600
! call MPI_Irecv(clwc(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
! &MPI_COMM_WORLD,reqs(j),mp_ierr)
! if (mp_ierr /= 0) goto 600
! j=j+1
! call MPI_Irecv(ciwc(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
! &MPI_COMM_WORLD,reqs(j),mp_ierr)
! if (mp_ierr /= 0) goto 600
end if
......@@ -2088,7 +2091,9 @@ contains
implicit none
integer, parameter :: li=1, ui=3 ! wfmem indices (i.e, operate on entire field)
integer :: li=1, ui=2 ! wfmem indices (i.e, operate on entire field)
if (.not.lmp_sync) ui=3
! Variables transferred at initialization only
......
src/par_mod.f90
View file @ d6a09777
......@@ -75,7 +75,7 @@ module par_mod
real,parameter :: karman=0.40, href=15., convke=2.0
real,parameter :: hmixmin=100., hmixmax=4500., turbmesoscale=0.16
real,parameter :: d_trop=50., d_strat=0.1
real,parameter :: rho_water=1000 !ZHG 2015 [kg/m3]
! karman Karman's constant
! href [m] Reference height for dry deposition
! konvke Relative share of kinetic energy used for parcel lifting
......@@ -225,8 +225,8 @@ module par_mod
! Sabine Eckhardt: change of landuse inventary numclass=13
! ---------
integer,parameter :: maxwf=50000, maxtable=1000, numclass=13, ni=11
!integer,parameter :: numwfmem=2 ! Serial version/MPI with 2 fields
integer,parameter :: numwfmem=3 ! MPI with 3 fields
integer,parameter :: numwfmem=2 ! Serial version/MPI with 2 fields
!integer,parameter :: numwfmem=3 ! MPI with 3 fields
! maxwf maximum number of wind fields to be used for simulation
! maxtable Maximum number of chemical species that can be
......
src/pathnames
View file @ d6a09777
../options/
../output/
/
/xnilu_wrk/flex_wrk/WIND_FIELDS/AVAILABLE_ECMWF_OPER_fields_global/
/xnilu_wrk/flex_wrk/WIND_FIELDS/AVAILABLE_ECMWF_OPER_fields_global
============================================
/xnilu_wrk/flex_wrk/zhg/ECMWF_CLWC/Availables
src/readwind.f90
View file @ d6a09777
......@@ -105,9 +105,9 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
hflswitch=.false.
strswitch=.false.
!hg
!ZHG test the grib fields that have lcwc without using them
! readcloud=.false.
!hg end
levdiff2=nlev_ec-nwz+1
iumax=0
iwmax=0
......@@ -189,12 +189,13 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
isec1(6)=132 ! indicatorOfParameter
elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
isec1(6)=133 ! indicatorOfParameter
!hg
!ZHG ! READ CLOUD FIELD : assume these to be added together to one variable
elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
isec1(6)=246 ! indicatorOfParameter
elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
isec1(6)=247 ! indicatorOfParameter
!hg end
! ICE AND WATER IS ADDED TOGETHER IN NEW WINDFIELDS
! elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
! isec1(6)=247 ! indicatorOfParameter
!ZHG end
elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
isec1(6)=134 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
......@@ -225,9 +226,11 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
isec1(6)=146 ! indicatorOfParameter
elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
isec1(6)=176 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
! elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS --wrong
elseif ((parCat.eq.2).and.(parNum.eq.38) .or. parId .eq. 180) then ! EWSS --correct
isec1(6)=180 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
! elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS --wrong
elseif ((parCat.eq.2).and.(parNum.eq.37) .or. parId .eq. 181) then ! NSSS --correct
isec1(6)=181 ! indicatorOfParameter
elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
isec1(6)=129 ! indicatorOfParameter
......@@ -355,17 +358,17 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
zsec4(nxfield*(ny-j-1)+i+1)
if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
!hg READING CLOUD FIELDS ASWELL
if(isec1(6).eq.246) then !! CLWC Cloud liquid water content
!ZHG READING CLOUD FIELDS ASWELL
if(isec1(6).eq.246) then !! CLWC Cloud liquid water content [kg/kg]
clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
readclouds = .true.
!write(*,*) 'found water!'
!if (clwch(i,j,nlev_ec-k+2,n) .gt. 0) write(*,*) 'readwind: found water!', clwch(i,j,nlev_ec-k+2,n)
endif
if(isec1(6).eq.247) then !! CIWC Cloud ice water content
ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
! if(isec1(6).eq.247) then !! CIWC Cloud ice water content
! ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
!write(*,*) 'found ice!'
endif
!hg end
! endif
!ZHG end
end do
end do
......@@ -421,10 +424,10 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
!hg
!ZHG
call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,2,n)
!hg end
!ZHG end
endif
......
src/readwind_mpi.f90
View file @ d6a09777
......@@ -199,12 +199,13 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
isec1(6)=132 ! indicatorOfParameter
elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
isec1(6)=133 ! indicatorOfParameter
!hg
!ZHG ! READ CLOUD FIELD : assume these to be added together to one variable
elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
isec1(6)=246 ! indicatorOfParameter
elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
isec1(6)=247 ! indicatorOfParameter
!hg end
! ICE AND WATER IS ADDED TOGETHER IN NEW WINDFIELDS
! elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
! isec1(6)=247 ! indicatorOfParameter
!ZHG end
elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
isec1(6)=134 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
......@@ -367,18 +368,17 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
zsec4(nxfield*(ny-j-1)+i+1)
if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
!hg READING CLOUD FIELDS ASWELL
if(isec1(6).eq.246) then !! CLWC Cloud liquid water content
!ZHG READING CLOUD FIELDS ASWELL
if(isec1(6).eq.246) then !! CLWC Cloud liquid water content [kg/kg]
clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
readclouds = .true.
!write(*,*) 'found water!'
!if (clwch(i,j,nlev_ec-k+2,n) .gt. 0) write(*,*) 'readwind: found water!', clwch(i,j,nlev_ec-k+2,n)
endif
if(isec1(6).eq.247) then !! CIWC Cloud ice water content
ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
! if(isec1(6).eq.247) then !! CIWC Cloud ice water content
! ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
!write(*,*) 'found ice!'
endif
!hg end
! endif
!ZHG end
end do
end do
......
src/timemanager.f90
View file @ d6a09777
......@@ -135,6 +135,14 @@ subroutine timemanager
! Loop over the whole modelling period in time steps of mintime seconds
!**********************************************************************
!ZHG 2015
!CGZ-lifetime: set lifetime to 0
checklifetime(:,:)=0
species_lifetime(:,:)=0
print*, 'Initialized lifetime'
!CGZ-lifetime: set lifetime to 0
if (verbosity.gt.0) then
write (*,*) 'timemanager> starting simulation'
......@@ -410,6 +418,15 @@ subroutine timemanager
if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
if (iflux.eq.1) call fluxoutput(itime)
write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
!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
!ZHG
!CGZ-lifetime: output species lifetime
!write(*,46) float(itime)/3600,itime,numpart
45 format(i9,' SECONDS SIMULATED: ',i8, ' PARTICLES: Uncertainty: ',3f7.3)
46 format(' Simulated ',f7.1,' hours (',i9,' s), ',i8, ' particles')
......@@ -574,6 +591,16 @@ subroutine timemanager
if (xmass(npoint(j),ks).gt.0.) &
xmassfract=max(xmassfract,real(npart(npoint(j)))* &
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
!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
!CGZ-lifetime: Check mass fraction left/save lifetime
!ZHG 2015
else
xmassfract=1.
endif
......
src/timemanager_mpi.f90
View file @ d6a09777
......@@ -138,6 +138,25 @@ subroutine timemanager
! Loop over the whole modelling period in time steps of mintime seconds
!**********************************************************************
! itime=0
if (lroot) then
! write(*,45) itime,numpart*mp_partgroup_np,gridtotalunc,wetgridtotalunc,drygridtotalunc
write(*,46) float(itime)/3600,itime,numpart*mp_partgroup_np
if (verbosity.gt.0) then
write (*,*) 'timemanager> starting simulation'
end if
end if ! (lroot)
!CGZ-lifetime: set lifetime to 0
checklifetime(:,:)=0
species_lifetime(:,:)=0
print*, 'Initialized lifetime'
!CGZ-lifetime: set lifetime to 0
do itime=0,ideltas,lsynctime
! Computation of wet deposition, OH reaction and mass transfer
......@@ -543,13 +562,17 @@ subroutine timemanager
& mp_comm_used, mp_ierr)
endif
if (lroot) then
!CGZ-lifetime: output species lifetime
! 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
write(*,45) itime,numpart_tot_mpi,gridtotalunc,&
&wetgridtotalunc,drygridtotalunc
if (verbosity.gt.0) then
write (*,*) 'timemanager> starting simulation'
end if
end if
! if (verbosity.gt.0) then
! write (*,*) 'timemanager> starting simulation'
! end if
! end if
45 format(i13,' SECONDS SIMULATED: ',i13, ' PARTICLES: Uncertainty: ',3f7.3)
46 format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
......@@ -728,15 +751,29 @@ subroutine timemanager
if (mdomainfill.eq.0) then
if (xmass(npoint(j),ks).gt.0.) &
if (xmass(npoint(j),ks).gt.0.)then
xmassfract=max(xmassfract,real(npart(npoint(j)))* &
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
!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
!CGZ-lifetime: Check mass fraction left/save lifetime
endif
else
xmassfract=1.
endif
end do
if (xmassfract.lt.0.0001) then ! terminate all particles carrying less mass
if (xmassfract.lt.0.00005 .and. sum(real(npart(npoint(j)))*xmass1(j,:)).lt.1.0) then ! terminate all particles carrying less mass
! print*,'terminated particle ',j,' for small mass (', sum(real(npart(npoint(j)))* &
! xmass1(j,:)), ' of ', sum(xmass(npoint(j),:)),')'
itra1(j)=-999999999
endif
......@@ -757,6 +794,7 @@ subroutine timemanager
if (linit_cond.ge.1) &
call initial_cond_calc(itime+lsynctime,j)
itra1(j)=-999999999
!print*, 'terminated particle ',j,'for age'
endif
endif
......@@ -827,7 +865,7 @@ subroutine timemanager
deallocate(drygridunc,wetgridunc)
endif
deallocate(gridunc)
deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass)
deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass, checklifetime)
deallocate(ireleasestart,ireleaseend,npart,kindz)
deallocate(xmasssave)
if (nested_output.eq.1) then
......
src/verttransform.f90
View file @ d6a09777
......@@ -96,13 +96,19 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
logical :: init = .true.
!hg
!ZHG SEP 2014 tests
integer :: cloud_ver,cloud_min, cloud_max
integer ::teller(5), convpteller=0, lspteller=0
real :: cloud_col_wat, cloud_water
!hg temporary variables for testing
!ZHG 2015 temporary variables for testing
real :: rcw(0:nxmax-1,0:nymax-1)
real :: rpc(0:nxmax-1,0:nymax-1)
!hg
character(len=60) :: zhgpath='/xnilu_wrk/flex_wrk/zhg/'
character(len=60) :: fnameA,fnameB,fnameC,fnameD,fnameE,fnameF,fnameG,fnameH
CHARACTER(LEN=3) :: aspec
integer :: virr=0
real :: tot_cloud_h
!ZHG
!*************************************************************************
! If verttransform is called the first time, initialize heights of the *
......@@ -560,12 +566,75 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
endif
!***********************************************************************************
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
! cloud scavenging are assigned with numbers 2-5 (following the old metod).
! Distinction is done for lsp and convp though they are treated the same in regards
! 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'
! if (.not.readclouds) write(*,*) 'using cloud water from parameterization'
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)