text:
module feval
  use bdf
  use bl_types
  implicit none
  integer, parameter :: neq = 1
contains
  subroutine f(neq, npt, y, t, ydot, upar)
    integer,    intent(in)    :: neq, npt
    real(dp_t), intent(in)    :: y(neq,npt), t(npt)
    real(dp_t), intent(  out) :: ydot(neq,npt)
    real(dp_t), intent(inout), optional :: upar(:,:)
  end subroutine f
  subroutine J(neq, npt, y, t, pd, upar)
    integer,    intent(in)    :: neq, npt
    real(dp_t), intent(in)    :: y(neq,npt), t(npt)
    real(dp_t), intent(  out) :: pd(neq,neq,npt)
    real(dp_t), intent(inout), optional :: upar(:,:)
  end subroutine J
end module feval

program test
  use bdf
  use feval

  double precision, parameter :: tol = 1.d-12
  integer,          parameter :: NPT = 1

  type(bdf_ts)  :: ts
  double precision :: rtol(neq), atol(neq)
  double precision :: v(0:6)

  rtol = 1
  atol = 1

  call bdf_ts_build(ts, neq, 1, rtol, atol, max_order=6)

  print *, "====> order 1"
  ts%k = 1
  call random_number(ts%h(:,NPT))
  call bdf_update(ts)
  print *, 'l', ts%l(0:1,NPT)
  call assert(all(ts%l(0:1,NPT) == [ 1.d0, 1.d0 ]), "error in l")

  ts%h = 1
  call bdf_update(ts)
  print *, 'tq', ts%tq
  call assert(all(abs(ts%tq(:,NPT) - [ 1.d0, 0.5d0, 0.2222222222222d0, 2.0d0 ]) < tol), "error in tq")

  print *, "====> order 2"
  ts%k = 2
  call random_number(ts%h(:,NPT))
  call bdf_update(ts)
  print *, 'l', ts%l(0:2,NPT)
  call assert(all(abs(ts%l(0:2,NPT) - [ 1.d0, 1.5d0, 0.5d0 ]) < tol), "error in l")

  ts%h = 1
  call bdf_update(ts)
  print *, 'tq', ts%tq
  call assert(all(abs(ts%tq(:,NPT) - [ 1.0d0, 0.222222222222222d0, 0.13636363636363638d0, 6.0d0 ]) < tol), "error in tq")

  print *, "====> order 3"
  ts%k = 3
  ts%h = 1
  call bdf_update(ts)
  print *, 'l', ts%l(0:3,NPT)
  v(0:3) = [ 1.d0, 1.8333333333333d0, 1.d0, 0.1666666666666d0 ]
  call assert(all(abs(ts%l(0:3,NPT) - v(0:3)) < tol), "error in l")

  ts%h(1,NPT) = 4
  call bdf_update(ts)
  print *, 'l', ts%l(0:3,NPT)
  v(0:3) = [ 1.d0, 1.8333333333333d0, 0.96d0, 0.12666666666666d0 ]
  call assert(all(abs(ts%l(0:3,NPT) - v(0:3)) < tol), "error in l")

  ts%h = 1
  call bdf_update(ts)
  print *, 'tq', ts%tq
  call assert(all(abs(ts%tq(:,NPT) - [ 1.3333333333333d0, 0.13636363636363638d0, 9.6d-2, 24.0d0 ]) < tol), "error in tq")

  call bdf_ts_destroy(ts)

contains

  subroutine assert(cond, message)
    logical, intent(in) :: cond
    character(*), intent(in) :: message

    if (.not. cond) then
       print *, "ASSERTION FAILED: ", message
       call abort()
    end if
  end subroutine assert

end program test

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO


text:
!> \file sw_adtg.f90
!! \BRIEF
!> Module with sw_adtg function - compute adiabatic temp. gradient from S,T,P
MODULE msw_adtg
CONTAINS
!>  Function to calculate adiabatic temperature gradient as per UNESCO 1983 routines.
FUNCTION sw_adtg  (s,t,p)

  !     ==================================================================
  !     Calculates adiabatic temperature gradient as per UNESCO 1983 routines.
  !     Armin Koehl akoehl@ucsd.edu
  !     ==================================================================
  USE msingledouble
  IMPLICIT NONE
  !> salinity [psu (PSU-78)]
  REAL(kind=r8) :: s
  !> temperature [degree C (IPTS-68)]
  REAL(kind=r8) :: t
  !> pressure [db]
  REAL(kind=r8) :: p

  REAL(kind=r8) :: a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2
  REAL(kind=r8) :: sref

  REAL(kind=r8) :: sw_adtg

  sref = 35.d0
  a0 =  3.5803d-5
  a1 = +8.5258d-6
  a2 = -6.836d-8
  a3 =  6.6228d-10

  b0 = +1.8932d-6
  b1 = -4.2393d-8

  c0 = +1.8741d-8
  c1 = -6.7795d-10
  c2 = +8.733d-12
  c3 = -5.4481d-14

  d0 = -1.1351d-10
  d1 =  2.7759d-12

  e0 = -4.6206d-13
  e1 = +1.8676d-14
  e2 = -2.1687d-16

  sw_adtg =  a0 + (a1 + (a2 + a3*T)*T)*T &
       + (b0 + b1*T)*(S-sref) &
       + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P &
       + (  e0 + (e1 + e2*T)*T )*P*P

END FUNCTION sw_adtg

!>  Function to calculate adiabatic temperature gradient as per UNESCO 1983 routines.
!! and derivative with respect to temperature and salinity
FUNCTION sw_adtg_DNAD  (s,t,p)

  !     ==================================================================
  !     Calculates adiabatic temperature gradient as per UNESCO 1983 routines.
  !     Armin Koehl akoehl@ucsd.edu
  !     ==================================================================
  USE msingledouble
  USE Dual_Num_Auto_Diff
  IMPLICIT NONE
  !> salinity [psu (PSU-78)]
  TYPE(DUAL_NUM) :: s
  !> temperature [degree C (IPTS-68)]
  TYPE(DUAL_NUM) :: t
  !> pressure [db]
  TYPE(DUAL_NUM) :: p

  REAL(kind=r8) :: a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2
  REAL(kind=r8) :: sref

  TYPE(DUAL_NUM) :: sw_adtg_DNAD

  sref = 35.d0
  a0 =  3.5803d-5
  a1 = +8.5258d-6
  a2 = -6.836d-8
  a3 =  6.6228d-10

  b0 = +1.8932d-6
  b1 = -4.2393d-8

  c0 = +1.8741d-8
  c1 = -6.7795d-10
  c2 = +8.733d-12
  c3 = -5.4481d-14

  d0 = -1.1351d-10
  d1 =  2.7759d-12

  e0 = -4.6206d-13
  e1 = +1.8676d-14
  e2 = -2.1687d-16

  sw_adtg_DNAD =  a0 + (a1 + (a2 + a3*T)*T)*T &
       + (b0 + b1*T)*(S-sref) &
       + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P &
       + (  e0 + (e1 + e2*T)*T )*P*P

END FUNCTION sw_adtg_DNAD
END MODULE msw_adtg

max_stars_repo_name: tomaslovato/mocsy

max_issues_repo_name: tomaslovato/mocsy

max_forks_repo_name: tomaslovato/mocsy


text:
      program t
      integer::i,j,k,l,m
      save i
 100  save j,k
      save::l,m
      print *,'ok'
      end program t

ext: f

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
!-----------------------------------------------------------------------
! Fortran interface to CUDA utility Library
! written by Viktor K. Decyk, UCLA
      module gpulib2_h
      implicit none
!
      interface
         subroutine gpu_setgbsize(nblock)
         implicit none
         integer :: nblock
         end subroutine
      end interface
!
      interface
         integer function getmmcc()
         implicit none
         end function
      end interface
!
      interface
         subroutine gpu_fallocate(gp_f,nsize,irc)
         implicit none
         integer :: nsize, irc
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_iallocate(gp_i,nsize,irc)
         implicit none
         integer :: nsize, irc
         integer, dimension(2) :: gp_i
         end subroutine
      end interface
!
      interface
         subroutine gpu_callocate(gp_c,nsize,irc)
         implicit none
         integer :: nsize, irc
         integer, dimension(2) :: gp_c
         end subroutine
      end interface
!
      interface
         subroutine gpu_deallocate(gp_d,irc)
         implicit none
         integer :: irc
         integer, dimension(2) :: gp_d
         end subroutine
      end interface
!
      interface
         subroutine hpl_f1allocate(hp_f,nx,irc)
         implicit none
         integer :: nx, irc
         real, dimension(:), pointer :: hp_f
         end subroutine
      end interface
!
      interface
         subroutine hpl_f2allocate(hp_f,nx,ny,irc)
         implicit none
         integer :: nx, ny, irc
         real, dimension(:,:), pointer :: hp_f
         end subroutine
      end interface
!
      interface
         subroutine hpl_c2allocate(hp_f,nx,ny,irc)
         implicit none
         integer :: nx, ny, irc
         complex, dimension(:,:), pointer :: hp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_fcopyin(f,gp_f,nsize)
         implicit none
         integer :: nsize
         real, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_fcopyout(f,gp_f,nsize)
         implicit none
         integer :: nsize
         real, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_icopyin(f,gp_f,nsize)
         implicit none
         integer :: nsize
         integer, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_icopyout(f,gp_f,nsize)
         implicit none
         integer :: nsize
         integer, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_ccopyin(f,gp_f,nsize)
         implicit none
         integer :: nsize
         complex, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_ccopyout(f,gp_f,nsize)
         implicit none
         integer :: nsize
         complex, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_initstream(nstream)
         implicit none
         integer :: nstream
         end subroutine
      end interface
!
      interface
         subroutine gpu_delstream(nstream)
         implicit none
         integer :: nstream
         end subroutine
      end interface
!
      interface
         subroutine waitstream(nstream)
         implicit none
         integer :: nstream
         end subroutine
      end interface
!
      interface
         subroutine gpu_cascopyin(f,gp_f,noff,nsize,nstream)
         implicit none
         integer :: noff, nsize, nstream
         complex, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_cascopyout(f,gp_f,noff,nsize,nstream)
         implicit none
         integer :: noff, nsize, nstream
         complex, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_zfmem(gp_f,nsize)
         implicit none
         integer :: nsize
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_set_cache_size(nscache)
         implicit none
         integer :: nscache
         end subroutine
      end interface
!
      interface
         subroutine emptykernel()
         implicit none
         end subroutine
      end interface
!
      interface
         subroutine init_cu(dev,irc)
         implicit none
         integer :: dev, irc
         end subroutine
      end interface
!
      interface
         subroutine end_cu()
         implicit none
         end subroutine
      end interface
!
      end module

max_stars_repo_name: gcasabona/cuda

max_issues_repo_name: gcasabona/cuda

max_forks_repo_name: gcasabona/cuda


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

program location_test

! Simple test program to exercise annulus location module.

use     types_mod, only : r8
use utilities_mod, only : open_file, close_file, initialize_utilities, &
                          finalize_utilities
use  location_mod

implicit none

! version controlled file description for error handling, do not edit
character(len=256), parameter :: source   = &
   "$URL$"
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate  = "$Date$"

type(location_type) :: loc0(6), loc1, loc2, loc3, loc4, locA(7), locB(5)
integer             :: iunit, iunit1, iunit2, i
character(len=102)  :: testbuf

call initialize_utilities('location_test')

print *, ''
print *, 'reading in 6 locations to use in distance tests:'
do i = 1, 6
   call interactive_location(loc0(i))
enddo

print *, ''
print *, 'testing distances between pairs of points:'
do i = 1, 5
   loc1 = loc0(i)
   loc2 = loc0(i+1)
   print *, 'prev location is ', get_location(loc1)
   print *, 'this location is ', get_location(loc2)
   print *, 'distance prev to this is ', get_dist(loc1, loc2)
   print *, 'distance this to prev is ', get_dist(loc2, loc1)
end do

print *, ''
print *, 'reading in 7 locations for file read/write tests:'

! read in 7 different locs for the tests below
do i = 1, 7
   call interactive_location(locA(i))
enddo

loc1 = locA(1)
loc2 = locA(2)

print *, ''
print *, 'location 1 is ', get_location(loc1)
print *, 'location 2 is ', get_location(loc2)
print *, 'distance   is ', get_dist(loc1, loc2)

print *, ''
print *, 'testing ascii file read/write:'

! Open an output file
iunit = open_file('location_test_file.txt', form='formatted', action='write')

! Write these locations to the file
call write_location(iunit, loc1)
call write_location(iunit, loc2)

do i = 1, 5
   call write_location(iunit, locA(i))
enddo

call close_file(iunit)

loc3 = set_location_missing()
loc4 = set_location_missing()

! Now read them back in and compute the distances between them
iunit1 = open_file('location_test_file.txt', form='formatted', action='read')

loc3 = read_location(iunit1)
loc4 = read_location(iunit1)

print *, ''
print *, 'read-back location 1 is ', get_location(loc3)
print *, 'read-back location 2 is ', get_location(loc4)
print *, 'distance   is ', get_dist(loc3, loc4)

print *, ''
print *, 'testing binary file read/write:'

! Open an output file
iunit = open_file('location_test_file.bin', form='unformatted', action='write')

! Write this location to the file
call write_location(iunit, loc1, 'unformatted')
call write_location(iunit, loc2, 'unformatted')

do i = 1, 5
   call write_location(iunit, locA(i), 'unformatted')
enddo

close(iunit)

loc3 = set_location_missing()
loc4 = set_location_missing()

! Now read them back in and compute the distances between them
iunit2 = open_file('location_test_file.bin', form='unformatted', action='read')

loc3 = read_location(iunit2, 'unformatted')
loc4 = read_location(iunit2, 'unformatted')

print *, ''
print *, 'read-back location 1 is ', get_location(loc3)
print *, 'read-back location 2 is ', get_location(loc4)
print *, 'distance   is ', get_dist(loc3, loc4)

! testing write to a character buffer
print *, ''
write(*,*) 'testing write to a char buffer'

print *, ''
print *, 'raw locations: '
do i = 1, 5
   print *, get_location(locA(i))
enddo
print *, ''
print *, 'locations formatted to a char buffer:'
do i = 1, 5
   call write_location(0, locA(i), 'formatted', testbuf)
   print *, 'string length: ', len_trim(testbuf)
enddo
do i = 1, 5
   call write_location(0, locA(i), 'formatted', testbuf)
   print *, trim(testbuf)
enddo

print *, ''
print *, 'ascii readback: '
do i = 1, 5
   locB(i) = set_location_missing()
   locB(i) = read_location(iunit1, 'formatted')
   print *, get_location(locB(i))
enddo
do i = 1, 5
   call write_location(0, locB(i), 'formatted', testbuf)
   print *, trim(testbuf)
enddo

! Now read them back in and compute the distances between them
open(iunit, file = 'location_test_file.bin', form='unformatted', action='read')

print *, ''
print *, 'binary readback: '
do i = 1, 5
   locB(i) = set_location_missing()
   locB(i) = read_location(iunit2, 'unformatted')
   print *, get_location(locB(i))
enddo
do i = 1, 5
   call write_location(0, locB(i), 'formatted', testbuf)
   print *, trim(testbuf)
enddo

call close_file(iunit1)
call close_file(iunit2)

call finalize_utilities()

end program location_test

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
      SUBROUTINE DIAT(NI,NJ,XI,XJ,DI)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
***************************************************************************
*
*   DIAT CALCULATES THE DI-ATOMIC OVERLAP INTEGRALS BETWEEN ATOMS
*        CENTERED AT XI AND XJ.
*
*   ON INPUT NI  = ATOMIC NUMBER OF THE FIRST ATOM.
*            NJ  = ATOMIC NUMBER OF THE SECOND ATOM.
*            XI  = CARTESIAN COORDINATES OF THE FIRST ATOM.
*            XJ  = CARTESIAN COORDINATES OF THE SECOND ATOM.
*
*  ON OUTPUT DI  = DIATOMIC OVERLAP, IN A 9 * 9 MATRIX. LAYOUT OF
*                  ATOMIC ORBITALS IN DI IS
*                  1   2   3   4   5            6     7       8     9
*                  S   PX  PY  PZ  D(X**2-Y**2) D(XZ) D(Z**2) D(YZ) D(XY)
*
*   LIMITATIONS:  IN THIS FORMULATION, NI AND NJ MUST BE LESS THAN 54
*         EXPONENTS ARE ASSUMED TO BE PRESENT IN COMMON BLOCK EXPONT.
*
***************************************************************************
      INTEGER PQ1,A,PQ2,B,PQA,PQB,AA,BB,PVAL1,PVAL2,YETA
      LOGICAL FIRST
      COMMON /EXPONT/ EMUS(54),EMUP(54),EMUD(54)
      DIMENSION DI(9,9),S(3,3,3),UL1(3),UL2(3),C(3,5,5),NPQ(54)
     +          ,XI(3),XJ(3)
      DATA NPQ/1,1, 2,2,2,2,2,2,2,2, 3,3,3,3,3,3,3,3, 4,4,4,4,4,4,4,4,
     +4,4,4,4,4,4,4,4,4,4, 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5/
      DATA FIRST /.TRUE./
      X1=XI(1)
      X2=XJ(1)
      Y1=XI(2)
      Y2=XJ(2)
      Z1=XI(3)
      Z2=XJ(3)
      PQ1=NPQ(NI)
      PQA=PQ1
      PQ2=NPQ(NJ)
      PQB=PQ2
      UL1(1)=EMUS(NI)
      UL2(1)=EMUS(NJ)
      UL1(2)=EMUP(NI)
      UL2(2)=EMUP(NJ)
      DO 349 I=1,9
      DO 350 J=1,9
      DI(I,J)=0.0D0
 350  CONTINUE
 349  DI(I,I)=1.D0
      CALL COE(X1,Y1,Z1,X2,Y2,Z2,PQ1,PQ2,C,R)
          IF (PQ1.GT.3) GO TO 351
          A=PQ1-1
          GO TO 352
 351      A=2
 352      CONTINUE
          IF (PQ2.GT.3) GO TO 353
          B=PQ2-1
          GO TO 354
 353      B=2
 354      CONTINUE
          IA=A+1
          IB=B+1
       IF(PQ1.LT.4.AND.PQ2.LT.4) THEN
          CALL DIAT2(NI,EMUS(NI),EMUP(NI),R,NJ,EMUS(NJ),EMUP(NJ),S)
C      WRITE(6,'(3I4,F12.6)')(((I,J,K,S(I,J,K),I=1,2),J=1,2),K=1,2)
      ELSE
          DO 355 I=1,3
              DO 355 J=1,3
                  DO 355 K=1,3
                  S(I,J,K)=0.0D0
 355      CONTINUE
          DO 363 I=1,IA
              IB=B+1
              DO 363 J=1,IB
                  IF (A.LT.B) GO TO 357
                  NEWK=B
                  GO TO 358
 357              NEWK=A
 358              CONTINUE
                  NK1=NEWK+1
                  DO 363 K=1,NK1
                      IF(K.GT.I.OR.K.GT.J) GOTO 363
                      U1=UL1(I)
                      U2=UL2(J)
                      YETA=PQA+PQB+3
                      S(I,J,K)=SS(PQA,PQB,I,J,K,U1,U2,R,YETA,FIRST)
C      WRITE(6,'(3I4,F12.6)')I,J,K,S(I,J,K)
 363      CONTINUE
      ENDIF
      DO 369 I=1,IA
      DO 369 J=1,IB
      IF (I.EQ.1) GO TO 10
      IF (I.EQ.2) GO TO 11
      KMIN=1
      KMAX=5
      GO TO 12
 11   KMIN=2
      KMAX=4
      GO TO 12
 10   KMIN=3
      KMAX=3
 12   DO 369 K=KMIN,KMAX
      IF (J.EQ.1) GO TO 13
      IF (J.EQ.2) GO TO 14
      LMIN=1
      LMAX=5
      GO TO 15
 14   LMIN=2
      LMAX=4
      GO TO 15
 13   LMIN=3
      LMAX=3
 15   DO 369 L=LMIN,LMAX
      IF (J.EQ.2) GO TO 364
      AA=1
      GO TO 365
 364  AA=-1
 365  CONTINUE
      IF (J.GT.2) GO TO 366
      BB=1
      GO TO 367
 366  BB=-1
 367  CONTINUE
      CALL VAL(I,K,PVAL1)
      CALL VAL(J,L,PVAL2)
      DI((PVAL1+1),(PVAL2+1))=S(I,J,1)*C(I,K,3)*C(J,L,3)*AA+(C(I,K,4)*C(
     1J,L,4)+C(I,K,2)*C(J,L,2))*BB*S(I,J,2)+(C(I,K,5)*C(J,L,5)+C(I,K,1)*
     2C(J,L,1))*S(I,J,3)
 369  CONTINUE
      RETURN
      END
      SUBROUTINE VAL(L,M,PVAL)
      INTEGER L,M,PVAL
      IF (L.EQ.1) PVAL=0
      IF ((L.EQ.2).AND.(M.EQ.2)) PVAL=2
      IF ((L.EQ.2).AND.(M.EQ.3)) PVAL=3
      IF ((L.EQ.2).AND.(M.EQ.4)) PVAL=1
      IF (L.EQ.3) PVAL=12-L-M
      RETURN
      END
      FUNCTION SS(NA,NB,LA,LB,M,UC,UD,R1,YETA,FIRST)
      DOUBLE PRECISION UC,SS,UD,R1,R
      LOGICAL FIRST
      INTEGER A,C,D,PP,B,Q,YETA
      DOUBLE PRECISION ER,P,BA,EX,QUO,S1,S2,QQ,QEB,SUM,SUM1,X,
     1FA(14),BI(13,13),AFF(3,3,3),AF(20),BF(20),UA,UB,SA
      DATA AFF/27*0. D0/
      R=R1
      UA=UC
      UB=UD
      IF(UA.GT.0.D0) GOTO 88
      SA=0.D0
      GO TO 99
   88 IF(UB.GT.0.D0) GOTO 299
      SA=0.D0
      GO TO 99
 299  CONTINUE
      R=R/0.529167D0
      ER=R
      GO TO 304
  300 FA(1)=1.D0
      DO 301 I=1,13
      FA(I+1)=FA(I)*I
 301  CONTINUE
      FIRST=.FALSE.
      DO 302 I=1,13
      BI(I,1)=1.D0
      BI(I,I)=1.D0
 302  CONTINUE
      DO 303 I=1,12
      I1=I-1
      DO 303 J=1,I1
      BI(I+1,J+1)=BI(I,J+1)+BI(I,J)
 303  CONTINUE
      AFF(1,1,1)=1.D0
      AFF(2,1,1)=1.D0
      AFF(2,2,1)=1.D0
      AFF(3,1,1)=1.5D0
      AFF(3,2,1)=1.73205D0
      AFF(3,3,1)=1.224745D0
      AFF(3,1,3)=-0.5D0
      GO TO 305
  304 IF(FIRST)GOTO 300
 305  CONTINUE
      P=(UA+UB)*ER*0.5D0
      BA=(UA-UB)*ER*0.5D0
      EX=EXP(BA)
      QUO=1/P
      AF(1)=QUO*EXP(-P)
      NANB=NA+NB
      DO 306 N=1,19
      AF(N+1)=N*QUO*AF(N)+AF(1)
 306  CONTINUE
      NANB1=NANB+1
      DO 309 N=1,13
      IF(ABS(BA).LT.0.1D0) GOTO 308
      S1=0.D0
      S2=0.D0
      QQ=1.D0
      DO 307 I=1,N
      QQ=QQ*BA
      QEB=1.D0/(QQ*FA(N+1-I))
      S1=S1+QEB
      INI=(N-I+1)/2
      S2=S2+2.D0*((INI*2-(N-1)+I)-1.5D0)*QEB
 307  CONTINUE
      BF(N)=-FA(N)*(S1/EX+S2*EX)
      GO TO 309
 308  IN=N/2
      BF(N)=0.D0
      IF(N.NE.IN*2) BF(N)=2.D0/N
 309  CONTINUE
      SUM=0.D0
      LAM1=LA-M+1
      LBM1=LB-M+1
      DO 311 I=1,LAM1,2
      DO 311 J=1,LBM1,2
      A=NA+I-LA
      B=NB+J-LB
      C=LA-I-M+1
      D=LB-J-M+1
      SUM1=0.D0
      IA=A+1
      IB=B+1
      IC=C+1
      ID=D+1
      AB=A+B-1
      DO 310 K1=1,IA
      DO 310 K2=1,IB
      DO 310 K3=1,IC
      DO 310 K4=1,ID
      DO 310 K5=1,M
      DO 310 K6=1,M
      Q=AB-K1-K2+K3+K4+2*K5
      PP=K1+K2+K3+K4+2*K6-5
      JX=M+K2+K4+K5+K6-5
      IX=JX/2
      SUM1=SUM1+BI(IA,K1)*BI(IB,K2)*BI(IC,K3)*BI(ID,K4)*BI(M,K5)*BI(
     1M,K6)*2.D0*(IX*2-JX+0.5D0)*AF(Q)*BF(PP)
 310  CONTINUE
      SUM=SUM+SUM1*AFF(LA,M,I)*AFF(LB,M,J)
 311  CONTINUE
      X=R
      DO 312 I=1,NA
      X=X*R*UA
 312  CONTINUE
      DO 313 I=1,NB
      X=X*R*UB
 313  CONTINUE
      SA=SUM*X*SQRT(UA*UB/(FA(NA+NA+1)*FA(NB+NB+1))*((LA+LA-1
     1)*(LB+LB-1)))/(2.D0**M)
 99   CONTINUE
      SS=SA
      RETURN
      END
      SUBROUTINE COE(X1,Y1,Z1,X2,Y2,Z2,PQ1,PQ2,C,R)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER PQ1,PQ2,PQ,CO
      DIMENSION C(3,5,5)
      XY=(X2-X1)**2+(Y2-Y1)**2
      R=SQRT(XY+(Z2-Z1)**2)
      XY=SQRT(XY)
      IF (XY.EQ.0.D0) GO TO 800
      CA=(X2-X1)/XY
      CB=(Z2-Z1)/R
      SA=(Y2-Y1)/XY
      SB=XY/R
      GO TO 804
 800  IF (Z2-Z1) 801,802,803
 801  CA=-1.D0
      CB=-1.D0
      SA=0.D0
      SB=0.D0
      GO TO 804
  802 CA=0.D0
      CB=0.D0
      SA=0.D0
      SB=0.D0
      GO TO 804
  803 CA=1.D0
      CB=1.D0
      SA=0.D0
      SB=0.D0
 804  CONTINUE
      CO=0
      DO 805 I=1,5
      DO 805 J=1,5
      DO 805 K=1,3
      C(K,I,J)=0.D0
 805  CONTINUE
      IF (PQ1.GT.PQ2) GO TO 806
      PQ=PQ2
      GO TO 807
 806  PQ=PQ1
 807  CONTINUE
      C(1,3,3)=1.D0
      IF (PQ.LT.2) GO TO 808
      C(2,4,4)=CA*CB
      C(2,4,3)=CA*SB
      C(2,4,2)=-SA
      C(2,3,4)=-SB
      C(2,3,3)=CB
      C(2,3,2)=0.D0
      C(2,2,4)=SA*CB
      C(2,2,3)=SA*SB
      C(2,2,2)=CA
      IF (PQ.LT.3) GO TO 808
      C2A=2*CA*CA-1.D0
      C2B=2*CB*CB-1.D0
      S2A=2*SA*CA
      S2B=2*SB*CB
      C(3,5,5)=C2A*CB*CB+0.5D0*C2A*SB*SB
      C(3,5,4)=0.5D0*C2A*S2B
      C(3,5,3)=0.8660254037841D0*C2A*SB*SB
      C(3,5,2)=-S2A*SB
      C(3,5,1)=-S2A*CB
      C(3,4,5)=-0.5D0*CA*S2B
      C(3,4,4)=CA*C2B
      C(3,4,3)=0.8660254037841D0*CA*S2B
      C(3,4,2)=-SA*CB
      C(3,4,1)=SA*SB
      C(3,3,5)=0.5773502691894D0*SB*SB*1.5D0
      C(3,3,4)=-0.8660254037841D0*S2B
      C(3,3,3)=CB*CB-0.5D0*SB*SB
      C(3,2,5)=-0.5D0*SA*S2B
      C(3,2,4)=SA*C2B
      C(3,2,3)=0.8660254037841D0*SA*S2B
      C(3,2,2)=CA*CB
      C(3,2,1)=-CA*SB
      C(3,1,5)=S2A*CB*CB+0.5D0*S2A*SB*SB
      C(3,1,4)=0.5D0*S2A*S2B
      C(3,1,3)=0.8660254037841D0*S2A*SB*SB
      C(3,1,2)=C2A*SB
      C(3,1,1)=C2A*CB
 808  CONTINUE
      RETURN
      END

ext: for

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
      program wave
      implicit none

      integer i, k, ni, nk
      parameter (ni = 3, nk = 3)
      real dt,dx,v,a,b1,dp

      real p(0:ni,0:nk),p0(0:ni,0:nk)

! compute the constant coeffect of the stiffness matrix
      b1 = dt/dx*v
      a = b1*b1

! First inital conditiona and boundary conditions
! Left boundary
      p(0,0) = 0.0
! Time iteration -1 = space loop.
      do 1 i = 1, ni
         p(i,0) = p0(i*dx)
 1    continue
! Right boundary
      p(ni,0) = 0.0

! Second initial condition (with p1=0) and boundary conditions
! Left boundary
      p(0,1) = 0.0
! Time iteration 0 = space loop.
      do 2 i = 1, ni
         dp = p(i+1,0) - 2*p(i,0) + p(i-1,0)
         p(i,1) = p(i,0) + 0.5*a*dp
 2    continue
! Right boundary
      p(ni,1) = 0.0

! Evolution problem and boundary conditions

      do 3 k = 1, nk
! Left boundary
         p(0, k+1) = 0
! Time iteration k = space loop
         do 4 i = 1, ni
            dp = p(i+1,k) - 2*p(i,k) + p(i-1,k)
            p(i,k+1) = 2*p(i,k) - p(i,k-1) + a*dp
 4       continue
! Right boundary
      p(ni,k+1) = 0.0
 3    continue
      end

max_stars_repo_name: raehik/camfort

max_issues_repo_name: apthorpe/camfort

max_forks_repo_name: raehik/camfort


module calculate_parser

    use, intrinsic :: iso_c_binding, only: CalculateHandler => c_ptr
    implicit none
    private

    public :: CalculateParser, calculate_get_parser, calculate_get_default_parser

    type :: CalculateParser
        type(CalculateHandler), private :: handler
    contains
        procedure :: check
        final :: free
    end type

    interface
        module function calculate_get_parser() result (parser)
            type(CalculateParser) :: parser
        end function

        module function calculate_get_default_parser() result (parser)
            type(CalculateParser) :: parser
        end function

        module function check(self)
            class(CalculateParser), intent(in) :: self
            logical :: check
        end function

        module subroutine free(parser)
            type(CalculateParser), intent(inout) :: parser
        end subroutine
    end interface

end module


text:
C $Id$
************************************************************************
*                                                                      *
      subroutine ecp_up_k (m_min,m_max,k,j,l,n,ldQ,alpha,beta,
     &    gamma,Q_lo,Q_mid,Q_hi)
*                                                                      *
*   Perform upward recursion in k for k even, using 6.38:              *
*   2c Q^{k+2}_{m,m+j} = (k-j) Q^k_{m,m+j} + a Q^{k+1}_{m+1,m+j}       *
*                    + b Q^{k+1}_{m,m-1+j}                             *
*   = gamma [(gamma (k-j)/2 Q^k_{m,m+j} + alpha Q^{k+1}_{m+1,m+j}      *
*                    + beta Q^{k+1}_{m,m-1+j} ]                        *
*                                                                      *
*   Argument (status) - description                                    *
*                                                                      *
*   m_min (inp) - minimum value of m                                   *
*   m_max (inp) - maximum value of m                                   *
*   k (inp) - value of k                                               *
*   j (inp) - value of j                                               *
*   l (inp) - switches alpha and beta for different forms of recursion *
*             formula.                                                 *
*   n (inp) - number of functions                                      *
*   alpha (inp) -  a/2*sqrt(c)                                         *
*   beta (inp) -  b/2*sqrt(c)                                          *
*   gamma (inp) - 1/sqrt(c)                                            *
*   Q_mid (inp) - Q^{k+1}_{m+1,m}                                      *
*   Q_lo (inp) - Q^{k}_{m,m}                                           *
*   Q_hi (out) - Q^{k+2}_{m,m}                                         *
*                                                                      *
*   Written by K. G. Dyall                                             *
*                                                                      *
************************************************************************
      implicit none
      integer i,j,k,l,ldQ,m,m_min,m_max,n
      double precision two,fac,alpha(n),beta(n),gamma(n),
     &    Q_lo(ldQ,m_min:m_max),Q_mid(ldQ,m_min:m_max+1),
     &    Q_hi(ldQ,m_min:m_max)
*
      parameter (two = 2.0d00)
*
      fac = k-j
      fac = fac/two
C      write (6,'(/A,/)') 'upward_k_recursion'
C      write (6,*) 'k,j,m_min,m_max',k,j,m_min,m_max
      do m = m_min,m_max-j
C        write (6,*) 'm =',m
        do i = 1,n
          Q_hi(i,m) = gamma(i)*( gamma(i)*fac*Q_lo(i,m)
     &        + alpha(i)*Q_mid(i,m+l)+ beta(i)*Q_mid(i,m+1-l) )
C          write (6,'(1p4e20.10)') gamma(i),alpha(i),beta(i)
C          write (6,'(1p4e20.10)') Q_lo(i,m),Q_mid(i,m+l),Q_mid(i,m+1-l),
C     &        Q_hi(i,m)
        end do
      end do
*
      return
      end

max_stars_repo_name: dinisAbranches/nwchem

max_issues_repo_name: dinisAbranches/nwchem

max_forks_repo_name: dinisAbranches/nwchem


text:
*DECK CFFTB
      SUBROUTINE CFFTB (N, C, WSAVE)
C***BEGIN PROLOGUE  CFFTB
C***SUBSIDIARY
C***PURPOSE  Compute the unnormalized inverse of CFFTF.
C***LIBRARY   SLATEC (FFTPACK)
C***CATEGORY  J1A2
C***TYPE      COMPLEX (RFFTB-S, CFFTB-C)
C***KEYWORDS  FFTPACK, FOURIER TRANSFORM
C***AUTHOR  Swarztrauber, P. N., (NCAR)
C***DESCRIPTION
C
C  ********************************************************************
C  *   NOTICE   NOTICE   NOTICE   NOTICE   NOTICE   NOTICE   NOTICE   *
C  ********************************************************************
C  *                                                                  *
C  *   This routine uses non-standard Fortran 77 constructs and will  *
C  *   be removed from the library at a future date.  You are         *
C  *   requested to use CFFTB1.                                       *
C  *                                                                  *
C  ********************************************************************
C
C  Subroutine CFFTB computes the backward complex discrete Fourier
C  transform (the Fourier synthesis).  Equivalently, CFFTB computes
C  a complex periodic sequence from its Fourier coefficients.
C  The transform is defined below at output parameter C.
C
C  A call of CFFTF followed by a call of CFFTB will multiply the
C  sequence by N.
C
C  The array WSAVE which is used by subroutine CFFTB must be
C  initialized by calling subroutine CFFTI(N,WSAVE).
C
C  Input Parameters
C
C  N       the length of the complex sequence C.  The method is
C          more efficient when N is the product of small primes.
C
C  C       a complex array of length N which contains the sequence
C
C  WSAVE   a real work array which must be dimensioned at least 4*N+15
C          in the program that calls CFFTB.  The WSAVE array must be
C          initialized by calling subroutine CFFTI(N,WSAVE), and a
C          different WSAVE array must be used for each different
C          value of N.  This initialization does not have to be
C          repeated so long as N remains unchanged.  Thus subsequent
C          transforms can be obtained faster than the first.
C          The same WSAVE array can be used by CFFTF and CFFTB.
C
C  Output Parameters
C
C  C       For J=1,...,N
C
C              C(J)=the sum from K=1,...,N of
C
C                 C(K)*EXP(I*(J-1)*(K-1)*2*PI/N)
C
C                         where I=SQRT(-1)
C
C  WSAVE   contains initialization calculations which must not be
C          destroyed between calls of subroutine CFFTF or CFFTB
C
C***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
C                 Computations (G. Rodrigue, ed.), Academic Press,
C                 1982, pp. 51-83.
C***ROUTINES CALLED  CFFTB1
C***REVISION HISTORY  (YYMMDD)
C   790601  DATE WRITTEN
C   830401  Modified to use SLATEC library source file format.
C   860115  Modified by Ron Boisvert to adhere to Fortran 77 by
C           changing dummy array size declarations (1) to (*).
C   861211  REVISION DATE from Version 3.2
C   881128  Modified by Dick Valent to meet prologue standards.
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900131  Routine changed from user-callable to subsidiary
C           because of non-standard Fortran 77 arguments in the
C           call to CFFTB1.  (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  CFFTB
      COMPLEX C
      DIMENSION C(*), WSAVE(*)
C***FIRST EXECUTABLE STATEMENT  CFFTB
      IF (N .EQ. 1) RETURN
      IW1 = N+N+1
      IW2 = IW1+N+N
      CALL CFFTB1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END


!-------------------------------------------------------------------------------
!> module ATMOSPHERE / Physics Cloud Microphysics / SDM
!!
!! @par Description
!!          Coversion between atmospheric fluid variables
!!
!! - Reference
!!  - Shima et al., 2009:
!!    The super-droplet method for the numerical simulation of clouds and precipitation:
!!    A particle-based and probabilistic microphysics model coupled with a non-hydrostatic model.
!!    Quart. J. Roy. Meteorol. Soc., 135: 1307-1320
!!
!! @author Team SCALE
!!
!! @par History
!! @li      2014-07-11 (S.Shima) [new] Separated from scale_atmos_phy_mp_sdm.F90
!! @li      2014-07-14 (S.Shima) [rev] sdm_rhot_qtrc2cpexnr added
!! @li      2014-07-24 (Y.Sato)  [mod] Modify bugs accessing upper/lower boundary
!! @li      2014-07-25 (Y.Sato)  [mod] Modify a bug in sdm_rho_mom2uvw
!! @li      2018-06-05 (S.Shima) [fix] Exner function in sdm_rhot_qtrc2cpexnr
!!
!<
!-------------------------------------------------------------------------------
module m_sdm_fluidconv

  implicit none
  private
  public :: sdm_rhot_qtrc2p_t, sdm_rho_rhot2pt, sdm_rho_mom2uvw, sdm_rho_qtrc2rhod,sdm_rhot_qtrc2cpexnr

contains
  !-----------------------------------------------------------------------------
  subroutine sdm_rho_qtrc2rhod(dens,qtrc,rhod)
    use scale_precision
    use scale_grid_index, only: &
         IA,JA,KA
    use scale_tracer, only: &
         QAD => QA, QQS, QQE, I_QV
    ! Input variables
    real(RP), intent(in) :: dens(KA,IA,JA)        !! Density [kg/m3]
    real(RP), intent(in) :: qtrc(KA,IA,JA,QAD)    !! Ratio of mass of tracer to total mass[kg/kg]
    ! Output variables
    real(RP), intent(out) :: rhod(KA,IA,JA)  ! Density of dry air [kg/m3]

    real(RP) :: qdry(KA,IA,JA)
    integer :: i, j, k,iq  ! index
    !---------------------------------------------------------------------

    qdry(:,:,:) = 1.0_RP
    do j = 1, JA
    do i = 1, IA
    do k = 1, KA
!!$      do iq = QQS, QQE
!!$        qdry(k,i,j) = qdry(k,i,j) - qtrc(k,i,j,iq)
!!$      enddo
       ! rho = rho_d + rho_v  when the SDM is used
       qdry(k,i,j) = qdry(k,i,j) - qtrc(k,i,j,I_QV)
       rhod(k,i,j) = dens(k,i,j)*qdry(k,i,j)
    enddo
    enddo
    enddo

    return
  end subroutine sdm_rho_qtrc2rhod
  !-----------------------------------------------------------------------------
  subroutine sdm_rho_mom2uvw(dens,momx,momy,momz,u,v,w)
    use scale_precision
    use scale_grid_index, only: &
         IA,IS,IE,JA,JS,JE,KA,KS,KE
    ! Input variables
    real(RP), intent(in)  :: dens(KA,IA,JA) ! Density [kg/m3]
    real(RP), intent(in)  :: momx(KA,IA,JA) ! Momentum [kg/s/m2]
    real(RP), intent(in)  :: momy(KA,IA,JA)
    real(RP), intent(in)  :: momz(KA,IA,JA)
    ! Output variables
    real(RP), intent(out) :: u(KA,IA,JA)  ! wind velocity [m/s]
    real(RP), intent(out) :: v(KA,IA,JA)  !
    real(RP), intent(out) :: w(KA,IA,JA)  !

    integer :: i, j, k  ! index
    !---------------------------------------------------------------------

    do j = JS-1, JE
    do i = IS-1, IE
       do k = KS, KE
          u(k,i,j) = 2.0_RP * momx(k,i,j) / ( dens(k,i,j)+dens(k,i+1,j) )
          v(k,i,j) = 2.0_RP * momy(k,i,j) / ( dens(k,i,j)+dens(k,i,j+1) )
       enddo
       u(1:KS-1,i,j)  = u(KS,i,j)
       v(1:KS-1,i,j)  = v(KS,i,j)
       u(KE+1:KA,i,j)   = u(KE,i,j)
       v(KE+1:KA,i,j)   = v(KE,i,j)
    enddo
    enddo

    do j = JS-1, JE
    do i = IS-1, IE
       do k = KS, KE
          w(k,i,j) = 2.0_RP * momz(k,i,j) / ( dens(k,i,j)+dens(k+1,i,j) )
       enddo
       w(1:KS-1,i,j) = 0
       w(KE:KA,i,j) = 0
    enddo
    enddo

    return
  end subroutine sdm_rho_mom2uvw
  !-----------------------------------------------------------------------------
  subroutine sdm_rho_rhot2pt(dens,rhot,pt)
    use scale_precision
    use scale_grid_index, only: &
         IA,JA,KA
    ! Input variables
    real(RP), intent(in) :: dens(KA,IA,JA)        !! Density [kg/m3]
    real(RP), intent(in) :: rhot(KA,IA,JA)        !! DENS * POTT [K*kg/m3]
    ! Output variables
    real(RP), intent(out) :: pt(KA,IA,JA)         !  Potential temperature [K]

    integer :: i, j, k  ! index
    !---------------------------------------------------------------------

    do j = 1, JA
    do i = 1, IA
    do k = 1, KA
       pt(k,i,j) = RHOT(k,i,j) / DENS(k,i,j)
    enddo
    enddo
    enddo

    return
  end subroutine sdm_rho_rhot2pt
  !-----------------------------------------------------------------------------
  subroutine sdm_rhot_qtrc2p_t(rhot,qtrc,dens,p,t)
    use scale_precision
    use scale_grid_index, only: &
         IA,JA,KA
    use scale_const, only: &
         P00 => CONST_PRE00, &
         Rdry => CONST_Rdry, &
         Rvap => CONST_Rvap, &
         CPdry => CONST_CPdry
    use scale_atmos_thermodyn, only: &
         CPw => AQ_CP
    use scale_tracer, only: &
         QAD => QA, QQS, QQE, I_QV
    ! Input variables
    real(RP), intent(in) :: rhot(KA,IA,JA)        !! DENS * POTT [K*kg/m3]
    real(RP), intent(in) :: qtrc(KA,IA,JA,QAD)    !! Ratio of mass of tracer to total mass[kg/kg]
    real(RP), intent(in) :: dens(KA,IA,JA)        !! Density [kg/m3]
    ! Output variables
    real(RP), intent(out) :: p(KA,IA,JA)          !  Pressure [Pa]
    real(RP), intent(out) :: t(KA,IA,JA)          !  Temperature [K]

    real(RP) :: qdry(KA,IA,JA)
    real(RP) :: rtot(KA,IA,JA)
    real(RP) :: cptot(KA,IA,JA)
    real(RP) :: cpovcv(KA,IA,JA)
    integer :: i, j, k, iq  ! index
    !---------------------------------------------------------------------

    qdry(:,:,:) = 1.0_RP
    do j = 1, JA
    do i = 1, IA
    do k = 1, KA
!!$      do iq = QQS, QQE
!!$        qdry(k,i,j) = qdry(k,i,j) - qtrc(k,i,j,iq)
!!$      enddo
      ! rho = rho_d + rho_v  when the SDM is used
      qdry(k,i,j) = qdry(k,i,j) - qtrc(k,i,j,I_QV)
      rtot (k,i,j) = Rdry * qdry(k,i,j) + Rvap * qtrc(k,i,j,I_QV)
      cptot(k,i,j) = CPdry * qdry(k,i,j)
!!$      do iq = QQS, QQE
!!$        cptot(k,i,j) = cptot(k,i,j) + qtrc(k,i,j,iq) * CPw(iq)
!!$      enddo
      cptot(k,i,j) = cptot(k,i,j) + qtrc(k,i,j,I_QV) * CPw(I_QV)
      cpovcv(k,i,j) = cptot(k,i,j) / ( cptot(k,i,j) - rtot(k,i,j) )

      p(k,i,j) = P00 * ( rhot(k,i,j) * rtot(k,i,j) / P00 )**cpovcv(k,i,j)
      t(k,i,j) = (rhot(k,i,j)/dens(k,i,j)) * (p(k,i,j)/P00)**(rtot(k,i,j)/cptot(k,i,j))
    enddo
    enddo
    enddo

    return
  end subroutine sdm_rhot_qtrc2p_t
  !---------------------------------------------------------------------
  subroutine sdm_rhot_qtrc2cpexnr(rhot,qtrc,dens,cpexnr)
    use scale_precision
    use scale_grid_index, only: &
         IA,JA,KA
    use scale_const, only: &
         P00 => CONST_PRE00, &
         Rdry => CONST_Rdry, &
         Rvap => CONST_Rvap, &
         CPdry => CONST_CPdry
    use scale_atmos_thermodyn, only: &
         CPw => AQ_CP
    use scale_tracer, only: &
         I_QV, QAD => QA
    ! Input variables
    real(RP), intent(in) :: rhot(KA,IA,JA)        !! DENS * POTT [K*kg/m3]
    real(RP), intent(in) :: qtrc(KA,IA,JA,QAD)    !! Ratio of mass of tracer to total mass[kg/kg]
    real(RP), intent(in) :: dens(KA,IA,JA)        !! Density [kg/m3]
    ! Output variables
    real(RP), intent(out) :: cpexnr(KA,IA,JA)     !  cp*exner_function

    real(RP) :: qdry(KA,IA,JA)
    real(RP) :: rtot(KA,IA,JA)
    real(RP) :: cptot(KA,IA,JA)
    integer :: i, j, k, iq  ! index
    !---------------------------------------------------------------------

    qdry(:,:,:) = 1.0_RP
    do j = 1, JA
    do i = 1, IA
    do k = 1, KA
       qdry(k,i,j) = qdry(k,i,j) - qtrc(k,i,j,I_QV)
       rtot (k,i,j) = Rdry * qdry(k,i,j) + Rvap * qtrc(k,i,j,I_QV)
       cptot(k,i,j) = CPdry * qdry(k,i,j)
       cptot(k,i,j) = cptot(k,i,j) + qtrc(k,i,j,I_QV) * CPw(I_QV)
       cpexnr(k,i,j) = cptot(k,i,j)*(rhot(k,i,j)*rtot(k,i,j)/P00)**(rtot(k,i,j)/(cptot(k,i,j)-rtot(k,i,j)))
    enddo
    enddo
    enddo

    return
  end subroutine sdm_rhot_qtrc2cpexnr
end module m_sdm_fluidconv


text:
      program prog
      save
      parameter(num1=2)
      parameter(num2=num1**2)
      dimension nums(num2)
      write(6, *) nums
      end


!-----------------------------------------------------------------------
!     Make center of mass motion off
!-----------------------------------------------------------------------

SUBROUTINE vcom( Velocity )

  use atom_module, only: Natom,ki_atom,zn_atom
  use cpmd_variables, only: pmass,AMU
  use calc_kine_temp_module, only: calc_kine

  implicit none
  real(8),intent(INOUT) :: Velocity(3,Natom)
  real(8) rescale,pm,kine0,kine,MG,VG(3)
  integer i

  call calc_kine( Velocity, kine0 )

  MG=0.0d0
  VG=0.0d0

  do i=1,Natom
     pm = pmass( zn_atom(ki_atom(i)) )*AMU
     MG = MG + pm
     VG(1) = VG(1) + Velocity(1,i)*pm
     VG(2) = VG(2) + Velocity(2,i)*pm
     VG(3) = VG(3) + Velocity(3,i)*pm
  end do

  VG(1:3)=VG(1:3)/MG
  do i=1,Natom
     Velocity(1,i) = Velocity(1,i) - VG(1)
     Velocity(2,i) = Velocity(2,i) - VG(2)
     Velocity(3,i) = Velocity(3,i) - VG(3)
  end do

  call calc_kine( Velocity, kine )

  rescale = sqrt( kine0/kine )

  Velocity(:,:) = Velocity(:,:)*rescale

END SUBROUTINE vcom


text:
module graph_m
use constants_m

implicit none

type graph_t
  integer::mblk=0 ! number of blocks
  integer::mrow=0 ! number of rows

  integer,dimension(:),pointer::nblk=>null() ! number of block entries
  integer,dimension(:),pointer::blk=>null() ! block start index

  integer,dimension(:),pointer::nrow=>null() ! number of row entries
  integer,dimension(:),pointer::row=>null() ! row start index

  integer,dimension(:),pointer::col=>null()

  integer,dimension(:),pointer::ival=>null()
  real(rp),dimension(:),pointer::rval=>null()
end type
contains

! C
subroutine graph_to_dual ( &
  m_grph, grph_ia, grph_ja, &
  m_dual, dual_ia, dual_ja)

  integer,intent(in)::m_grph
  integer,intent(out)::m_dual
  integer,dimension(:),intent(in)::grph_ia, grph_ja
  integer,dimension(:),pointer::dual_ia, dual_ja

  integer::i, jj1, jjn, jj, j, jjd;
  integer::nnz_grph, nnz_dual;
  integer,dimension(:),allocatable::cnt;

  nnz_grph = grph_ia(m_grph + 1) - 1;

  m_dual = 0;
  do i = 1, m_grph;
    jj1 = grph_ia(i);
    jjn = grph_ia(i + 1) - 1;
    do jj = jj1, jjn;
      j = grph_ja(jj);
      m_dual = max (j, m_dual);
    end do
  end do

  allocate (cnt(m_dual));
  cnt = 0;

  do i = 1, m_grph;
    jj1 = grph_ia(i);
    jjn = grph_ia(i + 1) - 1;
    do jj = jj1, jjn;
      j = grph_ja(jj);
      if (j > 0) cnt(j) = cnt(j) + 1;
    end do
  end do

  allocate (dual_ia(m_dual+1));
  allocate (dual_ja(nnz_grph));

  dual_ia(1) = 1;
  do i = 1, m_dual;
    dual_ia(i + 1) = dual_ia(i) + cnt(i);
    cnt(i) = 0;
  end do

  do i = 1, m_grph;
    jj1 = grph_ia(i);
    jjn = grph_ia(i + 1) - 1;
    do jj = jj1, jjn;
      j = grph_ja(jj);
      if (j > 0) then
        jjd = dual_ia(j) + cnt(j);
        dual_ja(jjd) = i;
        cnt(j) = cnt(j) + 1;
      end if
    end do
  end do
end subroutine

! C
subroutine print_range (name, mrow, row)
  character(*),intent(in)::name
  integer,intent(in)::mrow
  integer,dimension(:),intent(in)::row

  integer::i,jj1,jjn,jj,j

  write (*,"(/a)",advance="yes") trim(name)//" graph:"
  do i = 1, mrow
    jj1 = row(i)
    jjn = row(i + 1) - 1
    write (*,"(a,i0,1x,i0,a,i0)",advance="no") &
         "row [", jj1, jjn, "] ",i
    write (*,"(a)",advance="yes") ""
  end do

  write (*,"(a)",advance="yes") ""
end subroutine

! C
subroutine print_graph ( &
  name, mrow, row, col, idx, val)

  character(*),intent(in)::name
  integer,intent(in)::mrow
  integer,dimension(:),intent(in)::row, col
  integer,dimension(:),intent(in),optional::idx
  real(rp),dimension(:),intent(in),optional::val

  integer::i,jj1,jjn,jj,j

  write (*,"(/a)",advance="yes") trim(name)//" graph:"
  do i = 1, mrow
    jj1 = row(i)
    jjn = row(i + 1) - 1
    if (present(idx)) then
      write (*,"(a,i0,1x,i0,a,i0,a,i0,a)",advance="no") &
        "row [", jj1, jjn, "] ",i," -> ",idx(i),", col("
    else
      write (*,"(a,i0,1x,i0,a,i0,a)",advance="no") &
        "row [", jj1, jjn, "] ",i,", col("
    end if
    do jj = jj1, jjn
      j = col(jj)
      if (present(val)) then
        write (*,"(1x,g12.4)",advance="no") val(jj)
      else
        write (*,"(1x,i0)",advance="no") j
      end if
    end do
    write (*,"(a)",advance="yes") ")"
  end do

!   write (*,"(a/)",advance="yes") "end of "//trim(name)//" graph"
  write (*,"(a)",advance="yes") ""
end subroutine

subroutine print_blockgraph ( &
  name, mblk, blk, row, col)

  character(*),intent(in)::name
  integer,intent(in)::mblk
  integer,dimension(:),intent(in)::blk,row, col

  integer::i,jj1,jjn,jj,j,kk1,kkn,kk,k

  write (*,"(/a)",advance="yes") trim(name)//" block graph:"

  do i = 1, mblk
!     ty = mesh%fmt_ty(i)
!     dms = mesh%fmt_dim(i)
    jj1 = blk(i)
    jjn = blk(i + 1) - 1
    write (*,"(/a,1x,i0)",advance="yes") "block",i
    j = 0
    do jj = jj1, jjn
      j = j + 1
      kk1 = row(jj)
      kkn = row(jj + 1) - 1
      write (*,"(a,1x,i0,a)",advance="no") "row",j,", col("
      do kk = kk1, kkn
        k = col(kk)
        write (*,"(1x,i0)",advance="no") k
      end do
      write (*,"(a)",advance="yes") ")"
    end do
  end do

!   write (*,"(a/)",advance="yes") "end of "//trim(name)//" block graph"
  write (*,"(a)",advance="yes") ""
end subroutine
end module

max_stars_repo_name: DominicJones/pcns

max_issues_repo_name: DominicJones/pcns

max_forks_repo_name: DominicJones/pcns


text:
Love em or I Hate Davis hate em, Davis Drivers are an unusual breed that make driving in Davis an experience.

Some go slow.  Some go fast.
They suck because:
  Nobody seems to pay attention to the traffic signal. Seriously. Almost every time the driver in front of me would sit at a green light for ten seconds.
  They slow down to nearly a complete stop before making a turn.
  Theyre soft.  They let bikers get away with running stop signs and crap instead of running them over, jumping out of the car, and yelling HOW DO YOU LIKE ME NOW?! at the bleeding excyclist.
  And when they do hit you, they panic and ask if youre okay and call an ambulance instead of letting natural selection take its course. Way to keep the gene pool weak, you hippies.
  They use Highway 113 as a surface road, and drive the same speed as they would on a surface road.  The speed limit is 65, not 35, jackass!
  They put too many bumper stickers on their cars. Screw you, coexisting is for the weak.
  Theyre great at going the wrong way in singledirection aisles in parking lots.
  They almost manage to drive right down the dead center of any street that doesnt have a center line.
  Many drive those electric or hybrid vehicles that are too quiet to be safe.
  They come to a stop on turns which GIVE THEM A WHOLE LANE TO THEMSELVES before making them merge
  They dont look for motorcycles motorcycles or scooters and therefore frequently violate rightofway.
Theyre awesome because:
  http://i55.tinypic.com/erk50l.jpg They dont give a frak.
  They signal
  They (usually) understand rights of way
Davis Driver Conversations
On rare occasions it’s possible to communicate with a Davis driver by means other than thought projection or yelling.

  Users/Dozer was minding his own business driving his vehicle one day, when he came across a Davis driver who had outrageously run a stop sign. He was “shocked”.  After a while he came to post a negative comment on the wiki.  You see, the Davis driver had the Vehicles With Signs name of the company on his vehicle.  Dozer was http://daviswiki.org/Friendly_Cab?actiondiff&version2113&version1112 hot.

 The manager of the company saw the comment and was on it.  He http://daviswiki.org/Friendly_Cab?actiondiff&version2114&version1113 replied.

 The magic of the wiki was starting to calm the troubled waters.  Still shocked, Dozerhttp://daviswiki.org/Friendly_Cab?actiondiff&version2116&version1115 responded.

 Memory of the incident had started to return.  It could be http://daviswiki.org/Friendly_Cab?actiondiff&version2120&version1119 explained.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Describe Users/albertaramayo here.

20100729 23:15:54 nbsp Hey Albert, Welcome to the Wiki, Ive noticed that the updating and smoothing out of the Davis College Democrats DCD page hasnt quite gone as planned, but Im sure if we work with the community we can reach some good compromises. Because weve posted information on the site, its now a part of the whole, and DCD no longer owns that writing. It might also be helpful to look at the Welcome to the Wiki/Business Owner section... Users/MaxMikalonis


text:
      PROGRAM xgammln
C     driver for routine gammln
      INTEGER i,nval
      REAL actual,calc,gammln,x
      CHARACTER text*14
      open(7,file='FNCVAL.DAT',status='OLD')
10    read(7,'(a)') text
      if (text.ne.'Gamma Function') goto 10
      read(7,*) nval
      write(*,*) 'Log of gamma function:'
      write(*,'(1x,t11,a1,t24,a6,t40,a10)')
     *     'X','Actual','GAMMLN(X)'
      do 11 i=1,nval
        read(7,*) x,actual
        if (x.gt.0.0) then
          if (x.ge.1.0) then
            calc=gammln(x)
          else
            calc=gammln(x+1.0)-log(x)
          endif
          write(*,'(f12.2,2f18.6)') x,log(actual),calc
        endif
11    continue
      close(7)
      END

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xgammln.for


text:
cdis
cdis    Open Source License/Disclaimer, Forecast Systems Laboratory
cdis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
cdis
cdis    This software is distributed under the Open Source Definition,
cdis    which may be found at http://www.opensource.org/osd.html.
cdis
cdis    In particular, redistribution and use in source and binary forms,
cdis    with or without modification, are permitted provided that the
cdis    following conditions are met:
cdis
cdis    - Redistributions of source code must retain this notice, this
cdis    list of conditions and the following disclaimer.
cdis
cdis    - Redistributions in binary form must provide access to this
cdis    notice, this list of conditions and the following disclaimer, and
cdis    the underlying source code.
cdis
cdis    - All modifications to this software must be clearly documented,
cdis    and are solely the responsibility of the agent making the
cdis    modifications.
cdis
cdis    - If significant modifications or enhancements are made to this
cdis    software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
cdis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
cdis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
cdis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
cdis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
cdis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
cdis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
cdis
cdis
cdis
cdis
cdis
        subroutine plot_obs(k_level,l_ask_questions,asc9_tim
     1    ,i_radar_start,i_radar_end,namelist_parms,plot_parms
     1    ,imax,jmax,kmax,n_plotted,grid_ra_ref,grid_ra_vel,lat,lon
     1    ,topo,grid_spacing_m,mode)

!       Steve A         Nov  1989       Original Version
!       Steve A         Nov  1991       Adjustable Dimensions

        include 'trigd.inc'
        include 'lapsplot.inc'

        common /zoom/ zoom

        real grid_ra_ref(imax,jmax,kmax)
        real grid_ra_vel(imax,jmax,kmax)
        real lat(imax,jmax)
        real lon(imax,jmax)
        real topo(imax,jmax)

        real aspect_a(imax,jmax)          ! local

        integer n_plotted(imax,jmax)

        real riob_a(1000000)
        real rjob_a(1000000)
        logical l_plot, l_plot_thisob

        character*150 directory
        character*31 ext,c3_obsext
        character*10  units_2d
        character*125 comment_2d
        character*13 filename13
        character*3 var_2d
        character*1 c1_plottype
        character*8 c8_obstype

        character*13 fileradar
        character*24 atime_24
        character*9 asc9_tim,asc9_tim_radar
        logical l_ask_questions

        data init/0/

        Real
     1  azimuth_deg,
     1  X,
     1  Y,
     1  spd_kt,
     1  SPEED_ms,
     1  DIR,
     1  DU,
     1  PROJROT,
     1  range_km,
     1  pix_per_km,
     1  mspkt

        data
     1  PROJROT/0./,
     1  mspkt/.518/

        Character
     1  c_obs_type*1,c_map*1,c_mode*1,c_anl*1,c_radial*1,c_pin*1,
     1  c_sfcpro*1,STRING,wx*25

        common /plotobs/ c_obs_type,c_map,c_mode,c_anl,c_radial

        integer               vert_rad_pirep,
     1                          vert_rad_sao,
     1                          vert_rad_meso,
     1                          vert_rad_prof

!       *** set up constants

        call get_r_missing_data(r_missing_data, istatus)
        if(istatus .ne. 1)then
            write(6,*)' Error in plot_obs'
        endif

        call get_vert_rads (    vert_rad_pirep,
     1                          vert_rad_sao,
     1                          vert_rad_meso,
     1                          vert_rad_prof,
     1                          istatus)
        if(istatus .ne. 1)then
            write(6,*)' Error in plot_obs'
        endif

!       mode (1: call from 'wd'/'ob' options, 2: call from 'rd'/'vi')

        grid_scale = 3.1
        pix_per_km = 1.65

        x0 = 170.
        y0 = 245.
        pix_per_km = 1.37
        dubase = 2.4

        i_low_offset  = -nint(pix_per_km * 1.51 * grid_scale)
        i_high_offset = i_low_offset
     1                + nint(pix_per_km * 3.02 * grid_scale)
        dusmall = dubase * pix_per_km/1.65

        jmax_ref = 209
        size_factor = (float(jmax_ref-1) / 300.)

        if(imax .gt. jmax)then ! horizontally oriented box
            size_factor = size_factor / (float(imax)/float(jmax))
        endif

        size_factor = size_factor * plot_parms%obs_size

        write(6,*)' Subroutine plot_obs: size_factor = ',size_factor

        size_prof = 3.   * size_factor / zoom
        size_pirep = 3.0 * size_factor / zoom
        size_maps = 2.   * size_factor / zoom
        size_vad = 2.    * size_factor / zoom
        size_anl = 0.7   * size_factor / zoom
        size_suw = 1.    * size_factor / zoom
        size_radar = 1.  * size_factor
        size_meso = 2.   * size_factor / zoom

        aspect = 1.0 ! initialize to default value

        if(namelist_parms%l_sphere)then
            do i = 1,imax
            do j = 1,jmax
                arg = cosd(lat(i,j))
                ratio_log = nint(log(arg) / log(0.5))
                projfrac = 0.5**ratio_log

                interval_i = nint(float(interval) / projfrac)

                if(arg .gt. 0.)then
                    aspect_a(i,j) = 1.0 / arg
                else
                    aspect_a(i,j) = 1.0
                endif

            enddo ! j
            enddo ! i

        else
            aspect_a = 1.0

        endif

        LUN_IN = 5

        call getset(mxa,mxb,mya,myb,umin,umax,vmin,vmax,ltype)

        write(6,*)' Subroutine plot_obs at ',asc9_tim,': Level '
     1           ,k_level

        call setusv_dum(2hIN,201)

!       Plot Radar Obs   *********************************************************************

        write(6,95)
95      format(' Want Radial Velocities     [y,n,a,<RET>=y] ? ',$)
        if(l_ask_questions)read(lun_in,211)c_radial

        if(c_radial(1:1) .eq. 'n')goto205

        write(6,*)' Skipping call to ht_of_level to get retheight'
!       retheight = ht_of_level(k_level)

        c1_plottype = 'y'

        do i_radar = i_radar_start,i_radar_end

          write(6,*)' mode / radar = ',mode,i_radar

          if(mode .eq. 1)then
            call cv_asc_i4time(asc9_tim,i4time_needed)

            if(i_radar .le. 9)then
                write(ext,151)i_radar
151             format('v0',i1)
            else
                write(ext,152)i_radar
152             format('v',i2)
            endif

            var_2d = 'VEL'

            write(6,*)' calling get_laps_3dgrid'

            call get_laps_3dgrid(i4time_needed,1200,i4time_found
     1          ,imax,jmax,kmax,ext,var_2d
     1          ,units_2d,comment_2d,grid_ra_vel,istatus)

            if(istatus.ne.1)goto221

            write(6,*)comment_2d

          endif

          do k_grid = 1,kmax

            if(k_grid .gt. k_level .and. k_level .gt. 0)goto221

            do j_grid = 1,jmax
            do i_grid = 1,imax

              if(k_level .eq. 0)then
                k_sfc =
     1              nint(height_to_zcoord(topo(i_grid,j_grid),istatus))
              else
                k_sfc = -99
              endif

              if(k_grid  .eq. k_level                   .OR.
     1           k_level .eq. 0 .and. k_sfc .eq. k_grid      )then

                if(grid_ra_vel(i_grid,j_grid,k_grid)
     1                          .ne. r_missing_data)then

                  n_plotted(i_grid,j_grid) = n_plotted(i_grid,j_grid)+1

                  call plot_vr(i_grid,j_grid
     1              ,grid_ra_vel(i_grid,j_grid,k_grid),imax,jmax
     1              ,c1_plottype,n_plotted(i_grid,j_grid))

                endif

              endif ! at the right level to plot

            enddo ! i_grid
            enddo ! j_grid

          enddo ! k_grid

221       continue

        enddo ! i_radar

        if(mode .eq. 2)return

!******  Derived Radar Obs ********************************

205     call cv_asc_i4time(asc9_tim,i4time)

        write(6,210)
210     format(' Derived radar obs  [r, Ret = None]',30x,'? ',$)
        if(l_ask_questions)read(lun_in,211)c_obs_type
211     format(a1)

!       if(c_obs_type(1:1) .ne. ' ')then
        if(c_obs_type(1:1) .eq. 'r')then

        write(6,*)
        write(6,*)' Derived Radar Obs'

        call setusv_dum(2hIN,5) ! Orange

        lun = 61

        do i_radar = i_radar_start,i_radar_end

          if(i_radar .le. 9)then
            write(ext,251)i_radar
 251        format('d0',i1)
          else
            write(ext,252)i_radar
 252        format('d',i2)
          endif

          call get_directory(ext,directory,len_dir)
          open(lun
     1        ,file=directory(1:len_dir)//filename13(i4time,ext(1:3))
     1        ,status='old',err=1300)

1211      read(61,*,end=1300)ri,rj,k,dir,speed_ms
!421      format(1x,f6.3,f8.3,i2,2f6.1)
          ri = ri + 1.
          rj = rj + 1.
          k_ob = k + 1.

          k_sfc =
     1    nint(height_to_zcoord(topo(nint(ri),nint(rj)),istatus))

          if(k_ob .eq. k_level .or.
     1       k_level .eq. 0 .and. k_sfc .eq. k_ob)then

!               write(6,421)alat,alon,k,dir,speed_ms
1322            format(' old ',2f8.3,2i5)
1323            format(' new ',2f8.3,10x,2i5)

                spd_kt = speed_ms / mspkt

                if(nint(ri) .ge. 1 .AND. nint(ri) .le. imax .AND.
     1             nint(rj) .ge. 1 .AND. nint(rj) .le. jmax      )then

                    aspect = aspect_a(nint(ri),nint(rj))

                else
                    aspect = 1.0

                endif

                write(6,*)nint(ri),nint(rj),dir,spd_kt,aspect

!               Note the 'false' passed as DXX winds are grid north
                call plot_windob(dir,spd_kt,ri,rj,lat,lon,imax,jmax
     1                          ,size_radar,aspect,'false')

          elseif(k_ob .gt. k_level)then
                goto1300

          endif ! k_ob .eq. k_level

          goto1211

1300      continue

        enddo ! i_radar

        endif

!       Plot SUW obs   ****************************************************************************************************************

        if(c_obs_type(1:1) .eq. 's')then

222         call setusv_dum(2hIN,201)

            lun = 32
            ext = 'suw'
            call get_directory(ext,directory,len_dir)
            open(lun
     1          ,file=directory(1:len_dir)//filename13(i4time,ext(1:3))
     1          ,status='old',err=240)

230         read(32,*,end=240)azimuth_deg,range_km,elev_deg,dir,speed_ms

            range_m = range_km * 1000.

            range_km = range_km * cosd(elev_deg)

            call radar_to_latlon(alat,alon,retheight
     1                      ,azimuth_deg,range_m,elev_deg
     1                      ,rlat_radar,rlon_radar,rheight_radar)

            k = nint(height_to_zcoord(retheight,istatus))

            if ( k .eq. k_level) then ! plot even range_km rings

              if ( range_km .gt. 20. )  then

!               call latlon_ram(alat,alon,x,y,x0,y0,pix_per_km)
!               call latlon_ram_laps(alat,alon,x,y,init,'p')

                spd_kt = SPEED_ms  / mspkt

c               write(6,112)elev_deg,k,range_km,azimuth_deg,dir,spd_kt
112             format(1x,f6.1,i4,2f7.0,4x,2f7.0,i4)

                if(nint(ri) .ge. 1 .AND. nint(ri) .le. imax .AND.
     1             nint(rj) .ge. 1 .AND. nint(rj) .le. jmax      )then

                    aspect = aspect_a(nint(ri),nint(rj))

                else
                    aspect = 1.0

                endif

                call plot_windob(dir,spd_kt,ri,rj,lat,lon
     1                          ,imax,jmax,size_suw,aspect,'true')
              endif

            end if ! Valid Wind

        goto230
240     continue

        endif

        dularge = dusmall * 2.

!       Plot VAD wind  *************************************************
        write(6,*)
        write(6,*)' Plotting VAD Wind'

        call setusv_dum(2hIN,229)

        lun = 21
        ext = 'vad'
        call get_directory(ext,directory,len_dir)
        open(lun,file=directory(1:len_dir)//filename13(i4time,ext(1:3))
     1  ,status='old',err=30)

        do while (.true.)
            read(21,*,end=30)el,dir,speed_ms

            retheight = el * 1000. + rheight_radar

            k = nint(height_to_zcoord(retheight,istatus))

            write(6,*)el,retheight,k

            if(speed_ms .lt. 90. .and. k .eq. k_level)then

                alat = rlat_radar
                alon = rlon_radar

!               call latlon_ram(alat,alon,x,y,x0,y0,pix_per_km)
!               call latlon_ram_laps(alat,alon,x,y,init,'p')

                spd_kt = SPEED_ms  / mspkt

                write(6,111)alat,alon,max(dir,-99.),spd_kt
111             format(1x,2f8.1,4x,f7.0,f7.0,i4,2f8.3)

                if(nint(ri) .ge. 1 .AND. nint(ri) .le. imax .AND.
     1             nint(rj) .ge. 1 .AND. nint(rj) .le. jmax      )then

                    aspect = aspect_a(nint(ri),nint(rj))

                else
                    aspect = 1.0

                endif

                call plot_windob(dir,spd_kt,ri,rj,lat,lon,imax,jmax
     1                          ,size_vad,aspect,'true')

            endif
        enddo

 30     continue

!       Plot Sfc/METAR winds  ***********************************************
        write(6,97)
97      format(' Plot SFC / Profiler Winds?     [y,n,<RET>=y] ? ',$)
        if(l_ask_questions)read(lun_in,211)c_sfcpro

        if(c_sfcpro(1:1) .ne. 'n')then

50          write(6,*)' Sfc/METAR Data'

            call setusv_dum(2hIN,14) ! RoyalBlue

            lun = 32
            ext = 'sag'
            call get_directory(ext,directory,len_dir)
            open(lun,file=directory(1:len_dir)//
     1           filename13(i4time,ext(1:3)),status='old',err=811)

55          read(32,*,end=60)ri,rj,rk,dir,speed_ms
            k = nint(rk)

            if(k_level .gt. 0)then
                wt_vert = weight_vertical(retheight,k_level,istatus)
            else
                wt_vert = 0.
            endif

            if(abs(k - k_level) .le. vert_rad_sao .or. k_level .eq. 0
     1                                                           )then

                if(nint(ri) .ge. 1 .AND. nint(ri) .le. imax .AND.
     1             nint(rj) .ge. 1 .AND. nint(rj) .le. jmax      )then

                    aspect = aspect_a(nint(ri),nint(rj))

                else
                    aspect = 1.0

                endif

                write(6,111)ri,rj,max(dir,-99.),speed_ms
     1                     ,k,wt_vert,aspect

                if(k .eq. k_level)then
                    call setusv_dum(2hIN,12) ! Aqua
                else
                    call setusv_dum(2hIN,15) ! Slate Blue
                endif

                spd_kt = SPEED_ms  / mspkt

                call plot_windob(dir,spd_kt,ri,rj,lat,lon,imax,jmax
     1                          ,size_meso,aspect,'true')

            endif ! k .eq. k_level

            goto55

60          continue

            close(32)

            call getset(mxa,mxb,mya,myb,umin,umax,vmin,vmax,ltype)

!           Plot Profile winds  ***********************************************
811         write(6,*)
            write(6,*)
     1          ' Profile Winds (e.g. Profiler, Tower, Raob, Dropsonde)'

            lun = 32
            ext = 'prg'
            call get_directory(ext,directory,len_dir)
            open(lun,file=directory(1:len_dir)//
     1           filename13(i4time,ext(1:3)),status='old',err=911)

            do while (.true.)
                read(32,*,end=41,err=39)ri,rj,rk,dir,speed_ms,c8_obstype
 33             format(1x,5f10.0,1x,a8)

                k = nint(rk)
!               write(6,*)k,alat,alon,retheight,dir,speed_ms
!               k_sfc = nint(height_to_zcoord(topo(i,j),istatus))
                k_sfc = 2

                if(k .eq. k_level
     1    .or. k_level .eq. 0 .and. k .eq. k_sfc
     1   .and. dir .ne. r_missing_data
     1   .and. speed_ms .ne. r_missing_data
     1                                                  )then

!                   call latlon_ram(alat,alon,x,y,x0,y0,pix_per_km)
!                   call latlon_ram_laps(alat,alon,x,y,init,'p')

                    spd_kt = SPEED_ms  / mspkt

                    if(nint(ri) .ge. 1 .AND. nint(ri) .le. imax .AND.
     1                 nint(rj) .ge. 1 .AND. nint(rj) .le. jmax  )then

                        aspect = aspect_a(nint(ri),nint(rj))

                    else
                        aspect = 1.0

                    endif

                    if(c8_obstype(1:4) .eq. 'RAOB')then
                        call setusv_dum(2hIN,3)  ! Red
                    else
                        call setusv_dum(2hIN,17) ! Lavender
                    endif

                    call plot_windob(dir,spd_kt,ri,rj,lat,lon,imax,jmax
     1                              ,size_prof,aspect,'true')
                    write(6,311,err=38)ri,rj,dir,spd_kt,aspect
     1                                ,c8_obstype
311                 format(1x,2f8.1,4x,f7.0,f7.0,f8.3,1x,a8)
38                  continue

                endif ! k .eq. k_level

39              continue

            enddo

41          continue

            close(32)

        endif ! plot sfc and profiler winds

!       Plot Pirep winds  ***********************************************
911     write(6,96)
 96     format(' Plot ACARS / Pireps Winds?     [y,w,n,<RET>=y] ? ',$)
        if(l_ask_questions)read(lun_in,211)c_pin
        write(6,*)
        if(c_pin(1:1) .eq. 'n')then
            write(6,*)' Plotting WISDOM / Cloud Drift Winds'
            write(6,*)' ACARS / Pireps Winds - not plotted'
        elseif(c_pin(1:1) .eq. 'w')then
            write(6,*)' Plotting WISDOM Balloon Winds'
            write(6,*)
     1           ' ACARS / Pireps / Cloud Drift Winds - not plotted'
        else
            write(6,*)' ACARS / Pireps / WISDOM / Cloud Drift Winds'
        endif

        lun = 32
        ext = 'pig'
        call get_directory(ext,directory,len_dir)
        open(lun,file=directory(1:len_dir)//filename13(i4time,ext(1:3))
     1  ,status='old',err=1011)

        iob = 0 ! initialize for obs density routine
        dist_plot = namelist_parms%dist_plot_ua ! threshold in grid points

        do while (.true.)
            read(32,*,end=51)ri,rj,rk,dir,speed_ms,c3_obsext

            k = nint(rk)

            if(abs(rk - float(k_level)) .le. float(vert_rad_pirep))then

                if(k .eq. k_level)then
                    if(c3_obsext .eq. 'pin')then
                        call setusv_dum(2hIN,12) ! Aqua
                    elseif(c3_obsext .eq. 'wis')then
                        call setusv_dum(2hIN,5)  ! Orange
                    else
                        call setusv_dum(2hIN,8)  ! Green-Yellow
                    endif
                else
                    if(c3_obsext .eq. 'pin')then
                        call setusv_dum(2hIN,15) ! Slate Blue
                    elseif(c3_obsext .eq. 'wis')then
                        call setusv_dum(2hIN,26) ! Pale Orange
                    else
                        call setusv_dum(2hIN,34) ! Grey
                    endif
                endif

!               call latlon_ram(alat,alon,x,y,x0,y0,pix_per_km)
!               call latlon_ram_laps(alat,alon,x,y,init,'p')

                spd_kt = SPEED_ms  / mspkt

                if(nint(ri) .ge. 1 .AND. nint(ri) .le. imax .AND.
     1             nint(rj) .ge. 1 .AND. nint(rj) .le. jmax      )then

                    aspect = aspect_a(nint(ri),nint(rj))

                else
                    aspect = 1.0

                endif

!               Check if ACARS/Pireps are included
                if(c3_obsext .ne. 'pin' .OR. c_pin(1:1) .eq. 'y')then
                    l_plot_thisob = .true.
                else
                    l_plot_thisob = .false.
                endif

!               Option for plotting WISDOM only
                if(c_pin(1:1) .eq. 'w')then
                    if(c3_obsext .eq. 'wis')then
                        l_plot_thisob = .true.
                    else
                        l_plot_thisob = .false.
                    endif
                endif

                if(l_plot_thisob)then
                  call check_ob_density(riob_a,rjob_a,1000000
     1                                 ,iob,ri,rj,dist_plot,l_plot)

                  if(l_plot)then ! passes density criteria
                    iob = iob + 1
                    riob_a(iob) = ri
                    rjob_a(iob) = rj

                    write(6,921)ri,rj,rk,max(dir,-99.),spd_kt,c3_obsext
921                 format(1x,3f8.1,4x,f7.0,f7.0,2x,a3)
                    call plot_windob(dir,spd_kt,ri,rj,lat,lon,imax,jmax
     1                          ,size_pirep,aspect,'true')
                  else
                    write(6,922)ri,rj,rk,max(dir,-99.),spd_kt,c3_obsext
922                 format(1x,3f8.1,4x,f7.0,f7.0,2x,a3,' not plotted')
                  endif

                endif ! ACARS/Pireps are included

            endif ! k .eq. k_level

        enddo

51      continue

        close(32)

!       Plot MODEL winds  ***********************************************
930     write(6,*)
        write(6,*)' MODEL Winds'

        call setusv_dum(2hIN,203)

        lun = 32
        ext = 'mpg'
        call get_directory(ext,directory,len_dir)
        open(lun,file=directory(1:len_dir)//filename13(i4time,ext(1:3))
     1  ,status='old',err=1011)

        do while (.true.)

            read(32,*,end=932)
     1          k_grid,dum,dum,i_grid,j_grid,u_grid,v_grid

            alat = lat(i_grid,j_grid)
            alon = lon(i_grid,j_grid)

            if(k_grid .eq. k_level)then

!               call latlon_ram(alat,alon,x,y,x0,y0,pix_per_km)
!               call latlon_ram_laps(alat,alon,x,y,init,'p')
                call uvgrid_to_disptrue(u_grid,
     1                          v_grid,
     1                          dir,
     1                          speed_ms,
     1                          alat,
     1                          alon)

                spd_kt = SPEED_ms / mspkt

                if(nint(ri) .ge. 1 .AND. nint(ri) .le. imax .AND.
     1             nint(rj) .ge. 1 .AND. nint(rj) .le. jmax      )then

                    aspect = aspect_a(nint(ri),nint(rj))

                else
                    aspect = 1.0

                endif

!               write(6,111)alat,alon,max(dir,-99.),spd_kt,aspect
                call plot_windob(dir,spd_kt,ri,rj,lat,lon,imax,jmax
     1                          ,size_maps,aspect,'true')

            endif ! k .eq. k_level

        enddo ! i

932     continue

        close(32)

1011    continue

        return
        end

        subroutine plot_temp_obs(k_level,i4time,imax,jmax,kmax
     1                          ,r_missing_data,lat,lon,topo,zoom
     1                          ,plot_parms)

        include 'lapsplot.inc'

        character*1 c_pin
        character*3 ext
        character*8 c8_obstype
        character*150 directory
        character*150 filename
        character*13 filename13

        logical l_found_file

!       Plot Temperature Obs  ***********************************************

!       size_temp = 8. * float(max(imax,jmax)) / 300.
        size_temp = 1.1 ! 3.33

        write(6,*)
        write(6,*)' Plot Temperature Obs, size_temp = ',size_temp

        write(6,96)
 96     format(' Plot ACARS Temps?     [y,n,<RET>=y] ? ',$)
        read(5,211)c_pin
 211    format(a1)

        lun = 32
        ext = 'tmg'
        call get_directory(ext,directory,len_dir)

        filename = directory(1:len_dir)//filename13(i4time,ext(1:3))
        call s_len(filename,len_file)

        open(lun,file=filename(1:len_file),status='old',err=31)
        l_found_file = .true.
        go to 32

 31     write(6,*)' Could not open ',filename(1:len_file)
        l_found_file = .false.
        go to 42

 32     td = r_missing_data
        p = r_missing_data
        dir = r_missing_data
        spd_kt = r_missing_data
        gust = r_missing_data

        nobs_temp = 0

        do while (.true.) ! Count the temperature obs
            read(32,*,end=41,err=36)ri,rj,rk,t_k,c8_obstype
 36         continue

            k = nint(rk)
            k_sfc = 2

            if( k .eq. k_level
     1                      .OR.
     1         (k_level .eq. 0 .and. k .eq. k_sfc)  )then

              if(t_k .ne. r_missing_data)then

                nobs_temp = nobs_temp + 1

              endif ! t_k .ne. r_missing_data

            endif ! k .eq. k_level

        enddo

41      continue

        rewind(32)

42      write(6,*)' Number of temperature obs = ',nobs_temp

        do while (l_found_file) ! Plot the temp obs
            read(32,*,end=141,err=150)ri,rj,rk,t_k,c8_obstype
150         continue

            k = nint(rk)
            k_sfc = 2

            if( k .eq. k_level
     1                      .OR.
     1         (k_level .eq. 0 .and. k .eq. k_sfc)  )then

              if(t_k .ne. r_missing_data)then

                t_c = t_k - 273.15

!               spd_kt = SPEED_ms  / mspkt

                iflag = 3

                if(c8_obstype(1:3) .eq. 'RAS')then      ! RASS
                    icol_in = 12 ! Aqua
                elseif(c8_obstype(1:3) .eq. 'WIS')then  ! WISDOM Balloons
                    icol_in = 12 ! Aqua
                elseif(c8_obstype(1:3) .eq. 'RAO')then  ! RAOB
                    icol_in = 7  ! Yellow
                elseif(c8_obstype(1:3) .eq. 'RAD')then  ! Radiometer
                    icol_in = 5  ! Orange
                elseif(c8_obstype(1:3) .eq. 'DRO')then  ! Dropsonde
                    icol_in = 14 ! Royal Blue
                elseif(c8_obstype(1:3) .eq. 'TOW')then  ! Tower
                    icol_in = 14 ! Royal Blue
                elseif(c8_obstype(1:2) .eq. 'SA')then   ! SATSND
                    icol_in = 17 ! Lavender
                elseif(c8_obstype(1:4) .eq. 'GOES')then ! GOES Satellite
                    icol_in = 17 ! Lavender
                elseif(c8_obstype(1:4) .eq. 'POES')then ! POES Satellite
                    icol_in = 16 ! Dark violet
                elseif(c8_obstype(1:5) .eq. 'METAR')then ! METAR
                    icol_in = 11 ! Green
                else                                    ! ACARS
                    icol_in = 3  ! Red
                endif

                call setusv_dum(2hIN,icol_in)

                iflag_cv = 0

                if(icol_in .ne. 3 .or. c_pin .eq. 'y')then ! filter out ACARS?
                  call plot_mesoob(dir,spd_kt,gust,t_c,td,p,ri,rj
     1                            ,lat,lon,imax,jmax,size_temp
     1                            ,zoom,nobs_temp
     1                            ,icol_in,du_loc,wx
     1                            ,iflag,iflag_cv,namelist_parms
     1                            ,plot_parms)

                  write(6,111,err=121)ri,rj,t_c,c8_obstype
111               format(1x,3f8.1,1x,a8)
121               continue

                else
                  write(6,*)' access restricted for ',c8_obstype

                endif

              endif ! t_k .ne. r_missing_data

            endif ! k .eq. k_level

        enddo

141     continue

        close(32)

        return
        end

        subroutine plot_td_obs(k_level,i4time,imax,jmax,kmax
     1                        ,r_missing_data,lat,lon,topo,zoom
     1                        ,namelist_parms,plot_parms
     1                        ,k_mb,mode,field_2d,i_overlay)

        include 'lapsplot.inc'
        include 'icolors.inc'

        character*3 ext, t1
        character*5 pw
        character*8 c8_obstype
        character*150 directory
        character*150 filename
        character*13 filename13
        character*100 c_label
        character*10 c_staname

        real field_2d(imax,jmax)

        logical l_found_file

!       Plot Dewpoint Obs  ***********************************************

!       size_temp = 8. * float(max(imax,jmax)) / 300.
        size_temp = 1.1 ! 3.33

        write(6,*)

        if(mode .eq. 1)then
            write(6,*)' Plot Dewpoint Obs, size_td = ',size_temp
            call get_border(imax,jmax,x_1,x_2,y_1,y_2)
            call set(x_1,x_2,y_1,y_2,1.,float(imax),1.,float(jmax),1)
        elseif(mode .eq. 2)then
            write(6,*)' Plot SH Obs (from dewpoint), size = ',size_temp
            call get_border(imax,jmax,x_1,x_2,y_1,y_2)
            call set(x_1,x_2,y_1,y_2,1.,float(imax),1.,float(jmax),1)
        elseif(mode .ge. 3)then
            write(6,*)' Plot PW (GPS) Obs, mode/size = ',mode,size_temp
!           c_label = 'Integrated Vapor Obs (cm)'
            c_label = 'IWV Obs, Grid-Obs (cm)'
            i_overlay = i_overlay + 1
            call setusv_dum(2hIN,icolors(i_overlay))
            call set(.00,1.0,.00,1.0,.00,1.0,.00,1.0,1)
            call make_fnam_lp(i4time,asc_tim_9,istatus)
            call write_label_lplot(imax,jmax,c_label,asc_tim_9
     1                            ,plot_parms,namelist_parms
     1                            ,i_overlay,'hsect')
            call get_border(imax,jmax,x_1,x_2,y_1,y_2)
            call set(x_1,x_2,y_1,y_2,1.,float(imax),1.,float(jmax),1)
        endif

        lun = 32
        ext = 'hmg'
        call get_directory(ext,directory,len_dir)

        filename = directory(1:len_dir)//filename13(i4time,ext(1:3))
        call s_len(filename,len_file)

        open(lun,file=filename(1:len_file),status='old',err=31)
        l_found_file = .true.
        go to 32

 31     write(6,*)' Could not open ',filename(1:len_file)
        l_found_file = .false.
        go to 42

 32     continue
        p = r_missing_data
        dir = r_missing_data
        spd_kt = r_missing_data
        gust = r_missing_data

        nobs_temp = 0

        do while (.true.) ! Count the dewpoint obs
            read(32,*,end=41,err=36)ri,rj,rk,t_k,c8_obstype,c_staname
 36         continue

            k = nint(rk)
            k_sfc = 2

            if( k .eq. k_level
     1                      .OR.
     1         (k_level .eq. 0 .and. k .eq. k_sfc)  )then

              if(t_k .ne. r_missing_data)then

                nobs_temp = nobs_temp + 1

              endif ! t_k .ne. r_missing_data

            endif ! k .eq. k_level

        enddo

41      continue

        rewind(32)

42      write(6,*)' Number of dewpoint obs = ',nobs_temp

        do while (l_found_file) ! Plot the dewpoint obs
            if(mode .lt. 3)then
              read(32,*,end=141,err=150)ri,rj,rk,t_k,c8_obstype
     1                                              ,c_staname
            else
              read(32,*,end=141,err=150)ri,rj,rk,t_k,c8_obstype
            endif
150         continue

            k = nint(rk)
            k_sfc = 2

            if( k .eq. k_level
     1                      .OR.
     1         (k_level .eq. 0 .and. k .eq. k_sfc)  )then

              if(t_k .ne. r_missing_data)then

                if(mode .eq. 1)then
                  t_c = t_k - 273.15
                elseif(mode .eq. 2)then
                  t_c = t_k
                elseif(mode .eq. 3)then
                  if(t_k .lt. 0.)goto 140
                  t_c = t_k
                  c_staname = ' '

                  if(trim(c8_obstype) .ne. 'GPSTPW')then
                    write(6,*)' Skipping: ',ri,rj,rk,c8_obstype
                    goto 140
                  else
                    write(6,*)' Keeping:  ',ri,rj,rk,c8_obstype
                  endif
                elseif(mode .eq. 4)then
                  if(t_k .lt. 0.)goto 140
                  field_value = field_2d(nint(ri),nint(rj)) * 100.
                  t_c = field_value - t_k
                  c_staname = ' '

                  if(trim(c8_obstype) .ne. 'GPSTPW')then
                    write(6,*)' Skipping: ',ri,rj,rk,c8_obstype
                    goto 140
                  else
                    write(6,*)' Keeping:  ',ri,rj,rk,c8_obstype,t_k
     1                                     ,field_value,t_c
                  endif
                endif

!               spd_kt = SPEED_ms  / mspkt

                iflag = 3

                if(c8_obstype(1:3) .eq. 'RAS')then      ! RASS
                    icol_in = 12 ! Aqua
                elseif(c8_obstype(1:3) .eq. 'RAO')then  ! RAOB
                    icol_in = 7  ! Yellow
                elseif(c8_obstype(1:3) .eq. 'RAD')then  ! Radiometer
                    icol_in = 5  ! Orange
                elseif(c8_obstype(1:3) .eq. 'DRO')then  ! Dropsonde
                    icol_in = 14 ! Royal Blue
                elseif(c8_obstype(1:2) .eq. 'SA')then   ! SATSND
                    icol_in = 17 ! Lavender
                elseif(c8_obstype(1:4) .eq. 'GOES')then ! GOES Satellite
                    icol_in = 17 ! Lavender
                elseif(c8_obstype(1:4) .eq. 'POES')then ! POES Satellite
                    icol_in = 16 ! Dark violet
                elseif(c8_obstype(1:5) .eq. 'METAR')then ! METAR
                    icol_in = 11 ! Green
                else                                    ! ACARS / GPS
                    icol_in = 3  ! Red
                endif

                if(mode .lt. 3)then
                    call setusv_dum(2hIN,icol_in)
                endif

                if(mode .eq. 2)then ! convert from td to q
                    svp = es(t_c)
                    t_c = (svp / float(k_mb)) * 1000.
                endif

                obs_size = plot_parms%contour_line_width

                zoom_max = 1.5
                if(zoom .lt. zoom_max)then
                    if(nobs_temp .gt. 30)then
                        zoom_eff = 1.0                 ! smaller obs
                    else
                        zoom_eff = zoom / zoom_max     ! larger obs
                    endif
                else
                    zoom_eff = zoom / zoom_max         ! larger obs
                endif

                zoom_eff = zoom_eff / obs_size

                du = float(jmax) / 252.
                du2 = (du / zoom_eff) * 0.6
                xsta = ri
                ysta = rj
                ANGD = 0.
                CNTR = 0.
                charsize = .0040 / zoom_eff

!               Plot station name if it exists
                if(c_staname .ne. 'NONAME')then
                    call s_len(c_staname,len_name)
                    CALL PCHIQU(xsta, ysta-du2*3.5,
     1                          c_staname(1:len_name),
     1                          charsize,ANGD,CNTR)
                endif

!               Plot ob location (with a plus sign)
                if(mode .ge. 3 .OR.
     1             c8_obstype(1:3) .eq. 'RAO')then
                    call line(xsta,ysta+du2*0.5,xsta,ysta-du2*0.5)
                    call line(xsta+du2*0.5,ysta,xsta-du2*0.5,ysta)
                endif

!               Plot ob
                if(mode .lt. 3)then
                  write(t1,100,err=101) nint(t_c)
 100              format(i3)
 101              call left_justify(t1)
                  call s_len(t1,len_t1)
                  if(c8_obstype(1:3) .eq. 'RAO')then
                      xloc = xsta + du2*4.0
                  else
                      xloc = xsta
                  endif
                  CALL PCMEQU(xloc,ysta,t1(1:len_t1),charsize,ANGD,CNTR)
                  write(6,102,err=103)ri,rj,t_c,t1,c8_obstype,c_staname
 102              format(1x,3f8.1,1x,a8,1x,a10)
 103              continue
                else ! GPS PW
                  write(pw,110,err=111) t_k
 110              format(f5.2)
 111              call left_justify(pw)
                  CALL PCHIQU(xsta+du2*7.0,ysta+du2*3.0,trim(pw)
     1                       ,charsize,ANGD,CNTR)

                  if(mode .eq. 3)then
                    write(6,121,err=122)ri,rj,t_c,pw,c8_obstype
     1                                              ,c_staname
 121                format(1x,3f8.1,1x,a8,1x,a10)
 122                continue

                  elseif(mode .eq. 4)then ! plot difference ob
                    write(pw,110,err=131) t_c
 131                call left_justify(pw)
                    CALL PCHIQU(xsta+du2*7.0,ysta-du2*3.0,trim(pw)
     1                       ,charsize,ANGD,CNTR)
                    write(6,132,err=133)ri,rj,t_k,t_c,c8_obstype
     1                                               ,c_staname
 132                format(1x,4f8.1,1x,a8,1x,a10)
 133                continue

                  endif
                endif

!               write(6,111,err=121)ri,rj,t_c,t1,c8_obstype,c_staname
!11             format(1x,3f8.1,1x,a8,1x,a10)
!21             continue

              endif ! t_k .ne. r_missing_data

            endif ! k .eq. k_level

140     enddo

141     continue

        close(32)

        return
        end

        subroutine plot_td_sndobs(k_level,i4time,imax,jmax,kmax
     1                        ,r_missing_data,lat,lon,topo,zoom
     1                        ,plot_parms)

        include 'lapsplot.inc'

        integer max_snd,max_snd_levels
        parameter (max_snd = 30000)
        parameter (max_snd_levels = 200)

        real heights_3d(imax,jmax,kmax)
        real pres_3d(imax,jmax,kmax)
        real lat(imax,jmax)
        real lon(imax,jmax)
        real topo(imax,jmax)

        real dum1_3d(imax,jmax,kmax) ! dummy variable

!       Arrays returned from read_tdsnd
        real lat_tdsnd(max_snd)
        real lon_tdsnd(max_snd)
        real bias_htlow(max_snd)
        real ob_pr_td(max_snd,kmax) ! Vertically interpolated TSND temp
        real inst_err_tdsnd(max_snd)
        character*5 c5_name(max_snd)
        character*8 c8_sndtype(max_snd)

        write(6,*)' Subroutine plot_td_sndobs'

        i4_window_raob_file = 3600
        ilaps_cycle_time = 3600

        call get_pres_3d(i4time,imax,jmax,kmax,pres_3d,istatus)

!       Set up standard atmosphere for heights
        do k = 1,kmax
        do i = 1,imax
        do j = 1,jmax
            heights_3d(i,j,k) = psatoz(pres_3d(i,j,k)*.01)
        enddo ! j
        enddo ! i
        enddo ! k

!       Read dewpoint obs from HMG file (generating this if needed)
        call  read_tdsnd(i4time,heights_3d,dum1_3d,               ! Input
     1                   pres_3d,                                 ! Input
     1                   lat_tdsnd,lon_tdsnd,                     ! Output
     1                   lat,lon,                                 ! Input
     1                   max_snd,max_snd_levels,                  ! Input
     1                   ob_pr_td,inst_err_tdsnd,                 ! Output
     1                   c5_name,c8_sndtype,                      ! Output
     1                   .true.,.false.,                          ! Input
     1                   i4_window_raob_file,                     ! Input
     1                   bias_htlow,                              ! Output
     1                   n_tdsnd,                                 ! Output
     1                   ilaps_cycle_time,                        ! Input
     1                   imax,jmax,kmax,                          ! Input
     1                   r_missing_data,                          ! Input
     1                   istatus)                                 ! Output

        write(6,*)' back from read_tdsnd, # soundings = ',n_tdsnd

        return
        end

        subroutine check_ob_density(riob_a,rjob_a,maxobs
     1                             ,iob,ri,rj,dist_plot,l_plot)

!       Determine whether to plot an ob based on its distance to the previously
!       plotted obs

        real riob_a(maxobs) ! I/O
        real rjob_a(maxobs) ! I/O
        real dist_plot    ! I    number of gridpoints allowed between obs
        integer iob       ! I/O
        logical l_plot    ! O

        l_plot = .false.

        dist_plot_sq = dist_plot**2

        if(iob .eq. 0)then
            l_plot = .true.
        else
            l_plot = .true.
            do i = 1,iob
                dist_sq = (ri - riob_a(i))**2 + (rj - rjob_a(i))**2
                if(dist_sq .lt. dist_plot_sq)then
                    l_plot = .false.
                endif
            enddo ! i
        endif

        if(l_plot)then ! add this ob into the arrays of plotted obs
            iob = iob + 1
            riob_a(iob) = ri
            rjob_a(iob) = rj
        endif

        return
        end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
!***********************************************************************
!  Integrated Water Flow Model (IWFM)
!  Copyright (C) 2005-2021
!  State of California, Department of Water Resources
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU General Public License
!  as published by the Free Software Foundation; either version 2
!  of the License, or (at your option) any later version.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  (http://www.gnu.org/copyleft/gpl.html)
!
!  You should have received a copy of the GNU General Public License
!  along with this program; if not, write to the Free Software
!  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
!
!  For tecnical support, e-mail: IWFMtechsupport@water.ca.gov
!***********************************************************************
!
! This module uses DATE_AND_TIME function instead of CPU_TIME to compute the
! elapsed run time.  This is because CPU_TIME does not produce correct results
! when the executable is run in parallel.  DATE_AND_TIME seems to be more robust
! in sequential and parallel executions.
!
! Can Dogrul ; 1/22/08
!
!***********************************************************************
MODULE ProgramTimer
  IMPLICIT NONE

  LOGICAL,SAVE ::TimerStopped = .FALSE.
  LOGICAL,SAVE ::TimerStarted = .FALSE.
  INTEGER,SAVE ::StartTimeValues(8),EndTimeValues(8)

  PRIVATE
  PUBLIC :: StartTimer    ,&
            StopTimer     ,&
            GetRunTime    ,&
            TimerStopped

CONTAINS

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

! *************************************************************
! *************************************************************
! **** START TIMER
! *************************************************************
! *************************************************************
  SUBROUTINE StartTimer()

    CALL DATE_AND_TIME(VALUES = StartTimeValues)
    TimerStopped = .FALSE.
    TimerStarted = .TRUE.

  END SUBROUTINE StartTimer

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

! *************************************************************
! *************************************************************
! **** STOP TIMER
! *************************************************************
! *************************************************************
  SUBROUTINE StopTimer()

    CALL DATE_AND_TIME(VALUES = EndTimeValues)
    TimerStopped = .TRUE.

  END SUBROUTINE StopTimer

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

! *************************************************************
! *************************************************************
! **** GET RUN TIME
! *************************************************************
! *************************************************************
  SUBROUTINE GetRunTime(Hour,Minute,Second)
    INTEGER,INTENT(OUT) :: Hour,Minute
    REAL(8),INTENT(OUT) :: Second

    !Local variables
    INTEGER :: iJulianStart,iJulianEnd,iYear,iMonth,iDay
    REAL(8) :: RunTime,JulianStart,JulianEnd,rHour,rMin,rSec,rmSec

    !Initialize
    Hour   = 0
    Minute = 0
    Second = 0.0

    !If timer hasn't been started at all run time cannot be computed; return
    IF (.NOT. TimerStarted) RETURN

    !Start date and time in Julian days
    rHour        = REAL(StartTimeValues(5),8)
    rMin         = REAL(StartTimeValues(6),8)
    rSec         = REAL(StartTimeValues(7),8)
    rmSec        = REAL(StartTimeValues(8),8)
    iDay         = StartTimeValues(3)
    iMonth       = StartTimeValues(2)
    iYear        = StartTimeValues(1)
    iJulianStart = iDay-32075+1461*(iYear+4800+(iMonth-14)/12)/4+367*(iMonth-2-(iMonth-14)/12*12)/12-3*((iYear+4900+(iMonth-14)/12)/100)/4
    JulianStart  = REAL(iJulianStart,8) + rHour/24.0 + rMin/1440.0 + rSec/86400.0 + rmSec/86400000.0

    !End date and time in Julian days
    rHour      = REAL(EndTimeValues(5),8)
    rMin       = REAL(EndTimeValues(6),8)
    rSec       = REAL(EndTimeValues(7),8)
    rmSec      = REAL(EndTimeValues(8),8)
    iDay       = EndTimeValues(3)
    iMonth     = EndTimeValues(2)
    iYear      = EndTimeValues(1)
    iJulianEnd = iDay-32075+1461*(iYear+4800+(iMonth-14)/12)/4+367*(iMonth-2-(iMonth-14)/12*12)/12-3*((iYear+4900+(iMonth-14)/12)/100)/4
    JulianEnd  = REAL(iJulianEnd,8) + rHour/24.0 + rMin/1440.0 + rSec/86400.0 + rmSec/86400000.0

    RunTime = (JulianEnd-JulianStart)*86400.0  !In seconds
    Hour    = INT(RunTime/3600.0)
    Minute  = INT((RunTime-Hour*3600.0)/60.0)
    Second  = RunTime-Hour*3600.0-Minute*60.0

  END SUBROUTINE GetRunTime

  ! -------------------------------------------------------------
  ! --- CONVERT AN INTEGER TO TEXT
  ! -------------------------------------------------------------
  FUNCTION IntToText(Number) RESULT(Text)
    INTEGER,PARAMETER           :: MaxTextLength=20
    INTEGER,INTENT(IN)          :: Number
    CHARACTER(LEN=MaxTextLength):: Text

    !Local variables
    CHARACTER(LEN=6) :: FormatStatement
    WRITE (FormatStatement,'(A,I3,A)') '(I',MaxTextLength,')'

    WRITE (Text,FMT=FormatStatement) Number
    Text=ADJUSTL(Text)

  END FUNCTION IntToText

END MODULE

max_stars_repo_name: giorgk/IWFM-SAFE

max_issues_repo_name: giorgk/IWFM-SAFE

max_forks_repo_name: giorgk/IWFM-SAFE