PROGRAM main
  USE, INTRINSIC :: ISO_FORTRAN_ENV, &
  & ONLY : DFP=>real64, I4B=>int32
  IMPLICIT NONE

  ABSTRACT INTERFACE
    PURE FUNCTION func( x ) RESULT( Ans )
    IMPORT :: DFP
    REAL( DFP ), INTENT( IN ) :: x
    REAL( DFP ) :: Ans
    END FUNCTION func
  END INTERFACE

  PROCEDURE( func ), POINTER :: f => NULL()
  INTEGER( I4B ) :: n

  WRITE( *, "(A)", ADVANCE="NO" ) &
    & "Enter the value of case (n) :: "
  READ( *, * ) n

  SELECT CASE( n )
  CASE( 1 )
    f => func_1
  CASE( 2 )
    f => func_2
  CASE( 3 )
    f =>  func_3
  CASE DEFAULT
    f => func_1
  END SELECT

  WRITE( *, * ) f( 2.0_DFP )

  CONTAINS
    PURE FUNCTION func_1( x ) RESULT( Ans )
      REAL( DFP ), INTENT( IN ) :: x
      REAL( DFP ) :: Ans
      Ans = SIN( x )
    END FUNCTION func_1

    PURE FUNCTION func_2( x ) RESULT( Ans )
      REAL( DFP ), INTENT( IN ) :: x
      REAL( DFP ) :: Ans
      Ans = COS( x )
    END FUNCTION func_2

    PURE FUNCTION func_3( x ) RESULT( Ans )
      REAL( DFP ), INTENT( IN ) :: x
      REAL( DFP ) :: Ans
      Ans = TAN( x )
    END FUNCTION func_3
END PROGRAM main


program main
  use iso_fortran_env, only : input_unit
  use domain_interface, only : domain_t
  use assertions_interface, only : assert
  use module_mp_driver, only: microphysics
  use timer_interface, only: timer_t
  implicit none

  if (this_image()==1) print *,"Number of images = ",num_images()

  block
    type(domain_t), save :: domain
    integer :: i,nz, ypos,xpos
    type(timer_t) :: timer

    if (this_image()==1) print *,this_image(),"domain%initialize_from_file('input-parameters.txt')"
    call domain%initialize_from_file('input-parameters.txt')

    if (this_image()==1) then
        nz = size(domain%pressure,2)
        print *, " Layer height       Pressure        Temperature      Water Vapor"
        print *, "     [m]              [hPa]             [K]            [kg/kg]"
        do i=nz,1,-4
            print *,domain%z(1,i,1), domain%pressure(1,i,1)/100, domain%temperature(1,i,1), domain%water_vapor%local(1,i,1)
        end do
    endif

    ypos = (ubound(domain%accumulated_precipitation,2)-lbound(domain%accumulated_precipitation,2))/2
    ypos = ypos + lbound(domain%accumulated_precipitation,2)

    ! initialize microphysics before starting the timer
    call microphysics(domain, dt = 20.0)
    if (this_image()==1) print*, ""
    if (this_image()==1) print*, "Beginning simulation..."
    sync all
    call timer%start()
    do i=1,200
        ! print *,"Microphysics"
        ! note should this be wrapped into the domain object(?)
        call microphysics(domain, dt = 20.0, halo=1)
        call domain%halo_send()
        call microphysics(domain, dt = 20.0, subset=1)

        call domain%halo_retrieve()

        call domain%advect(dt = 1.0)

        ! if (this_image()==(num_images()/2)) then
        !     print*, domain%accumulated_precipitation(::3,ypos)
        ! endif

    end do
    sync all
    call timer%stop()

    if (this_image()==1) then
        print *,"Model run time:",timer%as_string('(f8.3," seconds")')
    endif

    ypos = (domain%jde-domain%jds)/2 + domain%jds
    do i=1,num_images()
        sync all
        if (this_image()==i) then
            if ((ypos>=domain%jts).and.(ypos<=domain%jte)) then
                xpos = (domain%ite-domain%its)/2 + domain%its
                print*, this_image(), " : ", domain%accumulated_precipitation(domain%its:domain%ite:2,ypos)
            endif
        endif
    enddo

  end block

  if (this_image()==1) print *,"Test passed."

end program


text:
!==============================================================================!
  subroutine Surf_Mod_Handle_4_Points(surf, surf_v)
!------------------------------------------------------------------------------!
!   Surface intersects cell at four points                                     !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Surf_Type), target :: surf
  real                    :: surf_v(3)
!-----------------------------------[Locals]-----------------------------------!
  type(Vert_Type), pointer :: vert(:)
  type(Elem_Type), pointer :: elem(:)
  integer,         pointer :: nv, ne
  integer                  :: ver(4), loop
  real                     :: v_21(3), v_31(3), v_41(3)
  real                     :: tri_p_123(3), tri_p_134(3)
  integer                  :: permutations(4, 6)
!=============================================================================!

  ! Take aliases
  nv   => surf % n_verts
  ne   => surf % n_elems
  vert => surf % vert
  elem => surf % elem

  permutations = reshape((/ &
    0, 1, 2, 3,  &
    0, 2, 1, 3,  &
    0, 3, 1, 2,  &
    0, 1, 3, 2,  &
    0, 2, 3, 1,  &
    0, 3, 2, 1   &
  /), shape(permutations))

  !-------------------------------------------!
  !   Try to find a permutation which works   !
  !-------------------------------------------!
  do loop = 1, 6

    ver(1) = nv - permutations(1, loop)
    ver(2) = nv - permutations(2, loop)
    ver(3) = nv - permutations(3, loop)
    ver(4) = nv - permutations(4, loop)

    v_21(1) = vert(ver(2)) % x_n - vert(ver(1)) % x_n
    v_21(2) = vert(ver(2)) % y_n - vert(ver(1)) % y_n
    v_21(3) = vert(ver(2)) % z_n - vert(ver(1)) % z_n

    v_31(1) = vert(ver(3)) % x_n - vert(ver(1)) % x_n
    v_31(2) = vert(ver(3)) % y_n - vert(ver(1)) % y_n
    v_31(3) = vert(ver(3)) % z_n - vert(ver(1)) % z_n

    v_41(1) = vert(ver(4)) % x_n - vert(ver(1)) % x_n
    v_41(2) = vert(ver(4)) % y_n - vert(ver(1)) % y_n
    v_41(3) = vert(ver(4)) % z_n - vert(ver(1)) % z_n

    tri_p_123 = Math_Mod_Cross_Product(v_21, v_31)
    tri_p_134 = Math_Mod_Cross_Product(v_31, v_41)

    if(dot_product(surf_v, tri_p_123) > 0.0 .and.  &
       dot_product(surf_v, tri_p_134) > 0.0) then
      exit
    end if

  end do

  !--------------------------------------------------------------------!
  !   You've found a permutation which works, store new elements now   !
  !--------------------------------------------------------------------!

  ! One new element with three vertices
  ne = ne + 1
  elem(ne) % i = ver(1)
  elem(ne) % j = ver(2)
  elem(ne) % k = ver(3)

  ! Second new element with three vertices
  ne = ne + 1
  elem(ne) % i = ver(1)
  elem(ne) % j = ver(3)
  elem(ne) % k = ver(4)

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
#if 0
!>
!! @file
!! This file is part of SeisSol.
!!
!! @author Sebastian Rettenberger (sebastian.rettenberger AT tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger,_M.Sc.)
!!
!! @section LICENSE
!! Copyright (c) 2014, SeisSol Group
!! All rights reserved.
!!
!! Redistribution and use in source and binary forms, with or without
!! modification, are permitted provided that the following conditions are met:
!!
!! 1. Redistributions of source code must retain the above copyright notice,
!!    this list of conditions and the following disclaimer.
!!
!! 2. Redistributions in binary form must reproduce the above copyright notice,
!!    this list of conditions and the following disclaimer in the documentation
!!    and/or other materials provided with the distribution.
!!
!! 3. Neither the name of the copyright holder nor the names of its
!!    contributors may be used to endorse or promote products derived from this
!!    software without specific prior written permission.
!!
!! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
!! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
!! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
!! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
!! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!! POSSIBILITY OF SUCH DAMAGE.
!!
!! @section DESCRIPTION
!! Defines macros for instrumentation
#endif

#ifndef INSTRUMENTATION_FPP
#define INSTRUMENTATION_FPP

#if 0
! Manual instrumentation for Scalasca with epik.
! Override function calls if not compiled with EPIK.
#endif

#if defined(EPIK)
#if (!defined(__STDC__) || (defined(BG) && !defined(__bg__)))
#include "epik_user.inc"
#else
#include "epik_user.h"
#endif
#else
#define EPIK_FUNC_REG(str)
#define EPIK_FUNC_START()
#define EPIK_FUNC_END()
#define EPIK_USER_REG(id,str)
#define EPIK_USER_START(id)
#define EPIK_USER_END(id)
#define EPIK_TRACER(str)
#endif

#if 0
! empty score-p definitions without usage
#endif

#if defined(SCOREP_USER_ENABLE)

#if (!defined(__STDC__) || (defined(BG) && !defined(__bg__)))
#include <scorep/SCOREP_User.inc>
#else
#include <scorep/SCOREP_User.h>
#endif

#else

#if (!defined(__STDC__) || (defined(BG) && !defined(__bg__)))
#define SCOREP_USER_REGION_DEFINE( handle )
#define SCOREP_USER_REWIND_DEFINE( handle )
#define SCOREP_USER_OA_PHASE_BEGIN( handle, name, type  )
#define SCOREP_USER_OA_PHASE_END( handle )
#define SCOREP_USER_REWIND_POINT( handle, name )
#define SCOREP_USER_REGION_BEGIN( handle, name, type )
#define SCOREP_USER_REGION_INIT( handle, name, type )
#define SCOREP_USER_REGION_END( handle )
#define SCOREP_USER_REWIND_CHECK( handle, value )
#define SCOREP_USER_REGION_ENTER( handle )
#define SCOREP_USER_FUNC_DEFINE()
#define SCOREP_USER_FUNC_BEGIN( name )
#define SCOREP_USER_FUNC_END()
#define SCOREP_USER_GLOBAL_REGION_DEFINE( handle )
#define SCOREP_USER_GLOBAL_REGION_EXTERNAL( handle )
#define SCOREP_USER_PARAMETER_DEFINE( handle )
#define SCOREP_USER_PARAMETER_INT64( handle, name, value )
#define SCOREP_USER_PARAMETER_UINT64( handle, name, value )
#define SCOREP_USER_PARAMETER_STRING( handle, name, value )
#define SCOREP_USER_METRIC_GLOBAL( metricHandle )
#define SCOREP_USER_METRIC_EXTERNAL( metricHandle )
#define SCOREP_USER_METRIC_LOCAL( metricHandle )
#define SCOREP_USER_METRIC_INIT( metricHandle, name, unit, type, context )
#define SCOREP_USER_METRIC_INT64( metricHandle, value )
#define SCOREP_USER_METRIC_UINT64( metricHandle, value )
#define SCOREP_USER_METRIC_DOUBLE( metricHandle, value )
#define SCOREP_RECORDING_ON()
#define SCOREP_RECORDING_OFF()
#define SCOREP_RECORDING_IS_ON( isOn ) isOn = .false.
#else
#define SCOREP_USER_REGION( name, type )
#define SCOREP_USER_REGION_DEFINE( handle )
#define SCOREP_USER_OA_PHASE_BEGIN( handle, name, type  )
#define SCOREP_USER_OA_PHASE_END( handle )
#define SCOREP_USER_REGION_BEGIN( handle, name, type )
#define SCOREP_USER_REGION_INIT( handle, name, type )
#define SCOREP_USER_REGION_END( handle )
#define SCOREP_USER_REGION_ENTER( handle )
#define SCOREP_USER_FUNC_BEGIN()
#define SCOREP_USER_FUNC_END()
#define SCOREP_GLOBAL_REGION_DEFINE( handle )
#define SCOREP_GLOBAL_REGION_EXTERNAL( handle )
#define SCOREP_USER_PARAMETER_INT64( name, value )
#define SCOREP_USER_PARAMETER_UINT64( name, value )
#define SCOREP_USER_PARAMETER_STRING( name, value )
#define SCOREP_USER_METRIC_GLOBAL( metricHandle )
#define SCOREP_USER_METRIC_EXTERNAL( metricHandle )
#define SCOREP_USER_METRIC_LOCAL( metricHandle )
#define SCOREP_USER_METRIC_INIT( metricHandle, name, unit, type, context )
#define SCOREP_USER_METRIC_INT64( metricHandle, value )
#define SCOREP_USER_METRIC_UINT64( metricHandle, value )
#define SCOREP_USER_METRIC_DOUBLE( metricHandle, value )
#define SCOREP_RECORDING_ON()
#define SCOREP_RECORDING_OFF()
#define SCOREP_RECORDING_IS_ON() 0
#endif

#endif

#endif

ext: fpp

max_stars_repo_name: fabian-kutschera/SeisSol

max_issues_repo_name: fabian-kutschera/SeisSol

max_forks_repo_name: fabian-kutschera/SeisSol


!        Generated by TAPENADE     (INRIA, Ecuador team)
!  tapenade 3.x
!
!  Differentiation of lin_interpol in reverse (adjoint) mode (with options noISIZE i8):
!   gradient     of useful results: lin_interpol
!   with respect to varying inputs: da db
!>    @brief linear interpolation
!>    @param[in] da value A
!>    @param[in] db value B
!>    @param[in] dp relative position
!>    @param[in] dd (delta) cell size
!>    @return linear interpolated value
SUBROUTINE LIN_INTERPOL_AD0(da, da_ad, db, db_ad, dp, dd, &
& lin_interpol_ad)
  IMPLICIT NONE
!       da: value A; db: value B; dp: relative position; dd: (delta) cell size
  DOUBLE PRECISION :: da, db, dp, dd
  DOUBLE PRECISION :: da_ad, db_ad
  DOUBLE PRECISION :: lin_interpol
  DOUBLE PRECISION :: lin_interpol_ad
  da_ad = (dd-dp)*lin_interpol_ad/dd
  db_ad = dp*lin_interpol_ad/dd
END SUBROUTINE LIN_INTERPOL_AD0


text:
! { dg-do compile }
! PR20849 - An external symbol may not have a initializer.
! Contributed by Joost VandeVondele <jv244@cam.ac.uk>
REAL, EXTERNAL :: X=0 ! { dg-error "may not have an initializer" }
END

max_stars_repo_name: vidkidz/crossbridge

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED 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
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.

!> @brief Computes spatial derivatives of displacement.
!! @author Ilario Mazzieri
!> @date September, 2013
!> @version 1.0
!> @param[in] nn number of 1-D Legendre nodes
!> @param[in] ct GLL nodes
!> @param[in] ww GLL weights
!> @param[in] dd matrix of spectral derivatives
!> @param[in] ux x-displacement
!> @param[in] uy y-displacement
!> @param[in] uy z-displacement
!> @param[in] alfa11 costant for the bilinear map
!> @param[in] alfa12 costant for the bilinear map
!> @param[in] alfa13 costant for the bilinear map
!> @param[in] alfa21 costant for the bilinear map
!> @param[in] alfa22 costant for the bilinear map
!> @param[in] alfa23 costant for the bilinear map
!> @param[in] alfa31 costant for the bilinear map
!> @param[in] alfa32 costant for the bilinear map
!> @param[in] alfa33 costant for the bilinear map
!> @param[in] beta11 costant for the bilinear map
!> @param[in] beta12 costant for the bilinear map
!> @param[in] beta13 costant for the bilinear map
!> @param[in] beta21 costant for the bilinear map
!> @param[in] beta22 costant for the bilinear map
!> @param[in] beta23 costant for the bilinear map
!> @param[in] beta31 costant for the bilinear map
!> @param[in] beta32 costant for the bilinear map
!> @param[in] beta33 costant for the bilinear map
!> @param[in] gamma1  costant for the bilinear map
!> @param[in] gamma2  costant for the bilinear map
!> @param[in] gamma3  costant for the bilinear map
!> @param[in] delta1  costant for the bilinear map
!> @param[in] delta2  costant for the bilinear map
!> @param[in] delta3  costant for the bilinear map
!> @param[out] duxdx nodal values for spatial derivatives of the displacement
!> @param[out] duydx nodal values for spatial derivatives of the displacement
!> @param[out] duzdx nodal values for spatial derivatives of the displacement
!> @param[out] duxdy nodal values for spatial derivatives of the displacement
!> @param[out] duydy nodal values for spatial derivatives of the displacement
!> @param[out] duzdy nodal values for spatial derivatives of the displacement
!> @param[out] duxdz nodal values for spatial derivatives of the displacement
!> @param[out] duydz nodal values for spatial derivatives of the displacement
!> @param[out] duzdz nodal values for spatial derivatives of the displacement

      subroutine MAKE_STRAIN_TENSOR(nn,ct,ww,dd,&
                               alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,&
                               alfa31,alfa32,alfa33,beta11,beta12,beta13,&
                               beta21,beta22,beta23,beta31,beta32,beta33,&
                               gamma1,gamma2,gamma3,delta1,delta2,delta3,&
                               ux,uy,uz,&
                               duxdx,duydx,duzdx,&
                               duxdy,duydy,duzdy,&
                               duxdz,duydz,duzdz)

      implicit none

      integer*4 :: nn
      integer*4 :: i,j,k,p,q,r

      real*8 :: alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,alfa31,alfa32,alfa33
      real*8 :: beta11,beta12,beta13,beta21,beta22,beta23,beta31,beta32,beta33
      real*8 :: gamma1,gamma2,gamma3,delta1,delta2,delta3
      real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j
      real*8 :: t1ux,t2ux,t3ux,t1uy,t2uy,t3uy,t1uz,t2uz,t3uz

      real*8, dimension(nn) :: ct,ww
      real*8, dimension(nn,nn) :: dd

      real*8, dimension(nn,nn,nn) :: ux,uy,uz
      real*8, dimension(nn,nn,nn) :: duxdx,duxdy,duxdz,duydx,duydy,duydz,duzdx,duzdy,duzdz

!     STRAIN CALCULATION

      do r = 1,nn
         do q = 1,nn
            do p = 1,nn
               t1ux = 0.d0
               t1uy = 0.d0
               t1uz = 0.d0
               t2ux = 0.d0
               t2uy = 0.d0
               t2uz = 0.d0
               t3ux = 0.d0
               t3uy = 0.d0
               t3uz = 0.d0

               dxdx = alfa11 + beta12*ct(r) + beta13*ct(q) &
                    + gamma1*ct(q)*ct(r)
               dydx = alfa21 + beta22*ct(r) + beta23*ct(q) &
                    + gamma2*ct(q)*ct(r)
               dzdx = alfa31 + beta32*ct(r) + beta33*ct(q) &
                    + gamma3*ct(q)*ct(r)

               dxdy = alfa12 + beta11*ct(r) + beta13*ct(p) &
                    + gamma1*ct(r)*ct(p)
               dydy = alfa22 + beta21*ct(r) + beta23*ct(p) &
                    + gamma2*ct(r)*ct(p)
               dzdy = alfa32 + beta31*ct(r) + beta33*ct(p) &
                    + gamma3*ct(r)*ct(p)

               dxdz = alfa13 + beta11*ct(q) + beta12*ct(p) &
                    + gamma1*ct(p)*ct(q)
               dydz = alfa23 + beta21*ct(q) + beta22*ct(p) &
                    + gamma2*ct(p)*ct(q)
               dzdz = alfa33 + beta31*ct(q) + beta32*ct(p) &
                    + gamma3*ct(p)*ct(q)

               det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
                     - dydz * (dxdx*dzdy - dzdx*dxdy) &
                     + dzdz * (dxdx*dydy - dydx*dxdy)

               do i = 1,nn
                  t1ux = t1ux + ux(i,q,r) * dd(p,i)
                  t1uy = t1uy + uy(i,q,r) * dd(p,i)
                  t1uz = t1uz + uz(i,q,r) * dd(p,i)
               enddo

               do j = 1,nn
                  t2ux = t2ux + ux(p,j,r) * dd(q,j)
                  t2uy = t2uy + uy(p,j,r) * dd(q,j)
                  t2uz = t2uz + uz(p,j,r) * dd(q,j)
               enddo

               do k = 1,nn
                  t3ux = t3ux + ux(p,q,k) * dd(r,k)
                  t3uy = t3uy + uy(p,q,k) * dd(r,k)
                  t3uz = t3uz + uz(p,q,k) * dd(r,k)
               enddo

               duxdx(p,q,r) = 1.0d0 / det_j *(&
                    ((dydy*dzdz - dydz*dzdy) * t1ux) + &
                    ((dydz*dzdx - dydx*dzdz) * t2ux) + &
                    ((dydx*dzdy - dydy*dzdx) * t3ux))

               duydx(p,q,r) = 1.0d0 / det_j *(&
                    ((dydy*dzdz - dydz*dzdy) * t1uy) + &
                    ((dydz*dzdx - dydx*dzdz) * t2uy) + &
                    ((dydx*dzdy - dydy*dzdx) * t3uy))

               duzdx(p,q,r) = 1.0d0 / det_j *(&
                    ((dydy*dzdz - dydz*dzdy) * t1uz) + &
                    ((dydz*dzdx - dydx*dzdz) * t2uz) + &
                    ((dydx*dzdy - dydy*dzdx) * t3uz))

               duxdy(p,q,r) = 1.0d0 / det_j *(&
                    ((dzdy*dxdz - dzdz*dxdy) * t1ux) + &
                    ((dzdz*dxdx - dzdx*dxdz) * t2ux) + &
                    ((dzdx*dxdy - dzdy*dxdx) * t3ux))

               duydy(p,q,r) = 1.0d0 / det_j *(&
                    ((dzdy*dxdz - dzdz*dxdy) * t1uy) + &
                    ((dzdz*dxdx - dzdx*dxdz) * t2uy) + &
                    ((dzdx*dxdy - dzdy*dxdx) * t3uy))

               duzdy(p,q,r) = 1.0d0 / det_j *(&
                    ((dzdy*dxdz - dzdz*dxdy) * t1uz) + &
                    ((dzdz*dxdx - dzdx*dxdz) * t2uz) + &
                    ((dzdx*dxdy - dzdy*dxdx) * t3uz))

               duxdz(p,q,r) = 1.0d0 / det_j *(&
                    ((dxdy*dydz - dxdz*dydy) * t1ux) + &
                    ((dxdz*dydx - dxdx*dydz) * t2ux) + &
                    ((dxdx*dydy - dxdy*dydx) * t3ux))

               duydz(p,q,r) = 1.0d0 / det_j *(&
                    ((dxdy*dydz - dxdz*dydy) * t1uy) + &
                    ((dxdz*dydx - dxdx*dydz) * t2uy) + &
                    ((dxdx*dydy - dxdy*dydx) * t3uy))

               duzdz(p,q,r) = 1.0d0 / det_j *(&
                    ((dxdy*dydz - dxdz*dydy) * t1uz) + &
                    ((dxdz*dydx - dxdx*dydz) * t2uz) + &
                    ((dxdx*dydy - dxdy*dydx) * t3uz))

            enddo
         enddo
      enddo

      return

      end subroutine MAKE_STRAIN_TENSOR


program AdvectSlottedCylindersDriver

use NumberKindsModule
use OutputWriterModule
use LoggerModule
use ParticlesModule
use EdgesModule
use FacesModule
use PolyMesh2dModule
use FieldModule
use MPISetupModule
use SphereGeomModule
use RefinementModule
use SSRFPACKInterfaceModule
use SSRFPACKRemeshModule
use SphereTransportSolverModule
use SphereTracersModule
use SphereTransportModule
use SphereTransportVelocitiesModule, velFn => LauritzenEtalDeformationalVelocity

implicit none

include 'mpif.h'

! mesh variables
type(TransportMesh) :: sphere
type(RefineSetup) :: refinement
integer(kint) :: initNest
integer(kint) :: maxNest
integer(kint) :: faceKind
integer(kint) :: meshSeed
integer(kint) :: amrLimit
real(kreal), parameter :: radius = 1.0_kreal
logical(klog) :: continueAMR
logical(klog) :: AMR
integer(kint) :: nParticlesBefore
integer(kint) :: nParticlesAfter
namelist /meshDefine/ faceKind, initNest, amrLimit

! test case variables
real(kreal), allocatable, dimension(:) :: trMass
real(kreal) :: l2Err, l2Denom
real(kreal) :: lInfErr
real(kreal) :: qMinTrue, qMaxTrue
real(kreal) :: qMinComp, qMaxComp
real(kreal) :: qMinErr, qMaxErr
real(kreal), parameter :: qRange = 0.9_kreal
integer(kint), parameter :: nTracers = 3
integer(kint), dimension(3), parameter :: tracerDims = [1,1,1]

! remeshing variables
type(TransportRemesh) :: remesh
integer(kint) :: remeshInterval
integer(kint) :: remeshCounter
type(TransportMesh) :: tempSphere
logical(klog) :: useDirectRemesh

! timestepping
type(TransportSolver) :: solver
real(kreal) :: dt
real(kreal) :: t
real(kreal) :: tfinal
integer(kint) :: nTimesteps
integer(kint) :: timeJ
namelist /timestepping/ dt, tfinal, remeshInterval, useDirectRemesh

! i/o
character(len=MAX_STRING_LENGTH) :: outputDir
character(len=MAX_STRING_LENGTH) :: outputRoot
character(len=MAX_STRING_LENGTH) :: vtkFile
character(len=MAX_STRING_LENGTH) :: matlabFile
character(len=MAX_STRING_LENGTH) :: vtkRoot
character(len=MAX_STRING_LENGTH) :: meshString
integer(kint) :: frameOut
integer(kint) :: frameCounter
namelist /fileIO/ outputDir, outputRoot, frameOut

! computing environment / general
type(Logger) :: exeLog
character(len=MAX_STRING_LENGTH) :: logString
integer(kint), parameter :: logLevel = DEBUG_LOGGING_LEVEL
character(len=28) :: logKey = "AdvectSlottedC"
integer(kint) :: mpiErrCode
real(kreal) :: timeStart, timeEnd
integer(kint) :: i
real(kreal), dimension(3) :: vec

!--------------------------------
!	initialize : setup computing environment
!--------------------------------
call MPI_INIT(mpiErrCode)
call MPI_COMM_SIZE(MPI_COMM_WORLD, numProcs, mpiErrCode)
call MPI_COMM_RANK(MPI_COMM_WORLD, procRank, mpiErrCode)

call InitLogger(exeLog, procRank)

timeStart = MPI_WTIME()

call ReadNamelistFile( procRank )

!
!	initialize mesh and spatial fields
!
t = 0.0_kreal
call New( sphere, meshSeed, initNest, maxNest, amrLimit, radius, .FALSE.)
call AddTracers(sphere, nTracers, tracerDims)
sphere%tracers(1)%name = "slotC"
sphere%tracers(2)%name = "initialLatitude"
sphere%tracers(3)%name = "relError"
sphere%tracers(3)%N = sphere%mesh%particles%N
call SetInitialDensityOnMesh(sphere)
call SetTracerOnMesh( sphere, 1, SlottedCylinderTracer )

call SetVelocityOnMesh( sphere, velFn, t)
call SetDivergenceOnMesh(sphere)

do i = 1, sphere%mesh%particles%N
	vec = LagCoord(sphere%mesh%particles, i)
	call InsertScalarToField( sphere%tracers(2), Latitude(vec) )
enddo

! TO DO : AMR

!
!	Output initial data
!
frameCounter = 0
remeshCounter = 0
if ( procRank == 0 ) then
	call LogStats(sphere, exeLog)
	if (meshSeed == ICOS_TRI_SPHERE_SEED) then
		if ( initNest == maxNest) then
			write(meshString, '(A,I1,A)') '_icosTri', initNest, '_'
		else
			write(meshString, '(2(A,I1),A)') '_icosTriAMR', initNest, 'to', maxNest, '_'
		endif
	elseif (meshSeed == CUBED_SPHERE_SEED ) then
		if ( initNest == maxNest ) then
			write(meshString, '(A,I1,A)') '_cubedSphere', initNest, '_'
		else
			write(meshString, '(2(A,I1),A)') '_cubedSphereAMR', initNest, 'to', maxNest, '_'
		endif
	endif

	write(vtkRoot,'(4A)') trim(outputDir), '/vtkOut/', trim(outputRoot), trim(meshString)
	write(vtkFile,'(A,I0.4,A)') trim(vtkRoot), frameCounter, '.vtk'

	call OutputToVTK(sphere, vtkFile)
	frameCounter = frameCounter + 1

	call LogMessage(exeLog, TRACE_LOGGING_LEVEL, trim(logkey)//" t = ", t)
endif

!
!	initialize time stepping
!
call New(solver, sphere)
nTimesteps = floor( tfinal/dt )
allocate(trMass( nTimesteps + 1))
trMass(1) = TracerMass(sphere, 1)

qMinTrue = MinScalarVal(sphere%tracers(1))
qMaxTrue = MaxScalarVal(sphere%tracers(1))

!--------------------------------
!	run : evolve the problem in time
!--------------------------------

call LogMessage(exeLog, DEBUG_LOGGING_LEVEL, trim(logkey)//" ", "starting timestepping loop.")

do timeJ = 0, nTimesteps - 1

	if ( mod(timeJ+1, remeshInterval) == 0 ) then
		call LogMessage(exeLog, TRACE_LOGGING_LEVEL, trim(logkey)//" remesh triggered by remesh interval : ",&
			 remeshCounter)
		call New(remesh, sphere)

		call New(tempSphere, meshSeed, initNest, maxNest, amrLimit, radius, .FALSE.)
		call AddTracers(tempSphere, nTracers, tracerDims)
		tempSphere%tracers(1)%name = "slotC"
		tempSphere%tracers(2)%name = "initialLatitude"
		tempSphere%tracers(3)%name = "relError"

		call LagrangianRemeshTransportWithFunctions(remesh, sphere, tempSphere, .FALSE., velFn, t, &
			tracerFn1 = SlottedCylinderTracer, tracerFn2 = InitLatTracer )

		call Copy(sphere, tempSphere)
		remeshCounter = remeshCounter + 1

		call Delete(tempSphere)
		call Delete(remesh)

		call Delete(solver)
		call New(solver, sphere)
	endif

	call Timestep(solver, sphere, t, dt, velFn)

	t = real(timeJ +1, kreal) * dt
	sphere%mesh%t = t

	trMass(timeJ+2) = TracerMass(sphere, 1)

	if ( timeJ+1 == nTimesteps ) then
	!--------------------------------
	!	calculate error at each particle
	!--------------------------------
		do i = 1, sphere%mesh%particles%N
			sphere%tracers(3)%scalar(i) = abs(SlottedCylinderTracer( sphere%mesh%particles%x(i), &
				 sphere%mesh%particles%y(i), sphere%mesh%particles%z(i) ) - sphere%tracers(1)%scalar(i) ) / qRange
		enddo
	endif

	if ( procRank == 0 .AND. mod(timeJ+1, frameOut) == 0 ) then
		write(vtkFile,'(A,I0.4,A)') trim(vtkRoot), frameCounter, '.vtk'
		call OutputToVTK(sphere, vtkFile)
		frameCounter = frameCounter + 1
		call LogMessage(exelog, TRACE_LOGGING_LEVEL, trim(logKey)//" t = ", t)
	endif
enddo

!
!	write t = tfinal output
!
if ( procRank == 0 ) then
	write(matlabFile, '(5A)') trim(outputDir), '/', trim(outputRoot), trim(meshString), '.m'

	l2Err = 0.0_kreal
	l2Denom = 0.0_kreal
	do i = 1, sphere%mesh%particles%N
		if ( sphere%mesh%particles%isActive(i) ) then
			l2Err = l2Err + sphere%tracers(3)%scalar(i) * sphere%tracers(3)%scalar(i) * sphere%mesh%particles%area(i)
			l2Denom = l2Denom + SlottedCylinderTracer( sphere%mesh%particles%x(i), sphere%mesh%particles%y(i), &
				  sphere%mesh%particles%z(i) ) ** 2 * sphere%mesh%particles%area(i)
		endif
	enddo
	l2Err = l2Err / l2Denom
	lInfErr = MaxScalarVal(sphere%tracers(3))

	qMinComp = MinScalarVal(sphere%tracers(1))
	qMaxComp = MaxScalarVal(sphere%tracers(1))

	qMinErr = (qMinComp - qMinTrue) / qRange
	qMaxErr = (qMaxComp - qMaxTrue) / qRange

	call StartSection(exeLog, "Final Errors: "//meshString )
		call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "l2Err = ", l2Err )
		call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "lInfErr = ", lInfErr )
		call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "qMinErr = ", qMinErr )
		call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "qMaxErr = ", qMaxErr )
		call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "rel. tracer mass change = ", &
											maxval(abs(trMass - trMass(1))) / trMass(1) )
		call LogMessage(exeLog, TRACE_LOGGING_LEVEL, " ", " ")
	call EndSection(exeLog)

	open(unit=WRITE_UNIT_1, file=matlabFile, status='REPLACE', action='WRITE')
		write(WRITE_UNIT_1,'(A,F12.9,A,F12.6,A)') "t = 0:", dt,":", tfinal, ";"
		call WriteToMatlab(trMass, WRITE_UNIT_1, "trMass")
		call WriteToMatlab(l2Err, WRITE_UNIT_1, "l2Err")
		call WriteToMatlab(lInfErr, WRITE_UNIT_1, "lInfErr")
		call WriteToMatlab(qMinErr, WRITE_UNIT_1, "qMinErr")
		call WriteToMatlab(qMaxErr, WRITE_UNIT_1, "qMaxErr")
		call WriteToMatlab(maxval(abs(trMass - trMass(1))) / trMass(1), WRITE_UNIT_1, "tracerMassChange")
	close(WRITE_UNIT_1)
endif

!--------------------------------
!	finalize : clean up
!--------------------------------

timeEnd = MPI_WTIME()
write(logstring,'(A,F12.2,A)') "PROGRAM COMPLETE : elapsed time ", timeEnd - timeStart, " seconds."
call LogMessage(exeLog, TRACE_LOGGING_LEVEL, trim(logKey)//" ", logstring)

call Delete(solver)
call Delete(sphere)

deallocate(trMass)

call MPI_FINALIZE(mpiErrCode)

contains

subroutine ReadNamelistFile( rank )
	integer(kint), intent(in) :: rank
	!
	character(len=MAX_STRING_LENGTH) :: namelistFilename
	integer(kint), parameter :: initBcast_intSize = 7
	integer(kint), parameter :: initBcast_realSize = 2
	integer(kint), dimension(initBcast_intSize) :: bcastIntegers
	real(kreal), dimension(initBcast_realSize) :: bcastReals
	integer(kint) :: mpiErrCode, readStat

	if ( COMMAND_ARGUMENT_COUNT() /= 1 ) then
		call LogMessage(exeLog, ERROR_LOGGING_LEVEL, trim(logKey), " ERROR: expected namelist file as 1st argument.")
		stop
	endif

	if ( rank == 0 ) then
		call GET_COMMAND_ARGUMENT(1, namelistFilename)

		open(unit=READ_UNIT, file=namelistFilename, status='OLD', action='READ', iostat=readStat)
			if ( readStat /= 0 ) then
				call LogMessage(exeLog, ERROR_LOGGING_LEVEL, trim(logKey), " ERROR: cannot read namelist file.")
				stop
			endif

			read(READ_UNIT, nml=meshDefine)
			rewind(READ_UNIT)
			read(READ_UNIT, nml=timestepping)
			rewind(READ_UNIT)
			read(READ_UNIT, nml=fileIO)
		close(READ_UNIT)

		if ( faceKind == 3 ) then
			meshSeed = ICOS_TRI_SPHERE_SEED
		elseif ( faceKind == 4) then
			meshSeed = CUBED_SPHERE_SEED
		else
			call LogMessage(exeLog, WARNING_LOGGING_LEVEL, trim(logkey)//" ReadNamelistFile WARNING : ", &
				" invalid faceKind -- using triangles.")
			meshSeed = ICOS_TRI_SPHERE_SEED
		endif

		maxNest = initNest + amrLimit

!		namelist /meshDefine/ faceKind, initNest, maxNest, amrLimit, radius
!		namelist /timestepping/ dt, tfinal
!		namelist /fileIO/ outputDir, outputRoot, frameOut

		bcastIntegers(1) = meshSeed
		bcastIntegers(2) = initNest
		bcastIntegers(3) = maxNest
		bcastIntegers(4) = amrLimit
		bcastIntegers(5) = frameOut
		bcastIntegers(6) = remeshInterval
		if ( useDirectRemesh ) then
			bcastIntegers(7) = 1
		else
			bcastIntegers(7) = 0
		endif

		bcastReals(1) = dt
		bcastReals(2) = tfinal
	endif

	call MPI_BCAST(bcastIntegers, initBCAST_intSize, MPI_INTEGER, 0, MPI_COMM_WORLD, mpiErrCode)
	if ( mpiErrCode /= 0 ) then
		call LogMessage(exeLog, ERROR_LOGGING_LEVEL, trim(logKey)//" bcastIntegers, MPI_BCAST ERROR : ", mpiErrCode)
	endif

	call MPI_BCAST(bcastReals, initBCAST_realSize, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, mpiErrCode)
	if ( mpiErrCode /= 0 ) then
		call LogMessage(exeLog, ERROR_LOGGING_LEVEL, trim(logKey)//" bcastReals, MPI_BCAST ERROR : ", mpiErrCode)
	endif

	meshSeed = bcastIntegers(1)
	initNest = bcastIntegers(2)
	maxNest = bcastIntegers(3)
	amrLimit = bcastIntegers(4)
	frameOut = bcastIntegers(5)
	remeshInterval = bcastIntegers(6)
	if (bcastIntegers(7) > 0 ) then
		useDirectRemesh = .TRUE.
	else
		useDirectRemesh = .FALSE.
	endif

	dt = bcastReals(1)
	tfinal = bcastReals(2)
end subroutine

!> @brief Initializes a @ref Logger for this executable program.
!>
!> Output is controlled by message priority level and MPI rank.
!>
!> @param[in] log @ref Logger to initialize
!> @param[in] rank MPI rank
subroutine InitLogger(log, rank)
	type(Logger), intent(inout) :: log
	integer(kint), intent(in) :: rank
	if ( rank == 0 ) then
		call New(log, logLevel)
	else
		call New(log, ERROR_LOGGING_LEVEL)
	endif
	write(logKey,'(A,I0.2,A)') trim(logKey)//"_", rank, ":"
end subroutine

end program


text:
module EQS_EDROMO
! Description:
!    Contains the subroutine for computing the RHS of the EDromo equations in
!    double precision.
!
! Author:
!    Davide Amato
!    Space Dynamics Group - Technical University of Madrid
!    d.amato@upm.es
!
! ==============================================================================

use KINDS, only: ik,dk,qk
implicit none

contains

subroutine DEDROMO_RHS(neq,phi,z,zdot)
! Description:
!    Computes the value of the right-hand side of the equations of motion of the
!    EDromo formulation.
!
! ==============================================================================

! MODULES
use AUXILIARIES, only: inSoI
use SETTINGS,    only: flag_time_EDr
use TRANSFORM,   only: DINERT2ORB_EDROMO
use THIRDBODIES, only: DACC3B_ND
use CONSTANTS,   only: DU,TU,wEarth,smaEarth,muEarth,muSun

! VARIABLES
implicit none

! Arguments
integer(ik),intent(in)     ::  neq             ! Number of equations
real(dk),intent(in)        ::  phi             ! EDromo independent variable
real(dk),intent(in)        ::  z(1:neq)        ! EDromo state vector, ND
real(dk),intent(out)       ::  zdot(1:neq)     ! RHS of EoM's, ND

! Auxiliary quantities
real(dk)    ::  sph,cph
real(dk)    ::  rho,rmag,zeta,emme
real(dk)    ::  cnu,snu
real(dk)    ::  enne,L3,wz
real(dk)    ::  aux0,aux1,aux2,aux3,aux4

! State
real(dk)    ::  rV(1:3),t
real(dk)    ::  x_vec(1:3),y_vec(1:3)

! Potential
real(dk)    ::  Upot,dUdt,dUdr(1:3)

! Perturbations
real(dk)    ::  lon_Earth,r2(1:3)
real(dk)    ::  f(1:3),p(1:3)

! ==============================================================================

! INDEX OF STATE VECTORS

! z(1) = zeta1
! z(2) = zeta2
! z(3) = zeta3
! z(4) = zeta4
! z(5) = zeta5
! z(6) = zeta6
! z(7) = zeta7
! z(8) = zeta8 (physical time,
!                      constant time element,
!                      linear time element,
!                      depending on flag_time)

! ==============================================================================
! 01. COMPUTE AUXILIARY QUANTITIES (1)
! ==============================================================================

! Store trig functions
sph = sin(phi)
cph = cos(phi)

rho   = 1._dk - z(1)*cph - z(2)*sph
rmag  = z(3)*rho
zeta  = z(1)*sph - z(2)*cph
emme  = sqrt(1._dk - z(1)**2_ik - z(2)**2_ik)

cnu = (cph - z(1) + (zeta*z(2))/(emme + 1._dk))/rho
snu = (sph - z(2) - (zeta*z(1))/(emme + 1._dk))/rho

! ==============================================================================
! 02. COMPUTE POSITION IN INERTIAL FRAME AND POTENTIAL
! ==============================================================================

! Intermediate frame unit vectors
x_vec = 2._dk*[ .5_dk - z(5)**2_ik - z(6)**2_ik,  &
           &  z(4)*z(5) + z(6)*z(7),  &
           &  z(4)*z(6) - z(5)*z(7)]
y_vec = 2._dk*[  z(4)*z(5) - z(6)*z(7),  &
           & .5_dk - z(4)**2_ik - z(6)**2_ik,  &
           &  z(5)*z(6) + z(4)*z(7)]

! Position in inertial frame
rV = rmag*(x_vec*cnu + y_vec*snu)

! Potential
Upot = 0._dk; dUdt = 0._dk; dUdr = 0._dk

!! Velocity in inertial frame
!i_vec = x_vec*cnu + y_vec*snu
!j_vec = y_vec*cnu - x_vec*snu
!v_rad = zeta/(sqrt(z(3))*rho)
!v_tan = sqrt((1. - z(1)**2 - z(2)**2)/(z(3)*rho**2) - 2.*Upot)
!vV = v_rad*i_vec + v_tan*j_vec

! ==============================================================================
! 03. COMPUTE PERTURBING ACCELERATIONS
! ==============================================================================

p = 0._dk; f = 0._dk

! Get time
if ( flag_time_EDr == 0_ik ) then
    ! Physical time
    t = z(8)
elseif  ( flag_time_EDr == 1_ik ) then
    ! Constant time element
    t = z(8) - z(3)**1.5_dk*(zeta - phi)
elseif  ( flag_time_EDr == 2_ik ) then
    ! Linear time element
    t = z(8) - z(3)**1.5_dk*zeta
end if

lon_Earth = wEarth*(t/TU)
if (inSoI) then
    ! Third body: Sun
    r2 = -(smaEarth/DU)*[cos(lon_Earth),sin(lon_Earth),0._dk]
    p  = DACC3B_ND(rV,r2,muEarth,muSun)
else
    ! Third body: Earth
    r2 = (smaEarth/DU)*[cos(lon_Earth),sin(lon_Earth),0._dk]
    p  = DACC3B_ND(rV,r2,muSun,muEarth)
end if

! Inertial -> Orbital frame
p    = DINERT2ORB_EDROMO(p,z,cnu,snu)
!dUdr = DINERT2ORB_EDROMO(dUdr,z,cnu,snu)

f = p + dUdr

! ==============================================================================
! 04. COMPUTE AUXILIARY QUANTITIES (2)
! ==============================================================================

enne = sqrt(emme**2_ik - 2._dk*z(3)*rho**2_ik*Upot)

aux0 = (p(1)*zeta + p(2)*enne + dUdt*sqrt(z(3))*rho)
zdot(3) = 2._dk*z(3)**3_ik*aux0
L3 = zdot(3)/(2._dk*z(3))

aux1 = ((2._dk*Upot - f(1)*rmag)*(2._dk - rho + emme)*rmag)/(emme*(emme+1._dk))
aux2 = (L3*zeta*(rho-emme))/(emme*(emme+1._dk))
wz   = (enne - emme)/rho + aux1 + aux2

! ==============================================================================
! 05. COMPUTE RIGHT-HAND SIDE
! ==============================================================================

! In-plane
aux3 = (f(1)*rmag - 2._dk*Upot)*rmag

zdot(1) =  aux3*sph + L3*((1._dk+rho)*cph - z(1))
zdot(2) = -aux3*cph + L3*((1._dk+rho)*sph - z(2))

! Out-of-plane
aux4 = (f(3)*rmag**2_ik)/(2._dk*enne)

zdot(4) =  aux4*(z(7)*cnu - z(6)*snu) + 0.5_dk*wz*z(5)
zdot(5) =  aux4*(z(6)*cnu + z(7)*snu) - 0.5_dk*wz*z(4)
zdot(6) = -aux4*(z(5)*cnu - z(4)*snu) + 0.5_dk*wz*z(7)
zdot(7) = -aux4*(z(4)*cnu + z(5)*snu) - 0.5_dk*wz*z(6)

! Time / Time Element
if (flag_time_EDr == 0_ik) then
    zdot(8) = sqrt(z(3))*rmag
else if (flag_time_EDr == 1_ik) then   ! Constant Time Element
    zdot(8) = z(3)**1.5_dk * ( aux3 + (zeta - 1.5_dk*phi)*zdot(3)/z(3) )
else if (flag_time_EDr == 2_ik) then   ! Linear Time Element
    zdot(8) = z(3)**1.5_dk * ( 1._dk + aux3 + 2._dk*L3*zeta )
end if

end subroutine DEDROMO_RHS

end module EQS_EDROMO

max_stars_repo_name: dkamato/NAPLES

max_issues_repo_name: dkamato/NAPLES

max_forks_repo_name: dkamato/NAPLES


text:
      subroutine genif(nperms,p,wt,msq)
c----initial-final subtraction.
      implicit none
      include 'constants.f'
      include 'qcdcouple.f'
      include 'xmin.f'
      include 'debug.f'
      logical justjac

      integer i1,i2,i3,i4,i5,j,k,nperms
      double precision p(mxpart,4),u,z,dot,q(mxpart,4),
     . msq(-nf:nf,-nf:nf),omxmin
      double precision facq,si1i3,si1i4,si3i4,x,omx,wt0,
     . wt5_4,wt4,wt,jacbit
      parameter(wt0=one/eight/pisq)
      common/justjac/justjac
      integer,parameter:: j1(4)=(/1,1,2,2/)
      integer,parameter:: j2(4)=(/2,2,1,1/)
      integer,parameter:: j3(4)=(/3,3,3,3/)
      integer,parameter:: j4(4)=(/4,5,4,5/)
      integer,parameter:: j5(4)=(/5,4,5,4/)

      i1=j1(nperms)
      i2=j2(nperms)
      i3=j3(nperms)
      i4=j4(nperms)
      i5=j5(nperms)

c first of calculate the variable's which one started with
c---nb all incoming
      si1i3=two*dot(p,i1,i3)
      si1i4=two*dot(p,i1,i4)
      si3i4=two*dot(p,i3,i4)
c----note momenta all incoming
      x=one+si3i4/(si1i3+si1i4)
      u=si1i3/(si1i3+si1i4)
      z=one-u
      omx=one-x
      omxmin=one-xmin
c -xx(i1)
c---at this stage the p are momenta including radiation

      do j=1,4
      q(i1,j)=x*p(i1,j)
      q(i2,j)=p(i2,j)
      q(i3,j)=p(i3,j)
      q(i4,j)=p(i4,j)+p(i3,j)+omx*p(i1,j)
      q(i5,j)=p(i5,j)
      q(6,j)=p(6,j)
      q(7,j)=p(7,j)
      enddo

      jacbit=four*sqrt(omx)/(half/sqrt(z)+half/sqrt(u))
      wt5_4=-wt0*dot(q,i1,i4)/x**2*omxmin*jacbit
      call wt4gen(q,wt4)

c---calculate total weight
      wt=wt5_4*wt4

      if (debug) then
      write(6,*) 'x in genif',x
      write(6,*) 'omxmin in genif',omxmin
      write(6,*) 'i1 in genif',i1
      write(6,*) 'i2 in genif',i2
      write(6,*) 'i3 in genif',i3
      write(6,*) 'i4 in genif',i4
      write(6,*) 'i5 in genif',i5
      write(6,*) 'wt5_4 in genif',wt5_4
      write(6,*) 'wt4 in genif',wt4
      write(6,*) 'wt in genif',wt
      endif

      if (justjac) return

      call qqb_wbb(q,msq)

c----initial wrt to final
      facq=-gsq/(x*si1i3)*(two/(one-x+u)-one-x)
     &     -gsq/(x*si3i4)*(two/(one-x+u)-one-z)

      do j=-nf,nf
      do k=-nf,nf
      msq(j,k)=0d0
      if     ((j .gt. 0) .and. (k .lt. 0)) then
      msq(j,k)=facq*msq(j,k)
      elseif ((j .lt. 0) .and. (k .gt. 0)) then
      msq(j,k)=facq*msq(j,k)
      endif

      enddo
      enddo

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
Originally founded at UC Berkeley, CampusCred is a Groupon knockoff that caters to only those with a ucdavis.edu email addresses.  NonUCD staff Townies are out of luck.  Their focus is to introduce students to local businesses and eateries within their community.  For those not familiar with the Groupon business model, they basically sell coupons to local businesses at a discount of 50% or more.  For example, their first deal was a $20 credit for KetMoRee that you could buy for $10.  Each deal is available for a few days before being replaced with a new one.  Check the fine print section beneath each offer  there may be limitations.

CampusCred uses the tagline only for students by students! Although it is run and managed locally by Davis students (juniors Sherry Wang and Users/NSheel Nikhil Sheel are mentioned as local directors in an http://theaggie.org/article/2011/01/18/newcompanyplanstohelpstudentssave Aggie article), it is owned by nonstudents.

CampusCred has nothing to do with ASUCD aside from paying for the right to use their logo.  For a while the logo was displayed prominently on the webpage, as though the site were run by ASUCD.  At best it was confusing; and it was probably deliberately shadiness factor deceptive.  The process by which they first got permission to use the logo led them down a bumpy road, as well.  They originally signed an agreement with Mark Champagne who was unaware he had had no authority to license the logo. When it was pointed out that Mark lacked the proper authority to license the logo, they received the proper license from the senate (see File(sr07.w11.pdf)).  Meanwhile, Occupy Davis is protesting corporate money in politics.  Hmmm...

When they first started up, they also sent out an unsolicited spam email to many ucdavis.edu email addresses.  Presumably, they got the mailing list from ASUCD, which is pretty shady in itself.  The list seemed to only include undergrads, as some faculty, alumni, and graduate students have reported not getting the email.  They also distributed fliers in many classrooms when they were first getting started.

Those associated with the organization constantly try to whitewash this page as well.

 Media
  http://theaggie.org/article/2011/03/14/aggiesenatebriefs
  http://theaggie.org/article/2011/01/18/newcompanyplanstohelpstudentssave New Company Plans to Help Students Save, California Aggie, 20110118
  http://blog.baremagazine.org/2010/09/campuscred.html Bare Magazine, story about the Berkeley group.
  http://theaggie.org/article/2011/02/10/asucdviolatesitsownbylawwithcampuscredagreement ASUCD violates its own bylaw with CampusCred agreement

20110129 15:31:42 nbsp As I posted on LeoThoms page, neither Jack, nor Joey signed an agreement, plus the Senate has to approve one anyways. Users/EliYani
  Hey Eli! Thanks for your comment. I apologize if I was not clear in our communications. CampusCred and ASUCD, more specifically, Mr. Mark Champagne, have a signed working agreement allowing CampusCred to use ASUCDs logo. Again, if you would like to view the plain English contact, please send me your availability and I would be more than happy to accommodate your schedule for a facetoface meeting. Thank you.

 I still think without formal Senate approval, you cant use the logo, even if Mark signed the contract.  Im not even sure if ASUCD as a notforprofit can endorse a forprofit business.

 Regardless, even if you think it, there is an agreement signed and if any further questions can be answered by the written agreement signed by CampusCred and Mr. Mark Champagne representing ASUCD.

20110129 19:39:37 nbsp So its like Groupon but for colleges? Sounds like a good idea. I wish this had been around when I was attending.  Users/ARWENNHOLD

20110131 19:04:35 nbsp Was not a fan when UCD made a deal with US Bank and I tend to avoid most businesses that try to get through to students using some weird contracts with the school (Sodexo), so a business that gets started in Davis with a potential scandal with one of my favorite groups on campus (ASUCD)? Sounds awesome (yes sarcasm).  Users/hankim

20110201 19:51:26 nbsp Some feedback on the CampusCred website: the placement of UC Davis and the ASUCD logo up in the header, without any explanation or anything, makes it look as though its a program run by the ASUCD.  If its unintentional, its simply a matter of confusion.  If its intentional, its disingenuous.  I know a couple of people have voiced this concern here on the Wiki, and it came up with another person in an offwiki conversation recently.  To the several representatives whove contributed to this page: itd be great if you could take a closer look at the logos and figure out a way to use them that doesnt give readers a false impression. Users/TomGarberson

20110221 21:46:28 nbsp We have removed the ASUCD logo pending the passing of the bill despite our 90day grace period. We are not trying to misrepresent anything and all the confusion will be settled until it is passed. Thanks. Users/Nsheel
 Now are you going to explain your attempt to spam the student body?  Users/WilliamLewis

20110425 23:10:12 nbsp I remember walking into one of my classes in Olson Hall last quarter and seeing fliers for this taped to the back of every other chair. Although theyre not the only group spamming stuff on campus, I find it curiously shady. Kind of like those students who show up on the first day of classes and make a short presentation about that textbook website...  Users/CecilioPadilla
  I hate those people so much. I want to go around erasing every chalkboard just to spite those asshats. Users/hankim

20110429 14:09:11 nbsp They seem to get butthurt easily over negative comments about them. Users/hankim

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C  PROGRAM CONTOUR

	INCLUDE 'VICMAIN_FOR'
	SUBROUTINE MAIN44
	IMPLICIT NONE

C  PROGRAM CONTOUR
C  PURPOSE ---
C
C	Create a graphics file of vectors representing contours
C	derived from an input "elevation" image.
C
C  INPUT ---
C
C	DIM	Dimension of graphics file (2D or 3D)
C	ZSTART	Starting elevation for contour generation
C	ZEND	Ending elevation for contour generation
C	CONTINT	Contour interval
C	SMOOTH	Smoothing level
C	ZERO	Switch to INCLude or OMIT contouring next to 0
C		valued pixels
C
C  OUTPUT ---
C	A 2D or 3D graphics file is generated
C
C  RESTRICTIONS ---
C
C
C  SUBROUTINES CALLED ---
C	XVCLOSE		Close a file from within VICAR/TAE
C	XVOPEN		Open a file from within VICAR/TAE
C	XVP		Parameter acquisition routine
C	XVSIZE		Get window size of incomming image
C	XVUNIT		Get the unit number
C
C  COMMENTS ---
C
C  HISTORY:
C	1-95   MAC   CRI   MSTP S/W CONVERSION (VICAR PORTING)
C
C
C  MODE DECLARATIONS ---
	INTEGER LINE, SAMP, COUNT, SL, SS, NL, NS, NUMPTS
	INTEGER DIM, STATUS, UNITIN, I, NLI, NSI
	INTEGER LL(10000), BF(10000), CONTINV, ZMIN, ZMAX
	CHARACTER*4 ZERO
	LOGICAL FRINGE

C  COMMON STATEMENTS ---
C	None
C
C  LOCAL VARIABLE DESCRIPTIONS ---
C	None
C
C-----------*** BEGINNING OF EXECUTABLE CODE ***-----------------

C		+=================+
C		| INITIALIZATIONS |
C		+=================+

        CALL IFMESSAGE('CONTOUR version 02-JAN-95')
	CALL XVUNIT (UNITIN, 'INP', 1, STATUS,' ')
	CALL XVOPEN(UNITIN, STATUS,'OPEN_ACT','SA','IO_ACT','SA',
     *            'U_FORMAT','FULL',' ')
	CALL XVSIZE (SL,SS,NL,NS,NLI,NSI)

C		Extract the contouring parameters

	CALL XVP ('DIM',DIM,COUNT)
	CALL XVP ('ZSTART',ZMIN,COUNT)
	CALL XVP ('ZEND',ZMAX, COUNT)
	CALL XVP ('CONTINT',CONTINV,COUNT)
	CALL XVP ('NUMPTS',NUMPTS,COUNT)
	CALL XVP ('ZERO',ZERO,COUNT)
	IF (ZERO(1:4).EQ.'OMIT') THEN
	    FRINGE = .FALSE.
	ELSE
	    FRINGE = .TRUE.
	END IF

C		Open the output graphics file

	CALL WRGR(1,1,DIM) ! Open a graphics output file

C	   ************ Begin EXECUTION *************

C------------------------ C O N T O U R ----------------------------

C	Program Contour

	CALL XVREAD (UNITIN,LL,STATUS,'LINE',SL,
     *		     'SAMP',SS,'NSAMPS',NS,' ')

	DO LINE = SL+1, NL+SL-1
	    CALL XVREAD (UNITIN,BF,STATUS,'LINE',LINE,'SAMP',SS,
     &                   'NSAMPS',NS,' ')
	    DO SAMP = 1, NS-1
		IF (FRINGE) THEN
		    CALL CNTRGRID (LL(SAMP),LL(SAMP+1),
     *		    BF(SAMP),BF(SAMP+1),CONTINV,ZMIN,ZMAX,
     *		    NUMPTS,FLOAT(LINE-1),FLOAT(SAMP+SS-1))
		ELSE
		    IF (.NOT.((LL(SAMP).EQ.0).OR.(LL(SAMP+1).EQ.0).OR.
     *		    (BF(SAMP).EQ.0).OR.(BF(SAMP+1).EQ.0))) THEN
			CALL CNTRGRID (LL(SAMP),LL(SAMP+1),
     *			BF(SAMP),BF(SAMP+1),CONTINV,ZMIN,ZMAX,
     *			NUMPTS,FLOAT(LINE-1),FLOAT(SAMP+SS-1))
		    END IF
		END IF
	    END DO
	    DO I = 1, NS
		LL(I)=BF(I)
	    END DO
	END DO

	CALL XVCLOSE(UNITIN,STATUS,' ')
	CALL CLGR(1)

	RETURN
	END

C ------------- C O N T O U R G R I D (CNTRGRID) --------------------

	SUBROUTINE CNTRGRID (F00,F01,F10,F11,CONTINV,ZMIN,
     *		ZMAX,NUMPTS,XOFF,YOFF)

	IMPLICIT NONE

	INTEGER F00, F01, F10, F11, NUMPTS
	INTEGER NUMPAIRS, BRANCH, POINT, I, CONTINV, ZMIN, ZMAX
	REAL XOFF, YOFF, K, X0, Y0, R, X, Y, TEMP, PHI, PHI1, PHI2
	REAL A, B, C, D, L, XEDGE(4), YEDGE(4)

	A = F10 - F00
	B = F11 - F10 - F01 + F00
	C = F01 - F00
	D = F00

	L = ZMIN + 0.5

	DO WHILE (L.LE.ZMAX)

	    CALL GETECRDS (NUMPAIRS, XEDGE, YEDGE, F00, F01, F10, F11, L)

	    IF (NUMPAIRS .GT. 0) THEN
		IF (B .EQ. 0) THEN
		    CALL PUTGR (1,XEDGE(1)+XOFF,YEDGE(1)+YOFF,L)
		    CALL PUTGR (1,XEDGE(2)+XOFF,YEDGE(2)+YOFF,L)
		    CALL PUTGR (1,0,0,0)
		ELSE
		    X0 = - C/B
		    Y0 = - A/B
		    K = 2*(A*C + B*(L-D)) / (B*B)
		    IF ((NUMPAIRS .EQ. 2) .AND. (K .GT. 0)) THEN
			TEMP = YEDGE(2)
			YEDGE(2) = YEDGE(4)
			YEDGE(4) = TEMP
			TEMP = XEDGE(2)
			XEDGE(2) = XEDGE(4)
			XEDGE(4) = TEMP
		    END IF

		    DO BRANCH = 1, NUMPAIRS
			I = 2*BRANCH - 1
			PHI1 = ATAN2 (YEDGE(I)-Y0, XEDGE(I)-X0)
			PHI2 = ATAN2 (YEDGE(I+1)-Y0, XEDGE(I+1)-X0)

			DO POINT = 0, NUMPTS-1
			    PHI = ((PHI2-PHI1)/(NUMPTS-1))*POINT + PHI1
			    IF (PHI .EQ. 0) PHI = 1.0E-8
			    R = SQRT(ABS( K/SIN(2*PHI) ))
			    X = R*COS(PHI) + X0
			    Y = R*SIN(PHI) + Y0
			    CALL PUTGR (1,X+XOFF,Y+YOFF,L)
			END DO
			CALL PUTGR (1,0.0,0.0,0.0)
		    END DO
	    	END IF
	    END IF
	    L = L + CONTINV
	END DO
	RETURN
	END

C ------------ G E T E D G E C O O R D S (GETECRDS) ------------------

	SUBROUTINE GETECRDS (NUMPAIRS, X, Y, F00, F01, F10, F11, L)

	IMPLICIT NONE

	INTEGER NUMPAIRS, J, F00, F01, F10, F11
	REAL X(4), Y(4), EDGE, L

	J = 0
	EDGE = (L-F00) / (F10-F00+0.0001)
	IF ((EDGE .GE. 0) .AND. (EDGE .LE. 1)) THEN
	    J = J + 1
	    X(J) = EDGE
	    Y(J) = 0
	END IF
	EDGE = (L-F10) / (F11-F10+0.0001)
	IF ((EDGE .GE. 0) .AND. (EDGE .LE. 1)) THEN
	    J = J + 1
	    X(J) = 1
	    Y(J) = EDGE
	END IF
	EDGE = (L-F01) / (F11-F01+0.0001)
	IF ((EDGE .GE. 0) .AND. (EDGE .LE. 1)) THEN
	    J = J + 1
	    X(J) = EDGE
	    Y(J) = 1
	END IF
	EDGE = (L-F00) / (F01-F00+0.0001)
	IF ((EDGE .GE. 0) .AND. (EDGE .LE. 1)) THEN
	    J = J + 1
	    X(J) = 0
	    Y(J) = EDGE
	END IF
	NUMPAIRS = J/2

	RETURN
	END

C  ------------------ E N D   O F   S O U R C E ----------------------

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
*..............................................................

!-------------------------------------------------------------------------
!         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
!BOP
!
! !ROUTINE: ODS_GetI() --- Reads data from file to native integers
!
! !DESCRIPTION:
!    \label{ODS:GetI}
!     Reads the contents of an integer variable for a full synoptic
!     time from an opened ODS file.
!
!     Note: For a list of error codes, see Table~\ref{tab:errors}.
!
! !INTERFACE:
      subroutine ODS_GetI ( id, VarName, julian_day, syn_hour,
     .                          nval,    values,     ierr )

!
! !INPUT PARAMETERS:
      implicit        NONE
      integer         id              ! ODS file handle as returned
                                      !  from ODS_Create() or
                                      !  ODS_Open().
      character * (*) VarName         ! name of variable
      integer         julian_day      ! Julian Day.  Use the function
                                      !  ODS_Julian to obtain this
                                      !  number.
      integer         syn_hour        ! hour of synoptic time since
                                      !  0:00 GMT (e.g., 0, 6, 12, 18)

! !INPUT/OUTPUT PARAMETERS:
      integer         nval            ! on input: maximum number of
                                      !  values (usually determined
                                      !  by the space allocated for
                                      !  storage)
                                      ! on output: the number of
                                      !  values obtained from the
                                      !  file
                                      ! note: If the input value is
                                      !  smaller than the number
                                      !  of values to be read, then
                                      !  the routine exits with a
                                      !  system error message and
                                      !  code of ODS_DimErr ( = -3 ).
                                      !  Also, nval is set to the
                                      !  minimum value required in
                                      !  order for the routine to
                                      !  execute successfully.

! !OUTPUT PARAMETERS:
      integer         values ( nval ) ! Values read from file
      integer         ierr            ! Error code. If non-zero, an
                                      !  error has occurred. For a list
                                      !  of possible values, see the
                                      !  description section of this
                                      !  prologue.
!
! !SEE ALSO:
!     ODS_PutC()    write a list of character  strings to   an ODS file
!     ODS_PutI()    write a list of integer    values  to   an ODS file
!     ODS_PutR()    write a list of real       values  to   an ODS file
!     ODS_GetC()    get   a list of characters strings from an ODS file
!     ODS_GetR()    get   a list of real       values  from an ODS file
!     ODS_NGet()    get the number of values for a given julian day
!                    and synoptic hour
!     ODS_Julian()  convert the "calendar" date to the Julian day
!     ods_hdf.h     include file defining internal constants
!                    and global variables, and setting up data
!                    structures
!
! !REVISION HISTORY:
!     02Oct1995   da Silva    Specification.
!     16May1996   C. Redder   Original code
!     13Apr1998   C. Redder   Updated documentation to reflect changes
!                             made in other routines to create ODS
!                             version 2.00
!     19Nov1999   C. Redder   Added a latex label in and moved the
!                             subroutine statement into the prologue.
!                             Modified the comments for the return status
!                             code.
!     16Feb2000   Todling     Rename stdio.h to ods_stdio.h
!     30Nov2001   C. Redder   Fixed bug that prevents last batch from
!                             being processed twice when nval = 0
!     16Sep2002   C. Redder   Updated prologue and added error handling
!                             for call to NetCDF routine, NCDINQ
!EOP
!-------------------------------------------------------------------------

      include 'netcdf.inc'
      include 'ods_hdf.h'
      include 'ods_stdio.h'

*     Functions referenced
*     --------------------
      character * ( 10 )
     .         ODS_VarType ! recommended FORTRAN variable type for data

*     Arguments for NetCDF library routines
*     -------------------------------------
      integer  nc_id       ! file id
      integer  varid       ! variable id
      integer  start  ( 2 )! starting file index for hyperslab
      integer  count  ( 2 )! length of edges in hyperslab
      character * ( MAXNCNAM )
     .         Var_Name    ! variable name ( returned from NetCDF
                           !   routine NCVINQ )
      integer  Var_Type    ! variable type
      integer  NVDims      ! number of dimensions
      integer  DimIDs ( 2 )! dimension ids
      integer  DimID       ! dimension id
      integer  NVAtts      ! number of attributes
      character * ( MAXNCNAM )
     .         DimName     ! variable name ( returned from NetCDF
                           !   routine NCDINQ )

*     Indices defining the set of observation data
*     corresponding to a synoptic time in the ODS file
*     ------------------------------------------------
      integer  synbeg      ! starting index for the data block
      integer  synlen      ! number of observation reports in
                           !   the data block

*     Parameters for defining the observation
*     data batches to be written to the ODS file
*     ------------------------------------------
      integer  batch_len   ! number of values in each batch
      integer  first_batch ! index for the first batch
      integer  last_batch  ! index for the last batch
      integer  first_ival  ! data value index of the first
                           !   value in the first batch
      integer  first_nval  ! number of values in the first batch
      integer  last_nval   ! number of values in the last batch
      integer  ibatch      ! batch number index for do loop

*     Other variables
*     ---------------
      integer  ival        ! index value for the input array values
      integer  nval_max    ! maximum number of values
      character * ( 10 )
     .         VarType     ! recommended FORTRAN variable type for
                           !   data

*     Set ierr code to valid
*     ----------------------
      ierr    = NCNoErr

*     Check to determine if the file handle id is valid
*     -------------------------------------------------
      if ( id            .lt. 1      .or.
     .     id            .gt. id_max .or.
     .     IOMode ( id ) .eq. CLOSED ) then
         write ( stderr, 901 )
         ierr = NCEBadID
         return
      end if

*     Check to determine if variable is appropriate for this routine
*     --------------------------------------------------------------
      VarType    = ODS_VarType ( VarName )
      if ( VarType ( :7  ) .ne. 'integer'     .and.
     .     VarType ( :10 ) .ne. 'indefinite' ) then
         write ( stderr, * )
         write ( stderr, 902 ) VarName
         ierr = ODS_BadInt
         return

      end if

*     Get pointer data from tables stored in common
*     ---------------------------------------------
      call ODS_GetP    ( id,        'read',
     .                   julian_day, syn_hour,
     .                   synbeg,     synlen,  ierr )
      if ( ierr .ne. NCNoErr ) return

*     Get NetCDF file and variable ids
*     --------------------------------
      nc_id = ncid ( id )
      varid = NCVID ( nc_id, VarName, ierr )
      if ( ierr .ne. NCNoErr ) return

*     Set nval and nval_max
*     ---------------------
      nval_max = nval
      nval     = synlen

*     Return with an error message if the number of
*     observations to be read exceeds the user-defined maximum
*     --------------------------------------------------------
      if ( nval .gt. nval_max ) then
         write ( stderr, * )
         write ( stderr, 903 ) nval_max, nval
         ierr = ODS_DimErr
         return

      end if

*     Get the length of each batch of observation data
*     in an ODS file
*     ------------------------------------------------
      call NCVINQ ( nc_id,    varid,
     .              Var_Name, Var_Type,
     .              NVDims,   DimIDs,
     .              NVatts,   ierr )
      if ( ierr .ne. NCNoErr ) return
      DimID       = DimIDs ( 1 )
      call NCDINQ ( nc_id,    DimID,
     .              DimName,  batch_len, ierr )
      if ( ierr .ne. NCNoErr ) return

*     The remaining code uses the ODS IO routines to read
*     the data in batches.  Each batch contains a subset
*     of the values to be read.  The maximum size of
*     each subset is batchlen as defined in the header
*     file, ods_hdf.h.  The values are being processed in
*     this manner in order to prevent the execution of
*     the NetCDF library routines from being
*     prohibitively slow.
*     ---------------------------------------------------

*     Determine the indices of the first and last batches
*     of data to be read from the ODS file
*     ---------------------------------------------------
      first_batch = ( synbeg - 1 ) / batch_len + 1
      last_batch  = ( synbeg + synlen - 2 ) / batch_len + 1

*     Read the first batch of values
*     ------------------------------
      first_ival  =   mod ( synbeg - 1, batch_len ) + 1
      first_nval  =   min ( synlen, batch_len - first_ival + 1 )
      start ( 1 ) =   first_ival
      start ( 2 ) =   first_batch
      count ( 1 ) =   first_nval
      count ( 2 ) =   1
      call ODS_NCVGTI ( nc_id, varid, start, count, values, ierr )
      if ( ierr .ne. NCNoErr ) return

*     Read data not in the first and last batch of values
*     ---------------------------------------------------
      ival = first_nval + 1
      do 10, ibatch = first_batch + 1, last_batch - 1
         start ( 1 ) = 1
         start ( 2 ) = ibatch
         count ( 1 ) = batch_len
         count ( 2 ) = 1
         call ODS_NCVGTI ( nc_id, varid, start, count,
     .                     values ( ival ), ierr )
         if ( ierr .ne. NCNoErr ) return
         ival        = ival + batch_len
 10   continue

*     Read the last batch of values, if necessary
*     -------------------------------------------
      start ( 1 ) =   1
      start ( 2 ) =   last_batch
      count ( 1 ) =   mod ( synbeg + synlen - 2, batch_len ) + 1
      count ( 2 ) =   1
      if ( last_batch .gt. first_batch )
     .   call ODS_NCVGTI ( nc_id, varid, start, count,
     .                     values ( ival ), ierr )
      if ( ierr .ne. NCNoErr ) return

      return
*     ------

 901  format ( /, ' ODS_GetI: File handle id number does not ',
     .         /, '           correspond to an opened ODS file' )
 902  format ( /, ' ODS_GetI: Inappropriate variable for this ',
     .         /, '           routine.  Variable name is: ', a )
 903  format ( /, ' ODS_GetI: The number of values to be read ',
     .         /, '           exceeds the user defined maximum',
     .         /, '           of ', i10, '. Number of values ',
     .         /, '           to be read = ', i10,'.' )

      end

max_stars_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
      SUBROUTINE TG01WD( N, M, P, A, LDA, E, LDE, B, LDB, C, LDC,
     $                   Q, LDQ, Z, LDZ, ALPHAR, ALPHAI, BETA, DWORK,
     $                   LDWORK, INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To reduce the pair (A,E) to a real generalized Schur form
C     by using an orthogonal equivalence transformation
C     (A,E) <-- (Q'*A*Z,Q'*E*Z) and to apply the transformation
C     to the matrices B and C: B <-- Q'*B and C <-- C*Z.
C
C     ARGUMENTS
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the original state-space representation,
C             i.e., the order of the matrices A and E.  N >= 0.
C
C     M       (input) INTEGER
C             The number of system inputs, or of columns of B.  M >= 0.
C
C     P       (input) INTEGER
C             The number of system outputs, or of rows of C.  P >= 0.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, the leading N-by-N part of this array must
C             contain the original state dynamics matrix A.
C             On exit, the leading N-by-N part of this array contains
C             the matrix Q' * A * Z in an upper quasi-triangular form.
C             The elements below the first subdiagonal are set to zero.
C
C     LDA     INTEGER
C             The leading dimension of array A.  LDA >= MAX(1,N).
C
C     E       (input/output) DOUBLE PRECISION array, dimension (LDE,N)
C             On entry, the leading N-by-N part of this array must
C             contain the original descriptor matrix E.
C             On exit, the leading N-by-N part of this array contains
C             the matrix Q' * E * Z in an upper triangular form.
C             The elements below the diagonal are set to zero.
C
C     LDE     INTEGER
C             The leading dimension of array E.  LDE >= MAX(1,N).
C
C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
C             On entry, the leading N-by-M part of this array must
C             contain the input matrix B.
C             On exit, the leading N-by-M part of this array contains
C             the transformed input matrix Q' * B.
C
C     LDB     INTEGER
C             The leading dimension of array B.  LDB >= MAX(1,N).
C
C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
C             On entry, the leading P-by-N part of this array must
C             contain the output matrix C.
C             On exit, the leading P-by-N part of this array contains
C             the transformed output matrix C * Z.
C
C     LDC     INTEGER
C             The leading dimension of array C.  LDC >= MAX(1,P).
C
C     Q       (output) DOUBLE PRECISION array, dimension (LDQ,N)
C             The leading N-by-N part of this array contains the left
C             orthogonal transformation matrix used to reduce (A,E) to
C             the real generalized Schur form.
C             The columns of Q are the left generalized Schur vectors
C             of the pair (A,E).
C
C     LDQ     INTEGER
C             The leading dimension of array Q.  LDQ >= max(1,N).
C
C     Z       (output) DOUBLE PRECISION array, dimension (LDZ,N)
C             The leading N-by-N part of this array contains the right
C             orthogonal transformation matrix used to reduce (A,E) to
C             the real generalized Schur form.
C             The columns of Z are the right generalized Schur vectors
C             of the pair (A,E).
C
C     LDZ     INTEGER
C             The leading dimension of array Z.  LDZ >= max(1,N).
C
C     ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
C     ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
C     BETA    (output) DOUBLE PRECISION array, dimension (N)
C             On exit, if INFO = 0, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j),
C             j=1,...,N, will be the generalized eigenvalues.
C             ALPHAR(j) + ALPHAI(j)*i, and BETA(j), j=1,...,N, are the
C             diagonals of the complex Schur form that would result if
C             the 2-by-2 diagonal blocks of the real Schur form of
C             (A,E) were further reduced to triangular form using
C             2-by-2 complex unitary transformations.
C             If ALPHAI(j) is zero, then the j-th eigenvalue is real;
C             if positive, then the j-th and (j+1)-st eigenvalues are a
C             complex conjugate pair, with ALPHAI(j+1) negative.
C
C     Workspace
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, DWORK(1) returns the optimal value
C             of LDWORK.
C
C     LDWORK  INTEGER
C             The dimension of working array DWORK.  LDWORK >= 8*N+16.
C             For optimum performance LDWORK should be larger.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value;
C             > 0:  if INFO = i, the QZ algorithm failed to compute
C                   the generalized real Schur form; elements i+1:N of
C                   ALPHAR, ALPHAI, and BETA should be correct.
C
C     METHOD
C
C     The pair (A,E) is reduced to a real generalized Schur form using
C     an orthogonal equivalence transformation (A,E) <-- (Q'*A*Z,Q'*E*Z)
C     and the transformation is applied to the matrices B and C:
C     B <-- Q'*B and C <-- C*Z.
C
C     NUMERICAL ASPECTS
C                                     3
C     The algorithm requires about 25N  floating point operations.
C
C     CONTRIBUTOR
C
C     A. Varga, German Aerospace Center,
C     DLR Oberpfaffenhofen, July 2000.
C
C     REVISIONS
C
C     V. Sima, Research Institute for Informatics, Bucharest, Mar. 2001.
C
C     KEYWORDS
C
C     Orthogonal transformation, generalized real Schur form, similarity
C     transformation.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      INTEGER           INFO, LDA, LDB, LDC, LDE, LDQ, LDWORK, LDZ,
     $                  M, N, P
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*),
     $                  BETA(*),  C(LDC,*),  DWORK(*),  E(LDE,*),
     $                  Q(LDQ,*), Z(LDZ,*)
C     .. Local Scalars ..
      LOGICAL           BLAS3, BLOCK
      INTEGER           BL, CHUNK, I, J, MAXWRK, SDIM
C     .. Local Arrays ..
      LOGICAL           BWORK(1)
C     .. External Functions ..
      LOGICAL           LSAME, DELCTG
      EXTERNAL          LSAME, DELCTG
C     .. External Subroutines ..
      EXTERNAL          DCOPY, DGEMM, DGEMV, DGGES, DLACPY, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         DBLE, INT, MAX, MIN
C
C     .. Executable Statements ..
C
      INFO = 0
C
C     Check the scalar input parameters.
C
      IF( N.LT.0 ) THEN
         INFO = -1
      ELSE IF( M.LT.0 ) THEN
         INFO = -2
      ELSE IF( P.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -9
      ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
         INFO = -11
      ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
         INFO = -13
      ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN
         INFO = -15
      ELSE IF( LDWORK.LT.8*N+16 ) THEN
         INFO = -20
      END IF
C
      IF( INFO.NE.0 ) THEN
C
C        Error return.
C
         CALL XERBLA( 'TG01WD', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF( N.EQ.0 ) THEN
         DWORK(1) = ONE
         RETURN
      END IF
C
C     Reduce (A,E) to real generalized Schur form using an orthogonal
C     equivalence transformation (A,E) <-- (Q'*A*Z,Q'*E*Z), accumulate
C     the transformations in Q and Z, and compute the generalized
C     eigenvalues of the pair (A,E) in (ALPHAR, ALPHAI, BETA).
C
C     Workspace:  need   8*N+16;
C                 prefer larger.
C
      CALL DGGES( 'Vectors', 'Vectors', 'Not ordered', DELCTG, N,
     $            A, LDA, E, LDE, SDIM, ALPHAR, ALPHAI, BETA, Q, LDQ,
     $            Z, LDZ, DWORK, LDWORK, BWORK, INFO )
      IF( INFO.NE.0 )
     $   RETURN
      MAXWRK = INT( DWORK(1) )
C
C     Apply the transformation: B <-- Q'*B. Use BLAS 3, if enough space.
C
      CHUNK = LDWORK / N
      BLOCK = M.GT.1
      BLAS3 = CHUNK.GE.M .AND. BLOCK
C
      IF( BLAS3 ) THEN
C
C        Enough workspace for a fast BLAS 3 algorithm.
C
         CALL DLACPY( 'Full', N, M, B, LDB, DWORK, N )
         CALL DGEMM( 'Transpose', 'No transpose', N, M, N, ONE, Q, LDQ,
     $               DWORK, N, ZERO, B, LDB )
C
      ELSE IF ( BLOCK ) THEN
C
C        Use as many columns of B as possible.
C
         DO 10 J = 1, M, CHUNK
            BL = MIN( M-J+1, CHUNK )
            CALL DLACPY( 'Full', N, BL, B(1,J), LDB, DWORK, N )
            CALL DGEMM( 'Transpose', 'NoTranspose', N, BL, N, ONE, Q,
     $                  LDQ, DWORK, N, ZERO, B(1,J), LDB )
   10    CONTINUE
C
      ELSE
C
C        Use a BLAS 2 algorithm. Here, M <= 1.
C
         IF ( M.GT.0 ) THEN
            CALL DCOPY( N, B, 1, DWORK, 1 )
            CALL DGEMV( 'Transpose', N, N, ONE, Q, LDQ, DWORK, 1, ZERO,
     $                  B, 1 )
         END IF
      END IF
      MAXWRK = MAX( MAXWRK, N*M )
C
C     Apply the transformation: C <-- C*Z.  Use BLAS 3, if enough space.
C
      BLOCK = P.GT.1
      BLAS3 = CHUNK.GE.P .AND. BLOCK
C
      IF ( BLAS3 ) THEN
         CALL DLACPY( 'Full', P, N, C, LDC, DWORK, P )
         CALL DGEMM( 'No transpose', 'No transpose', P, N, N, ONE,
     $               DWORK, P, Z, LDZ, ZERO, C, LDC )
C
      ELSE IF ( BLOCK ) THEN
C
C        Use as many rows of C as possible.
C
         DO 20 I = 1, P, CHUNK
            BL = MIN( P-I+1, CHUNK )
            CALL DLACPY( 'Full', BL, N, C(I,1), LDC, DWORK, BL )
            CALL DGEMM( 'NoTranspose', 'NoTranspose', BL, N, N, ONE,
     $                  DWORK, BL, Z, LDZ, ZERO, C(I,1), LDC )
   20    CONTINUE
C
      ELSE
C
C        Use a BLAS 2 algorithm. Here, P <= 1.
C
         IF ( P.GT.0 ) THEN
            CALL DCOPY( N, C, LDC, DWORK, 1 )
            CALL DGEMV( 'Transpose', N, N, ONE, Z, LDZ, DWORK, 1, ZERO,
     $                  C, LDC )
         END IF
C
      END IF
      MAXWRK = MAX( MAXWRK, P*N )
C
      DWORK(1) = DBLE( MAXWRK )
C
      RETURN
C *** Last line of TG01WD ***
      END

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


! { dg-do compile }
! { dg-options "-fopenmp -fcray-pointer" }
!
! PR fortran/43985

subroutine pete(A)
  real(8) :: A
  print *, 'pete got ',A
  if (A /= 3.0) STOP 1
end subroutine pete

       subroutine bob()
         implicit none
         real(8) peted
         pointer (ipeted, peted(*))
         integer(4) sz
         ipeted = malloc(5*8)
         peted(1:5) = [ 1.,2.,3.,4.,5.]
         sz = 3
!$omp parallel default(shared)
         call pete(peted(sz))
!$omp end parallel
         return
       end subroutine bob

call bob()
end


text:
C	$Id: wmlabs.f,v 1.13 2008-07-27 00:17:36 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE WMLABS(X,Y,SYMTYP)
C
C  Plot the symbol indicated by the string SYMTYP.  The symbol
C  is centered at (X,Y), where X and Y are world coordinates.
C  Legal values for SYMTYP are:
C
C         'HI'
C         'LOW'
C         'ARROW'
C         'DOT'
C         'CLOUD'
C         'ICE'
C         'IS' (Intermittent showers)
C         'IT' (Sunny, chance of T-storms)
C         'MC' (Mostly cloudy)
C         'MS' (Mostly sunny)
C         'RAIN'
C         'RS' (Rain and snow)
C         'SNOWFLAKES'
C         'SUN'
C         'THUNDERSTORM'
C         'VECTOR'
C         'WIND'
C
      CHARACTER*(*) SYMTYP
      CHARACTER*1 FCODE
      CHARACTER*5 CHRDIST
C
      include 'wmcomn.h'
C
C  Reference arrow.
C
      PARAMETER (IADIM=8,D2RAD=.017453293)
      DIMENSION XX(IADIM),YY(IADIM),ARROWX(IADIM),ARROWY(IADIM)
      REAL OLDWN(4),OLDVP(4),OCLIP(4)
      DIMENSION TMPX(IADIM),TMPY(IADIM)
      DATA XX/-1.000, -0.330, -0.360,  0.000,
     +        -0.360, -0.330, -1.000, -1.000  /
      DATA YY/-0.030, -0.030, -0.120,  0.000,
     +         0.120,  0.030,  0.030, -0.030  /
C
C  Array for vector draws.
C
      PARAMETER (IVSDM=5)
      DIMENSION VECTSX(IVSDM),VECTSY(IVSDM)
C
C  Save the current text function code and set it to a colon for
C  this subroutine.
C
      CALL PCGETC('FC',FCODE)
      CALL PCSETC('FC',':')
C
C  Save the current line and fill colors and set them to ICOLOR.
C
      CALL GQFAIS(IER,IFAISO)
      CALL GQFACI(IER,IFCLRO)
      CALL GQPLCI(IER,ILCLRO)
      CALL GSFAIS(1)
      CALL GSFACI(ICOLOR)
      CALL GSPLCI(ICOLOR)
C
C  Convert X and Y to NDC and work in NDC space.
C
      CALL WMW2NX(1,X,XNDC)
      CALL WMW2NY(1,Y,YNDC)
C
C  Save normalization transformation number 1 along with its
C  window and viewport and define a new transformation 1
C  that has as its window and viewport the current clip
C  rectangle.
C
      CALL GQCNTN(IER,NTRO)
      CALL GQCLIP(IER,ICLIP,OCLIP)
      CALL GQNT(1,IER,OLDWN,OLDVP)
C
      CALL GSWN(1,OCLIP(1),OCLIP(2),OCLIP(3),OCLIP(4))
      CALL GSVP(1,OCLIP(1),OCLIP(2),OCLIP(3),OCLIP(4))
      CALL GSELNT(1)
C
C  Turn on text extent computation for the remainder of the subroutine.
C
      CALL PCSETI('TE',1)
C
      SIZEL = WSIZES
C
      IF (SYMTYP(1:1).EQ.'H' .OR. SYMTYP(1:1).EQ.'h') THEN
C
C  Draw a high symbol.
C
        RADIUS = 1.10*SIZEL
        OFFSET = 0.15*RADIUS
        CALL NGGETI('CT',ICTYPO)
        CALL NGSETI('CT',0)
C
C  Shadow.
C
        IF (IHIGC1 .GE. 0) THEN
          CALL NGDOTS(XNDC-OFFSET,YNDC-OFFSET,1,2.*RADIUS,IHIGC1)
        ENDIF
C
C  Character background.
C
        IF (IHIGC2 .GE. 0) THEN
          CALL NGDOTS(XNDC,YNDC,1,2.*RADIUS,IHIGC2)
        ENDIF
        CALL NGSETI('CT',1)
C
C  Circumscribed circle.
C
        IF (IHIGC4 .GE. 0) THEN
          CALL NGDOTS(XNDC,YNDC,1,2.*RADIUS,IHIGC4)
        ENDIF
        CALL NGSETI('CT',ICTYPO)
        CALL PCGETI('CC - character color',ICLRO)
        CALL PCGETI('FN - font name',IFNTO)
        CALL PCSETI('CC - character color',IHIGC3)
        CALL PCSETI('FN - font name',22)
C
C  "H" character.
C
        IF (IHIGC3 .GE. 0) THEN
          CALL PLCHHQ(XNDC,YNDC,'H',SIZEL,0.,0.)
        ENDIF
        CALL PCSETI('CC - character color',ICLRO)
        CALL PCSETI('FN - font name',IFNTO)
      ELSE IF (SYMTYP(1:1).EQ.'L' .OR. SYMTYP(1:1).EQ.'l') THEN
C
C  Draw a low symbol.
C
        RADIUS = 1.0*SIZEL
        OFFSET = 0.20*RADIUS
        CALL NGGETI('CT',ICTYPO)
        CALL NGSETI('CT',0)
C
C  Shadow.
C
        IF (ILOWC1 .GE. 0) THEN
          CALL NGDOTS(XNDC,YNDC,1,2.*RADIUS+OFFSET,ILOWC1)
        ENDIF
C
C  Character background.
C
        IF (ILOWC2 .GE. 0) THEN
          CALL NGDOTS(XNDC,YNDC,1,2.*RADIUS,ILOWC2)
        ENDIF
        CALL NGSETI('CT',1)
C
C  Circumscribed circle.
C
        IF (ILOWC4 .GE. 0) THEN
          CALL NGDOTS(XNDC,YNDC,1,2.*RADIUS,ILOWC4)
        ENDIF
        CALL NGSETI('CT',ICTYPO)
        CALL PCGETI('CC - character color',ICLRO)
        CALL PCGETI('FN - font name',IFNTO)
        CALL PCSETI('CC - character color',ILOWC3)
        CALL PCSETI('FN - font name',22)
C
C  "L" character.
C
        IF (ILOWC3 .GE. 0) THEN
          CALL PLCHHQ (XNDC,YNDC,':F22:L',SIZEL,0.,0.)
        ENDIF
        CALL PCSETI('CC - CHARACTER COLOR',ICLRO)
        CALL PCSETI('FN - font name',IFNTO)
      ELSE IF (SYMTYP(1:1).EQ.'A' .OR. SYMTYP(1:1).EQ.'a') THEN
C
C  Draw an arrow of length ARWSIZ in direction ARWDIR with the
C  point of the arrow at (X,Y).
C
C  Scale the arrow.
C
        SS = ARWSIZ
        DO 10 I=1,IADIM
          IF (I.EQ.1 .OR. I.EQ.7 .OR. I.EQ.8) THEN
            ARROWX(I) = SS*ARWLEN*XX(I)
          ELSE
            ARROWX(I) = SS*XX(I)
          ENDIF
          ARROWY(I) = SS*YY(I)
   10   CONTINUE
C
C  Rotate the arrow.
C
        COSANG = COS(D2RAD*ARWDIR)
        SINANG = SIN(D2RAD*ARWDIR)
        DO 20 I=1,IADIM
          TMPX(I) = ARROWX(I)*COSANG-ARROWY(I)*SINANG
          TMPY(I) = ARROWX(I)*SINANG+ARROWY(I)*COSANG
   20   CONTINUE
        DO 30 I=1,IADIM
          ARROWX(I) = TMPX(I)
          ARROWY(I) = TMPY(I)
   30   CONTINUE
C
C  Translate the arrow.
C
        DO 40 I=1,IADIM
          ARROWX(I) = XNDC+ARROWX(I)
          ARROWY(I) = YNDC+ARROWY(I)
   40   CONTINUE
C
C  Draw a shadow if the arrow shadow color is non-negative.
C
        IF (IARSHC .GE. 0) THEN
          DO 50 I=1,IADIM
            TMPX(I) = ARROWX(I)-.0019*ARWSIZ/0.035
            TMPY(I) = ARROWY(I)-.0019*ARWSIZ/0.035
   50     CONTINUE
          CALL GSFACI(IARSHC)
          CALL GFA(IADIM,TMPX,TMPY)
        ENDIF
C
C  Draw the arrow.
C
        CALL GSFACI(IAROWC)
        CALL GFA(IADIM,ARROWX,ARROWY)
C
C  Draw an outline around the arrow if the outline color is non-negative.
C
        IF (IAROUC .GE. 0) THEN
          CALL GSPLCI(IAROUC)
          CALL GPL(IADIM,ARROWX,ARROWY)
        ENDIF
      ELSE IF (SYMTYP(1:1).EQ.'V' .OR. SYMTYP(1:1).EQ.'v') THEN
C
C  Draw a vector scaled by a size based on the ratio between
C  VCNREF and VCUREF, in direction VCWDIR with the
C  tail of the vector at (X,Y).  Scale the vector arrow head
C  by VCHSIZ.
C
C
C  Create a base vector of the correct size that will be
C  rotated and translated appropirately.
C
        SCALE = VCNREF/VCUREF
        HEADANG = VHDANG*D2RAD
        VECTSX(1) = 0.
        VECTSX(2) = SCALE*VCSIZE
        VECTSX(3) = VECTSX(2)-COS(HEADANG)*VCHSIZ
        VECTSX(4) = VECTSX(2)
        VECTSX(5) = VECTSX(3)
C
        VECTSY(1) = 0.
        VECTSY(2) = 0.
        VECTSY(3) = SIN(HEADANG)*VCHSIZ
        VECTSY(4) = 0.
        VECTSY(5) = -VECTSY(3)
C
C  Rotate.
C
        COSANG = COS(D2RAD*VCWDIR)
        SINANG = SIN(D2RAD*VCWDIR)
        DO 110 I=1,IVSDM
          TMPX(I) = VECTSX(I)*COSANG-VECTSY(I)*SINANG
          TMPY(I) = VECTSX(I)*SINANG+VECTSY(I)*COSANG
  110   CONTINUE
        DO 111 I=1,IVSDM
          VECTSX(I) = TMPX(I)
          VECTSY(I) = TMPY(I)
  111   CONTINUE
C
C  Translate the vector.
C
        DO 120 I=1,IVSDM
          VECTSX(I) = VECTSX(I)+XNDC
          VECTSY(I) = VECTSY(I)+YNDC
  120   CONTINUE
C
C  Draw the vector.
C
        CALL GSPLCI(VCCOLR)
        CALL GSLWSC(VCLWID)
        CALL GPL(IVSDM,VECTSX,VECTSY)

      ELSE IF (SYMTYP(1:1).EQ.'D' .OR. SYMTYP(1:1).EQ.'d') THEN
C
C  Draw a dot to mark a city location.
C
        RADIUS = 0.5*CDOTSZ
        CALL NGGETI('CT',ICTYPO)
        CALL NGSETI('CT',0)
        CALL NGDOTS(XNDC,YNDC,1,3.5*RADIUS,IDOTBG)
        CALL NGDOTS(XNDC,YNDC,1,2.0*RADIUS,IDOTCO)
        CALL NGSETI('CT',ICTYPO)
      ELSE IF (SYMTYP(1:2).EQ.'SU' .OR. SYMTYP(1:2).EQ.'su' .OR.
     +         SYMTYP(1:2).EQ.'MC' .OR. SYMTYP(1:2).EQ.'mc' .OR.
     +         SYMTYP(1:2).EQ.'IS' .OR. SYMTYP(1:2).EQ.'is' .OR.
     +         SYMTYP(1:2).EQ.'IT' .OR. SYMTYP(1:2).EQ.'it' .OR.
     +         SYMTYP(1:2).EQ.'MS' .OR. SYMTYP(1:2).EQ.'ms') THEN
C
C  Draw a sun as a daily weather icon.
C
        CALL PCGETI('CC - character color',ICC)
C
C  Cloud first if mostly sunny or intermittant showers.
C
        IF(SYMTYP(1:2).EQ.'MS' .OR. SYMTYP(1:2).EQ.'ms') THEN
          CXOFF =  0.02
          CYOFF = -0.02
          IF (ICLDC3 .GE. 0) THEN
C
            CALL PCGETI('SF - shadow flag',ISF)
            CALL PCGETI('SC - shadow color',ISC)
            CALL PCGETR('SX - shadow X offset',OSX)
            CALL PCGETR('SY - shadow Y offset',OSY)
C
            CALL PCSETI('SF - shadow flag',1)
            CALL PCSETI('SC - shadow color',ICLDC3)
            CALL PCSETI('CC - character color',ICLDC3)
            CALL PCSETR('SX - shadow X offset', -.016)
            CALL PCSETR('SY - shadow Y offset', -.02)
C
            CALL PLCHHQ(XNDC+CXOFF,YNDC+CYOFF,':F37:k',10.*SIZEL,0.,0.)
C
            CALL PCSETI('SF - shadow flag',ISF)
            CALL PCSETI('SC - shadow color',ISC)
            CALL PCSETR('SX - shadow X offset',OSX)
            CALL PCSETR('SY - shadow Y offset',OSY)
          ENDIF
          IF (ICLDC1 .GE. 0) THEN
            CALL PCSETI('CC',ICLDC1)
            CALL PLCHHQ(XNDC+CXOFF,YNDC+CYOFF,':F37:k',10.*SIZEL,0.,0.)
          ENDIF
          IF (ICLDC2 .GE. 0) THEN
            CALL PCSETI('CC',ICLDC2)
            CALL PLCHHQ(XNDC+CXOFF,YNDC+CYOFF,':F137:k',10.*SIZEL,0.,0.)
          ENDIF
        ENDIF
C
C  Draw sun.
C
        SUNSIZ = 7.5*SIZEL
        SXOFF = 0.0
        SYOFF = 0.0
        IF (SYMTYP(1:2).EQ.'MS' .OR. SYMTYP(1:2).EQ.'ms' .OR.
     +      SYMTYP(1:2).EQ.'IT' .OR. SYMTYP(1:2).EQ.'it' .OR.
     +      SYMTYP(1:2).EQ.'IS' .OR. SYMTYP(1:2).EQ.'is' .OR.
     +      SYMTYP(1:2).EQ.'MC' .OR. SYMTYP(1:2).EQ.'mc') THEN
          SXOFF = -0.02
          SYOFF =  0.02
        ENDIF
        IF (ISUNC4 .GE. 0) THEN
C
          CALL PCGETI('SF - shadow flag',ISF)
          CALL PCGETI('SC - shadow color',ISC)
          CALL PCGETR('SX - shadow X offset',OSX)
          CALL PCGETR('SY - shadow Y offset',OSY)
C
          CALL PCSETI('SF - shadow flag',1)
          CALL PCSETI('SC - shadow color',ISUNC4)
          CALL PCSETI('CC - character color',ISUNC4)
          CALL PCSETR('SX - shadow X offset', -.016)
          CALL PCSETR('SY - shadow Y offset', -.02)
C
          CALL PLCHHQ(XNDC+SXOFF,YNDC+SYOFF,':F37:l',SUNSIZ,0.,0.)
C
          CALL PCSETI('SF - shadow flag',ISF)
          CALL PCSETI('SC - shadow color',ISC)
          CALL PCSETR('SX - shadow X offset',OSX)
          CALL PCSETR('SY - shadow Y offset',OSY)
        ENDIF
        IF (ISUNC2 .GE. 0) THEN
          CALL PCSETI('CC',ISUNC2)
          CALL PLCHHQ(XNDC+SXOFF,YNDC+SYOFF,':F37:l',SUNSIZ,0.,0.)
          IF (ISUNC3 .GE. 0) THEN
            CALL PCSETI('CC',ISUNC3)
            CALL PLCHHQ(XNDC+SXOFF,YNDC+SYOFF,':F137:l',SUNSIZ,0.,0.)
          ENDIF
        ENDIF
        IF (ISUNC1 .GE. 0) THEN
          CALL PCSETI('CC',ISUNC1)
          CALL PLCHHQ(XNDC+SXOFF,YNDC+SYOFF,':F37:m',SUNSIZ,0.,0.)
          IF (ISUNC3 .GE. 0) THEN
            CALL PCSETI('CC',ISUNC3)
            CALL PLCHHQ(XNDC+SXOFF,YNDC+SYOFF,':F137:m',SUNSIZ,0.,0.)
          ENDIF
        ENDIF
C
C  Add showers.
C
        IF (SYMTYP(1:2).EQ.'IS' .OR. SYMTYP(1:2).EQ.'is' .OR.
     +      SYMTYP(1:2).EQ.'IT' .OR. SYMTYP(1:2).EQ.'it') THEN
          CALL PCSETI('CC - character color',ICOLOR)
          CALL PLCHHQ(XNDC+CXOFF,YNDC+CYOFF,':F37:o',10.*SIZEL,0.,0.)
        ENDIF
C
C  Add thunder bolt.
C
        IF (SYMTYP(1:2).EQ.'IT' .OR. SYMTYP(1:2).EQ.'it') THEN
          OXNDC = XNDC+0.019*50.*SIZEL+CXOFF
          OYNDC = YNDC-0.06*50.*SIZEL+CYOFF
          TSIZ  = 6.*SIZEL
          CALL PCSETR('ZX',.8)
          CALL PCSETR('ZY',1.3)
          IF (ILTNC3 .GE. 0) THEN
            CALL PCGETI('SF - shadow flag',ISF)
            CALL PCGETI('SC - shadow color',ISC)
            CALL PCGETR('SX - shadow X offset',OSX)
            CALL PCGETR('SY - shadow Y offset',OSY)
C
            CALL PCSETI('SF - shadow flag',1)
            CALL PCSETI('SC - shadow color',ILTNC3)
            CALL PCSETI('CC - character color',ILTNC3)
            CALL PCSETR('SX - shadow X offset', -.015)
            CALL PCSETR('SY - shadow Y offset', -.04)
C
            CALL PLCHHQ(OXNDC,OYNDC,':F37:f',TSIZ,0.,0.)
C
            CALL PCSETI('SF - shadow flag',ISF)
            CALL PCSETI('SC - shadow color',ISC)
            CALL PCSETR('SX - shadow X offset',OSX)
            CALL PCSETR('SY - shadow Y offset',OSY)
          ENDIF
          IF (ILTNC1 .GE. 0) THEN
            CALL PCSETI('CC',ILTNC1)
            CALL PLCHHQ(OXNDC,OYNDC,':F37:f',TSIZ,0.,0.)
          ENDIF
          IF (ILTNC2 .GE. 0) THEN
            CALL PCSETI('CC',ILTNC2)
            CALL PLCHHQ(OXNDC,OYNDC,':F137:f',TSIZ,0.,0.)
          ENDIF
          CALL PCSETR('ZX',1.)
          CALL PCSETR('ZY',1.)
        ENDIF
C
C  Cloud last if mostly cloudy.
C
        IF(SYMTYP(1:2).EQ.'MC' .OR. SYMTYP(1:2).EQ.'mc' .OR.
     +     SYMTYP(1:2).EQ.'IT' .OR. SYMTYP(1:2).EQ.'it' .OR.
     +     SYMTYP(1:2).EQ.'IS' .OR. SYMTYP(1:2).EQ.'is') THEN
          CXOFF =  0.02
          CYOFF = -0.02
          IF (ICLDC3 .GE. 0) THEN
C
            CALL PCGETI('SF - shadow flag',ISF)
            CALL PCGETI('SC - shadow color',ISC)
            CALL PCGETR('SX - shadow X offset',OSX)
            CALL PCGETR('SY - shadow Y offset',OSY)
C
            CALL PCSETI('SF - shadow flag',1)
            CALL PCSETI('SC - shadow color',ICLDC3)
            CALL PCSETI('CC - character color',ICLDC3)
            CALL PCSETR('SX - shadow X offset', -.016)
            CALL PCSETR('SY - shadow Y offset', -.02)
C
            CALL PLCHHQ(XNDC+CXOFF,YNDC+CYOFF,':F37:k',10.*SIZEL,0.,0.)
C
            CALL PCSETI('SF - shadow flag',ISF)
            CALL PCSETI('SC - shadow color',ISC)
            CALL PCSETR('SX - shadow X offset',OSX)
            CALL PCSETR('SY - shadow Y offset',OSY)
          ENDIF
          IF (ICLDC1 .GE. 0) THEN
            CALL PCSETI('CC',ICLDC1)
            CALL PLCHHQ(XNDC+CXOFF,YNDC+CYOFF,':F37:k',10.*SIZEL,0.,0.)
          ENDIF
          IF (ICLDC2 .GE. 0) THEN
            CALL PCSETI('CC',ICLDC2)
            CALL PLCHHQ(XNDC+CXOFF,YNDC+CYOFF,':F137:k',10.*SIZEL,0.,0.)
          ENDIF
        ENDIF
        CALL PCSETI('CC - character color',ICC)
      ELSE IF (SYMTYP(1:1).EQ.'C' .OR. SYMTYP(1:1).EQ.'c'   .OR.
     +         SYMTYP(1:1).EQ.'R' .OR. SYMTYP(1:1).EQ.'r'   .OR.
     +         SYMTYP(1:2).EQ.'RS' .OR. SYMTYP(1:2).EQ.'rs' .OR.
     +         SYMTYP(1:2).EQ.'MC' .OR. SYMTYP(1:2).EQ.'mc' .OR.
     +         SYMTYP(1:2).EQ.'SN' .OR. SYMTYP(1:2).EQ.'sn') THEN
C
C  Draw a cloud as a daily weather icon.
C
        CALL PCGETI('CC - character color',ICC)
        IF (ICLDC3 .GE. 0) THEN
C
          CALL PCGETI('SF - shadow flag',ISF)
          CALL PCGETI('SC - shadow color',ISC)
          CALL PCGETR('SX - shadow X offset',OSX)
          CALL PCGETR('SY - shadow Y offset',OSY)
C
          CALL PCSETI('SF - shadow flag',1)
          CALL PCSETI('SC - shadow color',ICLDC3)
          CALL PCSETI('CC - character color',ICLDC3)
          CALL PCSETR('SX - shadow X offset', -.016)
          CALL PCSETR('SY - shadow Y offset', -.02)
C
          CALL PLCHHQ(XNDC,YNDC,':F37:k',10.*SIZEL,0.,0.)
C
          CALL PCSETI('SF - shadow flag',ISF)
          CALL PCSETI('SC - shadow color',ISC)
          CALL PCSETR('SX - shadow X offset',OSX)
          CALL PCSETR('SY - shadow Y offset',OSY)
        ENDIF
        IF (ICLDC1 .GE. 0) THEN
          CALL PCSETI('CC',ICLDC1)
          CALL PLCHHQ(XNDC,YNDC,':F37:k',10.*SIZEL,0.,0.)
        ENDIF
        IF (ICLDC2 .GE. 0) THEN
          CALL PCSETI('CC',ICLDC2)
          CALL PLCHHQ(XNDC,YNDC,':F137:k',10.*SIZEL,0.,0.)
        ENDIF
        IF (SYMTYP(1:1).EQ.'R'  .OR. SYMTYP(1:1).EQ.'r' .OR.
     +      SYMTYP(1:2).EQ.'RS' .OR. SYMTYP(1:2).EQ.'rs') THEN
          CALL PCSETI('CC - character color',ICOLOR)
          CALL PLCHHQ(XNDC,YNDC,':F37:o',10.*SIZEL,0.,0.)
        ENDIF
        IF (SYMTYP(1:2).EQ.'SN' .OR. SYMTYP(1:2).EQ.'sn' .OR.
     +      SYMTYP(1:2).EQ.'RS' .OR. SYMTYP(1:2).EQ.'rs') THEN
C
C  Draw snowflakes as part of daily weather icon.
C
          DSIZ = 1.3*SIZEL
          CALL PCGETI('CC',ISNFC1)
          CALL PCSETI('CC',ICOLOR)
          DO 60 I=1,3
            XP = XNDC+0.8*DSIZ+2.7*DSIZ*REAL(I-2)
            YP = YNDC-2.40*DSIZ
            CALL PLCHHQ(XP,YP,':F37:q',DSIZ,0.,0.)
   60     CONTINUE
          DO 70 I=1,2
            XP = XNDC+0.8*DSIZ+2.7*DSIZ*(REAL(I)-1.5)
            YP = YNDC-3.50*DSIZ
            CALL PLCHHQ(XP,YP,':F37:q',DSIZ,0.,0.)
   70     CONTINUE
          CALL PCSETI('CC',ISNFC1)
        ENDIF
        CALL PCSETI('CC - character color',ICC)
      ELSE IF (SYMTYP(1:1).EQ.'W' .OR. SYMTYP(1:1).EQ.'w') THEN
          CALL PCGETI('CC',ISNFC1)
          CALL PCSETI('CC',ICOLOR)
          CALL PLCHHQ(XNDC,YNDC,':F37:r',10.*SIZEL,0.,0.)
          CALL PCSETI('CC',ISNFC1)
      ELSE IF (SYMTYP(1:2).EQ.'IC' .OR. SYMTYP(1:2).EQ.'ic') THEN
C
C  Ice. Use the ARROW arrays as scratch space to store a region.
C
        DELX = 5.2*SIZEL
        DELY = 3.0*SIZEL
        ARROWX(1) = XNDC-DELX
        ARROWY(1) = YNDC+DELY
        ARROWX(2) = XNDC
        ARROWY(2) = YNDC+DELY
        ARROWX(3) = XNDC+DELX
        ARROWY(3) = YNDC+DELY
        ARROWX(4) = XNDC+DELX
        ARROWY(4) = YNDC-DELY
        ARROWX(5) = XNDC
        ARROWY(5) = YNDC-DELY
        ARROWX(6) = XNDC-DELX
        ARROWY(6) = YNDC-DELY
        ARROWX(7) = ARROWX(1)
        ARROWY(7) = ARROWY(1)
        CALL WMRGWT(7,ARROWX,ARROWY,37,105)
      ELSE IF (SYMTYP(1:1).EQ.'T' .OR. SYMTYP(1:1).EQ.'t') THEN
C
C  Draw a thunderstorm symbol as a daily weather icon.
C
        CALL PCGETI('CC - character color',ICC)
        CALL PCSETI('CC - character color',ICOLOR)
        CALL PLCHHQ(XNDC,YNDC,':F37:o',10.*SIZEL,0.,0.)
C
        OXNDC = XNDC+0.019*50.*SIZEL
        OYNDC = YNDC-0.06*50.*SIZEL
        TSIZ  = 6.*SIZEL
        CALL PCSETR('ZX',.8)
        CALL PCSETR('ZY',1.3)
        IF (ILTNC3 .GE. 0) THEN
          CALL PCGETI('SF - shadow flag',ISF)
          CALL PCGETI('SC - shadow color',ISC)
          CALL PCGETR('SX - shadow X offset',OSX)
          CALL PCGETR('SY - shadow Y offset',OSY)
C
          CALL PCSETI('SF - shadow flag',1)
          CALL PCSETI('SC - shadow color',ILTNC3)
          CALL PCSETI('CC - character color',ILTNC3)
          CALL PCSETR('SX - shadow X offset', -.015)
          CALL PCSETR('SY - shadow Y offset', -.04)
C
          CALL PLCHHQ(OXNDC,OYNDC,':F37:f',TSIZ,0.,0.)
C
          CALL PCSETI('SF - shadow flag',ISF)
          CALL PCSETI('SC - shadow color',ISC)
          CALL PCSETR('SX - shadow X offset',OSX)
          CALL PCSETR('SY - shadow Y offset',OSY)
        ENDIF
        IF (ILTNC1 .GE. 0) THEN
          CALL PCSETI('CC',ILTNC1)
          CALL PLCHHQ(OXNDC,OYNDC,':F37:f',TSIZ,0.,0.)
        ENDIF
        IF (ILTNC2 .GE. 0) THEN
          CALL PCSETI('CC',ILTNC2)
          CALL PLCHHQ(OXNDC,OYNDC,':F137:f',TSIZ,0.,0.)
        ENDIF
        CALL PCSETR('ZX',1.)
        CALL PCSETR('ZY',1.)
C
        IF (ICLDC3 .GE. 0) THEN
          CALL PCGETI('SF - shadow flag',ISF)
          CALL PCGETI('SC - shadow color',ISC)
          CALL PCGETR('SX - shadow X offset',OSX)
          CALL PCGETR('SY - shadow Y offset',OSY)
C
          CALL PCSETI('SF - shadow flag',1)
          CALL PCSETI('SC - shadow color',ICLDC3)
          CALL PCSETI('CC - character color',ICLDC3)
          CALL PCSETR('SX - shadow X offset', -.016)
          CALL PCSETR('SY - shadow Y offset', -.02)
C
          CALL PLCHHQ(XNDC,YNDC,':F37:k',10.*SIZEL,0.,0.)
C
          CALL PCSETI('SF - shadow flag',ISF)
          CALL PCSETI('SC - shadow color',ISC)
          CALL PCSETR('SX - shadow X offset',OSX)
          CALL PCSETR('SY - shadow Y offset',OSY)
        ENDIF
        IF (ICLDC1 .GE. 0) THEN
          CALL PCSETI('CC',ICLDC1)
          CALL PLCHHQ(XNDC,YNDC,':F37:k',10.*SIZEL,0.,0.)
        ENDIF
        IF (ICLDC2 .GE. 0) THEN
          CALL PCSETI('CC',ICLDC2)
          CALL PLCHHQ(XNDC,YNDC,':F137:k',10.*SIZEL,0.,0.)
        ENDIF
        CALL PCSETI('CC - character color',ICC)
      ENDIF
C
C  Restore original attrributes.
C
      CALL GSFAIS(IFAISO)
      CALL GSFACI(IFCLRO)
      CALL GSPLCI(ILCLRO)
      CALL GSELNT(NTRO)
      CALL GSWN(1,OLDWN(1),OLDWN(2),OLDWN(3),OLDWN(4))
      CALL GSVP(1,OLDVP(1),OLDVP(2),OLDVP(3),OLDVP(4))
      CALL PCSETI('TE',0)
      CALL PCSETC('FC',FCODE)
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl