Commit 7d7a1d71 authored by Sabine's avatar Sabine
Browse files

introduced the timestep variable for the watercycle

parent e7452f4a
......@@ -309,6 +309,7 @@ subroutine boundcond_domainfill(itime,loutend)
itrasplit(ipart)=itra1(ipart)+ldirect*itsplit
xmass1(ipart,1)=xmassperparticle
if (WATERCYCLE) then
firsttimestep(ipart)=1
call calculate_watercycle(ipart,itime)
endif
if (mdomainfill.eq.2) xmass1(ipart,1)= &
......
......@@ -38,8 +38,8 @@ subroutine calculate_watercycle(partnumber,itime)
integer :: ix,jy,ixp,jyp,ind,indz,indzp,il,m,indexh,itage,nage
integer :: interp_time
integer :: partnumber,forparticle ! -1 inititalize q, partnumber real do deposition
real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,topo
real centerposx_real, centerposy_real, centerposx, centerposy, diff, xl(2), yl(2), prec
real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4
real centerposx_real, centerposy_real, centerposx, centerposy, diff, xl(2), yl(2)
real qv1(2),qvprof(2),qvi,dz,dz1,dz2
......@@ -75,13 +75,12 @@ subroutine calculate_watercycle(partnumber,itime)
i=j ! loop over all particle
else
i=partnumber ! just for one particle
firsttimestep(i)=.true.
endif
! Take only valid particles, itime invalid when particle released!
!*****************************************************************
if ((itra1(i).eq.itime).or.(forparticle.eq.1)) then
if (itra1(i).eq.itime) then ! trajectory itra1 is not -999999
!*****************************************************************************
! Interpolate several variables (PV, specific humidity, etc.) to particle position
......@@ -109,14 +108,6 @@ subroutine calculate_watercycle(partnumber,itime)
p3=rddx*ddy
p4=ddx*ddy
! Topography
!***********
topo=p1*oro(ix ,jy) &
+ p2*oro(ixp,jy) &
+ p3*oro(ix ,jyp) &
+ p4*oro(ixp,jyp)
! Potential vorticity, specific humidity, temperature, and density
!*****************************************************************
......@@ -150,49 +141,33 @@ subroutine calculate_watercycle(partnumber,itime)
end do
qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz
if (firsttimestep(i).gt.1) then ! not the initialisation stage, it is assumed a value is already saved
! Calculate humidity differences
!*******************************
prec=p1*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh)) &
+p2*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh)) &
+p3*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh)) &
+p4*(lsprec(ix ,jy ,1,indexh)+convprec(ix ,jy ,1,indexh))
if (partnumber.eq.-1) then ! not the initialisation stage
! 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
!**************************************************************
xl(1)=xlon0+xtra1_q(i)*dx
xl(2)=xlon0+xtra1(i)*dx
yl(1)=ylat0+ytra1_q(i)*dy
yl(2)=ylat0+ytra1(i)*dy
endif
xl(1)=xlon0+xtra1_q(i)*dx
xl(2)=xlon0+xtra1(i)*dx
yl(1)=ylat0+ytra1_q(i)*dy
yl(2)=ylat0+ytra1(i)*dy
! save the old values - for particle after release - particlenumber = -1
!********************
xtra1_q(i)=xtra1(i)
ytra1_q(i)=ytra1(i)
val_q(i)=qvi
if (partnumber.eq.-1) then ! not the initialisation stage
! Do this only for the 2nd timestep i
! Do this only for the 2nd timestep
! 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
if ( ((diff.ge.0).or.(e_minus_p(ix,jy)/24.gt.-2.0)) .and. firsttimestep(i)) then
! if ( ((diff.ge.0)) .and. firsttimestep(i)) then
if ( ((diff.ge.0).or.(e_minus_p(ix,jy).gt.-2.0)) .and. (firsttimestep(i).eq.2)) then
itra1(i)=-999999999
! write(*,*) 'Trajectory terminated: ',i,diff,e_minus_p(ix,jy),firsttimestep(i)
firsttimestep(i)=1
! write(*,*) 'terminated: ',i,diff,e_minus_p(ix,jy),firsttimestep(i) &
! ,partnumber,val_q(i),qvi,itime,itra1(i),itage
else
! write(*,*) 'Trajectory accounted: ',i,diff,e_minus_p(ix,jy),firsttimestep(i),&
! val_q(i),qvi,xmass1(i,1)
prec_q(i)=prec
! write(*,*) 'accounted: ',i,diff,e_minus_p(ix,jy),firsttimestep(i) &
! ,partnumber,val_q(i),qvi,itime,itra1(i),itage
call centerofmass(xl,yl,2,centerposx_real,centerposy_real)
centerposx=(centerposx_real-xlon0)/dx
centerposy=(centerposy_real-ylat0)/dy
......@@ -202,10 +177,20 @@ subroutine calculate_watercycle(partnumber,itime)
call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
endif
endif
firsttimestep(i)=.false.
endif
endif ! firsttimestep gt 1
endif
! save the old values
!********************
xtra1_q(i)=xtra1(i)
ytra1_q(i)=ytra1(i)
val_q(i)=qvi
! write(*,*) 'saved: ',i,xtra1_q(i),ytra1_q(i),val_q(i),firsttimestep(i)
firsttimestep(i)=firsttimestep(i)+1 ! this is a valid particle keep in the simulation
endif ! for valid trajectories
end do
......
......@@ -679,8 +679,8 @@ module com_mod
real, allocatable, dimension(:,:) :: xscav_frac1
real, allocatable, dimension(:,:) :: e_minus_p ! E-P field for the watercycle
real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q, prec_q
logical, allocatable, dimension(:) :: firsttimestep ! to save the previous q value and position for the WATERCYCLE mode
real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q
integer, allocatable, dimension(:) :: firsttimestep ! to save the previous q value and position for the WATERCYCLE mode
! eso: Moved from timemanager
real, allocatable, dimension(:) :: uap,ucp,uzp,us,vs,ws
......
......@@ -278,6 +278,7 @@ subroutine init_domainfill
itsplit
xmass1(numpart+jj,1)=colmass(ix,jy)/real(ncolumn)
if (WATERCYCLE) then
firsttimestep(numpart+jj)=1
call calculate_watercycle(numpart+jj,0)
endif
......
......@@ -353,7 +353,6 @@ subroutine readcommand
allocate(xtra1_q(maxpart))
allocate(ytra1_q(maxpart))
allocate(val_q(maxpart))
allocate(prec_q(maxpart))
allocate(firsttimestep(maxpart))
WATERCYCLE=.true.
mdomainfill=1
......
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