source: branches/flexpart91_hasod/src_parallel/initial_cond_calc.f90 @ 8

!**********************************************************************
! 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 initial_cond_calc(itime,i)
  !                               i   i
  !*****************************************************************************
  !                                                                            *
  !     Calculation of the sensitivity to initial conditions for BW runs       *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     15 January 2010                                                        *
  !                                                                            *
  !*****************************************************************************

  use unc_mod, only: init_cond 
  use outg_mod, only: outheight
  use com_mod, only: linit_cond, dx, dy, xmass1, ztra1, xtra1, ytra1, npoint, &
      height, rho, itra1, xoutshift, youtshift, nz, numzgrid, numygrid, numxgrid, &
      mdomainfill, nspec, dxout, dyout, ioutputforeachrelease

  implicit none

  integer :: itime,i,ix,jy,ixp,jyp,kz,ks
  integer :: il,ind,indz,indzp,nrelpointer
  real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
  real :: ddx,ddy
  real :: rhoprof(2),rhoi,xl,yl,wx,wy,w


  ! For forward simulations, make a loop over the number of species;
  ! for backward simulations, make an additional loop over the release points
  !**************************************************************************


  if (itra1(i).ne.itime) return

  ! Depending on output option, calculate air density or set it to 1
  ! linit_cond: 1=mass unit, 2=mass mixing ratio unit
  !*****************************************************************


  if (linit_cond.eq.1) then     ! mass unit

    ix=int(xtra1(i))
    jy=int(ytra1(i))
    ixp=ix+1
    jyp=jy+1
    ddx=xtra1(i)-real(ix)
    ddy=ytra1(i)-real(jy)
    rddx=1.-ddx
    rddy=1.-ddy
    p1=rddx*rddy
    p2=ddx*rddy
    p3=rddx*ddy
    p4=ddx*ddy

    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)

  ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
  !*****************************************************************************
    do ind=indz,indzp
      rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
           +p2*rho(ixp,jy ,ind,2) &
           +p3*rho(ix ,jyp,ind,2) &
           +p4*rho(ixp,jyp,ind,2)
    end do
    rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
  elseif (linit_cond.eq.2) then    ! mass mixing ratio unit
    rhoi=1.
  endif


  !****************************************************************************
  ! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy
  !****************************************************************************


  ! For backward simulations, look from which release point the particle comes from
  ! For domain-filling trajectory option, npoint contains a consecutive particle
  ! number, not the release point information. Therefore, nrelpointer is set to 1
  ! for the domain-filling option.
  !*****************************************************************************

  if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then
    nrelpointer=1
  else
    nrelpointer=npoint(i)
  endif

  do kz=1,numzgrid                ! determine height of cell
    if (outheight(kz).gt.ztra1(i)) goto 21
  end do
21   continue
  if (kz.le.numzgrid) then           ! inside output domain


    xl=(xtra1(i)*dx+xoutshift)/dxout
    yl=(ytra1(i)*dy+youtshift)/dyout
    ix=int(xl)
    if (xl.lt.0.) ix=ix-1
    jy=int(yl)
    if (yl.lt.0.) jy=jy-1


  ! If a particle is close to the domain boundary, do not use the kernel either
  !****************************************************************************

    if ((xl.lt.0.5).or.(yl.lt.0.5).or. &
         (xl.gt.real(numxgrid-1)-0.5).or. &
         (yl.gt.real(numygrid-1)-0.5)) then             ! no kernel, direct attribution to grid cell
      if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
           (jy.le.numygrid-1)) then
!$OMP CRITICAL
        do ks=1,nspec
          init_cond(ix,jy,kz,ks,nrelpointer)= &
               init_cond(ix,jy,kz,ks,nrelpointer)+ &
               xmass1(i,ks)/rhoi
        end do
!$OMP END CRITICAL
      endif

    else                                 ! attribution via uniform kernel

      ddx=xl-real(ix)                   ! distance to left cell border
      ddy=yl-real(jy)                   ! distance to lower cell border
      if (ddx.gt.0.5) then
        ixp=ix+1
        wx=1.5-ddx
      else
        ixp=ix-1
        wx=0.5+ddx
      endif

      if (ddy.gt.0.5) then
        jyp=jy+1
        wy=1.5-ddy
      else
        jyp=jy-1
        wy=0.5+ddy
      endif


  ! Determine mass fractions for four grid points
  !**********************************************

!$OMP CRITICAL 
      if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
          w=wx*wy
          do ks=1,nspec
            init_cond(ix,jy,kz,ks,nrelpointer)= &
                 init_cond(ix,jy,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
          end do
        endif

        if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
          w=wx*(1.-wy)
          do ks=1,nspec
            init_cond(ix,jyp,kz,ks,nrelpointer)= &
                 init_cond(ix,jyp,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
          end do
        endif
      endif


      if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
        if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
          w=(1.-wx)*(1.-wy)
          do ks=1,nspec
            init_cond(ixp,jyp,kz,ks,nrelpointer)= &
                 init_cond(ixp,jyp,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
          end do
        endif

        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
          w=(1.-wx)*wy
          do ks=1,nspec
            init_cond(ixp,jy,kz,ks,nrelpointer)= &
                 init_cond(ixp,jy,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
          end do
        endif
      endif
!$OMP END CRITICAL

    endif
  endif

end subroutine initial_cond_calc