Commit 574a08de authored by Sabine's avatar Sabine
Browse files

this version works with 3hourly output, very similar to Stohl paper, a test...

this version works with 3hourly output, very similar to Stohl paper, a test field waterfield included to check particle distribution, has to be removed later
parent 95b08353
......@@ -44,7 +44,7 @@ subroutine calculate_watercycle(partnumber,itime)
real e_minus_p1(2), e_minus_pi
real precip_reference(360,180),eminusp_saved
integer ixr,jyr,idumx,idumy
integer ixr,jyr,idumx,idumy,istep
......@@ -63,6 +63,8 @@ subroutine calculate_watercycle(partnumber,itime)
if (partnumber.eq.-1) then
forparticle=numpart
totalparticle=0
istep=int(itime/loutstep)
write(*,*) 'Model step:',istep
! open(43,file='precip_ref_input.dat',status='old')
! write(*,*) 'Reading eminusp from file'
! do 75 jyr=1,90
......@@ -90,13 +92,12 @@ 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)
! write(*,*) 'particle: ',i,xtra1_q(i),xtra1(i),itime,itra1(i),itramem(i),status_q(i)
! Take only valid particles, itime invalid when particle released!
!*****************************************************************
if (itra1(i).eq.itime) then ! trajectory itra1 is not -999999
! Determine age class of the particle - nage is used for the kernel
!******************************************************************
itage=abs(itra1(i)-itramem(i))
......@@ -104,9 +105,7 @@ subroutine calculate_watercycle(partnumber,itime)
if (itage.lt.lage(nage)) goto 33
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
!*****************************************
......@@ -138,7 +137,6 @@ subroutine calculate_watercycle(partnumber,itime)
ddy=ytra1(i)-real(jy)
endif
ixp=ix+1
jyp=jy+1
if (jy.eq.numygrid) jyp=numygrid
......@@ -187,6 +185,8 @@ subroutine calculate_watercycle(partnumber,itime)
end do
qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz
if ( (mod(itage,loutstep).eq.0).and.(itage.gt.0) ) then !trajectory is loutstep after simulation start, not a new particle
do m=1,2
indexh=memind(m)
......@@ -205,13 +205,10 @@ subroutine calculate_watercycle(partnumber,itime)
e_minus_pi=(e_minus_p1(1)*dt2+dt1*e_minus_p1(2))*dtt
eminusp_saved=e_minus_pi
! write(*,*) 'all: ',i,diff,e_minus_pi &
! ,partnumber,val_q(i),qvi,itime,itra1(i),itage
! 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
!**************************************************************
......@@ -228,47 +225,52 @@ subroutine calculate_watercycle(partnumber,itime)
! 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)
if ( ((diff.ge.0).or.((e_minus_pi/24).gt.-2.0)) .and. (status_q(i).eq.1) ) then
! if (status_q(i).eq.-100) then ! never! all trajectories.
! write(*,*) 'terminated: ',i,diff,e_minus_pi,itage,itramem(i),status_q(i),qvi,val_q(i),xmass1(i,1)
status_q(i)=-9 ! this particle is invalid
itra1(i)=-999999999
else
totalparticle=totalparticle+1
call centerofmass(xl,yl,2,centerposx_real,centerposy_real)
centerposx=(centerposx_real-xlon0)/dx
centerposy=(centerposy_real-ylat0)/dy
! write(*,*) 'accounted: ',i,itime,itra1(i),itage,diff,itramem(i),status_q(i)
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
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
if (istep.lt.100) then
waterfieldp(istep,ix+1,jy+1)=waterfieldp(istep,ix+1,jy+1)+xmass1(i,1)
endif
waterfielde(ix+1,jy+1)=waterfielde(ix+1,jy+1)+xmass1(i,1)
xtra1_q(i)=xtra1(i)
ytra1_q(i)=ytra1(i)
val_q(i)=qvi
status_q(i)=2
endif
endif ! at louttimestep interval
if (status_q(i).eq.-1) then ! first time, a value is saved
! write(*,*) 'saved: ',i,xtra1(i),ytra1(i),qvi,itime,itage,itramem(i),status_q(i)
xtra1_q(i)=xtra1(i)
ytra1_q(i)=ytra1(i)
val_q(i)=qvi
status_q(i)=1
endif
endif ! for valid trajectories - itra1.eq.itime
end do
if (partnumber.eq.-1) then
if ((partnumber.eq.-1).and.(totalparticle.gt.0)) then
write (*,*) 'calculated watercycle, ',totalparticle,numpart,partnumber,forparticle,itime
endif
......
......@@ -342,7 +342,7 @@ module com_mod
real :: xlanduse(0:nxmax-1,0:nymax-1,numclass)
real :: waterfieldp(360,180),waterfielde(360,180)
real :: waterfieldp(100,360,180),waterfielde(360,180)
! oro [m] orography of the ECMWF model
! excessoro excess orography mother domain
......
......@@ -105,7 +105,7 @@ subroutine timemanager(metdata_format)
implicit none
integer :: metdata_format
integer :: metdata_format,iistep
integer :: j,ks,kp,l,n,itime=0,nstop,nstop1
! integer :: ksp
integer :: loutnext,loutstart,loutend
......@@ -761,9 +761,10 @@ subroutine timemanager(metdata_format)
if (WATERCYCLE) then
open(43,file='waterfield.dat')
write(*,*) 'Writing waterfield file'
do 76 iistep=1,100
do 76 jy=1,180
do 76 ix=1,360
76 write(43,'(2i4,2e15.3)') ix,jy,waterfielde(ix,jy),waterfieldp(ix,jy)
76 write(43,'(2i4,2e15.3)') ix,jy,waterfieldp(iistep,ix,jy),waterfielde(ix,jy)
close(43)
close(89) ! Euler rain, water budget file
......
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