text:
C     ==================================================================
      SUBROUTINE SFT1D( A, B, NTOT, N, NSPAN, ISN )
C     ==================================================================
C
C --- This routine calculates a one-dimensional complex multivariate
C --- Fast Fourier Transform (FFT) using a generalisation of the
C --- Cooley-Tukey algorithm due to SINGLETON,
C --- IEEE Trans Audio & Electroacoustics AU-17, 93-103 (1969).
C --- The routine works on all array sizes which factorise into PRIMES
C --- UP TO 23.
C
C --- The parameters are as follows:
C
C ---    A     a REAL array which contains the real components of the
C ---          data to be transformed.
C
C ---    B     a REAL array which contains the imaginary components of
C ---          the data to be transformed.
C
C ---    NTOT  the total number of complex values (see below).
C
C ---    N     the dimension of the current variable (see below).
C
C ---    NSPAN NSPAN/N defines the spacing of consecutive data values
C ---          while indexing the current variable (see below).
C
C ---    ISN   the sign of the transform kernel:
C ---                 ISN = -1 ==> exponent is negative (reverse).
C ---                 ISN = +1 ==> exponent is positive (forward).
C
C --- The meaning of NTOT, N and NSPAN can be understood from the
C --- following examples. (It is inportant to understand that Fortran
C --- arrays have element A(2,J,...) next to A(1,J,...) .)
C
C --- A single-variate transform with A(N) & B(N) requires one call:
C ---             CALL FFT( A, B, N, N, N, ISN )
C
C --- A bi-variate transform with A(N1,N2) & B(N1,N2) requires 2 calls:
C ---             CALL FFT( A, B, N1*N2, N1, N1,    ISN )
C ---             CALL FFT( A, B, N1*N2, N2, N1*N2, ISN )
C
C --- A tri-variate transform with A(N1,N2,N3) & B(N1,N2,N3) requires
C --- three calls:
C ---             CALL FFT( A, B, N1*N2*N3, N1, N1,       ISN )
C ---             CALL FFT( A, B, N1*N2*N3, N2, N1*N2,    ISN )
C ---             CALL FFT( A, B, N1*N2*N3, N3, N1*N2*N3, ISN )
C
C --- Note, however, that this is a very inefficient method of
C --- calculating multivariate (particularly 2-D) transforms.
C
C --- The transforms of COMPLEX arrays can also be calculated.
C --- If C is a complex array, the correct call is:
C ---             CALL FFT( C, C(2), NTOT, N, NSPAN, +/-2 )
C --- Note that the magnitude of ISN is two.
C --- The zero frequency is found in A(1),B(1) for A 1-D transform, etc.
C
C --- The error message "Array bounds exceeded"
C --- means one of two things:
C ---              the maximum prime factor is greater than 23.
C ---              the number of prime factors in N is greater than 209.
C
C --- This routine is modified from the program listing in Singletons's
C --- paper. The modifications are by Adrian Clark, Physics Dept.,
C --- Queen Elizabeth College, Campden Hill Road, London W8 7AH.
C
      DIMENSION A(1), B(1)
      DIMENSION NFAC(11), NP(209)
C
C --- The array storage allows a maximum prime of 23.
C
      DIMENSION AT(23), CK(23), BT(23), SK(23)
      EQUIVALENCE (I,II)
C
C --- The following two numbers should be equal to the array dimenions.
C
      MAXF=23
      MAXP=209
      IF(N .LT. 2) RETURN
      INC=ISN
      RAD=8.0*ATAN(1.0)
      S72=RAD/5.0
      C72=COS(S72)
      S72=SIN(S72)
      S120=SQRT(0.75)
      IF(ISN .GE. 0) GO TO 10
      S72=-S72
      S120=-S120
      RAD=-RAD
      INC=-INC
10    NT=INC*NTOT
      KS=INC*NSPAN
      KSPAN=KS
      NN=NT-INC
      JC=KS/N
      RADF=RAD*FLOAT(JC)*0.5
      I=0
      JF=0
C
C --- Determine the factors of N.
C
      M=0
      K=N
      GO TO 20
15    M=M+1
      NFAC(M)=4
      K=K/16
20    IF(K-(K/16)*16 .EQ. 0) GO TO 15
      J=3
      JJ=9
      GO TO 30
25    M=M+1
      NFAC(M)=J
      K=K/JJ
30    IF(MOD(K,JJ) .EQ. 0) GO TO 25
      J=J+2
      JJ=J**2
      IF(JJ .LE. K) GO TO 30
      IF(K .GT. 4) GO TO 40
      KT=M
      NFAC(M+1)=K
      IF(K .NE. 1) M=M+1
      GO TO 80
40    IF(K-(K/4)*4 .NE. 0) GO TO 50
      M=M+1
      NFAC(M)=2
      K=K/4
50    KT=M
      J=2
60    IF(MOD(K,J) .NE. 0) GO TO 70
      M=M+1
      NFAC(M)=J
      K=K/J
70    J=((J+1)/2)*2+1
      IF(J .LE. K) GO TO 60
80    IF(KT .EQ. 0) GO TO 100
      J=KT
90    M=M+1
      NFAC(M)=NFAC(J)
      J=J-1
      IF(J .NE. 0) GO TO 90
C
C --- Compute the Fourier transform.
C
100   SD=RADF/FLOAT(KSPAN)
      CD=2.0*SIN(SD)**2
      SD=SIN(SD+SD)
      KK=1
      I=I+1
      IF(NFAC(I) .NE. 2) GO TO 400
C
C --- Radix-2 transform.
C
      KSPAN=KSPAN/2
      K1=KSPAN+2
210   K2=KK+KSPAN
      AK=A(K2)
      BK=B(K2)
      A(K2)=A(KK)-AK
      B(K2)=B(KK)-BK
      A(KK)=A(KK)+AK
      B(KK)=B(KK)+BK
      KK=K2+KSPAN
      IF(KK .LE. NN) GO TO 210
      KK=KK-NN
      IF(KK .LE. JC) GO TO 210
      IF(KK .GT. KSPAN) GO TO 800
220   C1=1.0-CD
      S1=SD
230   K2=KK+KSPAN
      AK=A(KK)-A(K2)
      BK=B(KK)-B(K2)
      A(KK)=A(KK)+A(K2)
      B(KK)=B(KK)+B(K2)
      A(K2)=C1*AK-S1*BK
      B(K2)=S1*AK+C1*BK
      KK=K2+KSPAN
      IF(KK .LT. NT) GO TO 230
      K2=KK-NT
      C1=-C1
      KK=K1-K2
      IF(KK .GT. K2) GO TO 230
      AK=C1-(CD*C1+SD*S1)
      S1=(SD*C1-CD*S1)+S1
C
C --- The following three statements compensate for truncation errors.
C --- If rounded arithmetic is used, replace them with
C --- C1=AK
C
      C1=0.5/(AK**2+S1**2)+0.5
      S1=C1*S1
      C1=C1*AK
      KK=KK+JC
      IF(KK .LT. K2) GO TO 230
      K1=K1+INC+INC
      KK=(K1-KSPAN)/2+JC
      IF(KK .LE. JC+JC) GO TO 220
      GO TO 100
C
C --- Radix-3 transform.
C
320   K1=KK+KSPAN
      K2=K1+KSPAN
      AK=A(KK)
      BK=B(KK)
      AJ=A(K1)+A(K2)
      BJ=B(K1)+B(K2)
      A(KK)=AK+AJ
      B(KK)=BK+BJ
      AK=-0.5*AJ+AK
      BK=-0.5*BJ+BK
      AJ=(A(K1)-A(K2))*S120
      BJ=(B(K1)-B(K2))*S120
      A(K1)=AK-BJ
      B(K1)=BK+AJ
      A(K2)=AK+BJ
      B(K2)=BK-AJ
      KK=K2+KSPAN
      IF(KK .LT. NN) GO TO 320
      KK=KK-NN
      IF(KK .LE. KSPAN) GO TO 320
      GO TO 700
C
C --- Radix-4 transform.
C
400   IF(NFAC(I) .NE. 4) GO TO 600
      KSPNN=KSPAN
      KSPAN=KSPAN/4
410   C1=1.0
      S1=0
420   K1=KK+KSPAN
      K2=K1+KSPAN
      K3=K2+KSPAN
      AKP=A(KK)+A(K2)
      AKM=A(KK)-A(K2)
      AJP=A(K1)+A(K3)
      AJM=A(K1)-A(K3)
      A(KK)=AKP+AJP
      AJP=AKP-AJP
      BKP=B(KK)+B(K2)
      BKM=B(KK)-B(K2)
      BJP=B(K1)+B(K3)
      BJM=B(K1)-B(K3)
      B(KK)=BKP+BJP
      BJP=BKP-BJP
      IF(ISN .LT. 0) GO TO 450
      AKP=AKM-BJM
      AKM=AKM+BJM
      BKP=BKM+AJM
      BKM=BKM-AJM
      IF(S1 .EQ. 0.0) GO TO 460
430   A(K1)=AKP*C1-BKP*S1
      B(K1)=AKP*S1+BKP*C1
      A(K2)=AJP*C2-BJP*S2
      B(K2)=AJP*S2+BJP*C2
      A(K3)=AKM*C3-BKM*S3
      B(K3)=AKM*S3+BKM*C3
      KK=K3+KSPAN
      IF(KK .LE. NT) GO TO 420
440   C2=C1-(CD*C1+SD*S1)
      S1=(SD*C1-CD*S1)+S1
C
C --- The following three statements compensate for truncation errors.
C --- If rounede arithmetic is used, substitute:
C --- C1=C2
C
      C1=0.5/(C2**2+S1**2)+0.5
      S1=C1*S1
      C1=C1*C2
      C2=C1**2-S1**2
      S2=2.0*C1*S1
      C3=C2*C1-S2*S1
      S3=C2*S1+S2*C1
      KK=KK-NT+JC
      IF(KK .LE. KSPAN) GO TO 420
      KK=KK-KSPAN+INC
      IF(KK .LE. JC) GO TO 410
      IF(KSPAN .EQ. JC) GO TO 800
      GO TO 100
450   AKP=AKM+BJM
      AKM=AKM-BJM
      BKP=BKM-AJM
      BKM=BKM+AJM
      IF(S1 .NE. 0.0) GO TO 430
460   A(K1)=AKP
      B(K1)=BKP
      A(K2)=AJP
      B(K2)=BJP
      A(K3)=AKM
      B(K3)=BKM
      KK=K3+KSPAN
      IF(KK .LE. NT) GO TO 420
      GO TO 440
C
C --- Radix-5 transform.
C
510   C2=C72**2-S72**2
      S2=2.0*C72*S72
520   K1=KK+KSPAN
      K2=K1+KSPAN
      K3=K2+KSPAN
      K4=K3+KSPAN
      AKP=A(K1)+A(K4)
      AKM=A(K1)-A(K4)
      BKP=B(K1)+B(K4)
      BKM=B(K1)-B(K4)
      BKP=B(K1)+B(K4)
      BKM=B(K1)-B(K4)
      AJP=A(K2)+A(K3)
      AJM=A(K2)-A(K3)
      BJP=B(K2)+B(K3)
      BJM=B(K2)-B(K3)
      AA=A(KK)
      BB=B(KK)
      A(KK)=AA+AKP+AJP
      B(KK)=BB+BKP+BJP
      AK=AKP*C72+AJP*C2+AA
      BK=BKP*C72+BJP*C2+BB
      AJ=AKM*S72+AJM*S2
      BJ=BKM*S72+BJM*S2
      A(K1)=AK-BJ
      A(K4)=AK+BJ
      B(K1)=BK+AJ
      B(K4)=BK-AJ
      AK=AKP*C2+AJP*C72+AA
      BK=BKP*C2+BJP*C72+BB
      AJ=AKM*S2-AJM*S72
      BJ=BKM*S2-BJM*S72
      A(K2)=AK-BJ
      A(K3)=AK+BJ
      B(K2)=BK+AJ
      B(K3)=BK-AJ
      KK=K4+KSPAN
      IF(KK .LT. NN) GO TO 520
      KK=KK-NN
      IF(KK .LE. KSPAN) GO TO 520
      GO TO 700
C
C --- Transform for odd radices.
C
600   K=NFAC(I)
      KSPNN=KSPAN
      KSPAN=KSPAN/K
      IF(K .EQ. 3) GO TO 320
      IF(K .EQ. 5) GO TO 510
      IF(K .EQ. JF) GO TO 640
      JF=K
      S1=RAD/FLOAT(K)
      C1=COS(S1)
      S1=SIN(S1)
      IF(JF .GT. MAXF) GO TO 998
      CK(JF)=1.0
      SK(JF)=0.0
      J=1
630   CK(J)=CK(K)*C1+SK(K)*S1
      SK(J)=CK(K)*S1-SK(K)*C1
      K=K-1
      CK(K)=CK(J)
      SK(K)=-SK(J)
      J=J+1
      IF(J .LT. K) GO TO 630
640   K1=KK
      K2=KK+KSPNN
      AA=A(KK)
      BB=B(KK)
      AK=AA
      BK=BB
      J=1
      K1=K1+KSPAN
650   K2=K2-KSPAN
      J=J+1
      AT(J)=A(K1)+A(K2)
      AK=AT(J)+AK
      BT(J)=B(K1)+B(K2)
      BK=BT(J)+BK
      J=J+1
      AT(J)=A(K1)-A(K2)
      BT(J)=B(K1)-B(K2)
      K1=K1+KSPAN
      IF(K1 .LT. K2) GO TO 650
      A(KK)=AK
      B(KK)=BK
      K1=KK
      K2=KK+KSPNN
      J=1
660   K1=K1+KSPAN
      K2=K2-KSPAN
      JJ=J
      AK=AA
      BK=BB
      AJ=0.0
      BJ=0.0
      K=1
670   K=K+1
      AK=AT(K)*CK(JJ)+AK
      BK=BT(K)*CK(JJ)+BK
      K=K+1
      AJ=AT(K)*SK(JJ)+AJ
      BJ=BT(K)*SK(JJ)+BJ
      JJ=JJ+J
      IF(JJ .GT. JF) JJ=JJ-JF
      IF(K .LT. JF) GO TO 670
      K=JF-J
      A(K1)=AK-BJ
      B(K1)=BK+AJ
      A(K2)=AK+BJ
      B(K2)=BK-AJ
      J=J+1
      IF(J .LT. K) GO TO 660
      KK=KK+KSPNN
      IF(KK .LE. NN) GO TO 640
      KK=KK-NN
      IF(KK .LE. KSPAN) GO TO 640
C
C --- Multiply by rotation factor (not radices 2 & 4 - they do it
C --- themselves).
C
700   IF(I .EQ. M) GO TO 800
      KK=JC+1
710   C2=1.0-CD
      S1=SD
720   C1=C2
      S2=S1
      KK=KK+KSPAN
730   AK=A(KK)
      A(KK)=C2*AK-S2*B(KK)
      B(KK)=S2*AK+C2*B(KK)
      KK=KK+KSPNN
      IF(KK .LE. NT) GO TO 730
      AK=S1*S2
      S2=S1*C2+C1*S2
      C2=C1*C2-AK
      KK=KK-NT+KSPAN
      IF(KK .LE. KSPNN) GO TO 730
      C2=C1-(CD*C1+SD*S1)
      S1=S1+(SD*C1-CD*S1)
C
C --- The following three statements compensate for truncation errors.
C --- If rounded arithmetic is used, they may be deleted.
C
      C1=0.5/(C2**2+S1**2)+0.5
      S1=C1*S1
      C2=C1*C2
      KK=KK-KSPNN+JC
      IF(KK .LE. KSPAN) GO TO 720
      KK=KK-KSPAN+JC+INC
      IF(KK .LE. JC+JC) GO TO 710
      GO TO 100
C
C --- Permute the results to normal order...this is done in 2 stages.
C --- (i) permutation for square factors of N.
C
800   NP(1)=KS
      IF(KT .EQ. 0) GO TO 890
      K=KT+KT+1
      IF(M .LT. K) K=K-1
      J=1
      NP(K+1)=JC
810   NP(J+1)=NP(J)/NFAC(J)
      NP(K)=NP(K+1)*NFAC(J)
      J=J+1
      K=K-1
      IF(J .LT. K) GO TO 810
      K3=NP(K+1)
      KSPAN=NP(2)
      KK=JC+1
      J=1
      K2=KSPAN+1
      IF(N .NE. NTOT) GO TO 850
C
C --- Permutation for single-variate transform.
C
820   AK=A(KK)
      A(KK)=A(K2)
      A(K2)=AK
      BK=B(KK)
      B(KK)=B(K2)
      B(K2)=BK
      KK=KK+INC
      K2=KSPAN+K2
      IF(K2 .LT. KS) GO TO 820
830   K2=K2-NP(J)
      J=J+1
      K2=NP(J+1)+K2
      IF(K2 .GT. NP(J)) GO TO 830
      J=1
840   IF(KK .LT. K2) GO TO 820
      KK=KK+INC
      K2=KSPAN+K2
      IF(K2 .LT. KS) GO TO 840
      IF(KK .LT. KS) GO TO 830
      JC=K3
      GO TO 890
C
C --- Permutation for multi-variate transform.
C
850   K=KK+JC
860   AK=A(KK)
      A(KK)=A(K2)
      A(K2)=AK
      BK=B(KK)
      B(KK)=B(K2)
      B(K2)=BK
      KK=KK+INC
      K2=K2+INC
      IF(KK .LT. K) GO TO 860
      KK=KK+KS-JC
      K2=K2+KS-JC
      IF(KK .LT. NT) GO TO 850
      K2=K2-NT+KSPAN
      KK=KK-NT+JC
      IF(K2 .LT. KS) GO TO 850
870   K2=K2-NP(J)
      J=J+1
      K2=NP(J+1)+K2
      IF(K2 .GT. NP(J)) GO TO 870
      J=1
880   IF(KK .LT. K2) GO TO 850
      KK=KK+JC
      K2=KSPAN+K2
      IF(K2 .LT. KS) GO TO 880
      IF(KK .LT. KS) GO TO 870
      JC=K3
890   IF(2*KT+1 .GE. M) RETURN
      KSPNN=NP(KT+1)
C
C --- (ii) permutation for square-free factors of N.
C
      J=M-KT
      NFAC(J+1)=1
900   NFAC(J)=NFAC(J)*NFAC(J+1)
      J=J-1
      IF(J .NE. KT) GO TO 900
      KT=KT+1
      NN=NFAC(KT)-1
      IF(NN .GT. MAXP) GO TO 998
      JJ=0
      J=0
      GO TO 906
902   JJ=JJ-K2
      K2=KK
      K=K+1
      KK=NFAC(K)
904   JJ=KK+JJ
      IF(JJ .GE. K2) GO TO 902
      NP(J)=JJ
906   K2=NFAC(KT)
      K=KT+1
      KK=NFAC(K)
      J=J+1
      IF(J .LE. NN) GO TO 904
C
C --- Determine the permutation cycles of length greater than unity.
C
      J=0
      GO TO 914
910   K=KK
      KK=NP(K)
      NP(K)=-KK
      IF(KK .NE. J) GO TO 910
      K3=KK
914   J=J+1
      KK=NP(J)
      IF(KK .LT. 0) GO TO 914
      IF(KK .NE. J) GO TO 910
      NP(J)=-J
      IF(J .NE. NN) GO TO 914
      MAXF=INC*MAXF
C
C --- Re-order A & B, following the permutation cycles.
C
      GO TO 950
924   J=J-1
      IF(NP(J) .LT. 0) GO TO 924
      JJ=JC
926   KSPAN=JJ
      IF(JJ .GT. MAXF) KSPAN=MAXF
      JJ=JJ-KSPAN
      K=NP(J)
      KK=JC*K+II+JJ
      K1=KK+KSPAN
      K2=0
928   K2=K2+1
      AT(K2)=A(K1)
      BT(K2)=B(K1)
      K1=K1-INC
      IF(K1 .NE. KK) GO TO 928
932   K1=KK+KSPAN
      K2=K1-JC*(K+NP(K))
      K=-NP(K)
936   A(K1)=A(K2)
      B(K1)=B(K2)
      K1=K1-INC
      K2=K2-INC
      IF(K1 .NE. KK) GO TO 936
      KK=K2
      IF(K .NE. J) GO TO 932
      K1=KK+KSPAN
      K2=0
940   K2=K2+1
      A(K1)=AT(K2)
      B(K1)=BT(K2)
      K1=K1-INC
      IF(K1 .NE. KK) GO TO 940
      IF(JJ .NE. 0) GO TO 926
      IF(J .NE. 1) GO TO 924
950   J=K3+1
      NT=NT-KSPNN
      II=NT-INC+1
      IF(NT .GE. 0) GO TO 924
      RETURN
C
C --- Error abort: insufficient array storage.
C
998   ISN=0
      WRITE(6,999)
      STOP
999   FORMAT(36H FFT:FATAL - ARRAY BOUNDS EXCEEDED.   )
      END

ext: f


text:
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE CSA3LD(NI,XI,UI,KNOTS,NO,XO,YO,ZO,UO,NWRK,WORK,IER)
      DOUBLE PRECISION XI
      DOUBLE PRECISION UI
      DOUBLE PRECISION XO
      DOUBLE PRECISION YO
      DOUBLE PRECISION ZO
      DOUBLE PRECISION UO
      DOUBLE PRECISION WORK
      DOUBLE PRECISION WTS
      DOUBLE PRECISION SSMTH
C
C  Alternate entry where the output points are specified as a list.
C
      DIMENSION XI(3,NI),UI(NI),KNOTS(3),XO(NO),YO(NO),ZO(NO),UO(NO),
     +          WORK(NWRK),NDERIV(3)
      DATA NDERIV/0,0,0/
C
C  Check on the number of knots.
C
      NTOT = KNOTS(1)*KNOTS(2)*KNOTS(3)
C
      DO 20 I = 1,3
          IF (KNOTS(I).LT.4) THEN
              CALL CFAERR(202,
     +' CSA3LD - must have at least four knots in every coordinate direc
     +tion',69)
              IER = 202
              RETURN
          END IF
   20 CONTINUE
C
C  Check on the size of the workspace.
C
      IF (NWRK.LT.NTOT* (NTOT+3)) THEN
          CALL CFAERR(203,' CSA3LD - workspace too small',29)
          IER = 203
          RETURN
      END IF
C
C  Invoke the expanded function.
C
      WTS = -1.D0
      SSMTH = 0.D0
      CALL CSA3LXD(NI,XI,UI,WTS,KNOTS,SSMTH,NDERIV,NO,XO,YO,ZO,UO,NWRK,
     +             WORK,IER)
      IF (IERR.NE.0) RETURN
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      INCLUDE 'VICMAIN_FOR'
C**********************************************************************
      SUBROUTINE MAIN44
C
C     DEMLOGA:
C
C    THIS PROGRAM WILL LOG USGS DEM FILES THAT ARE IN THE FORMAT DESCRIBED
C	IN THE USGS NATIONAL MAPPING DIVISION OPEN-FILE REPORT 86-004
C	"STANDARDS FOR DIGITAL ELEVATION MODELS"
C
C    THIS FORMAT DIFFERS FROM THE ONE KNOWN TO VDEMLOG; IT INVOLVES
C	UNLABELLED TAPES WITH FIXED LENGTH BLOCKSIZE OF 4096 BYTES.
C	ALL DATA ARE STORED AS ASCII CHARACTERS IN 1024 BYTE LOGICAL
C	RECORDS.
C
C	WRITTEN BY RON ALLEY, NOV 1986
C	MODIFIED TO INCLUDE BLOCKING FACTOR PARAMETER, AUG 1987
C	converted to unix, Ron Alley, November 1994
C	modified to handle <CR> delimited records,  Ron Alley,  April, 1997
C
	REAL*8 CORNER(2,4),ZMIN,ZMAX,ANGLE,YLOC,OFFSET
	INTEGER*2 LEVEL,IPATTERN,ISYSTEM,IZONE,IXY,IZ,NSIDES
	INTEGER*2 NROWS,NCOLS,NUMREC,NPTS,IBUF2(10000),IOUT(10000)
	INTEGER*2 IZERO/0/
	CHARACTER*32768 XBUF
	CHARACTER*60 NAME
	CHARACTER*132 PRT
	CHARACTER*11 UNITS(5)/'Radians','Feet','Meters','Arc-Seconds',
     +			      ' '/
C							  get blocking factor
	CALL XVPARM('BLOCK',IBLKFAC,ICNT,IDEF,0)
	IF (IBLKFAC .EQ. 0) THEN
	    NSIN = 1025
	ELSE
	    NSIN = 1024 * IBLKFAC
	ENDIF
C
C							     open input file
	CALL XVUNIT(INUNIT,'INP',1,ISTAT,' ')
	CALL XVOPEN(INUNIT,ISTAT,'OPEN_ACT','SA','IO_ACT','SA',
     +		    'U_NS',NSIN,'COND','NOLABELS',' ')
C
	CALL XVREAD(INUNIT,XBUF,ISTAT,' ')
	READ (XBUF,100) NAME,LEVEL,IPATTERN,ISYSTEM,IZONE,IXY,IZ,
     +		NSIDES,CORNER,ZMIN,ZMAX,ANGLE,DX,DY,DZ,NROWS,NCOLS
  100	FORMAT(A60,84X,4I6,360X,3I6,4(2D24.15),3D24.15,6X,3E12.6,2I6)
	IF ((IXY.GT.4) .OR. (IXY.LT.0)) IXY=4
	IF ((IZ.GT.2) .OR. (IZ.LT.1)) IXY=4
C							report to screen/log
	CALL XVMESSAGE(NAME,' ')
	WRITE (PRT,102) ZMIN,ZMAX,UNITS(IZ+1)
  102	FORMAT('Range =',F7.1,' to',F7.1,' ',A11)
	CALL XVMESSAGE(PRT,' ')
C
C					Check for the oddball cases that
C					this version won't handle.
C
	IF (IPATTERN.NE.1) THEN
	    CALL XVMESSAGE(
     +		' Random pattern logging not implemented.',' ')
	    CALL ABEND
	ENDIF
	IF (NSIDES.NE.4) THEN
	    CALL XVMESSAGE(
     +	 ' Logging of non-quadrilateral regions is not implemented.',' ')
	    CALL ABEND
	ENDIF
	IF (NROWS.NE.1) THEN
	    CALL XVMESSAGE(
     +	     ' Logging of multiple row profiles is not implemented.',' ')
	    CALL ABEND
	ENDIF
C
	XBASE = MIN(CORNER(1,1),CORNER(1,2),CORNER(1,3),CORNER(1,4))
	XTOP = MAX(CORNER(1,1),CORNER(1,2),CORNER(1,3),CORNER(1,4))
	YBASE = MIN(CORNER(2,1),CORNER(2,2),CORNER(2,3),CORNER(2,4))
	YTOP = MAX(CORNER(2,1),CORNER(2,2),CORNER(2,3),CORNER(2,4))
	NS = (YTOP-YBASE)/DY + 1
	NL = NCOLS
C								Open output file
	CALL XVUNIT(IOUTUNIT,'OUT',1,ISTAT,' ')
	CALL XVOPEN(IOUTUNIT,ISTAT,'OPEN_ACT','SA','IO_ACT','SA',
     +		    'OP','WRITE','U_NL',NL,'U_NS',NS,'U_FORMAT','HALF',
     +		    'O_FORMAT','HALF',' ')
C
C						Add location name to VICAR label
	CALL XLADD(IOUTUNIT,'HISTORY','LOCATION',NAME,ISTAT,
     +		   'FORMAT','STRING',' ')
C
C					    Add geographic coordinates to label
	IF (ISYSTEM.LE.1) THEN
	    IF (ISYSTEM.EQ.1) THEN
		WRITE (PRT,110) IZONE
  110		FORMAT(I4)
		CALL XLADD(IOUTUNIT,'HISTORY','UTM_ZONE',PRT,ISTAT,
     +			   'FORMAT','STRING',' ')
	    END IF
	    WRITE (PRT,120) YBASE,YTOP
  120	    FORMAT(2F14.1)
	    CALL XLADD(IOUTUNIT,'HISTORY','NORTHING',PRT,ISTAT,
     +		   'FORMAT','STRING',' ')
	    WRITE (PRT,120) XBASE,XTOP
	    CALL XLADD(IOUTUNIT,'HISTORY','EASTING',PRT,ISTAT,
     +		   'FORMAT','STRING',' ')
	ELSE IF (ISYSTEM.EQ.2) THEN
	    WRITE (PRT,130) IZONE
  130	    FORMAT(I6)
	    CALL XLADD(IOUTUNIT,'HISTORY','ST_CODE',PRT,ISTAT,
     +		   'FORMAT','STRING',' ')
	ELSE
	    WRITE (PRT,130) ISYSTEM
	    CALL XLADD(IOUTUNIT,'HISTORY','PROJECT',PRT,ISTAT,
     +		   'FORMAT','STRING',' ')
	END IF
C						     Put pixel spacing in label
	IF (DX .EQ. DY) THEN
	    WRITE (PRT,150) DX,UNITS(IXY+1)
  150	    FORMAT(F8.2,' ',A11)
	ELSE
	    WRITE (PRT,160) DX,DY,UNITS(IXY+1)
  160	    FORMAT(F8.2,' by',F8.2,' ',A11)
	END IF
	CALL XLADD(IOUTUNIT,'HISTORY','1_PIXEL',PRT,ISTAT,
     +			'FORMAT','STRING',' ')
C							Put elevation resolution
C							in label
	WRITE (PRT,170) DZ,UNITS(IZ+1)
  170	FORMAT(F8.2,' ',A11)
	CALL XLADD(IOUTUNIT,'HISTORY','1_DN',PRT,ISTAT,
     +			'FORMAT','STRING',' ')
C							Put elevation range in
C							label
	WRITE (PRT,180) ZMIN,ZMAX
  180	FORMAT(F8.2,' -',F8.2)
	CALL XLADD(IOUTUNIT,'HISTORY','RANGE',PRT,ISTAT,
     +			'FORMAT','STRING',' ')
C							Put quality level code
C							in label
	WRITE (PRT,190) LEVEL
  190   FORMAT(I2)
	CALL XLADD(IOUTUNIT,'HISTORY','LEVEL',PRT,ISTAT,
     +			'FORMAT','STRING',' ')
C							Put rotation angle in
C							label, if not 0
	IF (ANGLE.NE.0) THEN
	    WRITE (PRT,195) ANGLE
  195	    FORMAT(F7.2)
	    CALL XLADD(IOUTUNIT,'HISTORY','ROTATION',PRT,ISTAT,
     +			'FORMAT','STRING',' ')
	END IF
C							Build image
	N = 2
	DO LINE=1,NL
	    IF (N.GT.IBLKFAC) THEN
		CALL XVREAD(INUNIT,XBUF,ISTAT,' ')
		N = 1
	    END IF
	    READ(XBUF(1024*(N-1)+1:1024*N),200) NUMREC,NPTS,YLOC,OFFSET,
     +					(IBUF2(I), I=1,146)
  200	    FORMAT(6X,2I6,30X,2D24.15,48X,146I6)
	    N = N+1
	    NBUF = 146
	    IF (NUMREC.NE.LINE) THEN
		CALL XVMESSAGE(' Record sequence error.',' ')
		CALL ABEND
	    ENDIF
	    DO WHILE (NPTS.GT.NBUF)
		IF (N.GT.IBLKFAC) THEN
		    CALL XVREAD(INUNIT,XBUF,ISTAT,' ')
		    N = 1
		END IF
		READ (XBUF(1024*(N-1)+1:1024*N),210) (IBUF2(I),
     +						      I=NBUF+1,NBUF+170)
  210		FORMAT (170I6)
		N = N+1
		NBUF = NBUF+170
	    END DO
	    IF (OFFSET.NE.0.0) THEN
		DO I=1,NPTS
		    IBUF2(I) = IBUF2(I)+OFFSET+0.5
		END DO
	    END IF
	    NPAD = (YLOC-YBASE)/DY
	    IF (NPAD .NE. 0) THEN
		DO LOC=1,NPAD
		    IOUT(LOC) = IZERO
		END DO
	    END IF
	    DO LOC=1,NPTS
		LOC2 = LOC + NPAD
		IOUT(LOC2) = IBUF2(LOC)
	    END DO
	    NUM = NPAD+NPTS
	    IF (LOC2 .LT. NS) THEN
		DO LOC=LOC2+1,NS
		    IOUT(LOC) = IZERO
		END DO
	    END IF
	    CALL XVWRIT(IOUTUNIT,IOUT,ISTAT,' ')
	END DO
	RETURN
	END

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X   libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X   Portions of this code were written by
! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X   Copyright 2006-2010.
! H0 X
! H0 X   These portions of the source code are released under the GNU General
! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X   If you would like to license the source code under different terms,
! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X   Portions of this code were written by Noam Bernstein as part of
! H0 X   his employment for the U.S. Government, and are not subject
! H0 X   to copyright in the USA.
! H0 X
! H0 X
! H0 X   When using this software, please cite the following reference:
! H0 X
! H0 X   http://www.libatoms.org
! H0 X
! H0 X  Additional contributions by
! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

program test_CASTEP_MM_buffer_crack
use libatoms_module
use potential_module
implicit none
  type(Atoms) :: at, cluster_qm, cluster_mm
  type(Inoutput) :: io
  character(len=STRING_LENGTH) :: infile, init_args, calc_args
  real(dp) :: rqm, rmm, g_width
  type(Potential) :: pot
  type(Dictionary) :: cli
  type(Table) :: cluster_info
  integer, pointer :: idx(:)
  logical, pointer :: is_in_qm(:)
  real(dp), allocatable :: f(:,:)
  integer :: i
  real(dp) :: charge(128)

  call system_initialise()

  call initialise(cli)
  infile = ""
  call param_register(cli, "infile", "stdin", infile, help_string="No help yet.  This source file was $LastChangedBy$")
  init_args = ""
  call param_register(cli, "init_args", param_mandatory, init_args, help_string="No help yet.  This source file was $LastChangedBy$")
  calc_args = ""
  call param_register(cli, "calc_args", "", calc_args, help_string="No help yet.  This source file was $LastChangedBy$")
  call param_register(cli, "rqm", param_mandatory, rqm, help_string="No help yet.  This source file was $LastChangedBy$")
  call param_register(cli, "rmm", "0.0", rmm, help_string="No help yet.  This source file was $LastChangedBy$")
  call param_register(cli, "g_width", "0.25", rmm, help_string="No help yet.  This source file was $LastChangedBy$")
  if (.not. param_read_args(cli, ignore_unknown=.false.)) &
    call system_abort("Failed to parse CLI parameters")

  charge = 0.0_dp
  charge(6) = -1.25_dp
  charge(14) = 1.25_dp

  call finalise(cli)

  call read(at, infile)

  call initialise(pot, init_args)

  call create_hybrid_weights(at, "hysteretic_buffer hysteretic_buffer_inner_radius="//(rqm+rmm)//" hysteretic_buffer_outer_radius="//(rqm+rmm))
  cluster_info = create_cluster_info_from_hybrid_mark(at, "cluster_periodic_z")
  call carve_cluster(at, "randomise_buffer=F cluster_periodic_z" , cluster_info, cluster_mm)

  call create_hybrid_weights(cluster_mm, "hysteretic_buffer hysteretic_buffer_inner_radius="//rqm//" hysteretic_buffer_outer_radius="//rqm)
  cluster_info = create_cluster_info_from_hybrid_mark(cluster_mm, "cluster_periodic_z")
  call carve_cluster(cluster_mm, "randomise_buffer=F cluster_periodic_z" , cluster_info, cluster_qm)
  call finalise(cluster_info)

  call add_property(cluster_mm, "is_in_qm", .false.)
  if (.not. assign_pointer(cluster_mm, "is_in_qm", is_in_qm)) &
    call system_abort("Failed to assign pointer for is_in_qm in cluster_mm")
  if (.not. assign_pointer(cluster_qm, "index", idx)) &
    call system_abort("Failed to assign pointer for index in cluster_qm")

  call initialise(io, "extcharges", OUTPUT)
  do i=1, cluster_qm%N
    is_in_qm(idx(i)) = .true.
  end do
  do i=1, cluster_mm%N
    if (.not. is_in_qm(i)) then
      call print(charge(cluster_mm%Z(i))//" "//cluster_mm%pos(:,i)//" "//g_width, file=io)
    endif
  end do
  call finalise(io)

  allocate(f(3,cluster_qm%N))
  call calc(pot, cluster_qm, f=f)

  do i=1, cluster_qm%N
    call print(i// " " // cluster_qm%pos(:,i) // " " // f(:,i))
  end do

  call system_finalise()

end program test_CASTEP_MM_buffer_crack

max_stars_repo_name: albapa/QUIP

max_issues_repo_name: albapa/QUIP

max_forks_repo_name: albapa/QUIP


module netcdf_file_mod

use  ncd_define_mod
use   ncfile_mod, only: ncfile_type, write_ncfile, read_ncfile, &
                        copy_ncfile
use   ncaxis_mod, only: put_ncaxis, get_ncaxis_no, remove_ncaxis_data
use    ncvar_mod, only: put_ncvar, get_ncvar_no
use    ncatt_mod, only: put_ncatt, get_ncatt
use ncvarmap_mod, only: put_varmap, get_varmap, get_varmap_no

implicit none
private

public  define_file, define_axis, define_time_axis,         &
        define_variable, define_variable_tave, define_tave, &
        read_file_header, write_file_header, copy_file_header
public  get_axis_units, get_axis_long_name, get_axis_cart_name, &
        get_axis_edges_name, get_axis_direction, get_calendar_type,  &
        get_var_units, get_var_long_name,       &
        get_missing_value, get_valid_range

!---- overloaded get routines ----

interface get_axis_units
  module procedure get_axis_units_id, get_axis_units_ch
end interface

interface get_axis_long_name
  module procedure get_axis_long_name_id, get_axis_long_name_ch
end interface

interface get_axis_edges_name
  module procedure get_axis_edges_name_id, get_axis_edges_name_ch
end interface

interface get_axis_cart_name
  module procedure get_axis_cart_name_id, get_axis_cart_name_ch
end interface

interface get_axis_direction
  module procedure get_axis_direction_id, get_axis_direction_ch
end interface

interface get_calendar_type
  module procedure get_calendar_type_id, get_calendar_type_ch,  &
                   get_calendar_type_nn
end interface

interface get_var_units
  module procedure get_var_units_id, get_var_units_ch
end interface

interface get_var_long_name
  module procedure get_var_long_name_id, get_var_long_name_ch
end interface

interface get_missing_value
  module procedure get_missing_value_id, get_missing_value_ch
end interface

interface get_valid_range
  module procedure get_valid_range_id, get_valid_range_ch
end interface

!---- time averaging character strings -----

 character(len=32) :: t1_long_name = 'Start time for average period'
 character(len=32) :: t2_long_name = 'End time for average period'
 character(len=32) :: t3_long_name = 'Length of average period'
!character(len=32) :: t3_long_name = 'Number of items in avg period'

integer, parameter :: FATAL = 2

CONTAINS

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

     subroutine define_axis (File, axis_name, values, units,  &
                             cart_axis, long_name, direction)

     type(ncfile_type), intent(inout) :: File
     character(len=*) , intent(in)    :: axis_name
     real             , intent(in), dimension(:) :: values
     character(len=*) , intent(in), optional     :: units, cart_axis, &
                                                    long_name, direction

     real(float_type), dimension(size(values)) :: coords
     integer  id, iatt

     coords = values
     call put_ncaxis (File%Axis, File%ndim, axis_name, coords, id)

!    -------- generate axis attributes --------

     if (present(long_name)) then
            call put_ncatt (File%Axis(id)%Att, File%Axis(id)%natt,  &
                           'long_name', long_name, iatt)
     endif

     if (present(cart_axis)) then
            call put_ncatt (File%Axis(id)%Att, File%Axis(id)%natt,  &
                           'cartesian_axis', cart_axis, iatt)
     endif

     if (present(units)) then
            call put_ncatt (File%Axis(id)%Att, File%Axis(id)%natt,  &
                           'units', units, iatt)
     endif

     if (present(direction)) then
            call put_ncatt (File%Axis(id)%Att, File%Axis(id)%natt,  &
                           'positive', direction, iatt)
     endif

     end subroutine define_axis

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

     subroutine define_time_axis (File, axis_name, units, long_name,  &
                                  calendar, date)

     type(ncfile_type), intent(inout) :: File
     character(len=*) , intent(in)    :: axis_name
     character(len=*) , intent(in), optional :: units, long_name, &
                                                calendar
     integer          , intent(in), optional :: date(:)

     real(double_type), dimension(0) :: values
     character(len=maxc)             :: lunits
     integer  i, nc, id, iatt

     call put_ncaxis (File%Axis, File%ndim, axis_name, values, id)

     File%recdim = id

!    -------- generate axis attributes --------

     if (present(long_name)) then
            call put_ncatt (File%Axis(id)%Att, File%Axis(id)%natt,  &
                           'long_name', long_name, iatt)
     endif

            call put_ncatt (File%Axis(id)%Att, File%Axis(id)%natt,  &
                           'cartesian_axis', 'T', iatt)

     if (present(units)) then
            do i=1,maxc; lunits(i:i)=' '; enddo
            nc = len_trim(units)
            if (present(date)) then
                lunits = units // ' since ' // encode_date (date)
                nc = nc + 28
            else
                lunits = units
            endif
            call put_ncatt (File%Axis(id)%Att, File%Axis(id)%natt,  &
                           'units', lunits(1:nc), iatt)
     endif

     if (present(calendar)) then
            call put_ncatt (File%Axis(id)%Att, File%Axis(id)%natt,  &
                           'calendar_type', calendar, iatt)
     endif

     end subroutine define_time_axis

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

     function encode_date (date) result (label)

     integer, intent(in) :: date(:)
     character(len=21)   :: label
     integer                        :: i, nd, idate(6)
     character(len=2), dimension(6) :: adate

       nd = size(date)
       if (nd < 3) then
           print *, 'ERROR in encode_date: date too small'
           stop 111
       endif
       nd = min (nd, 6)
       idate = 0; idate(1:nd) = date(1:nd)

       do i=2,6
         if (idate(i) < 10) then
            write (adate(i), 10) idate(i)
         else
            write (adate(i), 20) idate(i)
         endif
       enddo

       write (label, 30) idate(1),(adate(i),i=2,6)

 10    format ('0', i1)
 20    format (i2)
 30    format (i4,'-',a2,'-',a2,'-',a2,':',a2,':',a2,'.0')

     end function encode_date

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

   subroutine define_variable (File, var_name, axes, units, long_name, &
                               miss_data, valid_range, pack )

     type(ncfile_type),  intent(inout)        :: File
     character(len=*), intent(in)           :: var_name, axes(:)
     character(len=*), intent(in), optional :: units, long_name
     real            , intent(in), optional :: valid_range(2), miss_data
     logical         , intent(in), optional :: pack

     real(float_type)    :: vrange(2), mdata, scalef, aoffset
     real                :: sdat, srng(2), scale, offset
     integer(short_type) :: sdata, srange(2)
     integer, dimension(max_dims) :: vdims
     integer  i, ndim, id, iatt, tid(3)
     logical   do_short

     ndim = size(axes)

        vdims = 0
     do i = 1, ndim
        vdims(i) = get_ncaxis_no (File%Axis(1:File%ndim), axes(i))
     enddo

!   ------ check for packed 16-bit integers -------
     do_short = .false.
     if (present(pack)) then
         if (present(valid_range)) then
             scale = abs(valid_range(2)-valid_range(1))/65500.
             if (scale <= 1.0) then
                 do_short = pack
                 offset   = 0.5*(valid_range(1)+valid_range(2))
             endif
         else
             scale = 1.0; offset = 0.0
             do_short = pack
         endif
     endif

!   ------- add new variable --------
    if (do_short) then
      call put_ncvar (File%Var, File%nvar, var_name, NCSHORT,  &
                      ndim, vdims, id)
    else
      call put_ncvar (File%Var, File%nvar, var_name, NCFLOAT,  &
                      ndim, vdims, id)
    endif

!    ---------- generate variable attributes ---------

     if (present(long_name)) then
            call put_ncatt (File%Var(id)%Att, File%Var(id)%natt,  &
                           'long_name', long_name, iatt)
     endif

     if (present(units)) then
            call put_ncatt (File%Var(id)%Att, File%Var(id)%natt,  &
                           'units', units, iatt)
     endif

     if (present(valid_range)) then
            if (do_short) then
               if (scale == 1.0 .and. offset == 0.0) then
                   srange = int(valid_range+sign(0.5,valid_range))
               else
                   srng   = (valid_range-offset)/scale
                   srange = int(srng+sign(0.5,srng))
               endif
               call put_ncatt (File%Var(id)%Att, File%Var(id)%natt,  &
                              'valid_range', srange, iatt)
            else
               vrange = valid_range
               call put_ncatt (File%Var(id)%Att, File%Var(id)%natt,  &
                              'valid_range', vrange, iatt)
            endif
     endif

     if (present(miss_data)) then
            if (do_short) then
               if (scale == 1.0 .and. offset == 0.0) then
                   sdata = int(miss_data+sign(0.5,miss_data))
               else
                   sdat  = (miss_data-offset)/scale
                   sdata = int(sdat+sign(0.5,sdat))
               endif
               call put_ncatt (File%Var(id)%Att, File%Var(id)%natt,  &
                              'missing_value', sdata, iatt)
            else
               mdata = miss_data
               call put_ncatt (File%Var(id)%Att, File%Var(id)%natt,  &
                              'missing_value', mdata, iatt)
            endif
     endif

     if (do_short .and. present(valid_range)) then
            scalef = scale; aoffset = offset
            call put_ncatt (File%Var(id)%Att, File%Var(id)%natt,  &
                           'scale_factor', scalef, iatt)
            call put_ncatt (File%Var(id)%Att, File%Var(id)%natt,  &
                           'add_offset',  aoffset, iatt)
     endif

!    ----------- save new variable mapping information ------

     tid = 0
     if (do_short) then
         call put_varmap (File%Varmap, File%nvarmap,   &
                          var_name, id, tid, do_short, scale, offset)
     else
         call put_varmap (File%Varmap, File%nvarmap, var_name, id, tid)
     endif

     end subroutine define_variable

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

     subroutine define_variable_tave (File, var_name, avg_name )

     type(ncfile_type), intent(inout)        :: File
     character(len=*) , intent(in)           :: var_name
     character(len=*) , intent(in), optional :: avg_name
     integer, dimension(max_dims) :: vdims
     integer  i, n, nc, idim, ndim, id, iatt, iatt_units, tid(3)
     character(len=maxn)      :: aname, avg_name_local
     character(len=maxc)      :: units
     character(len=3*maxn+15) :: time_string
     real     scale, offset
     logical  pack
     integer, dimension(max_dims) :: axes

!    ---------- get variable number ---------

      call get_varmap (File%Varmap, File%nvarmap, var_name,  &
                       id, pack, scale, offset)

      if (id == 0) then
          print *, ' ERROR in define_variable_tave: &
                   & requested variable does not exist.'
          stop 111
      endif

!    -------- make sure this variable has a record dimension ------
!    --------    and then get the record dimension units     ------

      if (File%recdim == 0) then
          print *, ' ERROR in define_variable_tave: &
                   & no record dimension for this file.'
          stop 111
      endif

      idim = 0; iatt_units = 0
      do i = 1, File%Var(id)%ndim
         if (File%Var(id)%vdims(i) == File%recdim) then
            idim = File%recdim
            do n=1,maxc; units(n:n)=' '; enddo
            call get_ncatt (File%Axis(idim)%Att, File%Axis(idim)%natt, &
                            'units', units, iatt_units)
            exit
         endif
      enddo

      if (idim == 0) then
          print *, ' ERROR in define_variable_tave: &
                & requested variable does not have a record dimension.'
          stop 111
      endif

      avg_name_local = var_name
      if (present(avg_name)) avg_name_local = avg_name

!    ----------- time averaging information ----------

      tid = 0
      n = len_trim(avg_name_local)
      nc = 3*n + 11
!del  nc = 3*n + 15
      time_string(1:nc) = avg_name_local(1:n) // '_T1,' //  &
                          avg_name_local(1:n) // '_T2,' //  &
                          avg_name_local(1:n) // '_DT'
!del                      avg_name_local(1:n) // '_NITEMS'
      nc = min(nc, maxc)

!     ----- add/replace time average attribute ----

      call put_ncatt (File%Var(id)%Att, File%Var(id)%natt,  &
                     'time_avg_info', time_string(1:nc), iatt)

      call define_tave (File, avg_name_local, tid)

!    ---------- add variable mapping information ----------

      call put_varmap (File%Varmap, File%nvarmap,   &
                       var_name, id, tid, pack, scale, offset)

     end subroutine define_variable_tave

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

  subroutine define_tave (File, avg_name, tid)

  type(ncfile_type), intent(inout) :: File
  character(len=*) , intent(in)    :: avg_name
  integer, optional, intent(out)   :: tid(3)

  character(len=maxn)      :: aname
  character(len=maxc)      :: units
  integer :: n, id(3), idim, iatt, iatt_units, vdims(max_dims)

     n = len_trim(avg_name)
     if (n == 0) return

!---- make sure there is a record dimension and then get the units -----

   if (File%recdim == 0) then
       print *, ' ERROR in define_tave: &
                & no record dimension for this file.'
       stop 111
   endif

      idim = File%recdim
      call get_ncatt (File%Axis(idim)%Att, File%Axis(idim)%natt, &
                         'units', units, iatt_units)

     vdims = 0; vdims(1) = File%recdim

     aname(1:n+3) = avg_name(1:n) // '_T1'
     id(1) = get_ncvar_no (File%Var(1:File%nvar), aname(1:n+3))
     if (id(1) == 0) then
            call put_ncvar (File%Var, File%nvar, aname(1:n+3),  &
                            NCDOUBLE, 1, vdims, id(1))
            if (iatt_units > 0) then
               call put_ncatt (File%Var(id(1))%Att,   &
                               File%Var(id(1))%natt,  &
                               'units', trim(units), iatt)
            endif
            call put_ncatt (File%Var(id(1))%Att, File%Var(id(1))%natt, &
                            'long_name', trim(t1_long_name), iatt)
     endif

     aname(1:n+3) = avg_name(1:n) // '_T2'
     id(2) = get_ncvar_no (File%Var(1:File%nvar), aname(1:n+3))
     if (id(2) == 0) then
            call put_ncvar (File%Var, File%nvar, aname(1:n+3),  &
                            NCDOUBLE, 1, vdims, id(2))
            if (iatt_units > 0) then
               call put_ncatt (File%Var(id(2))%Att,   &
                               File%Var(id(2))%natt,  &
                               'units', trim(units), iatt)
            endif
            call put_ncatt (File%Var(id(2))%Att, File%Var(id(2))%natt, &
                            'long_name', trim(t2_long_name), iatt)
     endif

     aname(1:n+3) = avg_name(1:n) // '_DT'
!del aname(1:n+7) = avg_name(1:n) // '_NITEMS'
     id(3) = get_ncvar_no (File%Var(1:File%nvar), aname(1:n+3))
     if (id(3) == 0) then
            call put_ncvar (File%Var, File%nvar, aname(1:n+3),  &
                            NCDOUBLE, 1, vdims, id(3))
            if (iatt_units > 0) then
               call put_ncatt (File%Var(id(3))%Att,   &
                               File%Var(id(3))%natt,  &
                               'units', trim(units), iatt)
!del                           'units', 'none', iatt)
            endif
            call put_ncatt (File%Var(id(3))%Att, File%Var(id(3))%natt, &
                            'long_name', trim(t3_long_name), iatt)
     endif

     if (present(tid)) tid = id

  end subroutine define_tave

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

     subroutine define_file (File, file_name, title, history, coards)

     type(ncfile_type), intent(inout)        :: File
     character(len=*),  intent(in)           :: file_name
     character(len=*),  intent(in), optional :: title, history
     logical,           intent(in), optional :: coards
     integer   iatt

        allocate (File%Axis  (max_ndim))
        allocate (File%Var   (max_nvar))
        allocate (File%Varmap(max_nvar))
        allocate (File%Att   (max_natt))

        File%id      = 0
        File%name    = file_name
        File%ndim    = 0
        File%nvar    = 0
        File%nvarmap = 0
        File%natt    = 0
        File%recdim  = 0

        File%define_mode = .false.

        if (present(title)) then
                 call put_ncatt (File%Att, File%natt,  &
                                'title', trim(title), iatt)
        endif

        if (present(history)) then
                 call put_ncatt (File%Att, File%natt,  &
                                'history', trim(history), iatt)
        endif

        if (present(coards)) then
            if (coards)  call put_ncatt (File%Att, File%natt,  &
                                        'Conventions', 'COARDS', iatt)
        endif

     end subroutine define_file

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

 function read_file_header (file_name) result (File)

 character(len=*), intent(in)  :: file_name
 type(ncfile_type)             :: File

 integer  i, no, nc, n1, tno, sno, ano, tavg_no(3)
 real(float_type) :: scalef, aoffset
 real             :: scale ,  offset
 character(len=maxc) :: string, tstring
 logical             :: pack

!---- read header information ----

      File = read_ncfile(file_name)

!----------------------------------------
!---- setup varmap_type information -----
!----------------------------------------

     allocate (File%Varmap(max_nvar))
     File%nvarmap = 0

     do i = 1, File%nvar
       nc = len_trim(File%Var(i)%name)
!      --- skip time averaging fields ---
       if (nc >= 3) then
         n1 = nc-2
         if ( File%Var(i)%name(nc-2:nc) == '_T1' ) cycle
         if ( File%Var(i)%name(nc-2:nc) == '_T2' ) cycle
         if ( File%Var(i)%name(nc-2:nc) == '_DT' ) cycle
       endif
       if (nc >= 7) then
         n1 = nc-6
         if ( File%Var(i)%name(nc-2:nc) == '_NITEMS' ) cycle
       endif

!      call get_ncatt (File%Var(i)%Att, File%Var(i)%natt, 'long_name', &
!                     string, no)
!      if (no > 0) then

            call get_ncatt (File%Var(i)%Att, File%Var(i)%natt,&
                            'time_avg_info', tstring, tno)
            tavg_no = 0
            if (tno > 0) call set_tavg_no (File, tstring, tavg_no)

            call get_ncatt (File%Var(i)%Att, File%Var(i)%natt, &
                            'scale_factor', scalef, sno)

            call get_ncatt (File%Var(i)%Att, File%Var(i)%natt, &
                            'add_offset', aoffset, ano)

            pack = .false.
            if (File%Var(i)%datatype == NCSHORT) pack = .true.

            if (sno > 0 .and. ano > 0) then
               scale = scalef; offset = aoffset; pack = .true.
               call put_varmap (File%Varmap, File%nvarmap,  &
                                File%Var(i)%name, i,        &
                                tavg_no, pack, scale, offset)
            else
               call put_varmap (File%Varmap, File%nvarmap,  &
                                File%Var(i)%name, i,  &
                                tavg_no, pack)
            endif
!      endif

     enddo

 end function read_file_header

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

 subroutine write_file_header (File)
 type(ncfile_type), intent(inout) :: File

       call write_ncfile (File)

 end subroutine write_file_header

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

 function copy_file_header (File_in, file_name_out, title, coards)  &
                    result (File_out)
 type(ncfile_type), intent(in)           :: File_in
character(len=*), intent(in), optional :: file_name_out, title
logical,          intent(in), optional :: coards
 type(ncfile_type)                       :: File_out
  integer  iatt, rdim

!!!!!   File_out = File_in
        call copy_ncfile (File_in, File_out)

!     ------ remove record dimension data ------
        call remove_ncaxis_data (File_out%Axis(File_out%recdim))

        if (present(file_name_out)) File_out%name = file_name_out

        if (present(title)) then
                 call put_ncatt (File_out%Att, File_out%natt,  &
                                'title', trim(title), iatt)
        endif

        if (present(coards)) then
            if (coards)  call put_ncatt (File_out%Att, File_out%natt,  &
                                        'Conventions', 'COARDS', iatt)
        endif

 end function copy_file_header

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

 subroutine set_tavg_no (File, string, tavg_no)

 type(ncfile_type) , intent(in)  :: File
 character(len=*), intent(in)  :: string
 integer         , intent(out) :: tavg_no(3)

 integer   is(3), ie(3), nc, i, k

      k = 1; is(k) = 1
      do i = 1, len_trim(string)
          if (string(i:i) /= ',') cycle
          ie(k) = i-1
          k = k+1
          if (k == 4) exit
          is(k) = i+1
      enddo
          ie(3) = len_trim(string)

      do k = 1, 3
          tavg_no(k) =  &
              get_ncvar_no (File%Var(1:File%nvar), string(is(k):ie(k)))
      enddo

 end subroutine set_tavg_no

!-----------------------------------------------------------------------
!#######################################################################
!-----------------------------------------------------------------------
!---------- LOGICAL FUNCTIONS TO RETRIEVE OPTIONAL ARGUMENTS -----------
!-----------------------------------------------------------------------
!---------------------- axis units -------------------------------------

 function get_axis_units_id (File, axis_num, units) result (flag)

 type(ncfile_type) , intent(in)  :: File
 integer         , intent(in)  :: axis_num
 character(len=*), intent(out) :: units
 logical                       :: flag
 integer  i, iatt, nc, ncs
 character(len=maxc) :: string

   if (axis_num < 1 .or. axis_num > File%ndim) then
        call error_handler ('get_axis_units_id',    &
                            'invalid axis index.', FATAL)
   endif

   call get_ncatt (File%Axis(axis_num)%Att, File%Axis(axis_num)%natt, &
                  'units', string, iatt)

   nc = len(units)
   do i=1,nc; units(i:i)=' '; enddo

   if (iatt > 0) then
       flag = .true.
       ncs = min(nc, len_trim(string))
       units(1:ncs) = string(1:ncs)
   else
       flag = .false.
   endif

 end function get_axis_units_id

!---------------------- axis long_name ---------------------------------

 function get_axis_long_name_id (File, axis_num, long_name)  &
                         result (flag)

 type(ncfile_type), intent(in)  :: File
 integer          , intent(in)  :: axis_num
 character(len=*) , intent(out) :: long_name
 logical                        :: flag
 integer  i, iatt, nc, ncs
 character(len=maxc) :: string

   if (axis_num < 1 .or. axis_num > File%ndim) then
        call error_handler ('get_axis_long_name_id',    &
                            'invalid axis index.', FATAL)
   endif

   call get_ncatt (File%Axis(axis_num)%Att, File%Axis(axis_num)%natt, &
                  'long_name', string, iatt)

   nc = len(long_name)
   do i=1,nc; long_name(i:i)=' '; enddo

   if (iatt > 0) then
       flag = .true.
       ncs = min(nc, len_trim(string))
       long_name(1:ncs) = string(1:ncs)
   else
       flag = .false.
   endif

 end function get_axis_long_name_id

!---------------------- axis long_name ---------------------------------

 function get_axis_edges_name_id (File, axis_num, edges_name)  &
                          result (flag)

 type(ncfile_type), intent(in)  :: File
 integer          , intent(in)  :: axis_num
 character(len=*) , intent(out) :: edges_name
 logical                        :: flag
 integer  i, iatt, nc, ncs
 character(len=maxc) :: string

   if (axis_num < 1 .or. axis_num > File%ndim) then
        call error_handler ('get_axis_edges_name_id',    &
                            'invalid axis index.', FATAL)
   endif

   call get_ncatt (File%Axis(axis_num)%Att, File%Axis(axis_num)%natt, &
                  'edges', string, iatt)

   nc = len(edges_name)
   do i=1,nc; edges_name(i:i)=' '; enddo

   if (iatt > 0) then
       flag = .true.
       ncs = min(nc, len_trim(string))
       edges_name(1:ncs) = string(1:ncs)
   else
       flag = .false.
   endif

 end function get_axis_edges_name_id

!---------------------- axis cart_name ---------------------------------

 function get_axis_cart_name_id (File, axis_num, cart_name)  &
                         result (flag)

 type(ncfile_type), intent(in)  :: File
 integer          , intent(in)  :: axis_num
 character(len=*) , intent(out) :: cart_name
 logical                        :: flag
 integer  i, iatt, nc, ncs
 character(len=maxc) :: string

   if (axis_num < 1 .or. axis_num > File%ndim) then
        call error_handler ('get_axis_cart_name_id',    &
                            'invalid axis index.', FATAL)
   endif

   call get_ncatt (File%Axis(axis_num)%Att, File%Axis(axis_num)%natt, &
                   'cartesian_axis', string, iatt)

   nc = len(cart_name)
   do i=1,nc; cart_name(i:i)=' '; enddo

   if (iatt > 0) then
       flag = .true.
       ncs = min(nc, len_trim(string))
       cart_name(1:ncs) = string(1:ncs)
   else
       flag = .false.
   endif

 end function get_axis_cart_name_id

!---------------------- axis (positive) direction ----------------------

 function get_axis_direction_id (File, axis_num, direction)  &
                         result (flag)

 type(ncfile_type), intent(in)  :: File
 integer          , intent(in)  :: axis_num
 character(len=*) , intent(out) :: direction
 logical                        :: flag
 integer  i, iatt, nc, ncs
 character(len=maxc) :: string

   if (axis_num < 1 .or. axis_num > File%ndim) then
        call error_handler ('get_axis_direction_id',    &
                            'invalid axis index.', FATAL)
   endif

   call get_ncatt (File%Axis(axis_num)%Att, File%Axis(axis_num)%natt, &
                  'positive', string, iatt)

   nc = len(direction)
   do i=1,nc; direction(i:i)=' '; enddo

   if (iatt > 0) then
       flag = .true.
       ncs = min(nc, len_trim(string))
       direction(1:ncs) = string(1:ncs)
   else
       flag = .false.
   endif

 end function get_axis_direction_id

!---------------------- time axis calendar type ------------------------

 function get_calendar_type_id (File, axis_num, calendar)  &
                        result (flag)

 type(ncfile_type), intent(in)  :: File
 integer          , intent(in)  :: axis_num
 character(len=*) , intent(out) :: calendar
 logical                        :: flag
 integer  i, iatt, nc, ncs
 character(len=maxc) :: string

   if (axis_num < 1 .or. axis_num > File%ndim) then
        call error_handler ('get_calendar_type_id',     &
                            'invalid axis index.', FATAL)
   endif

   call get_ncatt (File%Axis(axis_num)%Att, File%Axis(axis_num)%natt, &
                  'calendar_type', string, iatt)

   nc = len(calendar)
   do i=1,nc; calendar(i:i)=' '; enddo

   if (iatt > 0) then
       flag = .true.
       ncs = min(nc, len_trim(string))
       calendar(1:ncs) = string(1:ncs)
   else
       flag = .false.
   endif

 end function get_calendar_type_id

!-----------------------------------------------------------------------
!---------------------- variable units ---------------------------------

 function get_var_units_id (File, var_num, units) result (flag)

 type(ncfile_type) , intent(in)  :: File
 integer         , intent(in)  :: var_num
 character(len=*), intent(out) :: units
 logical                       :: flag
 integer  i, iatt, no, nc, ncs
 character(len=maxc) :: string

   if (var_num < 1 .or. var_num > File%nvarmap) then
        call error_handler ('get_var_units_id',     &
                            'invalid variable index.', FATAL)
   endif

   no = File%Varmap(var_num)%varno
   call get_ncatt (File%Var(no)%Att, File%Var(no)%natt, &
                  'units', string, iatt)

   nc = len(units)
   do i=1,nc; units(i:i)=' '; enddo

   if (iatt > 0) then
       flag = .true.
       ncs = min(nc, len_trim(string))
       units(1:ncs) = string(1:ncs)
   else
       flag = .false.
   endif

 end function get_var_units_id

!---------------------- variable long_name -----------------------------

 function get_var_long_name_id (File, var_num, long_name)  &
                        result (flag)

 type(ncfile_type), intent(in)  :: File
 integer          , intent(in)  :: var_num
 character(len=*) , intent(out) :: long_name
 logical                        :: flag
 integer  i, iatt, no, nc, ncs
 character(len=maxc) :: string

   if (var_num < 1 .or. var_num > File%nvarmap) then
        call error_handler ('get_var_long_name_id',  &
                            'invalid variable index.', FATAL)
   endif

   no = File%Varmap(var_num)%varno
   call get_ncatt (File%Var(no)%Att, File%Var(no)%natt, &
                  'long_name', string, iatt)

   nc = len(long_name)
   do i=1,nc; long_name(i:i)=' '; enddo

   if (iatt > 0) then
       flag = .true.
       ncs = min(nc, len_trim(string))
       long_name(1:ncs) = string(1:ncs)
   else
       flag = .false.
   endif

 end function get_var_long_name_id

!---------------- special value feature --------------------------------

 function get_missing_value_id (File, var_num, value) result (flag)

 type(ncfile_type), intent(in)  :: File
 integer          , intent(in)  :: var_num
 real             , intent(out) :: value
 logical                        :: flag
 integer   varno, iatt
 logical   pack
 real      scale, offset
 real(float_type)    :: fval
 integer(short_type) :: sval

    if (var_num < 1 .or. var_num > File%nvarmap) then
        call error_handler ('get_missing_value_id',  &
                            'invalid variable index.', FATAL)
    endif

    call get_varmap (File%Varmap, File%nvarmap,  &
                     File%Varmap(var_num)%name,  &
                     varno, pack, scale, offset)

    if (pack) then
       call get_ncatt (File%Var(varno)%Att, File%Var(varno)%natt,  &
                       'missing_value', sval, iatt)
       if (iatt > 0) then
          if (scale == 1.0 .and. offset == 0.0) then
              value = real(sval)
          else
              value = real(sval)*scale+offset
          endif
          flag = .true.
       else
          flag = .false.
          value = 0.0
       endif
    else
       call get_ncatt (File%Var(varno)%Att, File%Var(varno)%natt,  &
                       'missing_value', fval, iatt)
       if (iatt > 0) then
          value = fval
          flag = .true.
       else
          flag = .false.
          value = 0.0
       endif
    endif

 end function get_missing_value_id

!---------------- valid range ------------------------------------------

 function get_valid_range_id (File, var_num, vrange) result (flag)

 type(ncfile_type), intent(in)  :: File
 integer          , intent(in)  :: var_num
 real             , intent(out) :: vrange(2)
 logical                        :: flag
 integer   varno, iatt
 logical   pack
 real      scale, offset
 real(float_type)    :: fval(2)
 integer(short_type) :: sval(2)

   if (var_num < 1 .or. var_num > File%nvarmap) then
       call error_handler ('get_valid_range_id',  &
                           'invalid variable index', FATAL)
   endif

!     ---- get scale facor and add offset ----

    call get_varmap (File%Varmap, File%nvarmap,  &
                     File%Varmap(var_num)%name,  &
                     varno, pack, scale, offset)

!     ---- get valid range attribute ----
!     pack=.false.    ! for old version

    if (pack) then
       call get_ncatt (File%Var(varno)%Att, File%Var(varno)%natt,  &
                       'valid_range', sval, iatt)
       if (iatt > 0) then
          if (scale == 1.0 .and. offset == 0.0) then
              vrange = real(sval)
          else
              vrange = real(sval)*scale+offset
          endif
          flag = .true.
       else
          flag = .false.
          vrange = 0.0
       endif
    else
       call get_ncatt (File%Var(varno)%Att, File%Var(varno)%natt,  &
                       'valid_range', fval, iatt)
       if (iatt > 0) then
          vrange = fval
          flag = .true.
       else
          flag = .false.
          vrange = 0.0
       endif
    endif

 end function get_valid_range_id

!-----------------------------------------------------------------------
!#######################################################################
!-----------------------------------------------------------------------
!----- OVERLOADED (using char string name) LOGICAL FUNCTIONS -----
!-----------------------------------------------------------------------
!---------------------- axis units -------------------------------------

 function get_axis_units_ch (File, axis_name, units) result (flag)

 type(ncfile_type), intent(in)  :: File
 character(len=*),  intent(in)  :: axis_name
 character(len=*),  intent(out) :: units
 logical :: flag
 integer :: index

   index = get_ncaxis_no (File%Axis(1:File%ndim), axis_name)

   if (index == 0) call error_handler ('get_axis_units_ch',       &
                                       'invalid axis name.', FATAL)

   flag = get_axis_units_id (File, index, units)

 end function get_axis_units_ch

!---------------------- axis long_name ---------------------------------

 function get_axis_long_name_ch (File, axis_name, long_name)  &
                         result (flag)

 type(ncfile_type), intent(in)  :: File
 character(len=*) , intent(in)  :: axis_name
 character(len=*) , intent(out) :: long_name
 logical :: flag
 integer :: index

   index = get_ncaxis_no (File%Axis(1:File%ndim), axis_name)

   if (index == 0) call error_handler ('get_axis_long_name_ch',   &
                                       'invalid axis name.', FATAL)

   flag = get_axis_long_name_id (File, index, long_name)

 end function get_axis_long_name_ch

!---------------------- axis long_name ---------------------------------

 function get_axis_edges_name_ch (File, axis_name, edges_name)  &
                          result (flag)

 type(ncfile_type), intent(in)  :: File
 character(len=*) , intent(in)  :: axis_name
 character(len=*) , intent(out) :: edges_name
 logical :: flag
 integer :: index

   index = get_ncaxis_no (File%Axis(1:File%ndim), axis_name)

   if (index == 0) call error_handler ('get_axis_edges_name_ch',   &
                                       'invalid axis name.', FATAL)

   flag = get_axis_edges_name_id (File, index, edges_name)

 end function get_axis_edges_name_ch

!---------------------- axis cart_name ---------------------------------

 function get_axis_cart_name_ch (File, axis_name, cart_name)  &
                         result (flag)

 type(ncfile_type), intent(in)  :: File
 character(len=*) , intent(in)  :: axis_name
 character(len=*) , intent(out) :: cart_name
 logical :: flag
 integer :: index

   index = get_ncaxis_no (File%Axis(1:File%ndim), axis_name)

   if (index == 0) call error_handler ('get_axis_cart_name_ch',   &
                                       'invalid axis name.', FATAL)

   flag = get_axis_cart_name_id (File, index, cart_name)

 end function get_axis_cart_name_ch

!---------------------- axis (positive) direction ----------------------

 function get_axis_direction_ch (File, axis_name, direction)  &
                         result (flag)

 type(ncfile_type), intent(in)  :: File
 character(len=*) , intent(in)  :: axis_name
 character(len=*) , intent(out) :: direction
 logical :: flag
 integer :: index

   index = get_ncaxis_no (File%Axis(1:File%ndim), axis_name)

   if (index == 0) call error_handler ('get_axis_direction_ch',   &
                                       'invalid axis name.', FATAL)

   flag = get_axis_direction_id (File, index, direction)

 end function get_axis_direction_ch

!---------------------- time axis calendar type ------------------------

 function get_calendar_type_ch (File, axis_name, calendar)  &
                        result (flag)

 type(ncfile_type), intent(in)  :: File
 character(len=*) , intent(in)  :: axis_name
 character(len=*) , intent(out) :: calendar
 logical :: flag
 integer :: index

   index = get_ncaxis_no (File%Axis(1:File%ndim), axis_name)

   if (index == 0) call error_handler ('get_calendar_type_ch',    &
                                       'invalid axis name.', FATAL)

   flag = get_calendar_type_id (File, index, calendar)

 end function get_calendar_type_ch

!---------------------- time axis calendar type ------------------------

 function get_calendar_type_nn (File, calendar) result (flag)

 type(ncfile_type), intent(in)  :: File
 character(len=*) , intent(out) :: calendar
 logical :: flag
 integer :: index

   index = File%recdim

   if (index > 0) then
       flag = get_calendar_type_id (File, index, calendar)
   else
       flag = .false.
   endif

 end function get_calendar_type_nn

!-----------------------------------------------------------------------
!---------------------- variable units ---------------------------------

 function get_var_units_ch (File, var_name, units) result (flag)

 type(ncfile_type), intent(in)  :: File
 character(len=*) , intent(in)  :: var_name
 character(len=*) , intent(out) :: units
 logical :: flag
 integer :: index

   index = get_varmap_no (File%Varmap(1:File%nvarmap), var_name)

   if (index == 0) call error_handler ('get_var_units_ch',      &
                                 'invalid variable name.', FATAL)

   flag = get_var_units_id (File, index, units)

 end function get_var_units_ch

!---------------------- variable long_name -----------------------------

 function get_var_long_name_ch (File, var_name, long_name)  &
                        result (flag)

 type(ncfile_type), intent(in)  :: File
 character(len=*) , intent(in)  :: var_name
 character(len=*) , intent(out) :: long_name
 logical :: flag
 integer :: index

   index = get_varmap_no (File%Varmap(1:File%nvarmap), var_name)

   if (index == 0) call error_handler ('get_var_long_name_ch',  &
                                 'invalid variable name.', FATAL)

   flag = get_var_long_name_id (File, index, long_name)

 end function get_var_long_name_ch

!---------------- special value feature --------------------------------

 function get_missing_value_ch (File, var_name, value) result (flag)

 type(ncfile_type), intent(in)  :: File
 character(len=*) , intent(in)  :: var_name
 real             , intent(out) :: value
 logical :: flag
 integer :: index

   index = get_varmap_no (File%Varmap(1:File%nvarmap), var_name)

   if (index == 0) call error_handler ('get_missing_value_ch',  &
                                 'invalid variable name.', FATAL)

   flag = get_missing_value_id (File, index, value)

 end function get_missing_value_ch

!---------------- valid range ------------------------------------------

 function get_valid_range_ch (File, var_name, vrange) result (flag)

 type(ncfile_type), intent(in)  :: File
 character(len=*) , intent(in)  :: var_name
 real             , intent(out) :: vrange(2)
 logical :: flag
 integer :: index

   index = get_varmap_no (File%Varmap(1:File%nvarmap), var_name)

   if (index == 0) call error_handler ('get_valid_range_ch',    &
                                 'invalid variable name.', FATAL)

   flag = get_valid_range_id (File, index, vrange)

 end function get_valid_range_ch

!-----------------------------------------------------------------------
!#######################################################################

   subroutine error_handler (routine, message, level)
   character(len=*), intent(in) :: routine, message
   integer         , intent(in) :: level

      if ( level == FATAL ) then
          print *, 'ERROR in ',trim(routine),' in netcdf_file_mod'
          print *, trim(message)
          call abort
      else
          print *, 'MESSAGE from ',trim(routine),' in netcdf_file_mod'
          print *, trim(message)
      endif

   end subroutine error_handler

!#######################################################################

end module netcdf_file_mod


text:
*> \brief \b dladivquad performs complex division in real arithmetic, avoiding unnecessary overflow.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download dladivquad + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE dladivquad( A, B, C, D, P, Q )
*
*       .. Scalar Arguments ..
*       real*16   A, B, C, D, P, Q
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> dladivquad performs complex division in  real arithmetic
*>
*>                       a + i*b
*>            p + i*q = ---------
*>                       c + i*d
*>
*> The algorithm is due to Michael Baudin and Robert L. Smith
*> and can be found in the paper
*> "A Robust Complex Division in Scilab"
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] A
*> \verbatim
*>          A is real*16
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is real*16
*> \endverbatim
*>
*> \param[in] C
*> \verbatim
*>          C is real*16
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*>          D is real*16
*>          The scalars a, b, c, and d in the above expression.
*> \endverbatim
*>
*> \param[out] P
*> \verbatim
*>          P is real*16
*> \endverbatim
*>
*> \param[out] Q
*> \verbatim
*>          Q is real*16
*>          The scalars p and q in the above expression.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date January 2013
*
*> \ingroup auxOTHERauxiliary
*
*  =====================================================================
      SUBROUTINE dladivquad( A, B, C, D, P, Q )
*
*  -- LAPACK auxiliary routine (version 3.5.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     January 2013
*
*     .. Scalar Arguments ..
      real*16   A, B, C, D, P, Q
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      real*16   BS
      PARAMETER          ( BS = 2.0q0 )
      real*16   HALF
      PARAMETER          ( HALF = 0.5q0 )
      real*16   TWO
      PARAMETER          ( TWO = 2.0q0 )
*
*     .. Local Scalars ..
      real*16   AA, BB, CC, DD, AB, CD, S, OV, UN, BE, EPS
*     ..
*     .. External Functions ..
      real*16   dlamchquad
      EXTERNAL           dlamchquad
*     ..
*     .. External Subroutines ..
      EXTERNAL           dladivquad1
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX
*     ..
*     .. Executable Statements ..
*
      AA = A
      BB = B
      CC = C
      DD = D
      AB = MAX( ABS(A), ABS(B) )
      CD = MAX( ABS(C), ABS(D) )
      S = 1.0q0

      OV = dlamchquad( 'Overflow threshold' )
      UN = dlamchquad( 'Safe minimum' )
      EPS = dlamchquad( 'Epsilon' )
      BE = BS / (EPS*EPS)

      IF( AB >= HALF*OV ) THEN
         AA = HALF * AA
         BB = HALF * BB
         S  = TWO * S
      END IF
      IF( CD >= HALF*OV ) THEN
         CC = HALF * CC
         DD = HALF * DD
         S  = HALF * S
      END IF
      IF( AB <= UN*BS/EPS ) THEN
         AA = AA * BE
         BB = BB * BE
         S  = S / BE
      END IF
      IF( CD <= UN*BS/EPS ) THEN
         CC = CC * BE
         DD = DD * BE
         S  = S * BE
      END IF
      IF( ABS( D ).LE.ABS( C ) ) THEN
         CALL dladivquad1(AA, BB, CC, DD, P, Q)
      ELSE
         CALL dladivquad1(BB, AA, DD, CC, P, Q)
         Q = -Q
      END IF
      P = P * S
      Q = Q * S
*
      RETURN
*
*     End of dladivquad
*
      END

      SUBROUTINE dladivquad1( A, B, C, D, P, Q )
*
*  -- LAPACK auxiliary routine (version 3.5.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     January 2013
*
*     .. Scalar Arguments ..
      real*16   A, B, C, D, P, Q
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      real*16   ONE
      PARAMETER          ( ONE = 1.0q0 )
*
*     .. Local Scalars ..
      real*16   R, T
*     ..
*     .. External Functions ..
      real*16   dladivquad2
      EXTERNAL           dladivquad2
*     ..
*     .. Executable Statements ..
*
      R = D / C
      T = ONE / (C + D * R)
      P = dladivquad2(A, B, C, D, R, T)
      A = -A
      Q = dladivquad2(B, A, C, D, R, T)
*
      RETURN
*
*     End of dladivquad1
*
      END

      real*16 FUNCTION dladivquad2( A, B, C, D, R, T )
*
*  -- LAPACK auxiliary routine (version 3.5.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     January 2013
*
*     .. Scalar Arguments ..
      real*16   A, B, C, D, R, T
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      real*16   ZERO
      PARAMETER          ( ZERO = 0.0q0 )
*
*     .. Local Scalars ..
      real*16   BR
*     ..
*     .. Executable Statements ..
*
      IF( R.NE.ZERO ) THEN
         BR = B * R
         if( BR.NE.ZERO ) THEN
            dladivquad2 = (A + BR) * T
         ELSE
            dladivquad2 = A * T + (B * T) * R
         END IF
      ELSE
         dladivquad2 = (A + D * (B / C)) * T
      END IF
*
      RETURN
*
*     End of dladivquad12
*
      END


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! The ODE Jacobian of Chemical Model File
!
! Generated by KPP-2.1 symbolic chemistry Kinetics PreProcessor
!       (http://www.cs.vt.edu/~asandu/Software/KPP)
! KPP is distributed under GPL, the general public licence
!       (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
!     With important contributions from:
!        M. Damian, Villanova University, USA
!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
!
! File                 : cbm4_Jacobian.f90
! Time                 : Mon Feb 27 11:08:22 2017
! Working directory    : /data/ksetigui/setigui/WRFV3_171015/chem/KPP/mechanisms/cbm4
! Equation file        : cbm4.kpp
! Output root filename : cbm4
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

MODULE cbm4_Jacobian

  USE cbm4_Parameters
  USE cbm4_JacobianSP

  IMPLICIT NONE

CONTAINS

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! Jac_SP_Vec - function for sparse multiplication: sparse Jacobian times vector
!   Arguments :
!      JVS       - sparse Jacobian of variables
!      UV        - User vector for variables
!      JUV       - Jacobian times user vector
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE Jac_SP_Vec ( JVS, UV, JUV )

! JVS - sparse Jacobian of variables
  REAL(kind=dp) :: JVS(LU_NONZERO)
! UV - User vector for variables
  REAL(kind=dp) :: UV(NVAR)
! JUV - Jacobian times user vector
  REAL(kind=dp) :: JUV(NVAR)

  JUV(1) = JVS(1)*UV(1)+JVS(2)*UV(25)
  JUV(2) = JVS(3)*UV(2)+JVS(4)*UV(27)+JVS(5)*UV(28)
  JUV(3) = JVS(6)*UV(3)+JVS(7)*UV(26)+JVS(8)*UV(32)
  JUV(4) = JVS(9)*UV(4)+JVS(10)*UV(14)+JVS(11)*UV(26)+JVS(12)*UV(27)+JVS(13)*UV(30)
  JUV(5) = JVS(14)*UV(5)+JVS(15)*UV(27)
  JUV(6) = JVS(16)*UV(6)+JVS(17)*UV(26)+JVS(18)*UV(30)
  JUV(7) = JVS(19)*UV(7)+JVS(20)*UV(27)
  JUV(8) = JVS(21)*UV(8)+JVS(22)*UV(13)+JVS(23)*UV(20)+JVS(24)*UV(22)+JVS(25)*UV(23)+JVS(26)*UV(27)+JVS(27)*UV(29)&
             &+JVS(28)*UV(30)+JVS(29)*UV(31)
  JUV(9) = JVS(30)*UV(9)+JVS(31)*UV(26)+JVS(32)*UV(27)+JVS(33)*UV(31)
  JUV(10) = JVS(34)*UV(10)+JVS(35)*UV(26)+JVS(36)*UV(27)+JVS(37)*UV(28)
  JUV(11) = JVS(38)*UV(5)+JVS(39)*UV(7)+JVS(40)*UV(11)+JVS(41)*UV(27)+JVS(42)*UV(31)
  JUV(12) = JVS(43)*UV(6)+JVS(44)*UV(12)+JVS(45)*UV(14)+JVS(46)*UV(21)+JVS(47)*UV(24)+JVS(48)*UV(26)+JVS(49)*UV(27)&
              &+JVS(50)*UV(30)
  JUV(13) = JVS(51)*UV(13)+JVS(52)*UV(20)+JVS(53)*UV(26)+JVS(54)*UV(27)
  JUV(14) = JVS(55)*UV(5)+JVS(56)*UV(7)+JVS(57)*UV(11)+JVS(58)*UV(14)+JVS(59)*UV(27)+JVS(60)*UV(30)
  JUV(15) = JVS(62)*UV(7)+JVS(63)*UV(15)+JVS(64)*UV(19)+JVS(65)*UV(22)+JVS(66)*UV(25)+JVS(67)*UV(27)
  JUV(16) = JVS(68)*UV(15)+JVS(69)*UV(16)+JVS(70)*UV(17)+JVS(71)*UV(19)+JVS(72)*UV(21)+JVS(73)*UV(22)+JVS(74)*UV(23)&
              &+JVS(75)*UV(24)+JVS(76)*UV(25)+JVS(77)*UV(27)+JVS(78)*UV(29)+JVS(79)*UV(30)
  JUV(17) = JVS(80)*UV(17)+JVS(81)*UV(22)+JVS(82)*UV(25)+JVS(83)*UV(27)+JVS(84)*UV(29)
  JUV(18) = JVS(85)*UV(5)+JVS(86)*UV(7)+JVS(87)*UV(13)+JVS(88)*UV(14)+JVS(89)*UV(15)+JVS(90)*UV(17)+JVS(91)*UV(18)&
              &+JVS(92)*UV(19)+JVS(93)*UV(20)+JVS(94)*UV(22)+JVS(95)*UV(23)+JVS(96)*UV(24)+JVS(97)*UV(25)+JVS(99)*UV(27)&
              &+JVS(100)*UV(28)+JVS(101)*UV(29)+JVS(102)*UV(30)+JVS(103)*UV(31)+JVS(104)*UV(32)
  JUV(19) = JVS(105)*UV(11)+JVS(106)*UV(14)+JVS(107)*UV(19)+JVS(108)*UV(25)+JVS(109)*UV(27)+JVS(111)*UV(31)
  JUV(20) = JVS(112)*UV(7)+JVS(113)*UV(13)+JVS(114)*UV(20)+JVS(115)*UV(22)+JVS(116)*UV(23)+JVS(117)*UV(25)+JVS(119)&
              &*UV(27)+JVS(120)*UV(29)+JVS(121)*UV(30)
  JUV(21) = JVS(122)*UV(17)+JVS(123)*UV(19)+JVS(124)*UV(21)+JVS(125)*UV(22)+JVS(126)*UV(23)+JVS(127)*UV(24)+JVS(128)&
              &*UV(25)+JVS(129)*UV(27)+JVS(130)*UV(28)+JVS(131)*UV(29)+JVS(132)*UV(30)+JVS(133)*UV(31)+JVS(134)*UV(32)
  JUV(22) = JVS(135)*UV(22)+JVS(136)*UV(25)+JVS(137)*UV(27)+JVS(138)*UV(29)+JVS(139)*UV(30)
  JUV(23) = JVS(140)*UV(22)+JVS(141)*UV(23)+JVS(142)*UV(25)+JVS(143)*UV(27)+JVS(144)*UV(29)+JVS(145)*UV(30)
  JUV(24) = JVS(146)*UV(13)+JVS(147)*UV(17)+JVS(148)*UV(19)+JVS(149)*UV(20)+JVS(150)*UV(22)+JVS(151)*UV(23)+JVS(152)&
              &*UV(24)+JVS(153)*UV(25)+JVS(155)*UV(27)+JVS(156)*UV(29)+JVS(157)*UV(30)
  JUV(25) = JVS(159)*UV(17)+JVS(160)*UV(19)+JVS(161)*UV(22)+JVS(162)*UV(23)+JVS(163)*UV(25)+JVS(164)*UV(26)+JVS(165)&
              &*UV(27)+JVS(166)*UV(28)+JVS(167)*UV(29)+JVS(169)*UV(31)
  JUV(26) = JVS(170)*UV(3)+JVS(171)*UV(4)+JVS(172)*UV(6)+JVS(173)*UV(9)+JVS(174)*UV(10)+JVS(175)*UV(11)+JVS(176)*UV(13)&
              &+JVS(178)*UV(18)+JVS(182)*UV(23)+JVS(184)*UV(25)+JVS(185)*UV(26)+JVS(186)*UV(27)+JVS(187)*UV(28)+JVS(188)&
              &*UV(29)+JVS(189)*UV(30)+JVS(190)*UV(31)+JVS(191)*UV(32)
  JUV(27) = JVS(192)*UV(1)+JVS(193)*UV(2)+JVS(194)*UV(5)+JVS(195)*UV(7)+JVS(196)*UV(9)+JVS(197)*UV(10)+JVS(198)*UV(12)&
              &+JVS(199)*UV(14)+JVS(200)*UV(15)+JVS(201)*UV(16)+JVS(202)*UV(17)+JVS(203)*UV(19)+JVS(204)*UV(20)+JVS(205)&
              &*UV(21)+JVS(206)*UV(22)+JVS(207)*UV(23)+JVS(208)*UV(24)+JVS(209)*UV(25)+JVS(210)*UV(26)+JVS(211)*UV(27)&
              &+JVS(212)*UV(28)+JVS(213)*UV(29)+JVS(215)*UV(31)+JVS(216)*UV(32)
  JUV(28) = JVS(217)*UV(2)+JVS(218)*UV(5)+JVS(219)*UV(7)+JVS(220)*UV(10)+JVS(221)*UV(11)+JVS(222)*UV(13)+JVS(223)*UV(14)&
              &+JVS(224)*UV(15)+JVS(225)*UV(16)+JVS(226)*UV(17)+JVS(227)*UV(19)+JVS(228)*UV(20)+JVS(229)*UV(21)+JVS(230)&
              &*UV(22)+JVS(231)*UV(23)+JVS(232)*UV(24)+JVS(233)*UV(25)+JVS(234)*UV(26)+JVS(235)*UV(27)+JVS(236)*UV(28)&
              &+JVS(237)*UV(29)+JVS(238)*UV(30)+JVS(239)*UV(31)+JVS(240)*UV(32)
  JUV(29) = JVS(241)*UV(1)+JVS(242)*UV(17)+JVS(243)*UV(21)+JVS(244)*UV(22)+JVS(245)*UV(23)+JVS(246)*UV(24)+JVS(247)&
              &*UV(25)+JVS(248)*UV(26)+JVS(251)*UV(29)+JVS(252)*UV(30)+JVS(253)*UV(31)
  JUV(30) = JVS(255)*UV(6)+JVS(256)*UV(12)+JVS(257)*UV(14)+JVS(258)*UV(21)+JVS(259)*UV(22)+JVS(260)*UV(23)+JVS(261)&
              &*UV(24)+JVS(262)*UV(25)+JVS(263)*UV(26)+JVS(264)*UV(27)+JVS(266)*UV(29)+JVS(267)*UV(30)+JVS(268)*UV(31)
  JUV(31) = JVS(270)*UV(8)+JVS(271)*UV(9)+JVS(272)*UV(11)+JVS(274)*UV(18)+JVS(280)*UV(25)+JVS(281)*UV(26)+JVS(282)&
              &*UV(27)+JVS(283)*UV(28)+JVS(284)*UV(29)+JVS(285)*UV(30)+JVS(286)*UV(31)+JVS(287)*UV(32)
  JUV(32) = JVS(288)*UV(3)+JVS(289)*UV(15)+JVS(290)*UV(19)+JVS(291)*UV(22)+JVS(292)*UV(24)+JVS(293)*UV(25)+JVS(294)&
              &*UV(26)+JVS(295)*UV(27)+JVS(296)*UV(28)+JVS(297)*UV(29)+JVS(298)*UV(30)+JVS(299)*UV(31)+JVS(300)*UV(32)

END SUBROUTINE Jac_SP_Vec

! End of Jac_SP_Vec function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! JacTR_SP_Vec - sparse multiplication: sparse Jacobian transposed times vector
!   Arguments :
!      JVS       - sparse Jacobian of variables
!      UV        - User vector for variables
!      JTUV      - Jacobian transposed times user vector
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE JacTR_SP_Vec ( JVS, UV, JTUV )

! JVS - sparse Jacobian of variables
  REAL(kind=dp) :: JVS(LU_NONZERO)
! UV - User vector for variables
  REAL(kind=dp) :: UV(NVAR)
! JTUV - Jacobian transposed times user vector
  REAL(kind=dp) :: JTUV(NVAR)

  JTUV(1) = JVS(1)*UV(1)+JVS(192)*UV(27)+JVS(241)*UV(29)
  JTUV(2) = JVS(3)*UV(2)+JVS(193)*UV(27)+JVS(217)*UV(28)
  JTUV(3) = JVS(6)*UV(3)+JVS(170)*UV(26)+JVS(288)*UV(32)
  JTUV(4) = JVS(9)*UV(4)+JVS(171)*UV(26)
  JTUV(5) = JVS(14)*UV(5)+JVS(38)*UV(11)+JVS(55)*UV(14)+JVS(85)*UV(18)+JVS(194)*UV(27)+JVS(218)*UV(28)
  JTUV(6) = JVS(16)*UV(6)+JVS(43)*UV(12)+JVS(172)*UV(26)+JVS(255)*UV(30)
  JTUV(7) = JVS(19)*UV(7)+JVS(39)*UV(11)+JVS(56)*UV(14)+JVS(62)*UV(15)+JVS(86)*UV(18)+JVS(112)*UV(20)+JVS(195)*UV(27)&
              &+JVS(219)*UV(28)
  JTUV(8) = JVS(21)*UV(8)+JVS(270)*UV(31)
  JTUV(9) = JVS(30)*UV(9)+JVS(173)*UV(26)+JVS(196)*UV(27)+JVS(271)*UV(31)
  JTUV(10) = JVS(34)*UV(10)+JVS(174)*UV(26)+JVS(197)*UV(27)+JVS(220)*UV(28)
  JTUV(11) = JVS(40)*UV(11)+JVS(57)*UV(14)+JVS(105)*UV(19)+JVS(175)*UV(26)+JVS(221)*UV(28)+JVS(272)*UV(31)
  JTUV(12) = JVS(44)*UV(12)+JVS(198)*UV(27)+JVS(256)*UV(30)
  JTUV(13) = JVS(22)*UV(8)+JVS(51)*UV(13)+JVS(87)*UV(18)+JVS(113)*UV(20)+JVS(146)*UV(24)+JVS(176)*UV(26)+JVS(222)*UV(28)
  JTUV(14) = JVS(10)*UV(4)+JVS(45)*UV(12)+JVS(58)*UV(14)+JVS(88)*UV(18)+JVS(106)*UV(19)+JVS(199)*UV(27)+JVS(223)*UV(28)&
               &+JVS(257)*UV(30)
  JTUV(15) = JVS(63)*UV(15)+JVS(68)*UV(16)+JVS(89)*UV(18)+JVS(200)*UV(27)+JVS(224)*UV(28)+JVS(289)*UV(32)
  JTUV(16) = JVS(69)*UV(16)+JVS(201)*UV(27)+JVS(225)*UV(28)
  JTUV(17) = JVS(70)*UV(16)+JVS(80)*UV(17)+JVS(90)*UV(18)+JVS(122)*UV(21)+JVS(147)*UV(24)+JVS(159)*UV(25)+JVS(202)&
               &*UV(27)+JVS(226)*UV(28)+JVS(242)*UV(29)
  JTUV(18) = JVS(91)*UV(18)+JVS(178)*UV(26)+JVS(274)*UV(31)
  JTUV(19) = JVS(64)*UV(15)+JVS(71)*UV(16)+JVS(92)*UV(18)+JVS(107)*UV(19)+JVS(123)*UV(21)+JVS(148)*UV(24)+JVS(160)&
               &*UV(25)+JVS(203)*UV(27)+JVS(227)*UV(28)+JVS(290)*UV(32)
  JTUV(20) = JVS(23)*UV(8)+JVS(52)*UV(13)+JVS(93)*UV(18)+JVS(114)*UV(20)+JVS(149)*UV(24)+JVS(204)*UV(27)+JVS(228)*UV(28)
  JTUV(21) = JVS(46)*UV(12)+JVS(72)*UV(16)+JVS(124)*UV(21)+JVS(205)*UV(27)+JVS(229)*UV(28)+JVS(243)*UV(29)+JVS(258)&
               &*UV(30)
  JTUV(22) = JVS(24)*UV(8)+JVS(65)*UV(15)+JVS(73)*UV(16)+JVS(81)*UV(17)+JVS(94)*UV(18)+JVS(115)*UV(20)+JVS(125)*UV(21)&
               &+JVS(135)*UV(22)+JVS(140)*UV(23)+JVS(150)*UV(24)+JVS(161)*UV(25)+JVS(206)*UV(27)+JVS(230)*UV(28)+JVS(244)&
               &*UV(29)+JVS(259)*UV(30)+JVS(291)*UV(32)
  JTUV(23) = JVS(25)*UV(8)+JVS(74)*UV(16)+JVS(95)*UV(18)+JVS(116)*UV(20)+JVS(126)*UV(21)+JVS(141)*UV(23)+JVS(151)*UV(24)&
               &+JVS(162)*UV(25)+JVS(182)*UV(26)+JVS(207)*UV(27)+JVS(231)*UV(28)+JVS(245)*UV(29)+JVS(260)*UV(30)
  JTUV(24) = JVS(47)*UV(12)+JVS(75)*UV(16)+JVS(96)*UV(18)+JVS(127)*UV(21)+JVS(152)*UV(24)+JVS(208)*UV(27)+JVS(232)&
               &*UV(28)+JVS(246)*UV(29)+JVS(261)*UV(30)+JVS(292)*UV(32)
  JTUV(25) = JVS(2)*UV(1)+JVS(66)*UV(15)+JVS(76)*UV(16)+JVS(82)*UV(17)+JVS(97)*UV(18)+JVS(108)*UV(19)+JVS(117)*UV(20)&
               &+JVS(128)*UV(21)+JVS(136)*UV(22)+JVS(142)*UV(23)+JVS(153)*UV(24)+JVS(163)*UV(25)+JVS(184)*UV(26)+JVS(209)&
               &*UV(27)+JVS(233)*UV(28)+JVS(247)*UV(29)+JVS(262)*UV(30)+JVS(280)*UV(31)+JVS(293)*UV(32)
  JTUV(26) = JVS(7)*UV(3)+JVS(11)*UV(4)+JVS(17)*UV(6)+JVS(31)*UV(9)+JVS(35)*UV(10)+JVS(48)*UV(12)+JVS(53)*UV(13)&
               &+JVS(164)*UV(25)+JVS(185)*UV(26)+JVS(210)*UV(27)+JVS(234)*UV(28)+JVS(248)*UV(29)+JVS(263)*UV(30)+JVS(281)&
               &*UV(31)+JVS(294)*UV(32)
  JTUV(27) = JVS(4)*UV(2)+JVS(12)*UV(4)+JVS(15)*UV(5)+JVS(20)*UV(7)+JVS(26)*UV(8)+JVS(32)*UV(9)+JVS(36)*UV(10)+JVS(41)&
               &*UV(11)+JVS(49)*UV(12)+JVS(54)*UV(13)+JVS(59)*UV(14)+JVS(67)*UV(15)+JVS(77)*UV(16)+JVS(83)*UV(17)+JVS(99)&
               &*UV(18)+JVS(109)*UV(19)+JVS(119)*UV(20)+JVS(129)*UV(21)+JVS(137)*UV(22)+JVS(143)*UV(23)+JVS(155)*UV(24)&
               &+JVS(165)*UV(25)+JVS(186)*UV(26)+JVS(211)*UV(27)+JVS(235)*UV(28)+JVS(264)*UV(30)+JVS(282)*UV(31)+JVS(295)&
               &*UV(32)
  JTUV(28) = JVS(5)*UV(2)+JVS(37)*UV(10)+JVS(100)*UV(18)+JVS(130)*UV(21)+JVS(166)*UV(25)+JVS(187)*UV(26)+JVS(212)*UV(27)&
               &+JVS(236)*UV(28)+JVS(283)*UV(31)+JVS(296)*UV(32)
  JTUV(29) = JVS(27)*UV(8)+JVS(78)*UV(16)+JVS(84)*UV(17)+JVS(101)*UV(18)+JVS(120)*UV(20)+JVS(131)*UV(21)+JVS(138)*UV(22)&
               &+JVS(144)*UV(23)+JVS(156)*UV(24)+JVS(167)*UV(25)+JVS(188)*UV(26)+JVS(213)*UV(27)+JVS(237)*UV(28)+JVS(251)&
               &*UV(29)+JVS(266)*UV(30)+JVS(284)*UV(31)+JVS(297)*UV(32)
  JTUV(30) = JVS(13)*UV(4)+JVS(18)*UV(6)+JVS(28)*UV(8)+JVS(50)*UV(12)+JVS(60)*UV(14)+JVS(79)*UV(16)+JVS(102)*UV(18)&
               &+JVS(121)*UV(20)+JVS(132)*UV(21)+JVS(139)*UV(22)+JVS(145)*UV(23)+JVS(157)*UV(24)+JVS(189)*UV(26)+JVS(238)&
               &*UV(28)+JVS(252)*UV(29)+JVS(267)*UV(30)+JVS(285)*UV(31)+JVS(298)*UV(32)
  JTUV(31) = JVS(29)*UV(8)+JVS(33)*UV(9)+JVS(42)*UV(11)+JVS(103)*UV(18)+JVS(111)*UV(19)+JVS(133)*UV(21)+JVS(169)*UV(25)&
               &+JVS(190)*UV(26)+JVS(215)*UV(27)+JVS(239)*UV(28)+JVS(253)*UV(29)+JVS(268)*UV(30)+JVS(286)*UV(31)+JVS(299)&
               &*UV(32)
  JTUV(32) = JVS(8)*UV(3)+JVS(104)*UV(18)+JVS(134)*UV(21)+JVS(191)*UV(26)+JVS(216)*UV(27)+JVS(240)*UV(28)+JVS(287)&
               &*UV(31)+JVS(300)*UV(32)

END SUBROUTINE JacTR_SP_Vec

! End of JacTR_SP_Vec function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

END MODULE cbm4_Jacobian


text:
module test_aomp

implicit none

integer, parameter :: rstd = 8
integer :: nsize
REAL(rstd), allocatable :: a_dev(:,:,:), b_dev(:,:,:), c_dev(:,:,:)
logical :: touch_limit
REAL(rstd) :: limit
integer :: use_gpu = 1, count = 0, count_cpu = 0, count_gpu = 0
!$acc declare create(a_dev,b_dev,c_dev,nsize)
!!$omp declare target(a_dev,b_dev,c_dev,nsize)

contains
    subroutine init()
        integer i,j,k
        nsize = 100

        allocate(a_dev(nsize,nsize,nsize))
        allocate(b_dev(nsize,nsize,nsize))
        allocate(c_dev(nsize,nsize,nsize))
        do i=1,nsize
            do j=1,nsize
                do k=1,nsize
                    b_dev(i,j,k) = 0
                    c_dev(i,j,k) = 0
                end do
            end do
        end do
    end subroutine init

    subroutine print_output()
        integer i,j,k
        do i=1,nsize
            do j=1,nsize
                do k=1,nsize
                    !print *, "a(", i, ",", j, ",", k, ") = ", a_dev(i,j,k)
                end do
            end do
        end do
        print *, "touch limit: ", touch_limit
        print *, "count: ", count
    end subroutine print_output

    subroutine _compute_dev_cpu()
        integer i,j,k
!$omp parallel do reduction(.or.:touch_limit) reduction(+:count)
        do i=1,nsize
            do j=1,nsize
                do k=1,nsize
                    a_dev(i,j,k) = b_dev(i,j,k) * c_dev(i,j,k) * i * nsize*nsize + j * nsize + k
                    if (a_dev(i,j,k) > limit) then
                        touch_limit = touch_limit .or. .true.
                        count = count + 1
                    else
                        touch_limit = touch_limit .or. .false.
                        count = count + 0
                    end if
                end do
            end do
        end do
    end subroutine _compute_dev_cpu

    subroutine _compute_dev_gpu()
        integer i,j,k
!$omp target teams distribute parallel do reduction(.or.:touch_limit) reduction(+:count) if(target:use_gpu)
        do i=1,nsize
            do j=1,nsize
                do k=1,nsize
                    a_dev(i,j,k) = b_dev(i,j,k) * c_dev(i,j,k) * i * nsize*nsize + j * nsize + k
                    if (a_dev(i,j,k) > limit) then
                        touch_limit = touch_limit .or. .true.
                        count = count + 1
                    else
                        touch_limit = touch_limit .or. .false.
                        count = count + 0
                    end if
                end do
            end do
        end do
    end subroutine _compute_dev_gpu

    subroutine compute_dev()
!!$omp target update to(b_dev,c_dev,nsize,limit)
        CALL _compute_dev_cpu()
        print *, "count CPU: ", count
        count_cpu = count
        CALL _compute_dev_gpu()
        print *, "count GPU: ", count
        CALL _compute_dev_gpu()
        count_gpu = count
        if (count_cpu == count_gpu) then
            print *, "PASS"
        else
            print *, "FAIL"
        end if
!!$omp target update from(a_dev,count,touch_limit)
    end subroutine compute_dev

end module test_aomp

program test
   use test_aomp

   character(len=12) :: args

   limit= 10099
   touch_limit= .false.

   call get_command_argument(1, args)
   read(args,'(i)') use_gpu

   call init()
   call compute_dev()
end program test

max_stars_repo_name: aperdeus-amd/aomp

max_issues_repo_name: aperdeus-amd/aomp

max_forks_repo_name: aperdeus-amd/aomp


     SUBROUTINE LA_TEST_SGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, &
     &  WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
!
!  -- LAPACK95 interface driver routine (version 0.0) --
!     UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK
!     May 21, 1999
!
!  .. Use Statements ..
      USE LA_PRECISION, ONLY: WP => SP
      USE F95_LAPACK, ONLY: LA_GEESX
!  .. Implicit Statement ..
      IMPLICIT NONE
!  .. Scalar Arguments ..
      INTEGER, INTENT(IN) :: N, LDA, LDVS, LWORK, LIWORK
      INTEGER, INTENT(OUT) :: SDIM
      INTEGER, INTENT(INOUT) :: INFO
      CHARACTER*1, INTENT(IN) :: JOBVS, SORT, SENSE
      REAL(WP), INTENT(OUT) :: RCONDV, RCONDE
!   LOGICAL SELECT
      INTERFACE
      LOGICAL FUNCTION SELECT(WR, WI)
      USE LA_PRECISION, ONLY: WP => SP
      REAL(WP), INTENT(IN) :: WR, WI
      END FUNCTION SELECT
      END INTERFACE
!  .. Array Arguments ..
      REAL(WP), INTENT(INOUT) :: A(1:LDA,1:N)
      REAL(WP), INTENT(OUT) :: WR(1: N), WI(1: N), VS(1: LDVS, 1: N), &
     &  WORK (1: LWORK)
      LOGICAL :: BWORK(1: N)
      INTEGER, INTENT(OUT) :: IWORK(1: LIWORK)
      LOGICAL LSAME
!  .. Parameters ..
      CHARACTER(LEN=8),  PARAMETER :: SRNAME = 'LA_GEESX'
      CHARACTER(LEN=14), PARAMETER :: SRNAMT = 'LA_TEST_SGEESX'
!  .. Common blocks ..
      INTEGER :: INFOTC
      COMMON /LINFO95/ INFOTC
!  .. Local Scalars ..
      INTEGER :: I, J, IA1, IA2, ISDIM, IWI, IWR, IVS1, IVS2
      CHARACTER*1 :: IJOBVS
!  .. Local Arrays ..
      LOGICAL, SAVE :: CTEST = .TRUE., ETEST = .TRUE.
!  .. Executable Statements ..
      IA1 = N; IA2 = N; IWI = N; IWR = N;
      IVS1 = N; IVS2 = N; ISDIM = SDIM; IJOBVS = JOBVS

      I = INFO / 100; J = INFO - I*100
      SELECT CASE(I)
      CASE (1)
      IA2 = IA2-1
      CASE (2)
      IWR = IWR-1
      CASE(3)
      IWI = IWI - 1
      CASE (4)
      IVS1 = IVS1 - 1
      IJOBVS = 'V'
      CASE(:-1,5:)
      CALL UESTOP(SRNAMT)
      END SELECT
      SDIM = 0
      if (lsame(sense,'E')) then
        IF (LSAME (IJOBVS,'V')) THEN
          IF (LSAME(SORT, 'S')) THEN
            CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
     &        VS(1: IVS1, 1: IVS2), SELECT, SDIM, RCONDE, INFO=INFO)
          ELSE
            CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
     &        VS(1: IVS1, 1: IVS2), RCONDE=RCONDE, INFO=INFO)
          END IF
        ELSE
          IF (LSAME(SORT, 'S')) THEN
            CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
     &        SELECT=SELECT, SDIM=SDIM, RCONDE=RCONDE, INFO=INFO)
          ELSE
            CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
     &        RCONDE=RCONDE, INFO=INFO)
          END IF
        END IF
      else if (lsame (sense,'V')) then
        IF (LSAME (IJOBVS,'V')) THEN
          IF (LSAME(SORT, 'S')) THEN
            CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
     &        VS(1: IVS1, 1: IVS2), SELECT, SDIM, RCONDV=RCONDV, &
     &        INFO=INFO)
          ELSE
            CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
     &        VS(1: IVS1, 1: IVS2), RCONDV=RCONDV, INFO=INFO)
          END IF
        ELSE
          IF (LSAME(SORT, 'S')) THEN
            CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
     &        SELECT=SELECT, SDIM=SDIM, RCONDV=RCONDV, INFO=INFO)
          ELSE
            CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
     &        RCONDV=RCONDV, INFO=INFO)
          END IF
        END IF
      else
        if (lsame (sense,'B')) then
          IF (LSAME (IJOBVS,'V')) THEN
            IF (LSAME(SORT, 'S')) THEN
              CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
&               VS(1: IVS1, 1: IVS2), SELECT, SDIM, RCONDE, RCONDV, &
&               INFO=INFO)
            ELSE
              CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
&               VS(1: IVS1, 1: IVS2), RCONDE=RCONDE, RCONDV=RCONDV, &
&               INFO=INFO)
            END IF
          ELSE
            IF (LSAME(SORT, 'S')) THEN
              CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
&               SELECT=SELECT, SDIM=SDIM, RCONDE=RCONDE, RCONDV=RCONDV, &
&               INFO=INFO)
            ELSE
              CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
&               RCONDE=RCONDE, RCONDV=RCONDV, INFO=INFO)
            END IF
          END IF
        else
          IF (LSAME (IJOBVS,'V')) THEN
            IF (LSAME(SORT, 'S')) THEN
              CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
&               VS(1: IVS1, 1: IVS2), SELECT, SDIM, INFO=INFO)
            ELSE
              CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
&               VS(1: IVS1, 1: IVS2), INFO=INFO)
            END IF
          ELSE
            IF (LSAME(SORT, 'S')) THEN
              CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
&               SELECT=SELECT, SDIM=SDIM, INFO=INFO)
            ELSE
              CALL LA_GEESX( A(1:IA1, 1: IA2), WR(1:IWR), WI(1:IWI), &
&               INFO=INFO)
            END IF
          END IF
        end if
      endif
      CALL LA_AUX_AA01( I, CTEST, ETEST, SRNAMT )
END SUBROUTINE LA_TEST_SGEESX


! { dg-do compile }
! See PR fortran/36251.
module a
  implicit none
  integer :: i = 42
end module a

module b
  use a
  implicit none
  public a  ! { dg-error "attribute applied to" }
end module b

module d
  use a
  implicit none
  private a  ! { dg-error "attribute applied to" }
end module d


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!     subroutine rdsolvff_two: read in informations for QMDFF2
!
!     part of QMDFF
!
subroutine rdsolvff_two(n,xyz,at,q,imass,d,scalehb,scalexb,fname)
use qmdff
implicit none
integer::n,at(n),imass(n)
real(kind=8)::xyz(3,n),q(n),d
real(kind=8)::scalehb(94),scalexb(94)
character(len=*)::fname

integer::i,j,nn,nterm,idum
real(kind=8)::qdum,r,xx(10),dum1,dum2
character(len=80)::atmp

open(unit=1,file=fname)

read (1,'(a)') atmp
call readl(atmp,xx,nn)
n=idint(xx(1))
qdum=xx(2)
d=xx(3)

if (nn.eq.3.and.d.lt.0) then
   d=xx(3)
   write(*,*) 'taking solvent density from file ',d
end if
read (1,'(a)') atmp

do i=1,n
   read (1,*)at(i),xyz(1:3,i),q(i),imass(i)
end do
read (1,*) nbond_two,nangl_two,ntors_two,nhb_two,nnci_two
if (nbond_two.gt.ndim.or.nangl_two.gt.ndim.or.ntors_two.gt.ndim) &
    &  stop 'too many FF terms'
do i=1,nbond_two
   read (1,*)bond_two(1:2,i),vbond_two(1:3,i)
end do
do i=1,nangl_two
   read (1,*)angl_two(1,i),angl_two(2,i),angl_two(3,i), &
     &     vangl_two(1:2,i)
enddo
do i=1,ntors_two
   read (1,*)tors_two(1:6,i),vtors_two(1:2,i), &
       &   (vtors_two(3*(j-1)+3,i), &
       &    vtors_two(3*(j-1)+4,i), &
       &    vtors_two(3*(j-1)+5,i),j=1,tors_two(5,i))
end do
do i=1,ntors_two
   do j=1,tors_two(5,i)
      vtors_two(3*(j-1)+4,i)=vtors_two(3*(j-1)+4,i)*pi
   end do
end do
if(nhb_two.gt.0) read(1,'(6(3i4,2x))')hb_two(1:3,1:nhb_two)
do i=1,nhb_two
   if (at(hb_two(3,i)).eq.1) then
!
!     subroutine hbpara is totally FF-unspecific!
!
      call hbpara(10.0d0,5.0d0,q(hb_two(1,i)),dum1)
      vhb_two(1,i)=dum1*scalehb(at(hb_two(1,i)))
      call hbpara(10.0d0,5.0d0,q(hb_two(2,i)),dum1)
      vhb_two(2,i)=dum1*scalehb(at(hb_two(2,i)))
   else
      call hbpara(-6.5d0,1.0d0,q(hb_two(3,i)),dum1)
      vhb_two(1,i)=scalexb(at(hb_two(3,i)))*dum1
   end if
end do
read(1,'(8(3i4,2x))')nci_two(1:3,1:nnci_two)
close(1)

q=q*qdum

return
end subroutine rdsolvff_two

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense