Skip to content

GitLab

Commit 8100a04a authored by Sabine's avatar Sabine
Browse files

Bug: Release in domain filling corrected

parent a89b0606
Hide whitespace changes
Inline Side-by-side
src/boundcond_domainfill.f90
View file @ 8100a04a
......@@ -309,6 +309,8 @@ subroutine boundcond_domainfill(itime,loutend)
itrasplit(ipart)=itra1(ipart)+ldirect*itsplit
xmass1(ipart,1)=xmassperparticle
if (WATERCYCLE) then
status_q(ipart)=-1
! write(*,*) 'calculate watercycle in domainfilling',numactiveparticles,ipart
call calculate_watercycle(ipart,itime)
endif
if (mdomainfill.eq.2) xmass1(ipart,1)= &
......@@ -543,6 +545,11 @@ subroutine boundcond_domainfill(itime,loutend)
itramem(ipart)=itra1(ipart)
itrasplit(ipart)=itra1(ipart)+ldirect*itsplit
xmass1(ipart,1)=xmassperparticle
if (WATERCYCLE) then
status_q(ipart)=-1
! write(*,*) 'calculate watercycle in domainfilling 1',numactiveparticles,ipart
call calculate_watercycle(ipart,itime)
endif
if (mdomainfill.eq.2) xmass1(ipart,1)= &
xmass1(ipart,1)*pvpart*48./29.*ozonescale/10.**9
else
......
src/calculate_watercycle.f90
View file @ 8100a04a
......@@ -37,11 +37,15 @@ subroutine calculate_watercycle(partnumber,itime)
integer :: itime,i,j,k,jjjjmmdd,ihmmss,numshortout,numshortall
integer :: ix,jy,ixp,jyp,ind,indz,indzp,il,m,indexh,itage,nage
integer :: interp_time, ngrid
integer :: partnumber,forparticle ! -1 inititalize q, partnumber real do deposition
integer :: partnumber,forparticle,totalparticle ! -1 inititalize q, partnumber real do deposition
real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn
real centerposx_real, centerposy_real, centerposx, centerposy, diff, xl(2), yl(2)
real qv1(2),qvprof(2),qvi,dz,dz1,dz2
real e_minus_p1(2), e_minus_pi
real precip_reference(360,180),eminusp_saved
integer ixr,jyr,idumx,idumy
! Some variables needed for temporal interpolation
......@@ -51,11 +55,6 @@ subroutine calculate_watercycle(partnumber,itime)
dt2=real(memtime(2)-itime)
dtt=1./(dt1+dt2)
indexh=memind(2)
if (abs(memtime(1)-itime).lt.abs(memtime(2)-itime)) &
indexh=memind(1)
! if called from timemanager a loop over all particle is necessary
! if particle is released, the initialization has just to be done
! for the single particle
......@@ -63,7 +62,21 @@ subroutine calculate_watercycle(partnumber,itime)
if (partnumber.eq.-1) then
forparticle=numpart
write (*,*) 'calculate watercycle, ',numpart,partnumber,forparticle,itime
totalparticle=0
! open(43,file='precip_ref_input.dat',status='old')
! write(*,*) 'Reading eminusp from file'
! do 75 jyr=1,90
! do 75 ixr=1,360
!75 read(43,'(2i4,e15.3)') idumx,idumy,precip_reference(ixr,jyr)
! close(43)
! open(43,file='precip_ref_plus_eminusp.dat')
! write(*,*) 'Writing eminusp from file'
! do 76 jyr=1,90
! do 76 ixr=1,360
!6 write(43,'(2i4,3e15.3)') ixr,jyr,precip_reference(ixr,jyr), &
! e_minus_p(ixr,jyr,1), e_minus_p(ixr,jyr+90,1)
! close(43)
! write(*,*) 'Written eminusp from file'
else
forparticle=1
endif
......@@ -77,6 +90,7 @@ subroutine calculate_watercycle(partnumber,itime)
else
i=partnumber ! just for one particle
endif
waterfieldp(ix+1,jy+1)=waterfieldp(ix+1,jy+1)+xmass1(i,1)
! Take only valid particles, itime invalid when particle released!
!*****************************************************************
......@@ -91,6 +105,9 @@ subroutine calculate_watercycle(partnumber,itime)
end do
33 continue
if (((mod(itage,loutstep).eq.0).and.(status_q(i).ne.-9)).or.(status_q(i).eq.-1)) then
! at a loutstep timestep, and valid part, or first timestep after release
! Determine which nesting level to be used
!*****************************************
ngrid=0
......@@ -132,8 +149,8 @@ subroutine calculate_watercycle(partnumber,itime)
p3=rddx*ddy
p4=ddx*ddy
! Interpolate specific humidity
!******************************
! Interpolate specific humidity - code from partoutput.f90
!*********************************************************
do il=2,nz
if (height(il).gt.ztra1(i)) then
......@@ -186,18 +203,15 @@ subroutine calculate_watercycle(partnumber,itime)
endif
end do
e_minus_pi=(e_minus_p1(1)*dt2+dt1*e_minus_p1(2))*dtt
! e_minus_pi=e_minus_p1(memtime(1)) !test how the results change without interpolation - does not change a lot
eminusp_saved=e_minus_pi
! write(*,*) 'all: ',i,diff,e_minus_pi &
! write(*,*) 'all: ',i,diff,e_minus_pi &
! ,partnumber,val_q(i),qvi,itime,itra1(i),itage
if (itage.gt.0) then ! not the initialisation stage, it is assumed a value is already saved
! Calculate humidity differences
!*******************************
diff=(qvi-val_q(i))*xmass1(i,1)*ldirect
diff=(qvi-val_q(i))*xmass1(i,1)*ldirect
! Determine center of mass position on earth and average height
!**************************************************************
......@@ -211,43 +225,52 @@ subroutine calculate_watercycle(partnumber,itime)
! there has to be a deltaq which is negative and the column has to have precipitation
! e_minus_p units is mm/day .. convert to mm/h
! we cannot be sure 3600 is the second timestep, could also be 1800!
! if (((e_minus_pi)/24.gt.-2.0) .and. (itage.le.3600)) then
if (((diff.ge.0).or.((e_minus_pi/24).gt.-2.0)) .and. (itage.le.abs(loutstep))) then
! if (((diff.ge.0).or.((e_minus_pi).gt.-2.0)) .and. (itage.le.3600)) then
! if ((diff.ge.-1000000) .and. (itage.le.3600)) then
! write(*,*) 'terminated: ',i,diff,e_minus_pi &
! ,partnumber,val_q(i),qvi,itime,itra1(i),itage
itra1(i)=-999999999
else
! write(*,*) 'accounted: ',i,diff,e_minus_pi &
! ,partnumber,val_q(i),qvi,itime,itra1(i),itage,xtra1(i),ytra1(i),xmass1(i,1)
call centerofmass(xl,yl,2,centerposx_real,centerposy_real)
centerposx=(centerposx_real-xlon0)/dx
centerposy=(centerposy_real-ylat0)/dy
if (diff.gt.0) then
call drydepokernel(1,diff,centerposx,centerposy,nage,1)
if (nested_output.eq.1) &
call drydepokernel_nest(1,diff,centerposx,centerposy,nage,1)
else
call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
if (nested_output.eq.1) &
call wetdepokernel_nest(1,abs(diff),centerposx,centerposy,nage,1)
endif
endif
endif
! save the old values
!********************
xtra1_q(i)=xtra1(i)
ytra1_q(i)=ytra1(i)
val_q(i)=qvi
endif ! for valid trajectories
! e_minus_pi=precip_reference(ix+1,jy-90+2)*(-24.)
! First timestep and criteria fullfilled, keep it or throw it?
! if ( ((diff.ge.0).or.((e_minus_pi/24).gt.-2.0)) .and. (itage.eq.abs(loutstep)) ) then
! if ( (diff.ge.0) .and. (itage.eq.abs(loutstep)) ) then
! if ( (diff.ge.0) .and. (itage.eq.abs(loutstep)) ) then
! write(*,*) 'terminated :',i,diff,e_minus_pi,itage,status_q(i)
! status_q(i)=-9 ! this particle is invalid
! itra1(i)=-999999999
! else
if ((status_q(i).ne.-1).and.(itage.ge.abs(loutstep))) then
! write(*,*) 'accounted :',i,diff,e_minus_pi,itage,status_q(i)
totalparticle=totalparticle+1
call centerofmass(xl,yl,2,centerposx_real,centerposy_real)
centerposx=(centerposx_real-xlon0)/dx
centerposy=(centerposy_real-ylat0)/dy
if (diff.gt.0) then
call drydepokernel(1,diff,centerposx,centerposy,nage,1)
if (nested_output.eq.1) &
call drydepokernel_nest(1,diff,centerposx,centerposy,nage,1)
waterfielde(ix+1,jy+1)=waterfielde(ix+1,jy+1)+diff
else
call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
if (nested_output.eq.1) &
call wetdepokernel_nest(1,abs(diff),centerposx,centerposy,nage,1)
endif
! waterfieldp(ix+1,jy+1)=waterfieldp(ix+1,jy+1)+xmass1(i,1)
endif
! save the old values
! from timemanager ??????????????? if (mod(itime-loutstart,loutsample).eq.0) then
xtra1_q(i)=xtra1(i)
ytra1_q(i)=ytra1(i)
val_q(i)=qvi
! write(*,*) 'saved :',i,diff,e_minus_pi,itage,loutstep,status_q(i)
status_q(i)=1
! endif
endif ! at louttimestep interval
endif ! for valid trajectories - itra1.eq.itime
end do
if (partnumber.eq.-1) then
write (*,*) 'calculated watercycle, ',totalparticle,numpart,partnumber,forparticle,itime
endif
end subroutine calculate_watercycle
src/com_mod.f90
View file @ 8100a04a
......@@ -341,6 +341,9 @@ module com_mod
real :: lsm(0:nxmax-1,0:nymax-1)
real :: xlanduse(0:nxmax-1,0:nymax-1,numclass)
real :: waterfieldp(360,180),waterfielde(360,180)
! oro [m] orography of the ECMWF model
! excessoro excess orography mother domain
! lsm land sea mask of the ECMWF model
......@@ -681,7 +684,8 @@ module com_mod
real, allocatable, dimension(:,:,:) :: e_minus_p ! E-P field for the watercycle
real, allocatable, dimension(:,:,:,:) :: e_minus_p_nest ! E-P field for the watercycle
real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q
real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q
integer, allocatable, dimension(:) :: status_q
! eso: Moved from timemanager
real, allocatable, dimension(:) :: uap,ucp,uzp,us,vs,ws
......
src/init_domainfill.f90
View file @ 8100a04a
......@@ -278,6 +278,7 @@ subroutine init_domainfill
itsplit
xmass1(numpart+jj,1)=colmass(ix,jy)/real(ncolumn)
if (WATERCYCLE) then
status_q(numpart+jj)=-1
call calculate_watercycle(numpart+jj,0)
endif
......
src/readcommand.f90
View file @ 8100a04a
......@@ -353,6 +353,7 @@ subroutine readcommand
allocate(xtra1_q(maxpart))
allocate(ytra1_q(maxpart))
allocate(val_q(maxpart))
allocate(status_q(maxpart))
WATERCYCLE=.true.
mdomainfill=1
end select
......@@ -568,6 +569,13 @@ subroutine readcommand
stop
endif
if ((loutsample.gt.3600).and.(WATERCYCLE)) then
write(*,*) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### '
write(*,*) ' #### FOR WATERCYCLE RUNS HAS TO BE MAX 3600 #### '
write(*,*) ' #### CHANGE INPUT IN FILE COMMAND. #### '
stop
endif
if (loutsample.gt.loutaver) then
write(*,*) ' #### FLEXPART MODEL ERROR! SAMPLING TIME OF #### '
write(*,*) ' #### CONCENTRATION FIELD OUTPUT MUST NOT BE #### '
......
src/releaseparticles.f90
View file @ 8100a04a
......@@ -70,6 +70,7 @@ subroutine releaseparticles(itime)
write (*,*) 'releasepart: '
! Determine the actual date and time in Greenwich (i.e., UTC + correction for daylight savings time)
!*****************************************************************************
......
src/timemanager.f90
View file @ 8100a04a
......@@ -239,32 +239,45 @@ subroutine timemanager(metdata_format)
write (*,*) 'timemanager> Release particles'
endif
if (mdomainfill.ge.1) then
if (itime.eq.0) then
if (verbosity.gt.0) then
write (*,*) 'timemanager> call init_domainfill'
endif
call init_domainfill
else
if (abs(itime).le.abs(ireleasestart(1))) then
if (verbosity.gt.0) then
write (*,*) 'timemanager> call boundcond_domainfill'
endif
call boundcond_domainfill(itime,loutend)
if (WATERCYCLE) then ! works only backward for the moment
if (itime.eq.ireleaseend(1)) then !this has to be changed to +- timeintervall
if (verbosity.gt.0) then
write (*,*) 'timemanager> call boundcond_domainfill',itime,ireleasestart(1)
endif
call init_domainfill
endif
if ( (abs(itime).le.abs(ireleasestart(1))).and. &
(abs(itime).gt.abs( ireleaseend(1))) ) then
if (verbosity.gt.0) then
write (*,*) 'timemanager> call boundcond_domainfill'
endif
call boundcond_domainfill(itime,loutend)
endif
endif
else
if (verbosity.gt.0) then
print*,'call releaseparticles'
endif
call releaseparticles(itime)
if (verbosity.gt.1) then
CALL SYSTEM_CLOCK(count_clock)
WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
endif
if (mdomainfill.ge.1) then
if (itime.eq.0) then
if (verbosity.gt.0) then
write (*,*) 'timemanager> call init_domainfill'
endif
call init_domainfill
else
if (verbosity.gt.0) then
write (*,*) 'timemanager> call boundcond_domainfill'
endif
call boundcond_domainfill(itime,loutend)
endif
else
if (verbosity.gt.0) then
print*,'call releaseparticles'
endif
call releaseparticles(itime)
if (verbosity.gt.1) then
CALL SYSTEM_CLOCK(count_clock)
WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
endif
endif
endif
! Compute convective mixing for forward runs
! for backward runs it is done before next windfield is read in
!**************************************************************
......@@ -360,10 +373,11 @@ subroutine timemanager(metdata_format)
! Check whether concentrations are to be calculated
!**************************************************
if ((ldirect*itime.ge.ldirect*loutstart).and. &
(ldirect*itime.le.ldirect*loutend)) then ! add to grid
if (mod(itime-loutstart,loutsample).eq.0) then
if (WATERCYCLE) &
call calculate_watercycle(-1,itime)
! If we are exactly at the start or end of the concentration averaging interval,
! give only half the weight to this sample
......@@ -376,6 +390,7 @@ subroutine timemanager(metdata_format)
endif
outnum=outnum+weight
call conccalc(itime,weight)
endif
......@@ -447,8 +462,8 @@ subroutine timemanager(metdata_format)
46 format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
! if (ipout.ge.1) call partoutput(itime) ! dump particle positions
if (WATERCYCLE) &
call calculate_watercycle(-1,itime)
! if (WATERCYCLE) &
! call calculate_watercycle(-1,itime)
if (ipout.ge.1) call partoutput_short(itime) ! dump particle positions
loutnext=loutnext+loutstep
......@@ -745,6 +760,13 @@ subroutine timemanager(metdata_format)
if (linit_cond.ge.1) call initial_cond_output(itime) ! dump initial cond. field
if (WATERCYCLE) then
open(43,file='waterfield.dat')
write(*,*) 'Writing waterfield file'
do 76 jy=1,180
do 76 ix=1,360
76 write(43,'(2i4,2e15.3)') ix,jy,waterfielde(ix,jy),waterfieldp(ix,jy)
close(43)
close(89) ! Euler rain, water budget file
close(90) ! Euler rain, water budget file
endif
......
Supports Markdown
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