Ticket #93: coordtrafo.f90

!***********************************************************************
!* Copyright 2012,2013                                                *
!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
!*                                                                     *
!* This file is part of FLEXPART WRF                                   *
!*                                                                     *
!* 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 coordtrafo
!**********************************************************************
!                                                                     * 
! Note:  This is the FLEXPART_WRF version of subroutine coordtrafo.   *
!                                                                     * 
!             FLEXPART MODEL SUBROUTINE COORDTRAFO                    *
!                                                                     *
!**********************************************************************
!                                                                     * 
! AUTHOR:      G. WOTAWA                                              *
! DATE:        1994-02-07                                             *
! LAST UPDATE: 1996-05-18   A. STOHL                                  *
!                                                                     * 
! Dec 2005, R. Easter - changed names of "*lon0*" & "*lat0*" variables*
! July 2012, J Brioude: modification following flexpart 9             * 
!**********************************************************************
!                                                                     *
! DESCRIPTION: This subroutine transforms x and y coordinates of      *
! particle release points to grid coordinates.                        *
!                                                                     *
!**********************************************************************

  use point_mod
  use par_mod
  use com_mod

  implicit none

  integer :: i,j,k

      if (numpoint.eq.0) goto 30

! TRANSFORM X- AND Y- COORDINATES OF STARTING POINTS TO GRID COORDINATES
!***********************************************************************

      do i=1,numpoint
        xpoint1(i)=(xpoint1(i)-xmet0)/dx
        xpoint2(i)=(xpoint2(i)-xmet0)/dx
        ypoint1(i)=(ypoint1(i)-ymet0)/dy
      ypoint2(i)=(ypoint2(i)-ymet0)/dy
   end do

15    continue


! CHECK IF RELEASE POINTS ARE WITHIN DOMAIN
!******************************************

      do i=1,numpoint
      if ((ypoint1(i).lt.1.e-6).or.(ypoint1(i).ge.real(nymin1)-1.e-6)  &
      .or.(ypoint2(i).lt.1.e-6).or.(ypoint2(i).ge.real(nymin1)-1.e-6)  &
      .or.(xpoint1(i).lt.1.e-6).or.(xpoint1(i).ge.real(nxmin1)-1.e-6) &
      .or.(xpoint2(i).lt.1.e-6).or.(xpoint2(i).ge.real(nxmin1)-1.e-6)) then

!      if ((ypoint1(i).lt.1.e-6).or.(ypoint1(i).ge.real(nymin1)-1.e-6) &
!      .or.(ypoint2(i).lt.1.e-6).or.(ypoint2(i).ge.real(nymin1)-1.e-6) &
!      .or.((.not.xglobal).and.((xpoint1(i).lt.1.e-6).or. &
!      (xpoint1(i).ge.real(nxmin1)-1.e-6).or.(xpoint2(i).lt.1.e-6).or. &
!      (xpoint2(i).ge.real(nxmin1)-1.e-6)))) then

          write(*,*) ' NOTICE: RELEASE POINT OUT OF DOMAIN DETECTED.'
          write(*,*) ' IT IS REMOVED NOW ... '
          if (i.le.2000) then
          write(*,*) ' i, COMMENT: ',i,compoint(i)
          endif
          write(*,*) ' x,y,z,np,mass ',xpoint1(i),ypoint1(i),zpoint1(i),npart(i),xmass(i,1)
          write(*,*) ' REPLACEMENT x,y,z,np,mass ',xpoint1(i+1),ypoint1(i+1),zpoint1(i+1),npart(i+1),xmass(i+1,1)

          if (i.lt.numpoint) then
            do j=i+1,numpoint
              xpoint1(j-1)=xpoint1(j)
              ypoint1(j-1)=ypoint1(j)
              xpoint2(j-1)=xpoint2(j)
              ypoint2(j-1)=ypoint2(j)
              zpoint1(j-1)=zpoint1(j)
              zpoint2(j-1)=zpoint2(j)
              npart(j-1)=npart(j)
           kindz(j-1)=kindz(j)
          ireleasestart(j-1)=ireleasestart(j)
           ireleaseend(j-1)=ireleaseend(j)
              xpoint12(j-1) = xpoint12(j)
              xpoint22(j-1) = xpoint22(j)
              ypoint12(j-1) = ypoint12(j)
              ypoint22(j-1) = ypoint22(j)
              releases_swlon(j-1) = releases_swlon(j)
              releases_swlat(j-1) = releases_swlat(j)
              releases_nelon(j-1) = releases_nelon(j)
              releases_nelat(j-1) = releases_nelat(j)
           if (j.le.2000) compoint(j-1)=compoint(j)
           do k=1,nspec
             xmass(j-1,k)=xmass(j,k)
          end do
              rho_rel(j-1) = rho_rel(j)

         enddo
          endif

          numpoint=numpoint-1
          if (numpoint.gt.0) goto 15

        endif
  end do

30    if (numpoint.eq.0) then
        write(*,*) ' FLEXPART MODEL SUBROUTINE COORDTRAFO: ERROR ! '
        write(*,*) ' NO PARTICLE RELEASES ARE DEFINED!'
        write(*,*) ' CHECK FILE RELEASES...'
        stop
      endif

end subroutine coordtrafo