!********************************************************************** ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 * ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * ! * ! This file is part of FLEXPART. * ! * ! FLEXPART is free software: you can redistribute it and/or modify * ! it under the terms of the GNU General Public License as published by* ! the Free Software Foundation, either version 3 of the License, or * ! (at your option) any later version. * ! * ! FLEXPART is distributed in the hope that it will be useful, * ! but WITHOUT ANY WARRANTY; without even the implied warranty of * ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * ! GNU General Public License for more details. * ! * ! You should have received a copy of the GNU General Public License * ! along with FLEXPART. If not, see . * !********************************************************************** subroutine calculate_watercycle(partnumber,itime) !***************************************************************************** ! * ! calculate difference in q and dump it to the dry/wetdepofield * ! VAR: real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q * ! * !***************************************************************************** use par_mod use com_mod implicit none real(kind=dp) :: jul 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,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, threshold real precip_reference(360,180),eminusp_saved integer ixr,jyr,idumx,idumy,istep, watersynctime watersynctime=3600*3 ! was loutstep before, but for monthly run not possible ! watersynctime=3600 ! watersynctime=loutstep threshold = -.5 ! units mm/h if ((partnumber.eq.1).and.(itime.eq.0)) write(*,*) 'Watersynctime ',watersynctime,lsynctime, ' s, threshold: ',threshold, 'mm/h' ! Some variables needed for temporal interpolation !************************************************* dt1=real(itime-memtime(1)) dt2=real(memtime(2)-itime) dtt=1./(dt1+dt2) ! 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 !****************************************************************** if (partnumber.eq.-1) then forparticle=numpart totalparticle=0 istep=int(itime/watersynctime) ! write(*,*) 'Model step:',istep ! 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 ! Loop about all particles !************************* do j=1,forparticle if (partnumber.eq.-1) then i=j ! loop over all particle else i=partnumber ! just for one particle endif ! 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)) do nage=1,nageclass if (itage.lt.lage(nage)) goto 33 end do 33 continue ! Determine which nesting level to be used !***************************************** ngrid=0 do k=numbnests,1,-1 if ((xtra1(i).gt.xln(k)).and.(xtra1(i).lt.xrn(k)).and. & (ytra1(i).gt.yln(k)).and.(ytra1(i).lt.yrn(k))) then ngrid=k goto 23 endif end do 23 continue ! Determine nested grid coordinates !********************************** if (ngrid.gt.0) then xtn=(xtra1(i)-xln(ngrid))*xresoln(ngrid) ytn=(ytra1(i)-yln(ngrid))*yresoln(ngrid) ix=int(xtn) jy=int(ytn) ddy=ytn-real(jy) ddx=xtn-real(ix) else ix=int(xtra1(i)) jy=int(ytra1(i)) ddx=xtra1(i)-real(ix) ddy=ytra1(i)-real(jy) endif ixp=ix+1 jyp=jy+1 if (jy.eq.numygrid) jyp=numygrid rddx=1.-ddx rddy=1.-ddy p1=rddx*rddy p2=ddx*rddy p3=rddx*ddy p4=ddx*ddy ! Interpolate specific humidity - code from partoutput.f90 !********************************************************* do il=2,nz if (height(il).gt.ztra1(i)) then indz=il-1 indzp=il goto 6 endif end do 6 continue dz1=ztra1(i)-height(indz) dz2=height(indzp)-ztra1(i) dz=1./(dz1+dz2) do ind=indz,indzp do m=1,2 indexh=memind(m) if (ngrid.eq.0) then qv1(m)=p1*qv(ix ,jy ,ind,indexh) & +p2*qv(ixp,jy ,ind,indexh) & +p3*qv(ix ,jyp,ind,indexh) & +p4*qv(ixp,jyp,ind,indexh) else qv1(m)=p1*qvn(ix ,jy ,ind,indexh,ngrid) & +p2*qvn(ixp,jy ,ind,indexh,ngrid) & +p3*qvn(ix ,jyp,ind,indexh,ngrid) & +p4*qvn(ixp,jyp,ind,indexh,ngrid) endif end do qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt end do qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz if ( (mod(itage,watersynctime).eq.0).and.(itage.gt.0) ) then !trajectory is loutstep after simulation start, not a new particle do m=1,2 indexh=memind(m) if (ngrid.eq.0) then e_minus_p1(m)=p1*e_minus_p(ix ,jy ,indexh) & +p2*e_minus_p(ixp,jy ,indexh) & +p3*e_minus_p(ix ,jyp,indexh) & +p4*e_minus_p(ixp,jyp,indexh) else e_minus_p1(m)=p1*e_minus_p_nest(ix ,jy ,indexh,ngrid) & +p2*e_minus_p_nest(ixp,jy ,indexh,ngrid) & +p3*e_minus_p_nest(ix ,jyp,indexh,ngrid) & +p4*e_minus_p_nest(ixp,jyp,indexh,ngrid) endif end do e_minus_pi=(e_minus_p1(1)*dt2+dt1*e_minus_p1(2))*dtt eminusp_saved=e_minus_pi ! Calculate humidity differences !******************************* 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 ! 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 ! 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.threshold)) .and. (status_q(i).eq.1) ) then ! if ( ((diff.ge.0).or.((e_minus_pi/24).gt.-2.0)) .and. (status_q(i).eq.1) ) then ! if ( ((diff.ge.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) 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 ! 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).and.(totalparticle.gt.0)) then ! write (*,*) 'calculated watercycle, ',totalparticle,numpart,partnumber,forparticle,itime endif end subroutine calculate_watercycle