calculate_watercycle.f90

!**********************************************************************
! 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 <http://www.gnu.org/licenses/>.   *
!**********************************************************************

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