prismtrs_line_integral.F90

Go to the documentation of this file.

!------------------------------------------------------------------------
! Copyright 2006-2010, CERFACS, Toulouse, France.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
!
! !ROUTINE: PRISMTrs_line_integral
!
! !INTERFACE!
      subroutine prismtrs_line_integral(lcl_weights, & 
                             in_phi1, in_phi2, theta1, theta2, &
                             grid1_lon, grid2_lon, num_wts)

!-----------------------------------------------------------------------
!
!     intent(in):
!
!-----------------------------------------------------------------------
      USE PRISMDrv, dummy_INTERFACE => PRISMTrs_line_integral
      IMPLICIT NONE

      DOUBLE PRECISION, intent(in) ::  
           in_phi1, in_phi2,      ! longitude endpoints for the segment
           theta1, theta2,           ! latitude  endpoints for the segment
           grid1_lon,  ! reference coordinates for each
           grid2_lon  ! grid (to ensure correct 0,2pi interv.
      INTEGER, intent(in) :: num_wts ! Number of weights for 2D conservative remapping

!-----------------------------------------------------------------------
!
!     intent(out):
!
!-----------------------------------------------------------------------

      DOUBLE PRECISION, dimension(2*num_wts), intent(out) :: 
           lcl_weights   ! line integral contribution to lcl_weights

! !DESCRIPTION
!
! SUBROUTINE PRISMTrs_line_integral computes the line integral 
!     of the flux function that results in the interpolation lcl_weights.  
!     the line is defined by the input lat/lon of the endpoints.
!     This routine comes from the SCRIP 1.4 library available 
!     from Los Alamos National Laboratory.
!
! !REVISED HISTORY
!   Date      Programmer   Description
! ----------  ----------   -----------
! 20/03/2007  S. Valcke    Duplicated from the SCRIP library
!
! EOP
!----------------------------------------------------------------------
! $Id: prismtrs_line_integral.F90 2082 2009-10-21 13:31:19Z hanke $
! $Author: hanke $
!----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!     local variables
!
!-----------------------------------------------------------------------
  CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
     '$Id: prismtrs_line_integral.F90 2082 2009-10-21 13:31:19Z hanke $'
      DOUBLE PRECISION :: dphi, sinth1, sinth2, costh1, costh2, fac, 
                              phi1, phi2, phidiff1, phidiff2
      DOUBLE PRECISION :: f1, f2, fint

#ifdef DEBUGX
  PRINT *, '| | | | | | | Enter PRISMTrs_line_integral'
  call psmile_flushstd
#endif
!-----------------------------------------------------------------------
!
!     lcl_weights for the general case based on a trapezoidal approx to
!     the integrals.
!
!-----------------------------------------------------------------------

      sinth1 = SIN(theta1)
      sinth2 = SIN(theta2)
      costh1 = COS(theta1)
      costh2 = COS(theta2)

      dphi = in_phi1 - in_phi2
      if (dphi >  dble_pi) then
        dphi = dphi - dble_pi2
      else if (dphi < -dble_pi) then
        dphi = dphi + dble_pi2
      endif
      dphi = half*dphi

!-----------------------------------------------------------------------
!
!     the first weight is the area overlap integral. the second and
!     fourth are second-order latitude gradient lcl_weights.
!
!-----------------------------------------------------------------------
      lcl_weights(:) = zero
      lcl_weights(        1) = dphi*(sinth1 + sinth2)
      lcl_weights(num_wts+1) = dphi*(sinth1 + sinth2)
      lcl_weights(        2) = dphi*(costh1 + costh2 + (theta1*sinth1 + &
                                                    theta2*sinth2))
      lcl_weights(num_wts+2) = dphi*(costh1 + costh2 + (theta1*sinth1 + &
                                                    theta2*sinth2))

!-----------------------------------------------------------------------
!
!     the third and fifth lcl_weights are for the second-order phi gradient
!     component.  must be careful of longitude range.
!
!-----------------------------------------------------------------------
      f1 = half*(costh1*sinth1 + theta1)
      f2 = half*(costh2*sinth2 + theta2)

      phi1 = in_phi1 - grid1_lon
      if (phi1 >  dble_pi) then
        phi1 = phi1 - dble_pi2
      else if (phi1 < -dble_pi) then
        phi1 = phi1 + dble_pi2
      endif

      phi2 = in_phi2 - grid1_lon
      if (phi2 >  dble_pi) then
        phi2 = phi2 - dble_pi2
      else if (phi2 < -dble_pi) then
        phi2 = phi2 + dble_pi2
      endif

      if ((phi2-phi1) <  dble_pi .and. (phi2-phi1) > -dble_pi) then
        lcl_weights(3) = dphi*(phi1*f1 + phi2*f2)
      else
        if (phi1 > zero) then
          fac = dble_pi
        else
          fac = -dble_pi
        endif
        fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
        lcl_weights(3) = half*phi1*(phi1-fac)*f1 - &
                     half*phi2*(phi2+fac)*f2 + &
                     half*fac*(phi1+phi2)*fint
      endif

      phi1 = in_phi1 - grid2_lon
      if (phi1 >  dble_pi) then
        phi1 = phi1 - dble_pi2
      else if (phi1 < -dble_pi) then
        phi1 = phi1 + dble_pi2
      endif

      phi2 = in_phi2 - grid2_lon
      if (phi2 >  dble_pi) then
        phi2 = phi2 - dble_pi2
      else if (phi2 < -dble_pi) then
        phi2 = phi2 + dble_pi2
      endif

      if ((phi2-phi1) <  dble_pi .and. (phi2-phi1) > -dble_pi) then
        lcl_weights(num_wts+3) = dphi*(phi1*f1 + phi2*f2)
      else
        if (phi1 > zero) then
          fac = dble_pi
        else
          fac = -dble_pi
        endif
        fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
        lcl_weights(num_wts+3) = half*phi1*(phi1-fac)*f1 - &
                             half*phi2*(phi2+fac)*f2 + &
                             half*fac*(phi1+phi2)*fint
      endif

#ifdef DEBUGX
  PRINT *, '| | | | | | | Quit PRISMTrs_line_integral'
  call psmile_flushstd
#endif

!-----------------------------------------------------------------------

      end subroutine prismtrs_line_integral