text:
MODULE STMAS4D_CORE

  USE PRMTRS_STMAS
  USE GENERALTOOLS
  USE PREP_STMAS4D
  USE POST_STMAS4D
  USE COSTFUN_GRAD
  USE WCOMPT_GRADT,  ONLY : WCOMPGERNL

CONTAINS

SUBROUTINE MGANALYSS
!*************************************************
! MULTI-GRID ANALYSIS
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
!------------------------
  INTEGER       :: IT,FW
!                  FW=1: FROM COARSE GRID TO FINE GRID
!                  FW=0: FROM FINE GRID TO COARSE GRID
   INTEGER  :: I,J,K,T,S

  INTEGER  ::NUM,NO,N,O   !shuyuan 20100901
  INTEGER  :: NP(MAXDIMS)  !shuyuan
!------------------------

  PRINT*,'PREPROCSS'
  CALL PREPROCSS

  GRDLEVL=0
  IT=0
  FW=1
  NUM=0
    DO WHILE(.TRUE.)

      GRDLEVL=GRDLEVL+1
      PRINT*,'NUMBER',GRDLEVL,'LEVEL GRID IS IN PROCESSING'
      CALL RDINITBGD
      PRINT*,'PHYSCPSTN'
      CALL PHYSCPSTN
      PRINT*,'GETCOEFFT'
      CALL GETCOEFFT

      IF(NUMVARS.EQ.0)EXIT

      PRINT*,'MINIMIZER'
      CALL MINIMIZER_XIE
      IF(GRDLEVL.EQ.FNSTGRD .AND. IT.EQ.ITREPET)EXIT

      IF(IFREPET.EQ.0) THEN
        PRINT*,'COAS2FINE'
        CALL COAS2FINE_XIE		! USE XIE'S AS THE OLD HAS ERROR
      ELSE
        IF(GRDLEVL.EQ.FNSTGRD) FW=0
        IF(GRDLEVL.EQ.1) THEN
          FW=1
          IT=IT+1
        ENDIF
        IF(FW.EQ.1) THEN
          PRINT*,'COAS2FINE'
          CALL COAS2FINE_XIE		! USE XIE'S AS THE OLD HAS ERROR
        ELSE
          GRDLEVL=GRDLEVL-2
          PRINT*,'FINE2COAS'
          CALL FINE2COAS_XIE		! USE XIE'S AS THE OLD HAS ERROR
        ENDIF
      ENDIF

    ENDDO
!   CALL GETW
  PRINT*,'PSTPROCSS'
  CALL PSTPROCSS
  RETURN
END SUBROUTINE MGANALYSS

SUBROUTINE TMPMEMALC
!*************************************************
! MEMORY ALLOCATE FOR TMP ARRAY
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,S,ER
! --------------------
  ALLOCATE(TMPANALS(NTMPGRD(1),NTMPGRD(2),NTMPGRD(3),NTMPGRD(4),NUMSTAT),STAT=ER)
  IF(ER.NE.0)STOP 'TMPANALS ALLOCATE WRONG'
  DO T=1,NTMPGRD(4)
  DO K=1,NTMPGRD(3)
  DO J=1,NTMPGRD(2)
  DO I=1,NTMPGRD(1)
    DO S=1,NUMSTAT
      TMPANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)
    ENDDO
  ENDDO
  ENDDO
  ENDDO
  ENDDO
  RETURN
END SUBROUTINE TMPMEMALC

SUBROUTINE TMPMEMRLS
!*************************************************
! MEMORY RELEASE FOR TMP ARRAY
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
  DEALLOCATE(TMPANALS)
  RETURN
END SUBROUTINE TMPMEMRLS

SUBROUTINE GETCOEFFT
!*************************************************
! GET GRID POINT COEFFICENT FOR EVERY OBSERVATION
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,O,M,N,UV,RR
  INTEGER  :: LM(MAXDIMS)
  INTEGER  :: NP(MAXDIMS),NN(MAXDIMS)
  REAL     :: AC(NUMDIMS,NGPTOBS),OC(NUMDIMS),CO(NGPTOBS)
! --------------------
! GET INTERPOLATION COEFFICENT
  UV=0
  RR=0
  IF(NALLOBS.EQ.0) RETURN
  DO O=1,NALLOBS
    DO N=1,MAXDIMS
      NP(N)=1
    ENDDO
    DO N=1,NUMDIMS
      NP(N)=INT((OBSPOSTN(N,O)-ORIPSTN(N))/GRDSPAC(N))+1
      IF(NP(N).EQ.NUMGRID(N).AND.NUMGRID(N).NE.1)NP(N)=NUMGRID(N)-1
      IF(NUMGRID(N).EQ.1)NP(N)=1
    ENDDO
    DO N=1,MAXDIMS
      OBSIDXPC(N,O)=NP(N)
    ENDDO
    M=0
!====================================================
!    DO T=NP(4),MIN0(NP(4)+1,NUMGRID(4))
!    DO K=NP(3),MIN0(NP(3)+1,NUMGRID(3))
!    DO J=NP(2),MIN0(NP(2)+1,NUMGRID(2))
!    DO I=NP(1),MIN0(NP(1)+1,NUMGRID(1))
!      NN(1)=I
!      NN(2)=J
!      NN(3)=K
!      NN(4)=T
!      M=M+1
!      DO N=1,NUMDIMS
!        AC(N,M)=NN(N)*1.0
!      ENDDO
!    ENDDO
!    ENDDO
!    ENDDO
!    ENDDO
!    DO N=1,NUMDIMS
!      OC(N)=(OBSPOSTN(N,O)-ORIPSTN(N))/GRDSPAC(N)+1
!    ENDDO
!================== modified by zhongjie he ========
    DO N=1,MAXDIMS
      LM(N)=NP(N)+1
      IF(NUMDIMS.LT.N) LM(N)=NP(N)
    ENDDO
    DO T=NP(4),LM(4)
    DO K=NP(3),LM(3)
    DO J=NP(2),LM(2)
    DO I=NP(1),LM(1)
      NN(1)=MIN0(I,NUMGRID(1))
      NN(2)=MIN0(J,NUMGRID(2))
      NN(3)=MIN0(K,NUMGRID(3))
      NN(4)=MIN0(T,NUMGRID(4))
      M=M+1
      DO N=1,NUMDIMS
        AC(N,M)=NN(N)*1.0
      ENDDO
    ENDDO
    ENDDO
    ENDDO
    ENDDO
    DO N=1,NUMDIMS
      IF(NUMGRID(N).GE.2) THEN
        OC(N)=(OBSPOSTN(N,O)-ORIPSTN(N))/GRDSPAC(N)+1
      ELSE
        OC(N)=1
      ENDIF
    ENDDO
!====================================================
    ! CALL INTERPLTN_XIE(NUMDIMS,NGPTOBS,CO,AC,OC,3,NUMGRID(3),PPM)
    CALL INTERPLTN(NUMDIMS,NGPTOBS,CO,AC,OC)
    DO M=1,NGPTOBS
      OBSCOEFF(M,O)=CO(M)
    ENDDO
  ENDDO
  UV=NOBSTAT(U_CMPNNT)+NOBSTAT(V_CMPNNT)
  RR=NOBSTAT(NUMSTAT+1)
  OBSRADAR = 1.0
  IF(RR.NE.0)OBSRADAR=100.0*UV/FLOAT(RR)
  RR=NOBSTAT(NUMSTAT+2)
  OBS_SFMR = 1.0
  IF(RR.NE.0)OBS_SFMR=1.0/RR
!jhui
  RR=NOBSTAT(NUMSTAT+3)
  IF(RR.NE.0)OBSREF=1.0/RR

  RETURN
END SUBROUTINE GETCOEFFT

SUBROUTINE MINIMIZER_XIE
!*************************************************
! MINIMIZE THE COST FUNCTION
! HISTORY: AUGUST 2007, CODED by WEI LI.
!
!          MODIFIED BY YUANFU FOR USING A SINGLE
!          PRECISION LBFGS ROUTINE
!
!          DECEMBER 2013 MODIFIED BY YUANFU FOR
!          a) CHANGING AUTOMATIC ARRAYS TO ALLOCATABLES
!          b) SWITCHING TO LBFGSB 3.0 WITH YUANFU'S
!             MODIFICATION FOR SUPER LARGE MINIMIZATION
!          c) REMOVING INITIALIZATION OF IW and WA FOR
!             EFFICIENCY
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER ,PARAMETER :: MM=5
  CHARACTER(LEN=60) :: TA,CS
  LOGICAL           :: LS(4)
  INTEGER    :: I0,IC,IP,IT,ISBMN,N,O,S,T,K,J,I,NO,ER
  INTEGER    :: IS(44)
  INTEGER, ALLOCATABLE :: NB(:),IW(:)
  INTEGER    :: NN(MAXDIMS+1),NG(MAXDIMS+1),NC
  REAL       :: MN(NUMSTAT+3),MX(NUMSTAT+3)

  REAL, ALLOCATABLE :: LB(:), UB(:)
  REAL :: FA, PG, DS(29)
  REAL,ALLOCATABLE :: WA(:)

  REAL  ::temp,temp1,temp2,dif1,dif2,temp0
  INTEGER::  locx,locy,locz,loct,locv
! --------------------

  ! ALLOCATE WORKING ARRAY:
  ALLOCATE(WA(2*MM*NUMVARS+5*NUMVARS+12*MM*MM+12*MM), &
           LB(NUMVARS),UB(NUMVARS),NB(NUMVARS),IW(3*NUMVARS),STAT=ER)
  IF (ER .NE. 0) THEN
    PRINT*,'MINIMIZER: Cannot allocate enough memory for the working array'
    STOP
  ELSE
    PRINT*,'Successfully allocate memory for LBFGS!'
  ENDIF

  IP=1
  FA=1.0d+2
  PG=1.0d-10
  ISBMN = 1	! SUBSPACE MINIMIZATION
  DO N=1,NUMVARS
    NB(N)=0
    LB(N)=0.0
    UB(N)=0.0
  ENDDO

  ! BOUND: IF IFBOUND EQ 1, USE (MIN(OBS), MAX(OBS)) TO BOUND
  !        ANALYSIS
  IF(IFBOUND.EQ.1 .AND. NALLOBS.GE.1)THEN

    ! FIND THE MIN(OBS) AND MAX(OBS):

    ! BOUND AT LEAST ALLOW ZERO INCREMENT: BY YUANFU CORRECTING ZHONGJIE'S SETTING
    MN = 0.0
    MX = 0.0

    ! CONVENTIONAL OBS:
    DO N=1,MAXDIMS
      NG(N)=NUMGRID(N)
    ENDDO
    NG(MAXDIMS+1)=NUMSTAT+1
    O=0
    DO S=1,NUMSTAT
      DO NO=1,NOBSTAT(S)
        O=O+1
        MN(S)=MIN(MN(S),OBSVALUE(O))
        MX(S)=MAX(MX(S),OBSVALUE(O))
      ENDDO
    ENDDO

    ! RADAR AND SFMR:
    S=NUMSTAT+1
    DO NO=1,NOBSTAT(S)
      O=O+1
      MN(U_CMPNNT)=MIN(MN(U_CMPNNT),OBSVALUE(O))
      MX(U_CMPNNT)=MAX(MX(U_CMPNNT),OBSVALUE(O))
      MN(V_CMPNNT)=MIN(MN(V_CMPNNT),OBSVALUE(O))
      MX(V_CMPNNT)=MAX(MX(V_CMPNNT),OBSVALUE(O))
    ENDDO
    S=NUMSTAT+2
    DO NO=1,NOBSTAT(S)
      O=O+1
      MN(U_CMPNNT)=MIN(MN(U_CMPNNT),-1.*OBSVALUE(O))
      MX(U_CMPNNT)=MAX(MX(U_CMPNNT),OBSVALUE(O))
      MN(V_CMPNNT)=MIN(MN(V_CMPNNT),-1.*OBSVALUE(O))
      MX(V_CMPNNT)=MAX(MX(V_CMPNNT),OBSVALUE(O))
    ENDDO
!jhui
!----------------------------------
!changed by shuyuan 20100722
!      MN(6) =0.0
!      MX(6) =0.0

    S=NUMSTAT+3
    MN(S) =0.0
    MX(S) =0.0
    DO NO=1,NOBSTAT(S)
      O=O+1
      MN(S)=MIN(MN(S),OBSVALUE(O))
      MX(S)=MAX(MX(S),OBSVALUE(O))
    ENDDO

    DO S=6,7
     MN(S) =0.0
     MX(S) =0.0
     DO NO=1,NOBSTAT(S)
      O=O+1
      MN(S)=0.
      MX(S)=1000.!just for test
     ENDDO
    ENDDO
!--------------------------------------------------
    DO S=1,5!!NUMSTAT+1
      DO T=1,NUMGRID(4)
      DO K=1,NUMGRID(3)
      DO J=1,NUMGRID(2)
      DO I=1,NUMGRID(1)
        NN(1)=I
        NN(2)=J
        NN(3)=K
        NN(4)=T
        NN(MAXDIMS+1)=S
        CALL PSTN2NUMB(MAXDIMS+1,NN,NG,NC)
        IF(MN(S).LE.MX(S))THEN
          NB(NC)=2
          LB(NC)=MN(S)
          UB(NC)=MX(S)
        ELSE
          NB(NC)=0
          LB(NC)=0.0D0
          UB(NC)=0.0D0
        ENDIF
      ENDDO
      ENDDO
      ENDDO
      ENDDO
    ENDDO

    IF(IFBKGND.EQ.1)THEN
      DO N=1,NUMVARS
        IF(NB(N).EQ.2)THEN
          LB(N)=MIN(LB(N),0.0)
          UB(N)=MAX(UB(N),0.0)
        ENDIF
      ENDDO
    ENDIF

  ENDIF

! BOUND
! FOR 3D RADAR
  IF(W_CMPNNT.NE.0)THEN
    DO N=1,MAXDIMS
      NG(N)=NUMGRID(N)
    ENDDO
    NG(MAXDIMS+1)=NUMSTAT
    S=W_CMPNNT
    K=1
    DO T=1,NUMGRID(4)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      NN(1)=I
      NN(2)=J
      NN(3)=K
      NN(4)=T
      NN(MAXDIMS+1)=S
      CALL PSTN2NUMB(MAXDIMS+1,NN,NG,NC)
      NB(NC)=2
      LB(NC)=0.0
      UB(NC)=0.0
    ENDDO
    ENDDO
    ENDDO
  ENDIF

  TA='START'
  I0=0
  IC=0
  WRITE(*,9001)

  IF(W_CMPNNT.NE.0) THEN         ! BY ZHONGJIE HE
    CALL COSTFUNCT1
    CALL COSTGRADT1
  ELSE                           ! BY ZHONGJIE HE
    CALL WCOMPGERNL

    CALL COSTFUNCT2
    CALL COSTGRADT2
  ENDIF
print*,'minvalue of bk: ',minval(grdbkgnd(1:numgrid(1),1:numgrid(2),1:numgrid(3),1:numgrid(4),5))

  ! COUNT NUMBER OF ITERATIONS:
  IT = 0
  ITERLOOP: DO WHILE (.TRUE.)

    ! ADD LOW BOUND FOR SPECIFIC HUMIDITY:
    NG(1:4) = NUMGRID(1:4)
    NG(5) = NUMSTAT  ! THE NUMBER OF STATES IN THE 5TH DIMENSION
    DO T=1,NUMGRID(4)
      DO K=1,NUMGRID(3)
        DO J=1,NUMGRID(2)
          DO I=1,NUMGRID(1)
            NN(1)=I
            NN(2)=J
            NN(3)=K
            NN(4)=T
            NN(5)=5  ! POSITION OF SPECIFIC HUMIDITY

            CALL PSTN2NUMB(5,NN,NG,NC)

            NB(NC) = 0
            IF ((MAXGRID(1)-1)/2+1 .LE. NUMGRID(1)) NB(NC) = 1

            ! LOW BOUND FROM REFLECTIVITY DID NOT GET SCALED AND SO HERE IT DOES:
            LB(NC) = -GRDBKGND(I,J,K,T,5)+GRDBKGND(I,J,K,T,NUMSTAT+1)
            UB(NC) = -GRDBKGND(I,J,K,T,5)+GRDBKGND(I,J,K,T,NUMSTAT+2)
          ENDDO
        ENDDO
      ENDDO
    ENDDO
    ! ADD LOWER BOUNDS FOR RAIN AND SNOW:
    IF (NUMSTAT .GT. 5) THEN ! RAIN AND SNOW ARE 6TH AND 7TH:
    DO T=1,NUMGRID(4)
      DO K=1,NUMGRID(3)
        DO J=1,NUMGRID(2)
          DO I=1,NUMGRID(1)
            NN(1)=I
            NN(2)=J
            NN(3)=K
            NN(4)=T
            NN(5)=ROUR_CMPNNT  ! POSITION OF RAIN

            CALL PSTN2NUMB(5,NN,NG,NC)

            NB(NC) = 1
            LB(NC) = -GRDBKGND(I,J,K,T,ROUR_CMPNNT)

            NN(5)=ROUS_CMPNNT  ! POSITION OF SNOW

            CALL PSTN2NUMB(5,NN,NG,NC)

            NB(NC) = 1
            LB(NC) = -GRDBKGND(I,J,K,T,ROUS_CMPNNT)
          ENDDO
        ENDDO
      ENDDO
    ENDDO
    ENDIF

    CALL SETULB(NUMVARS,MM,GRDANALS,LB,UB,NB,COSTFUN,GRADINT,FA,PG,WA,IW, &
                TA,IP,CS,LS,IS,DS)

    IF(TA(1:2).EQ.'FG')THEN
      IF(W_CMPNNT.NE.0) THEN         ! BY ZHONGJIE HE
        CALL COSTFUNCT1
        CALL COSTGRADT1
      ELSE                           ! BY ZHONGJIE HE
        CALL WCOMPGERNL
        CALL COSTFUNCT2
        CALL COSTGRADT2
      ENDIF
      IC=IC+1
      CYCLE ITERLOOP
    ELSEIF(TA(1:5).EQ.'NEW_X')THEN
      IF(I0.EQ.0)THEN
        WRITE(*,1003)IS(30),IS(34),DS(13),COSTFUN
      ELSE
        WRITE(*,3001)IS(30),IS(34),DS(14),DS(13),COSTFUN
      ENDIF
      I0=I0+1
      IT = IT+1

      ! Exceeds maximum iterations:
      IF(GRDLEVL.LE.MIDGRID.AND.IT.GT.COSSTEP) EXIT ITERLOOP
      IF(GRDLEVL.GT.MIDGRID.AND.IT.GT.FINSTEP) EXIT ITERLOOP

      CYCLE ITERLOOP
    ELSE
      IF(IP.LE.-1.AND.TA(1:4).NE.'STOP')WRITE(6,*)TA
      EXIT ITERLOOP
    ENDIF
  END DO ITERLOOP
  1003 FORMAT (2(1x,i4),5x,'-',3x,1p,2(1x,d10.3))
  9001 FORMAT (/,3x,'it',3x,'nf',2x,'stepl',5x,'projg',8x,'f')
  3001 FORMAT (2(1x,i4),1p,2x,d8.1,1p,2(1x,d10.3))

  IF(W_CMPNNT.NE.0) THEN         ! BY ZHONGJIE HE
    CALL COSTFUNCT1
!    CALL COSTGRADT1
  ELSE                           ! BY ZHONGJIE HE
    CALL WCOMPGERNL
    CALL COSTFUNCT2

!    CALL COSTGRADT2
  ENDIF

  ! DEALLOCATE WORKING ARRAY:
  DEALLOCATE(WA,LB,UB,NB,IW,STAT=ER)
  IF (ER .NE. 0) THEN
    PRINT*,'MINIMIZER: Cannot deallocate enough memory for the working array'
    STOP
  ENDIF

  WRITE(8,*)'BOTTOM FUNCTION =',COSTFUN
  RETURN
END SUBROUTINE MINIMIZER_XIE

SUBROUTINE COAS2FINE_XIE
!*************************************************
! INTERPOLATION FROM COARSE GRID TO FINE GRID
! HISTORY: AUGUST 2007, CODED by WEI LI.
!		MODIFIED BY YUANFU FROM WEI LI'S
!		COAS2FINE,WHICH IS INCORRECTLY
!		INTERPOLATING COASE TO FINE GRID.
!
!	NOTE: THIS ROUTINE ASSUME FINE GRID HALVES
!		THE COASE RESOLUTION.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,S,N,I0,J0,K0,T0,SM
  INTEGER  :: IC(MAXDIMS)
! --------------------
  DO N=1,MAXDIMS
    NTMPGRD(N)=NUMGRID(N)
  ENDDO
  CALL TMPMEMALC
  CALL GRDMEMRLS
  DO N=1,MAXDIMS
    IC(N)=1
    IF(NTMPGRD(N).LT.MAXGRID(N)) THEN
!      IF(N.NE.3 .OR. GRDLEVL.GT.3) THEN         ! ADDED BY ZHONGJIE HE
        IC(N)=2
        NUMGRID(N)=2*NUMGRID(N)-1
!      ENDIF
    ENDIF
  ENDDO

  ! PRESSURE COORDINATE LEVELS:
  I = (MAXGRID(PRESSURE)-1)/(NUMGRID(PRESSURE)-1)
  DO K=1,NUMGRID(PRESSURE)
    PPM(K) = PP0((K-1)*I+1)
  ENDDO

  CALL GRDMEMALC
! GET THE GRID SPACING INFORMATION FOR NEW GRD ARRAY
  DO N=1,NUMDIMS
    IF(IC(N).EQ.2)GRDSPAC(N)=GRDSPAC(N)/2.0
  ENDDO

! NOTE FINE GRID DOUBLES RESOLUTION FROM COASE GRID:

! 1. PROJECT COASE GRID ONTO FINE GRID:
  DO T=1,NUMGRID(4),IC(4)
    T0 = (T-1)/IC(4)+1
  DO K=1,NUMGRID(3),IC(3)
    K0 = (K-1)/IC(3)+1
  DO J=1,NUMGRID(2),IC(2)
    J0 = (J-1)/IC(2)+1
  DO I=1,NUMGRID(1),IC(1)
    I0 = (I-1)/IC(1)+1

    GRDANALS(I,J,K,T,1:NUMSTAT) = TMPANALS(I0,J0,K0,T0,1:NUMSTAT)
  ENDDO
  ENDDO
  ENDDO
  ENDDO

  ! RELEASE TEMPORARY MEMORY:
  CALL TMPMEMRLS

! 2. X DIRECTION:
  IF (IC(1) .EQ. 2) THEN
    DO T=1,NUMGRID(4),IC(4)
    DO K=1,NUMGRID(3),IC(3)
    DO J=1,NUMGRID(2),IC(2)
    DO I=2,NUMGRID(1),IC(1)

      GRDANALS(I,J,K,T,1:NUMSTAT) = 0.5*(GRDANALS(I-1,J,K,T,1:NUMSTAT) &
                                        +GRDANALS(I+1,J,K,T,1:NUMSTAT))
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF

! 3. Y DIRECTION:
  IF (IC(2) .EQ. 2) THEN
    DO T=1,NUMGRID(4),IC(4)
    DO K=1,NUMGRID(3),IC(3)
    DO J=2,NUMGRID(2),IC(2)
    DO I=1,NUMGRID(1)

      GRDANALS(I,J,K,T,1:NUMSTAT) = 0.5*(GRDANALS(I,J-1,K,T,1:NUMSTAT) &
                                        +GRDANALS(I,J+1,K,T,1:NUMSTAT))
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF

! 4. Z DIRECTION:
  IF (IC(3) .EQ. 2) THEN
    DO T=1,NUMGRID(4),IC(4)
    DO K=2,NUMGRID(3),IC(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)

      GRDANALS(I,J,K,T,1:NUMSTAT) = 0.5*(GRDANALS(I,J,K-1,T,1:NUMSTAT) &
                                        +GRDANALS(I,J,K+1,T,1:NUMSTAT))
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF

! 5. T DIRECTION:
  IF (IC(4) .EQ. 2) THEN
    DO T=2,NUMGRID(4),IC(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)

      GRDANALS(I,J,K,T,1:NUMSTAT) = 0.5*(GRDANALS(I,J,K,T-1,1:NUMSTAT) &
                                        +GRDANALS(I,J,K,T+1,1:NUMSTAT))
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF

  ! MAKE SURE INTERPOLATED SH ANALYSIS POSITIVE:
  DO T=1,NUMGRID(4)
  DO K=1,NUMGRID(3)
  DO J=1,NUMGRID(2)
  DO I=1,NUMGRID(1)
    GRDANALS(I,J,K,T,HUMIDITY) = &
      MAX(-GRDBKGND(I,J,K,T,HUMIDITY),GRDANALS(I,J,K,T,HUMIDITY))
  ENDDO
  ENDDO
  ENDDO
  ENDDO

  RETURN
END SUBROUTINE COAS2FINE_XIE

SUBROUTINE COAS2FINE
!*************************************************
! INTERPOLATION FROM COARSE GRID TO FINE GRID
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,S,N,I0,J0,K0,T0,SM
  INTEGER  :: IC(MAXDIMS)
! --------------------
  DO N=1,MAXDIMS
    NTMPGRD(N)=NUMGRID(N)
  ENDDO
  CALL TMPMEMALC
  CALL GRDMEMRLS
  DO N=1,MAXDIMS
    IC(N)=1
    IF(NTMPGRD(N).LT.MAXGRID(N)) THEN
!      IF(N.NE.3 .OR. GRDLEVL.GT.3) THEN         ! ADDED BY ZHONGJIE HE
        IC(N)=2
        NUMGRID(N)=2*NUMGRID(N)-1
!      ENDIF
    ENDIF
  ENDDO
  CALL GRDMEMALC
! GET THE GRID SPACING INFORMATION FOR NEW GRD ARRAY
  DO N=1,NUMDIMS
    IF(IC(N).EQ.2)GRDSPAC(N)=GRDSPAC(N)/2.0
  ENDDO
! GET THE INFORMATION FOR NEW GRID ARRAY BY INTERPOLATION
!======================== MODIFIED BY ZHONGJIE HE
  DO T=1,NUMGRID(4),IC(4)
  DO K=1,NUMGRID(3),IC(3)
  DO J=1,NUMGRID(2),IC(2)
  DO I=1,NUMGRID(1),IC(1)
    I0=0.5*(I+1)
    J0=0.5*(J+1)
    K0=0.5*(K+1)
    T0=0.5*(T+1)
    IF(IC(1).EQ.1)I0=I
    IF(IC(2).EQ.1)J0=J
    IF(IC(3).EQ.1)K0=K
    IF(IC(4).EQ.1)T0=T

    DO S=1,NUMSTAT
      SM=0
      IF(IC(1).EQ.2) THEN
        IF(I0.EQ.1) THEN
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0+1,J0,K0,T0,S)
          SM=SM+1
        ELSEIF(I0.EQ.NTMPGRD(1)) THEN
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0-1,J0,K0,T0,S)
          SM=SM+1
        ELSE
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0-1,J0,K0,T0,S)+TMPANALS(I0+1,J0,K0,T0,S)
          SM=SM+2
        ENDIF
      ENDIF
      IF(IC(2).EQ.2) THEN
        IF(J0.EQ.1) THEN
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0,J0+1,K0,T0,S)
          SM=SM+1
        ELSEIF(J0.EQ.NTMPGRD(2)) THEN
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0,J0-1,K0,T0,S)
          SM=SM+1
        ELSE
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0,J0-1,K0,T0,S)+TMPANALS(I0,J0+1,K0,T0,S)
          SM=SM+2
        ENDIF
      ENDIF
      IF(IC(3).EQ.2) THEN
        IF(K0.EQ.1) THEN
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0,J0,K0+1,T0,S)
          SM=SM+1
        ELSEIF(K0.EQ.NTMPGRD(3)) THEN
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0,J0,K0-1,T0,S)
          SM=SM+1
        ELSE
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0,J0,K0-1,T0,S)+TMPANALS(I0,J0,K0+1,T0,S)
          SM=SM+2
        ENDIF
      ENDIF
      IF(IC(4).EQ.2) THEN
        IF(T0.EQ.1) THEN
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0,J0,K0,T0+1,S)
          SM=SM+1
        ELSEIF(T0.EQ.NTMPGRD(4)) THEN
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0,J0,K0,T0-1,S)
          SM=SM+1
        ELSE
          GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0,J0,K0,T0-1,S)+TMPANALS(I0,J0,K0,T0+1,S)
          SM=SM+2
        ENDIF
      ENDIF

      IF(SM.GE.1) THEN
        GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)/FLOAT(SM)
      ELSE
        GRDANALS(I,J,K,T,S)=TMPANALS(I0,J0,K0,T0,S)
      ENDIF
    ENDDO
  ENDDO
  ENDDO
  ENDDO
  ENDDO
!========================================================== END OF MODIFICATION BY ZHONGJIE HE
! X DIRECTION
  IF(IC(1).EQ.2)THEN
    DO T=1,NUMGRID(4)  ,IC(4)
    DO K=1,NUMGRID(3)  ,IC(3)
    DO J=1,NUMGRID(2)  ,IC(2)
    DO I=2,NUMGRID(1)-1,IC(1)
      DO S=1,NUMSTAT
        GRDANALS(I,J,K,T,S)=0.5*(GRDANALS(I-1,J,K,T,S)+GRDANALS(I+1,J,K,T,S))
      ENDDO
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
! Y DIRECTION
  IF(IC(2).EQ.2)THEN
    DO T=1,NUMGRID(4)  ,IC(4)
    DO K=1,NUMGRID(3)  ,IC(3)
    DO J=2,NUMGRID(2)-1,IC(2)
    DO I=1,NUMGRID(1)
      DO S=1,NUMSTAT
        GRDANALS(I,J,K,T,S)=0.5*(GRDANALS(I,J-1,K,T,S)+GRDANALS(I,J+1,K,T,S))
      ENDDO
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
! Z DIRECTION
  IF(IC(3).EQ.2)THEN
    DO T=1,NUMGRID(4)  ,IC(4)
    DO K=2,NUMGRID(3)-1,IC(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      DO S=1,NUMSTAT
        GRDANALS(I,J,K,T,S)=0.5*(GRDANALS(I,J,K-1,T,S)+GRDANALS(I,J,K+1,T,S))
      ENDDO
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
! T DIRECTION
  IF(IC(4).EQ.2)THEN
    DO T=2,NUMGRID(4)-1,IC(4)
    DO K=1,NUMGRID(3)
    DO J=1,NUMGRID(2)
    DO I=1,NUMGRID(1)
      DO S=1,NUMSTAT
        GRDANALS(I,J,K,T,S)=0.5*(GRDANALS(I,J,K,T-1,S)+GRDANALS(I,J,K,T+1,S))
      ENDDO
    ENDDO
    ENDDO
    ENDDO
    ENDDO
  ENDIF
  CALL TMPMEMRLS
  RETURN
END SUBROUTINE COAS2FINE

SUBROUTINE FINE2COAS_XIE
!*************************************************
! INTERPOLATION FROM COARSE GRID TO FINE GRID
! HISTORY: AUGUST 2007, CODED by WEI LI.
!		MODIFIED BY YUANFU FROM WEI LI'S
!		COAS2FINE,WHICH IS INCORRECTLY
!		INTERPOLATING COASE TO FINE GRID.
!		THUS YUANFU REWRITES THE FINE2COAS.
!
!	NOTE: THIS ROUTINE ASSUME FINE GRID HALVES
!		THE COASE RESOLUTION.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,S,N,I0,J0,K0,T0,SM
  INTEGER  :: IC(MAXDIMS)
! --------------------
  DO N=1,MAXDIMS
    NTMPGRD(N)=NUMGRID(N)
  ENDDO
  CALL TMPMEMALC		! TEMPORARY MEMORY AND COPY GRDANALS TO HERE
  CALL GRDMEMRLS
  DO N=1,MAXDIMS
    IC(N)=1
    IF(NTMPGRD(N).GT.INIGRID(N)) THEN
        IC(N)=2
        NUMGRID(N)=(NUMGRID(N)-1)/2+1			! ASSUME NUMGRID ODD NUMBER
    ENDIF
  ENDDO
  CALL GRDMEMALC
! GET THE GRID SPACING INFORMATION FOR NEW GRD ARRAY
  DO N=1,NUMDIMS
    IF(IC(N) .EQ. 2) GRDSPAC(N)=GRDSPAC(N)*2.0
  ENDDO

! NOTE FINE GRID DOUBLES RESOLUTION FROM COASE GRID:

! 1. PROJECT COASE GRID ONTO FINE GRID:
  DO T=1,NUMGRID(4)
    T0 = (T-1)*IC(4)+1
  DO K=1,NUMGRID(3)
    K0 = (K-1)*IC(3)+1
  DO J=1,NUMGRID(2)
    J0 = (J-1)*IC(2)+1
  DO I=1,NUMGRID(1)
    I0 = (I-1)*IC(1)+1

    GRDANALS(I,J,K,T,1:NUMSTAT) = TMPANALS(I0,J0,K0,T0,1:NUMSTAT)
  ENDDO
  ENDDO
  ENDDO
  ENDDO

  ! RELEASE TEMPORARY MEMORY:
  CALL TMPMEMRLS

  RETURN
END SUBROUTINE FINE2COAS_XIE

SUBROUTINE FINE2COAS
!*************************************************
! PROJECTION FROM FINE GRID TO COARSE GRID (NOT USED)
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,S,N,M,I0,J0,K0,T0,SM
  INTEGER  :: IC(MAXDIMS)
! --------------------
  DO N=1,MAXDIMS
    NTMPGRD(N)=NUMGRID(N)
  ENDDO
  CALL TMPMEMALC
  DO N=1,MAXDIMS
    IC(N)=1
    IF(NTMPGRD(N).GT.INIGRID(N))IC(N)=2
  ENDDO

!======================== MODIFIED BY ZHONGJIE HE
  DO T=1,NTMPGRD(4)
  DO K=1,NTMPGRD(3)
  DO J=1,NTMPGRD(2)
  DO I=1,NTMPGRD(1)
    DO S=1,NUMSTAT
      SM=0
      TMPANALS(I,J,K,T,S)=0
      IF(IC(1).EQ.2) THEN
        IF(I.EQ.1) THEN
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I+1,J,K,T,S)
          SM=SM+1
        ELSEIF(I.EQ.NTMPGRD(1)) THEN
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I-1,J,K,T,S)
          SM=SM+1
        ELSE
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I-1,J,K,T,S)+GRDANALS(I+1,J,K,T,S)
          SM=SM+2
        ENDIF
      ENDIF
      IF(IC(2).EQ.2) THEN
        IF(J.EQ.1) THEN
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I,J+1,K,T,S)
          SM=SM+1
        ELSEIF(J.EQ.NTMPGRD(2)) THEN
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I,J-1,K,T,S)
          SM=SM+1
        ELSE
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I,J-1,K,T,S)+GRDANALS(I,J+1,K,T,S)
          SM=SM+2
        ENDIF
      ENDIF
      IF(IC(3).EQ.2) THEN
        IF(K.EQ.1) THEN
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I,J,K+1,T,S)
          SM=SM+1
        ELSEIF(K.EQ.NTMPGRD(3)) THEN
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I,J,K-1,T,S)
          SM=SM+1
        ELSE
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I,J,K-1,T,S)+GRDANALS(I,J,K+1,T,S)
          SM=SM+2
        ENDIF
      ENDIF
      IF(IC(4).EQ.2) THEN
        IF(T.EQ.1) THEN
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I,J,K,T+1,S)
          SM=SM+1
        ELSEIF(T.EQ.NTMPGRD(4)) THEN
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I,J,K,T-1,S)
          SM=SM+1
        ELSE
          TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)+GRDANALS(I,J,K,T-1,S)+GRDANALS(I,J,K,T+1,S)
          SM=SM+2
        ENDIF
      ENDIF

      IF(SM.GE.1) THEN
        TMPANALS(I,J,K,T,S)=TMPANALS(I,J,K,T,S)/FLOAT(SM)
      ELSE
        TMPANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)
      ENDIF
    ENDDO
  ENDDO
  ENDDO
  ENDDO
  ENDDO
!========================================================== END OF MODIFICATION BY ZHONGJIE HE

  CALL GRDMEMRLS
  DO N=1,MAXDIMS
    IF(IC(N).EQ.2)NUMGRID(N)=0.5*(NUMGRID(N)+1)
  ENDDO
  CALL GRDMEMALC
  DO N=1,NUMDIMS
    IF(IC(N).EQ.2)GRDSPAC(N)=GRDSPAC(N)*2.0
  ENDDO
! ----
  IF(.FALSE.)THEN
! ----
  DO T=1,NTMPGRD(4),IC(4)
  DO K=1,NTMPGRD(3),IC(3)
  DO J=1,NTMPGRD(2),IC(2)
  DO I=1,NTMPGRD(1),IC(1)
    DO S=1,NUMSTAT
      GRDANALS(I,J,K,T,S)=0.0
      M=0
      DO T0=MAX0(T-1,1),MIN0(T+1,NTMPGRD(4))
      DO K0=MAX0(K-1,1),MIN0(K+1,NTMPGRD(3))
      DO J0=MAX0(J-1,1),MIN0(J+1,NTMPGRD(2))
      DO I0=MAX0(I-1,1),MIN0(I+1,NTMPGRD(1))
        IF(I0.EQ.I.AND.J0.EQ.J.AND.K0.EQ.K.AND.T0.EQ.T)CYCLE
        M=M+1
        GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)+TMPANALS(I0,J0,K0,T0,S)
      ENDDO
      ENDDO
      ENDDO
      ENDDO
      IF(M.GT.0)GRDANALS(I,J,K,T,S)=GRDANALS(I,J,K,T,S)/M
    ENDDO
  ENDDO
  ENDDO
  ENDDO
  ENDDO
! ----
  ELSE
! ----
  DO T=1,NTMPGRD(4),IC(4)
  DO K=1,NTMPGRD(3),IC(3)
  DO J=1,NTMPGRD(2),IC(2)
  DO I=1,NTMPGRD(1),IC(1)
    I0=0.5*(I+1)
    J0=0.5*(J+1)
    K0=0.5*(K+1)
    T0=0.5*(T+1)
    IF(IC(1).EQ.1)I0=I
    IF(IC(2).EQ.1)J0=J
    IF(IC(3).EQ.1)K0=K
    IF(IC(4).EQ.1)T0=T
    DO S=1,NUMSTAT
      GRDANALS(I0,J0,K0,T0,S)=TMPANALS(I,J,K,T,S)
    ENDDO
  ENDDO
  ENDDO
  ENDDO
  ENDDO
! ----
  ENDIF
! ----
  CALL TMPMEMRLS
  RETURN
END SUBROUTINE FINE2COAS

SUBROUTINE CHECK_F_G
!*************************************************
! CHECK WHETHER THE GRADIENT MATCH COST FUNCTION (AFFILIATE)
! HISTORY: AUGUST 2007, CODED by WEI LI.
!*************************************************
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T,S
  REAL     :: CONTROL(NUMGRID(1),NUMGRID(2),NUMGRID(3),NUMGRID(4),NUMSTAT)
  REAL     :: ED,F2,F1,GG
! --------------------
  ED=0.00001
  DO S=1,NUMSTAT
  DO T=1,NUMGRID(4)
  DO K=1,NUMGRID(3)
  DO J=1,NUMGRID(2)
  DO I=1,NUMGRID(1)
    CONTROL(I,J,K,T,S)=GRDANALS(I,J,K,T,S)
  ENDDO
  ENDDO
  ENDDO
  ENDDO
  ENDDO
  CALL WCOMPGERNL
  CALL COSTGRADT2
  DO S=1,NUMSTAT
  DO T=1,NUMGRID(4)
  DO K=1,NUMGRID(3)
  DO J=1,NUMGRID(2)
  DO I=1,NUMGRID(1)
    PRINT*,'       '
    PRINT*,'       '
    PRINT*,'       '
    PRINT*,'       '
    PRINT*,'       '
    GRDANALS(I,J,K,T,S)=CONTROL(I,J,K,T,S)+ED
    CALL COSTFUNCT2
    F2=COSTFUN
    GRDANALS(I,J,K,T,S)=CONTROL(I,J,K,T,S)-ED
    CALL COSTFUNCT2
    F1=COSTFUN
    GRDANALS(I,J,K,T,S)=CONTROL(I,J,K,T,S)
    GG=(F2-F1)/(2.0*ED)
!    PRINT*,F2,F1,GG
	IF(ABS(GRADINT(I,J,K,T,S)-GG).GT.1.0E-19) &
    PRINT*,I,J,K,T,S,GRADINT(I,J,K,T,S),GG,GRADINT(I,J,K,T,S)-GG
	IF(ABS(GRADINT(I,J,K,T,S)-GG).GT.1000.0)STOP
  ENDDO
  ENDDO
  ENDDO
  ENDDO
  ENDDO
  RETURN
END SUBROUTINE CHECK_F_G

SUBROUTINE GETW
  IMPLICIT NONE
! --------------------
  INTEGER  :: I,J,K,T
! --------------------
  CALL WCOMPGERNL      ! MODIFIED BY ZHONGJIE HE
  DO T=1,NUMGRID(4)
  DO K=1,NUMGRID(3)
  DO J=1,NUMGRID(2)
  DO I=1,NUMGRID(1)
    GRDANALS(I,J,K,T,3)=WWW(I,J,K,T) !*SCL(U_CMPNNT)
  ENDDO
  ENDDO
  ENDDO
  ENDDO
  RETURN
END SUBROUTINE GETW

END MODULE STMAS4D_CORE

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
MODULE gold_seaice_netcdf

  USE gold_seaice_lib
  IMPLICIT NONE

  PRIVATE
  PUBLIC :: ini_netcdf_sic
  PUBLIC :: write_netcdf_sic
  PUBLIC :: end_netcdf_sic
  PUBLIC :: nclon1, nclon2, nclon3, nclat1, nclat2, nclat3

  REAL, DIMENSION(:), ALLOCATABLE :: nclon1, nclon2, nclon3
  REAL, DIMENSION(:), ALLOCATABLE :: nclat1, nclat2, nclat3

  INTEGER, PARAMETER :: nall=100, nfiles=2
  INTEGER :: nco(nfiles), iddimo(nall, nfiles), idvaro(nall, nfiles)

CONTAINS

  ! Begin netCDF initialisation process
  SUBROUTINE ini_netcdf_sic(istep, imode)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: istep, imode

    REAL :: day, rtime
    INTEGER :: iyear, imonth

    day = istep * dtsic * tsc / 86400.0
    ! Small constant added to day calculation for round-off reasons
    ! (i.e. 365.2499 vs. 365.25)
    iyear = INT((day + 0.001) / yearlen)
    imonth = INT((day - iyear * yearlen) / (yearlen / 12.0)) + 1
    rtime = iyear + ((day - iyear * yearlen)) / yearlen
    IF (debug_loop) &
         & PRINT *, 'istep', istep, 'day', day, 'iyear', iyear, &
         & 'imonth', imonth, 'rtime', rtime

    CALL ini_netcdf_sic1(outdir_name, lenout, lout, imonth, rtime, &
         & nclon1, nclat1, nclon2, nclat2, nclon3, nclat3, maxi, maxj, imode)
  END SUBROUTINE ini_netcdf_sic

  ! Continue netCDF initialisation process
  SUBROUTINE ini_netcdf_sic1(dir_name, ilen, runid, imonth, rtime, &
       & alon1, alat1, alon2, alat2, alon3, alat3, mg, jgg, imode)
    USE writenc6
    IMPLICIT NONE

    CHARACTER(LEN=100), INTENT(IN) :: dir_name
    INTEGER, INTENT(IN) :: ilen
    CHARACTER(LEN=3), INTENT(IN) :: runid
    INTEGER, INTENT(IN) :: mg, jgg
    REAL, DIMENSION(mg), INTENT(IN) :: alon1, alon2, alon3
    REAL, DIMENSION(jgg), INTENT(IN) :: alat1, alat2, alat3
    REAL, INTENT(IN) :: rtime
    INTEGER, INTENT(IN) :: imonth, imode

    INTEGER, PARAMETER :: nmaxdims=4
    INTEGER :: ndim, nvar
    INTEGER, DIMENSION(nall) :: natts, nattsvar, vdims, ndims
    INTEGER, DIMENSION(nmaxdims,nall) :: vadims
    CHARACTER(LEN=200), DIMENSION(nall) :: dimname, varname
    CHARACTER(LEN=200), DIMENSION(2,nmaxdims,nall) :: attdimname, attvarname

    INTEGER :: itime, ifname1
    CHARACTER(LEN=200) :: fname1
    CHARACTER(LEN=10) :: cyear
    CHARACTER(LEN=2), DIMENSION(13) :: cmon = &
         & (/ '00','01','02','03','04','05','06', &
         &    '07','08','09','10','11','12' /)

    INTEGER, EXTERNAL :: lnsig

    itime = 1

    CALL setup_nc_sic(mg, jgg, itime, nmaxdims, ndim, nvar, natts, &
         & nattsvar, vdims, vadims, ndims, dimname, varname, &
         & attdimname, attvarname)
    IF (imode == 1) THEN
       fname1 = dir_name(1:ilen)//'gsic_'//runid(1:3)//'_rs_'// &
            & cyear//'_'//cmon(imonth)//'.nc'
    ELSE IF (imode == 2) THEN
       fname1 = dir_name(1:ilen)//'gsic_'//runid(1:3)//'_av_'// &
            & cyear//'_'//cmon(imonth)//'.nc'
    END IF
    ifname1 = lnsig(fname1)

    CALL ininc(fname1(1:ifname1), nmaxdims, ndim, nvar, natts, nattsvar, &
         & vdims, vadims, ndims, dimname, varname, attdimname, attvarname, &
         & nco(imode), iddimo(1,imode), idvaro(1,imode))

    ! Longitude coordinates (tracer, u-point, v-point)
    CALL writedim(nco(imode), iddimo(1,imode), alon1)
    CALL writedim(nco(imode), iddimo(2,imode), alon2)
    CALL writedim(nco(imode), iddimo(3,imode), alon3)

    ! Latitude coordinates (tracer, u-point, v-point)
    CALL writedim(nco(imode), iddimo(4,imode), alat1)
    CALL writedim(nco(imode), iddimo(5,imode), alat2)
    CALL writedim(nco(imode), iddimo(6,imode), alat3)

    ! Time
    CALL writedim(nco(imode), iddimo(7,imode), (/ REAL(rtime) /))
  END SUBROUTINE ini_netcdf_sic1

  ! Set up netCDF file's array names, units, descriptions, etc.
  SUBROUTINE setup_nc_sic(nlon, nlat, ntime, nmaxdims, &
       & ndim, nvar, natts, nattsvar, vdims, vadims, ndims, &
       & dimname, varname, attdimname, attvarname)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: nlon, nlat, ntime, nmaxdims
    INTEGER, INTENT(OUT) :: ndim, nvar
    INTEGER, DIMENSION(nall), INTENT(OUT) :: natts, nattsvar, vdims, ndims
    INTEGER, INTENT(OUT) :: vadims(nmaxdims,nall)
    CHARACTER(LEN=200), DIMENSION(nall), INTENT(OUT) :: dimname, varname
    CHARACTER(LEN=200), DIMENSION(2,nmaxdims,nall), INTENT(OUT) :: &
         & attdimname, attvarname
    INTEGER, EXTERNAL :: loc_dim

    ndim = 0
    nvar = 0

    ndim = ndim + 1
    dimname(ndim) = 'longitude'
    ndims(ndim) = nlon
    natts(ndim) = 2
    attdimname(1,1,ndim) = 'long_name'
    attdimname(2,1,ndim) = 'longitude'
    attdimname(1,2,ndim) = 'units'
    attdimname(2,2,ndim) = 'degrees_east'

    ndim = ndim + 1
    dimname(ndim) = 'longitude_1'
    ndims(ndim) = nlon
    natts(ndim) = 2
    attdimname(1,1,ndim) = 'long_name'
    attdimname(2,1,ndim) = 'longitude'
    attdimname(1,2,ndim) = 'units'
    attdimname(2,2,ndim) = 'degrees_east'

    ndim = ndim + 1
    dimname(ndim) = 'longitude_2'
    ndims(ndim) = nlon
    natts(ndim) = 2
    attdimname(1,1,ndim) = 'long_name'
    attdimname(2,1,ndim) = 'longitude'
    attdimname(1,2,ndim) = 'units'
    attdimname(2,2,ndim) = 'degrees_east'

    ndim = ndim + 1
    dimname(ndim) = 'latitude'
    ndims(ndim) = nlat
    natts(ndim) = 2
    attdimname(1,1,ndim) = 'long_name'
    attdimname(2,1,ndim) = 'latitude'
    attdimname(1,2,ndim) = 'units'
    attdimname(2,2,ndim) = 'degrees_north'

    ndim = ndim + 1
    dimname(ndim) = 'latitude_1'
    ndims(ndim) = nlat
    natts(ndim) = 2
    attdimname(1,1,ndim) = 'long_name'
    attdimname(2,1,ndim) = 'latitude'
    attdimname(1,2,ndim) = 'units'
    attdimname(2,2,ndim) = 'degrees_north'

    ndim = ndim + 1
    dimname(ndim) = 'latitude_2'
    ndims(ndim) = nlat
    natts(ndim) = 2
    attdimname(1,1,ndim) = 'long_name'
    attdimname(2,1,ndim) = 'latitude'
    attdimname(1,2,ndim) = 'units'
    attdimname(2,2,ndim) = 'degrees_north'

    ndim = ndim + 1
    dimname(ndim) = 'time'
    ndims(ndim) = ntime
    natts(ndim) = 2
    attdimname(1,1,ndim) = 'long_name'
    attdimname(2,1,ndim) = 'time'
    attdimname(1,2,ndim) = 'units'
    attdimname(2,2,ndim) = 'years'

    nvar = nvar + 1
    varname(nvar) = 'sic_height'
    vdims(nvar) = 3
    vadims(1,nvar) = loc_dim('longitude',dimname,nall)
    vadims(2,nvar) = loc_dim('latitude',dimname,nall)
    vadims(3,nvar) = loc_dim('time',dimname,nall)
    nattsvar(nvar) = 2
    attvarname(1,1,nvar) = 'long_name'
    attvarname(2,1,nvar) = 'Sea-ice height'
    attvarname(1,2,nvar) = 'units'
    attvarname(2,2,nvar) = 'm'

    nvar = nvar + 1
    varname(nvar) = 'sic_cover'
    vdims(nvar) = 3
    vadims(1,nvar) = loc_dim('longitude',dimname,nall)
    vadims(2,nvar) = loc_dim('latitude',dimname,nall)
    vadims(3,nvar) = loc_dim('time',dimname,nall)
    nattsvar(nvar) = 2
    attvarname(1,1,nvar) = 'long_name'
    attvarname(2,1,nvar) = 'Sea-ice fractional cover'
    attvarname(1,2,nvar) = 'units'
    attvarname(2,2,nvar) = 'dimensionless'

    nvar = nvar + 1
    varname(nvar) = 'sic_temp'
    vdims(nvar) = 3
    vadims(1,nvar) = loc_dim('longitude',dimname,nall)
    vadims(2,nvar) = loc_dim('latitude',dimname,nall)
    vadims(3,nvar) = loc_dim('time',dimname,nall)
    nattsvar(nvar) = 2
    attvarname(1,1,nvar) = 'long_name'
    attvarname(2,1,nvar) = 'Sea-ice temperature'
    attvarname(1,2,nvar) = 'units'
    attvarname(2,2,nvar) = 'C'

    nvar = nvar + 1
    varname(nvar) = 'sic_albedo'
    vdims(nvar) = 3
    vadims(1,nvar) = loc_dim('longitude',dimname,nall)
    vadims(2,nvar) = loc_dim('latitude',dimname,nall)
    vadims(3,nvar) = loc_dim('time',dimname,nall)
    nattsvar(nvar) = 2
    attvarname(1,1,nvar) = 'long_name'
    attvarname(2,1,nvar) = 'Sea-ice albedo'
    attvarname(1,2,nvar) = 'units'
    attvarname(2,2,nvar) = 'dimensionless'

    nvar = nvar + 1
    varname(nvar) = 'delta_height'
    vdims(nvar) = 3
    vadims(1,nvar) = loc_dim('longitude',dimname,nall)
    vadims(2,nvar) = loc_dim('latitude',dimname,nall)
    vadims(3,nvar) = loc_dim('time',dimname,nall)
    nattsvar(nvar) = 2
    attvarname(1,1,nvar) = 'long_name'
    attvarname(2,1,nvar) = 'Change in sea-ice height'
    attvarname(1,2,nvar) = 'units'
    attvarname(2,2,nvar) = 'm/s'

    nvar = nvar + 1
    varname(nvar) = 'delta_cover'
    vdims(nvar) = 3
    vadims(1,nvar) = loc_dim('longitude',dimname,nall)
    vadims(2,nvar) = loc_dim('latitude',dimname,nall)
    vadims(3,nvar) = loc_dim('time',dimname,nall)
    nattsvar(nvar) = 2
    attvarname(1,1,nvar) = 'long_name'
    attvarname(2,1,nvar) = 'Change in sea-ice fractional cover'
    attvarname(1,2,nvar) = 'units'
    attvarname(2,2,nvar) = 'dimensionless'

    nvar = nvar + 1
    varname(nvar) = 'sic_heat'
    vdims(nvar) = 3
    vadims(1,nvar) = loc_dim('longitude',dimname,nall)
    vadims(2,nvar) = loc_dim('latitude',dimname,nall)
    vadims(3,nvar) = loc_dim('time',dimname,nall)
    nattsvar(nvar) = 2
    attvarname(1,1,nvar) = 'long_name'
    attvarname(2,1,nvar) = 'Sea-ice heat flux'
    attvarname(1,2,nvar) = 'units'
    attvarname(2,2,nvar) = 'W/m2'

    nvar = nvar + 1
    varname(nvar) = 'sic_fw'
    vdims(nvar) = 3
    vadims(1,nvar) = loc_dim('longitude',dimname,nall)
    vadims(2,nvar) = loc_dim('latitude',dimname,nall)
    vadims(3,nvar) = loc_dim('time',dimname,nall)
    nattsvar(nvar) = 2
    attvarname(1,1,nvar) = 'long_name'
    attvarname(2,1,nvar) = 'Sea-ice freshwater flux'
    attvarname(1,2,nvar) = 'units'
    attvarname(2,2,nvar) = 'mm/s'
  END SUBROUTINE setup_nc_sic

  ! Writes data to netCDF file
  SUBROUTINE write_netcdf_sic(k1, varice, tice, albice, dtha, &
       & fx0_sic, fw_sic, work, maxi, maxj, imode)
    USE writenc6
    IMPLICIT NONE

    integer maxi,maxj,imode,k1(0:maxi+1,0:maxj+1)
    real varice(2,maxi,maxj),tice(maxi,maxj),albice(maxi,maxj), &
         & dtha(2,maxi,maxj),fx0_sic(maxi,maxj),fw_sic(maxi,maxj)
    real work(maxi, maxj)

    ! Sea-ice height (i.e. final argument = 1)
    call prep_netcdf_sic(varice(1,:,:), work, k1)
    call writevar(nco(imode), idvaro(1,imode), work)

    ! Sea-ice area (i.e. final argument = 2)
    call prep_netcdf_sic(varice(2,:,:), work, k1)
    call writevar(nco(imode), idvaro(2,imode), work)

    ! Sea-ice temperature
    call prep_netcdf_sic(tice, work, k1)
    call writevar(nco(imode), idvaro(3,imode), work)

    ! Sea-ice albedo
    call prep_netcdf_sic(albice, work, k1)
    call writevar(nco(imode), idvaro(4,imode), work)

    ! Change in sea-ice height
    call prep_netcdf_sic(dtha(1,:,:), work, k1)
    call writevar(nco(imode), idvaro(5,imode), work)

    ! Change in sea-ice fractional cover
    call prep_netcdf_sic(dtha(2,:,:), work, k1)
    call writevar(nco(imode), idvaro(6,imode), work)

    ! Sea-ice heat flux
    call prep_netcdf_sic(fx0_sic, work, k1)
    call writevar(nco(imode), idvaro(7,imode), work)

    ! Sea-ice freshwater flux
    call prep_netcdf_sic(fw_sic, work, k1)
    call writevar(nco(imode), idvaro(8,imode), work)
  END SUBROUTINE write_netcdf_sic

  ! Organises a two-dimensional array that's on the tracer grid
  SUBROUTINE prep_netcdf_sic(temper, temp1, iland)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: iland(0:maxi+1,0:maxj+1)
    REAL,    INTENT(IN) :: temper(maxi,maxj)
    REAL,    INTENT(OUT) :: temp1(maxi,maxj)

    WHERE (iland(1:maxi,1:maxj) >= 90)
       temp1 = -99999.0
    ELSEWHERE
       temp1 = temper
    END WHERE
  END SUBROUTINE prep_netcdf_sic

  ! Ends netCDF-writing process and closes netCDF file
  SUBROUTINE end_netcdf_sic(imode)
    USE writenc6
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: imode

    CALL closenc(nco(imode))
  END SUBROUTINE END_NETCDF_SIC

END MODULE gold_seaice_netcdf

max_stars_repo_name: derpycode/cgenie

max_issues_repo_name: derpycode/cgenie

max_forks_repo_name: derpycode/cgenie


text:
*DECK CPBDI
      SUBROUTINE CPBDI (ABD, LDA, N, M, DET)
C***BEGIN PROLOGUE  CPBDI
C***PURPOSE  Compute the determinant of a complex Hermitian positive
C            definite band matrix using the factors computed by CPBCO or
C            CPBFA.
C***LIBRARY   SLATEC (LINPACK)
C***CATEGORY  D3D2
C***TYPE      COMPLEX (SPBDI-S, DPBDI-D, CPBDI-C)
C***KEYWORDS  BANDED, DETERMINANT, INVERSE, LINEAR ALGEBRA, LINPACK,
C             MATRIX, POSITIVE DEFINITE
C***AUTHOR  Moler, C. B., (U. of New Mexico)
C***DESCRIPTION
C
C     CPBDI computes the determinant
C     of a complex Hermitian positive definite band matrix
C     using the factors computed by CPBCO or CPBFA.
C     If the inverse is needed, use CPBSL  N  times.
C
C     On Entry
C
C        ABD     COMPLEX(LDA, N)
C                the output from CPBCO or CPBFA.
C
C        LDA     INTEGER
C                the leading dimension of the array  ABD .
C
C        N       INTEGER
C                the order of the matrix  A .
C
C        M       INTEGER
C                the number of diagonals above the main diagonal.
C
C     On Return
C
C        DET     REAL(2)
C                determinant of original matrix in the form
C                determinant = DET(1) * 10.0**DET(2)
C                with  1.0 .LE. DET(1) .LT. 10.0
C                or  DET(1) .EQ. 0.0 .
C
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   780814  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  CPBDI
      INTEGER LDA,N,M
      COMPLEX ABD(LDA,*)
      REAL DET(2)
C
      REAL S
      INTEGER I
C***FIRST EXECUTABLE STATEMENT  CPBDI
C
C     COMPUTE DETERMINANT
C
      DET(1) = 1.0E0
      DET(2) = 0.0E0
      S = 10.0E0
      DO 50 I = 1, N
         DET(1) = REAL(ABD(M+1,I))**2*DET(1)
         IF (DET(1) .EQ. 0.0E0) GO TO 60
   10    IF (DET(1) .GE. 1.0E0) GO TO 20
            DET(1) = S*DET(1)
            DET(2) = DET(2) - 1.0E0
         GO TO 10
   20    CONTINUE
   30    IF (DET(1) .LT. S) GO TO 40
            DET(1) = DET(1)/S
            DET(2) = DET(2) + 1.0E0
         GO TO 30
   40    CONTINUE
   50 CONTINUE
   60 CONTINUE
      RETURN
      END

ext: f


text:
Miller Hall is one of the three buildings located in Segundo North. It houses the approximately 120 members of the Integrated Studies program.  The building was first opened for use in Fall 2004. It has a separate meeting area on the ground floor, which has been used for summer advising Freshman Seminars.

During the 20062007 school year, Miller Hall instead housed 171 Integrated Studies students, because of increased interest in the program and a larger freshman class. This meant that all of the double occupancy rooms became triple occupancy.

 Photos

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C$Procedure      CUTSTR ( Cut a long string into substrings )

      SUBROUTINE CUTSTR ( STRING, START, WIDTH, BREAKS, BEG, END )
      IMPLICIT NONE

C$ Abstract
C
C     Cut a long string into substrings, breaking at "good" points
C     whenever possible.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     STRING
C     UTILITY
C
C$ Declarations

      CHARACTER*(*)         STRING
      INTEGER               START
      INTEGER               WIDTH
      CHARACTER*(*)         BREAKS
      INTEGER               BEG
      INTEGER               END

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     STRING     I   Long string.
C     START      I   Nominal beginning of substring.
C     WIDTH      I   Maximum width of substrings.
C     BREAKS     I   Characters indicating good break points.
C     BEG        I   Beginning of substring.
C     END        O   End of substring.
C
C$ Detailed_Input
C
C     STRING      is an arbitrary character string. Typically, this
C                 is too wide to fit into an area of limited width:
C                 an element of a character array, for instance, or
C                 an area on a terminal screen.
C
C     START       is the nominal beginning of the next substring.
C                 (STRING(START:START) is the first character that
C                 can appear in the substring.) It is used to skip
C                 past substrings returned by previous calls.
C
C     WIDTH       is the width (in characters) of the limited area.
C                 Thus, it is the maximum width of the substrings
C                 to be returned.
C
C     BREAKS      is a collection of characters indicating preferred
C                 places to break the string into substrings: commas,
C                 colons, and periods, for instance. BREAKS is always
C                 treated as though it contains a space, whether it
C                 does or not. (That is, '+-=' is treated as ' +-='.)
C
C$ Detailed_Output
C
C     BEG,
C     END         are the endpoints of a substring no wider than
C                 WIDTH. Substrings always begin and end with non-blank
C                 characters. BEG is zero if no non-blank substring
C                 was found.
C
C$ Exceptions.
C
C     1) If STRING(START:) is blank or BEG is greater than the declared
C        length of STRING, both BEG and END are zero.
C
C     2) If START is less than one, the error 'SPICE(BEFOREBEGSTR)'
C        is signalled.
C
C     4) If WIDTH is less than one, the error 'SPICE(WIDTHTOOSMALL)'
C        is signalled.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C      This routine is useful primarily for displaying messages on
C      terminal screens or within log files. Since messages can run
C      to hundreds of characters, while most output devices cannot
C      handle more than 80 or 132 characters at a time, it is necessary
C      to break the strings. The friendliness of a message is enhanced
C      if these breaks occur at "natural" places within the message,
C      rather than at rigid intervals.
C
C      The most natural breaks occur before spaces. Slightly less
C      natural breaks occur at the characters
C
C         Comma             ,
C         Period            .
C         Semicolon         ;
C         Colon             :
C         Question          ?
C         Exclamation       !
C         End parenthesis   )
C         End bracket       ]
C         End brace         }
C         End angle         >
C
C      or before the characters
C
C         Begin parenthesis   (
C         Begin bracket       [
C         Begin brace         {
C         Begin angle         <
C
C      At any rate, breaks should occur between adjacent letters or
C      numeric characters only as a last resort.
C
C      In the absence of other instructions, CUTSTR tries to break:
C
C         1) before      space
C
C         2) at           , . ; : - ) ] } >
C            or before              ( [ { <
C
C         3) at           ' "                   (even occurrence)
C            or before    ' "                   (odd  occurrence)
C
C         4) at           ? ! = _ %
C            or before    \ ~ $ @ ^ * / | & +
C
C      before forcing a break at an aribitrary location.
C
C      You may override these rules by supplying a set of preferred
C      characters in BREAKS. Before applying the rules shown above,
C      CUTSTR will try to break AT these characters. (However, breaks
C      always occur BEFORE spaces.)
C
C$ Examples
C
C      CUTSTR might typically be used to display a long diagnostic
C      messages on a terminal screen. For example, suppose that the
C      following message has been returned by a subroutine.
C
C         'I believe you have made an significant error by requesting
C          that I send to the printer a file containing 250 megabytes
C          of text information. The system manager is likely to be
C          very unhappy with such a request. I suggest you reconsider
C          your intended action.'
C
C      and that this needs to be displayed on a 40-character monitor.
C      The following code fragment
C
C         WIDTH  = 40
C         BREAKS = ' ,.'
C
C         CALL CUTSTR ( MSSG, 1, WIDTH, BREAKS, BEG, END )
C
C         DO WHILE ( BEG .NE. 0 )
C            WRITE (6,*) MSSG (BEG:END)
C
C            START = END + 1
C            CALL CUTSTR ( MSSG, START, WIDTH, BREAKS, BEG, END )
C         END DO
C
C      would display something like
C
C         I believe you have made an significant
C         error by requesting that I send to the
C         printer a file containing 250 megabytes
C         of text information. The system manager
C         is likely to be very unhappy with such a
C         request.  I suggest you reconsider your
C         intended action.
C
C      On a more whimsical note, you could indent each successive lines
C      by three characters: the code fragment
C
C         WIDTH  = 40
C         BREAKS = ' ,.'
C         INDENT = 1
C
C         CALL CUTSTR ( MSSG, 1, WIDTH, BREAKS, BEG, END )
C
C         DO WHILE ( BEG .NE. 0 )
C            TEMP           = ' '
C            TEMP(INDENT: ) = MSSG(BEG:END)
C            WRITE (6,*) TEMP
C
C            INDENT = INDENT + 3
C            WIDTH  = MAX ( WIDTH-3, 9 )
C
C            START = END + 1
C            CALL CUTSTR ( MSSG, START, WIDTH, BREAKS, BEG, END )
C         END DO
C
C      would display something like
C
C         I believe you have made an significant
C            error by requesting that I send to
C               the printer a file containing 250
C                  megabytes of text information.
C                     The system manager is likely
C                        to be very unhappy with
C                           such a request. I
C                              suggest you
C                                 reconsider your
C                                    intended
C                                       action.
C
C     Note that both loops terminate when BEG is zero. This indicates
C     that no substring was found (and that none will be returned by
C     subsequent calls). If the string is full, the loop will terminate
C     normally when START becomes greater than the length of the string.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     W.L. Taber     (JPL)
C     I.M. Underwood (JPL)
C
C$ Version
C
C-     META/2 Configured Version 2.0.0, 9-MAY-1994 (WLT)
C
C         This is the configured version of the Command Loop
C         software as of May 9, 1994
C
C
C-     META/2 Configured Version 1.0.0, 3-MAY-1994 (WLT)
C
C         This is the configured version of META/2
C         software as of May 3, 1994
C
C
C     Version B1.0.0, 29-APR-1988 (WLT) (IMU)
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               EVEN
      INTEGER               OCCURS
      LOGICAL               RETURN

C
C     Local variables
C

      INTEGER               A
      INTEGER               B
      INTEGER               BLANK
      INTEGER               P

      INTEGER               LENGTH
      INTEGER               LONG

      CHARACTER*15          PUNCT
      CHARACTER*15          QUOTE
      CHARACTER*15          OTHER
      CHARACTER*15          DEF
      CHARACTER*15          QTYPE
      CHARACTER*15          PTYPE
      CHARACTER*15          OTYPE
      CHARACTER*15          DTYPE

      INTEGER               PASS

      CHARACTER*1           THIS
      CHARACTER*1           NEXT
      INTEGER               HERE
      INTEGER               THERE
C
C     Parameters used to simulate an enumerated type for
C     the various passes required to break the string at
C     good places.  Note that the order is important.
C     This forces the routine to try spaces first, user
C     supplied preferences next, etc.    It is also
C     critical that these be defined to be a sequence
C     of consecutive integers.
C
      INTEGER               SPACE
      PARAMETER           ( SPACE  = 1 )

      INTEGER               PREF
      PARAMETER           ( PREF   = SPACE  + 1 )

      INTEGER               PUNC
      PARAMETER           ( PUNC   = PREF   + 1 )

      INTEGER               QUOT
      PARAMETER           ( QUOT   = PUNC   + 1 )

      INTEGER               OTHR
      PARAMETER           ( OTHR   = QUOT   + 1 )

      INTEGER               FAIL
      PARAMETER           ( FAIL   = OTHR   + 1 )

      INTEGER               FINSH
      PARAMETER           ( FINSH  = FAIL   + 1 )
C
C     The ASCII character value for the backslash is needed for
C     uniformity of porting this routine (Some platforms treat the
C     backslah as a special character and so we can't just use
C     the character in strings.)
C
      INTEGER               BSLASH
      PARAMETER           ( BSLASH = 92 )

C
C     Standard error handling.
C
      IF ( RETURN() ) THEN
         RETURN
      END IF

C
C     Exceptions first. Is START outside the bounds of the string?
C
      LENGTH = LEN ( STRING )

      IF ( START .GT. LENGTH ) THEN
         BEG = 0
         END = 0
         RETURN

      ELSE IF ( START .LT. 1 ) THEN
         CALL CHKIN  ( 'CUTSTR' )
         CALL SIGERR ( 'SPICE(BEFOREBEGSTR)' )
         CALL CHKOUT ( 'CUTSTR' )
         RETURN
      END IF

C
C     Is the width reasonable?
C
      IF ( WIDTH .LT. 1 ) THEN
         CALL CHKIN  ( 'CUTSTR' )
         CALL SIGERR ( 'SPICE(WIDTHTOOSMALL)' )
         CALL CHKOUT ( 'CUTSTR' )
         RETURN
      END IF

C
C     Does the remainder of the string contain anything besides blanks?
C
      IF ( STRING(START: ) .EQ. ' ' ) THEN
         BEG   = 0
         END   = 0
         RETURN
      END IF

C
C     Obviously, we should try to get the longest possible substring.
C
C
      BEG   = START
      BLANK = ICHAR( ' ' )

      DO WHILE ( ICHAR(STRING(BEG:BEG)) .EQ. BLANK )
         BEG = BEG + 1
      END DO

      LONG = BEG + WIDTH - 1

C
C     The remainder of the substring may fit without a trim.
C     But drop trailing blanks anyway.
C
      IF ( LENGTH .LE. LONG ) THEN

         END = LENGTH

         DO WHILE ( ICHAR(STRING(END:END)) .EQ. BLANK )
            END = END - 1
         END DO

         RETURN
      END IF

C
C     Assign the default break characters. Each character in PUNCT,
C     QUOTE, or OTHER indicates a good place to break. The associated
C     type indicates whether the break should occur at or before the
C     the character:
C
C        Type         Break occurs
C        ----         ------------------------------------------------
C         A           At the character.
C         B           Before the character.
C         P           At an EVEN occurrence, or
C                     Before an ODD occurrence.
C
C
      PUNCT = ',.;:-)]}>([{<'
      PTYPE = 'AAAAAAAAABBBB'

      QUOTE = '"'''
      QTYPE = 'PP'

      OTHER = '?!~$@^=_%*/|&+' // CHAR(BSLASH)
      OTYPE = 'AABBBBAAABBBBB' // 'B'

C
C     We will do this in five passes. During the first pass, we will
C     try to break before a space. During the second pass, we will try
C     to break at one of the preferred characters. During the third,
C     fourth, and fifth passes, we will try to break at or before one
C     of the quotation, punctuation, or other default characters.
C
      PASS      = SPACE
      A         = ICHAR( 'A' )
      P         = ICHAR( 'P' )
      B         = ICHAR( 'B' )

      DO WHILE ( PASS .NE. FAIL )

         END = LONG

         DO WHILE ( END .GE. BEG )

            THIS = STRING(END   : END  )
            NEXT = STRING(END+1 : END+1)

C
C           Always break BEFORE a space.
C
            IF ( PASS .EQ. SPACE ) THEN

              IF ( ICHAR(NEXT) .EQ. BLANK ) THEN
                PASS = FINSH
              END IF

C
C           Always break AT a preferred character.
C
            ELSE IF ( PASS .EQ. PREF ) THEN

              IF ( INDEX ( BREAKS, THIS ) .GT. 0 ) THEN
                PASS = FINSH
              END IF

C
C           But with default characters, some break at, some
C           before, and some depend on the parity of strangers.
C
            ELSE

               HERE  = INDEX ( DEF, THIS )
               THERE = INDEX ( DEF, NEXT )

               IF ( HERE .GT. 0 ) THEN

                 IF ( ICHAR(DTYPE(HERE:HERE)) .EQ. A ) THEN
                   PASS = FINSH

                 ELSE IF ( ICHAR(DTYPE(HERE:HERE)) .EQ. P ) THEN
                   IF ( EVEN ( OCCURS ( STRING( :END), THIS ) ) ) THEN
                     PASS = FINSH
                   END IF
                 END IF
               END IF

               IF (      THERE .GT. 0
     .             .AND. PASS  .NE. FINSH ) THEN

                 IF ( ICHAR(DTYPE(THERE:THERE)) .EQ. B ) THEN
                   PASS = FINSH

                 ELSE IF ( ICHAR(DTYPE(THERE:THERE)) .EQ. P ) THEN
                   IF ( EVEN ( OCCURS ( STRING( :END), NEXT ) ) ) THEN
                     PASS = FINSH
                   END IF
                 END IF

               END IF

            END IF

C
C           If we've found a break point, remove any trailing blanks
C           before returning.
C
            IF ( PASS .EQ. FINSH ) THEN
               DO WHILE ( ICHAR(STRING(END:END)) .EQ. BLANK )
                  END = END - 1
               END DO

               RETURN
            ELSE
               END = END - 1
            END IF

         END DO

C
C        We may have to try another pass.
C
         PASS  = PASS + 1

C
C        In the final passes, only the character set changes.
C
         IF ( PASS .EQ. PUNC ) THEN
            DEF   = PUNCT
            DTYPE = PTYPE

         ELSE IF ( PASS .EQ. QUOT ) THEN
            DEF   = QUOTE
            DTYPE = QTYPE

         ELSE IF ( PASS .EQ. OTHR ) THEN
            DEF   = OTHER
            DTYPE = OTYPE
         END IF

      END DO

C
C     Looks like we'll have to do this the hard way.
C
      END = LONG

      RETURN

      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
program ewald_conv4

! Runs a convergence study on the test_ewald4 problem, the example 4.1 in [1].
!
! [1] Pask, J. E., Sukumar, N., & Mousavi, S. E. (2012). Linear scaling
! solution of the all-electron Coulomb problem in solids. International Journal
! for Multiscale Computational Engineering, 10(1), 83–99.
! doi:10.1615/IntJMultCompEng.2011002201

use types, only: dp
use constants, only: ang2bohr, kJmol2Ha
use ewald_sums, only: ewald_box, ewald_fft1, fred2fcart
use utils, only: assert, init_random
implicit none

integer :: natom, ntypat
real(dp), parameter :: alpha = 0.90980836237716_dp ! Madelung constant
real(dp), parameter :: stress0 = 1.213077816502878_dp

real(dp) :: stress(6)
real(dp), allocatable :: xred(:, :), forces(:, :), q(:)
real(dp) :: L, E_ewald, E_ewald_fft, E_madelung
integer :: i, Ng, u

natom = 2
ntypat = 1
allocate(xred(3, natom), forces(3, natom), q(natom))
xred(:, 1) = [0, 0, 0]
xred(:, 2) = [1, 1, 1] / 2._dp
q = [1, 1]
L = 2 / sqrt(3._dp)

call ewald_box(L, xred*L, q, E_ewald, forces, stress)
E_ewald = E_ewald / (natom/ntypat)
E_madelung = -2*alpha/L
print *, "a =", L/ang2bohr*100, "pm"
print *, "Madelung: ", E_madelung / kJmol2Ha, "kJ/mol =", E_madelung, "a.u."
print *, "Ewald:    ", E_ewald / kJmol2Ha, "kJ/mol =", E_ewald, "a.u."
print *, "Ewald error:     ", abs(E_ewald - E_madelung), "a.u."
call assert(abs(E_ewald - E_madelung) < 5e-14_dp)
call assert(all(abs(forces) < 1e-15_dp))
call assert(all(abs(stress - (-stress0/L)*[1, 1, 1, 0, 0, 0]) < 1e-15_dp))

open(newunit=u, file="ewald_conv4.txt", status="replace")

Ng = 1
do i = 1, 8
    Ng = Ng*2
    print *
    print *, "Ng =",Ng

    call ewald_fft1(L, xred*L, q, Ng, 0.25_dp*L, E_ewald_fft)
    E_ewald_fft = E_ewald_fft / (natom/ntypat)

    print *, "Ewald FFT:", E_ewald_fft / kJmol2Ha, "kJ/mol =",E_ewald_fft,"a.u."
    print *, "Ewald FFT error: ", abs(E_ewald_fft - E_madelung), "a.u."

    write(u, *) Ng, E_ewald_fft
end do

close(u)

end program

max_stars_repo_name: certik/hfsolver

max_issues_repo_name: certik/hfsolver

max_forks_repo_name: certik/hfsolver


text:
      SUBROUTINE CLAKF2( M, N, A, LDA, B, D, E, Z, LDZ )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            LDA, LDZ, M, N
*     ..
*     .. Array Arguments ..
      COMPLEX            A( LDA, * ), B( LDA, * ), D( LDA, * ),
     $                   E( LDA, * ), Z( LDZ, * )
*     ..
*
*  Purpose
*  =======
*
*  Form the 2*M*N by 2*M*N matrix
*
*         Z = [ kron(In, A)  -kron(B', Im) ]
*             [ kron(In, D)  -kron(E', Im) ],
*
*  where In is the identity matrix of size n and X' is the transpose
*  of X. kron(X, Y) is the Kronecker product between the matrices X
*  and Y.
*
*  Arguments
*  =========
*
*  M       (input) INTEGER
*          Size of matrix, must be >= 1.
*
*  N       (input) INTEGER
*          Size of matrix, must be >= 1.
*
*  A       (input) COMPLEX, dimension ( LDA, M )
*          The matrix A in the output matrix Z.
*
*  LDA     (input) INTEGER
*          The leading dimension of A, B, D, and E. ( LDA >= M+N )
*
*  B       (input) COMPLEX, dimension ( LDA, N )
*  D       (input) COMPLEX, dimension ( LDA, M )
*  E       (input) COMPLEX, dimension ( LDA, N )
*          The matrices used in forming the output matrix Z.
*
*  Z       (output) COMPLEX, dimension ( LDZ, 2*M*N )
*          The resultant Kronecker M*N*2 by M*N*2 matrix (see above.)
*
*  LDZ     (input) INTEGER
*          The leading dimension of Z. ( LDZ >= 2*M*N )
*
*  ====================================================================
*
*     .. Parameters ..
      COMPLEX            ZERO
      PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IK, J, JK, L, MN, MN2
*     ..
*     .. External Subroutines ..
      EXTERNAL           CLASET
*     ..
*     .. Executable Statements ..
*
*     Initialize Z
*
      MN = M*N
      MN2 = 2*MN
      CALL CLASET( 'Full', MN2, MN2, ZERO, ZERO, Z, LDZ )
*
      IK = 1
      DO 50 L = 1, N
*
*        form kron(In, A)
*
         DO 20 I = 1, M
            DO 10 J = 1, M
               Z( IK+I-1, IK+J-1 ) = A( I, J )
   10       CONTINUE
   20    CONTINUE
*
*        form kron(In, D)
*
         DO 40 I = 1, M
            DO 30 J = 1, M
               Z( IK+MN+I-1, IK+J-1 ) = D( I, J )
   30       CONTINUE
   40    CONTINUE
*
         IK = IK + M
   50 CONTINUE
*
      IK = 1
      DO 90 L = 1, N
         JK = MN + 1
*
         DO 80 J = 1, N
*
*           form -kron(B', Im)
*
            DO 60 I = 1, M
               Z( IK+I-1, JK+I-1 ) = -B( J, L )
   60       CONTINUE
*
*           form -kron(E', Im)
*
            DO 70 I = 1, M
               Z( IK+MN+I-1, JK+I-1 ) = -E( J, L )
   70       CONTINUE
*
            JK = JK + M
   80    CONTINUE
*
         IK = IK + M
   90 CONTINUE
*
      RETURN
*
*     End of CLAKF2
*
      END

max_stars_repo_name: nashmit/lapack-release

max_issues_repo_name: nashmit/lapack-release

max_forks_repo_name: nashmit/lapack-release


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details
C=======================================================================
      SUBROUTINE RDMAP (NDB, NUMEL, MAPEL, ISEOF)
C=======================================================================

C   --*** RDMAP *** (EXPLORE) Read database element order map
C   --
C   --RDMAP reads the element order map from the database.  An error
C   --message is displayed if the end of file is read.
C   --
C   --Parameters:
C   --   NDB - IN - the database file
C   --   NUMEL - IN - the number of elements
C   --   MAPEL - OUT - the element order map
C   --   ISEOF - IN/OUT - set true if end of file read

      include 'exodusII.inc'

      INTEGER MAPEL(*)
      LOGICAL ISEOF

      CHARACTER*80 ERRMSG

      CALL INIINT (NUMEL, 0, MAPEL)

C ... Don't warn about no map stored in file
      call exopts (0, ierr1)
      call exgenm(ndb, mapel, ierr)
      call exopts (EXVRBS, ierr1)
      if (ierr .lt. 0) then
        go to 100
      else if (ierr .ne. 0) then
        do 10 i=1, numel
          mapel(i) = i
 10     continue
      end if

      RETURN

  100 CONTINUE
      WRITE (ERRMSG, 10000) 'ELEMENT ORDER MAP'
      GOTO 110
  110 CONTINUE
      CALL WDBERR (IERR, ERRMSG)
      ISEOF = .TRUE.

      RETURN

10000  FORMAT (A)
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
program real32_vs_real64
    use, intrinsic :: iso_fortran_env, only : error_unit, I8 => INT64, &
        SP => REAL32, DP => REAL64
    implicit none
    integer(kind=I8) :: nr_values
    real(kind=SP), dimension(:), allocatable :: x_sp, y_sp
    real(kind=SP) :: result
    real :: start_time, end_time
    integer :: status, i

    call get_argument(nr_values)

    ! implicit loops
    allocate (x_sp(nr_values), y_sp(nr_values), stat=status)
    if (status /= 0) then
        write (unit=error_unit, fmt='(A)') 'error: can not allocate x_sp/y_sp'
        stop 1
    end if
    call random_number(x_sp)
    call random_number(y_sp)

    call cpu_time(start_time)
    result = sum(sqrt(x_sp**2 + y_sp**2))
    call cpu_time(end_time)
    print '(A, F12.6)', 'implicit loop: ', end_time - start_time
    print '(A, F15.6)', '    result: ', result
    deallocate (x_sp, y_sp)

    ! explicit loops
    allocate (x_sp(nr_values), y_sp(nr_values), stat=status)
    if (status /= 0) then
        write (unit=error_unit, fmt='(A)') 'error: can not allocate x_sp/y_sp'
        stop 1
    end if
    call random_number(x_sp)
    call random_number(y_sp)

    call cpu_time(start_time)
    result = 0.0_SP
    do i = 1, size(x_sp)
        result = result + sqrt(x_sp(i)**2 + y_sp(i)**2)
    end do
    call cpu_time(end_time)
    print '(A, F12.6)', 'explicit loop: ', end_time - start_time
    print '(A, F15.6)', '    result: ', result
    deallocate (x_sp, y_sp)

contains

    subroutine get_argument(nr_values)
        implicit none
        integer(Kind=I8), intent(out) :: nr_values
        character(len=1024) :: buffer, msg
        integer :: status

        if (command_argument_count() /= 1) then
            write (unit=error_unit, fmt='(A)') 'error: expect number of values'
            stop 2
        end if
        call get_command_argument(1, buffer)
        read (buffer, fmt=*, iomsg=msg, iostat=status) nr_values
        if (status /= 0) then
            write (unit=error_unit, fmt='(2A)') 'error: ', trim(msg)
            stop 3
        end if
    end subroutine get_argument

end program real32_vs_real64

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
      MODULE m_qsf

!      interface qsf
!        module procedure qsfReal, qsfComplex
!      end interface

      CONTAINS

      SUBROUTINE qsf(h,y,z,ndim,isave)
c     ..................................................................
c     subroutine qsf
c     purpose
c        to integrate an equidistant table of function values.
c     description of parameters
c         h     - the increment of arguement values.
c         y     - the input vector of function values.  y must be
c                 dimensioned ndim.
c         z     - if isave = 0,z is the final value of the integration.
c                 if isave not 0,z is the resulting vector of integrals
c                 for each y value. in this case it must be dimensioned
c                 ndim.
c         ndim  - the number of points in the table.
c         isave - determines whether a single value or multiple results
c                 will be returned
c                  isave =  0  - only the final value of the integra-
c                                tion will be returned.
c                  isave ne 0  - the results of integration at each
c                                point in the table will be returned.
c     remarks
c        no action in case ndim less than 3.
c     subroutines and function subprograms required
c        none
c     method
c        integration is done by means of simpsons rule together with
c        newtons 3/8 rule or a combination of these two rules. trunca-
c        tion error is of order h**5 if ndim greater than 3. if ndim = 3
c        truncation error is of order h**4.
c     reference
c        (1) f.b.hildebrand, introduction to numerical analysis,
c            mcgraw-hill, new york/toronto/london, 1956, pp.71-76.
c        (2) r. zurmuehl, praktische mathematic fuer ingenieure und
c            physiker, springer, berlin/goettingen/heidelberg, 1963,
c            pp. 214-221.
c     ..................................................................
C     .. Scalar Arguments ..

      IMPLICIT NONE
      REAL h
      INTEGER isave,ndim
C     ..
C     .. Array Arguments ..
      REAL y(ndim),z(*)
C     ..
C     .. Local Scalars ..
      REAL aux,aux1,aux2,ht,sum1,sum2,val
      INTEGER i
      LOGICAL save
C     ..
      save = isave .NE. 0
      ht = h / 3.0
      IF ( ndim == 5 ) THEN
        GOTO 140
      ELSEIF ( ndim < 5 ) THEN
        GOTO 130
      ENDIF
c     ndim is greater than 5. preparations of integration loop
      sum1 = y(2) + y(2)
      sum1 = sum1 + sum1
      sum1 = ht* (y(1)+sum1+y(3))
      aux1 = y(4) + y(4)
      aux1 = aux1 + aux1
      aux1 = sum1 + ht* (y(3)+aux1+y(5))
      aux2 = ht* (y(1)+3.875e0* (y(2)+y(5))+2.625e0* (y(3)+y(4))+y(6))
      IF (.NOT.save) GO TO 30
   20 z(1) = 0.e0
      sum2 = y(5) + y(5)
      sum2 = sum2 + sum2
      sum2 = aux2 - ht* (y(4)+sum2+y(6))
      aux = y(3) + y(3)
      aux = aux + aux
      z(2) = sum2 - ht* (y(2)+aux+y(4))
      z(3) = sum1
      z(4) = sum2
   30 CONTINUE
      IF ( ndim <= 6 ) THEN
         GOTO 70
      ENDIF
c     integration loop
   40 DO 60 i = 7,ndim,2
         sum1 = aux1
         sum2 = aux2
         aux1 = y(i-1) + y(i-1)
         aux1 = aux1 + aux1
         aux1 = sum1 + ht* (y(i-2)+aux1+y(i))
         IF (save) z(i-2) = sum1
         IF ( i >= ndim ) THEN
           GOTO 100
         ENDIF
         aux2 = y(i) + y(i)
         aux2 = aux2 + aux2
         aux2 = sum2 + ht* (y(i-1)+aux2+y(i+1))
         IF (save) z(i-1) = sum2
   60 CONTINUE
   70 IF (.NOT.save) GO TO 90
   80 z(ndim-1) = aux1
      z(ndim) = aux2
      RETURN
   90 z(1) = aux2
      RETURN
  100 IF (.NOT.save) GO TO 120
  110 z(ndim-1) = sum2
      z(ndim) = aux1
      RETURN
  120 z(1) = aux1
      RETURN
c     end of integration loop
  130 CONTINUE
      IF ( ndim == 3 ) THEN
        GOTO 210
      ELSEIF ( ndim < 3 ) THEN
        GOTO 230
      ENDIF
c     ndim is equal to 4 or 5
  140 sum2 = 1.125e0*ht* (y(1)+y(2)+y(2)+y(2)+y(3)+y(3)+y(3)+y(4))
      sum1 = y(2) + y(2)
      sum1 = sum1 + sum1
      sum1 = ht* (y(1)+sum1+y(3))
      IF (.NOT.save) GO TO 160
  150 z(1) = 0.e0
      aux1 = y(3) + y(3)
      aux1 = aux1 + aux1
      z(2) = sum2 - ht* (y(2)+aux1+y(4))
  160 CONTINUE
      IF ( ndim < 5 ) THEN
        GOTO 190
      ENDIF
  170 aux1 = y(4) + y(4)
      aux1 = aux1 + aux1
      val = sum1 + ht* (y(3)+aux1+y(5))
      IF (save) GO TO 180
      z(1) = val
      RETURN
  180 z(5) = val
      GO TO 200
  190 IF (save) GO TO 200
      z(1) = sum2
      RETURN
  200 z(3) = sum1
      z(4) = sum2
      RETURN
c     ndim is equal to 3
  210 sum2 = y(2) + y(2)
      sum2 = sum2 + sum2
      val = ht* (y(1)+sum2+y(3))
      IF (save) GO TO 220
      z(1) = val
      RETURN
  220 z(1) = 0.e0
      z(2) = ht* (1.25e0*y(1)+y(2)+y(2)-.25e0*y(3))
      z(3) = val
  230 RETURN
      END SUBROUTINE qsf

      subroutine qsfComplex( h, y, z, ndim, isave )

        implicit none

        REAL h
        INTEGER isave,ndim
        COMPLEX y(ndim),z(*)

        real zr(ndim), zi(ndim)

        call qsf( h,  real(y), zr, ndim, isave )
        call qsf( h, aimag(y), zi, ndim, isave )
        z(1:ndim) = cmplx( zr(1:ndim), zi(1:ndim) )

      end subroutine qsfComplex

      END MODULE m_qsf

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
!! R713 rel-op
!   is .EQ.
!   or .NE.
!   or .LT.
!   or .LE.
!   or .GT.
!   or .GE.
!   or ==
!   or /=
!   or <
!   or <=
!   or >
!   or >=

a = 1 .EQ. 2
a = 1==2
a = 1.ne.2
a = 1 /= 2
a = 1.LT.2
a = 1 < 2
a = 1 .LE. 2
a = 1<=2
a = 1 .gt. 2
a = 1>2
a = 1.ge.2
a = 1 >= 2

end

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
*=======================================================================
* File:    tscord.f
* Author:  Alex Stivala
* Created: July 2008
*
* Tableau (discrete) element matching score function
*
* $Id: tscord.F 2969 2009-11-22 01:29:39Z astivala $
*=======================================================================

      integer function tscord(x, y)

      implicit none
*
* tscord - Tableau (discrete) matching score function
*
*    Return the tableau matching score between two tableau entries
*    x and y.
*    The scores are all postiive to help ensure integer solutions
*    to relaxed QP.
*    The score is 2 if the tableau entries are equal, 1 if they are
*    equal in only one position, else 0.
*
*     .. Scalar Arguments ..
      character  x(*),y(*)
*     ..
*
* Arguments
* =========
*
* x (input) CHARACTER*2
*     two character tableau code
*
* y (input) CHARACTER*2
*     two character tableau code
*
*=======================================================================
*

*     ..
*     .. Executable Statements ..
*
      if (x(1) .eq. y(1)) then
        if (x(2) .eq. y(2)) then
            tscord = 2
        else
            tscord = 1
        endif
      elseif (x(2) .eq. y(2)) then
        tscord = 1
      else
        tscord = -2
      endif
      return
      end

max_stars_repo_name: stivalaa/qptabsearch

max_issues_repo_name: stivalaa/qptabsearch

max_forks_repo_name: stivalaa/qptabsearch


text:
Julie
Mom to a wild child, 3 year old boy
fiance to a pretty cool dude
leader of Yolo County Eyes for Others
residential counselor by day (and sometimes night)
loves gardening


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine read_opt3 ( lun     , lchar   , ndim1   , ndim2p  , ndim3   ,
     &                       ifact   , dtflg   , iwidth  , dtflg3  , vrsion  ,
     &                       ioutpt  , subject , datacoll, ierr    )

!     Deltares Software Centre

!>\File
!>               Reads time dependent variables, old style in new style work file

!     Global declarations

      use rd_token       ! for the reading of tokens
      use dlwq_data      ! for definition and storage of data
      use timers         ! performance timers
      implicit none

!     declaration of arguments

      integer               , intent(in   ) :: lun(*)       !< unit numbers used
      character(len=*)      , intent(inout) :: lchar(*)     !< filenames
      integer               , intent(in   ) :: ndim1        !< first dimension (items)
      integer               , intent(in   ) :: ndim2p       !< second dimension (values)
      integer               , intent(in   ) :: ndim3        !< number of scale factors
      integer               , intent(in   ) :: ifact        !< factor between timescales
      logical               , intent(in   ) :: dtflg        !< option "date"-format
      integer               , intent(in   ) :: iwidth       !< width of output file
      logical               , intent(in   ) :: dtflg3       !< 'date'-format (F;ddmmhhss,T;yydddhh)
      real                  , intent(in   ) :: vrsion       !< program version number
      integer               , intent(in   ) :: ioutpt       !< output file option
      integer               , intent(in   ) :: subject      !< subject of input(boundary,waste,function..)
      type(t_dlwqdatacoll)  , intent(inout) :: datacoll     !< all the data blocks from file
      integer               , intent(inout) :: ierr         !< Cumulative error count

!     local declarations

      type(t_dlwqdata)                      :: dlwqdata     ! one data block
      integer                               :: iopt3        ! function type option
      integer                               :: ndim2        ! local copy of ndim2p, not zero
      integer                               :: nvarnw       ! nr of items in specific block
      integer                               :: ndim1_left   ! nr of items left to read
      integer                               :: nobrk        ! nr of breakpoint, harmonics + 1
      integer                               :: nhar         ! nr of harmonics
      integer                               :: ntotal       ! number of items read
      integer                               :: ntt          ! number of values to be read
      integer               , pointer       :: ipnt(:)      ! item pointers, either location or parameter in one block
      integer               , pointer       :: ipnt_all(:)  ! item pointers, either location or parameter all
      integer               , pointer       :: times(:)     ! times / periods
      real                  , pointer       :: scale(:)     ! scale factors
      real                  , pointer       :: phase(:)     ! phase
      real                  , pointer       :: values(:,:,:)! avregae and amplitudes
      character(len=255)                    :: cdummy       ! dummy, but used somewhere with length 255 !!!
      real                                  :: adummy       ! dummy
      real                                  :: vrswrk       ! version wrk files
      integer                               :: idummy       ! dummy
      integer                               :: itype        ! type to be read
      integer                               :: i, i2        ! loop index
      integer                               :: idata        ! index
      integer                               :: lunuit       ! dummy unit number
      logical                               :: found        ! then found
      integer                               :: ierr2        ! local error
      integer                               :: ierr_alloc   ! local error allocating
      integer(4) :: ithndl = 0
      if (timon) call timstrt( "read_opt3", ithndl )

      ! some help variables

      ndim2  = max(1,ndim2p)
      ntotal = 0
      lunut  = lun(29)
      allocate(ipnt_all(ndim1))

      write ( lunut , 2000 )

      ! read data till every item is read

      do while ( ntotal .lt. ndim1 )

         ! read function type

         if ( gettoken( iopt3, ierr2 ) .ne. 0 ) goto 105
         write ( lunut , 2010 ) iopt3
         if ( iopt3 .lt. 1 .or. iopt3 .gt. 4 ) then
            write ( lunut , 2020 )
            goto 105
         endif

         ! read item pointer for next data block

         ndim1_left = ndim1 - ntotal
         call rdpoin ( ndim1 , iopt3  , ioutpt , ipnt_all(ntotal+1:ndim1), nvarnw,
     &                 ierr  )
         allocate(ipnt(nvarnw))
         ipnt = abs(ipnt_all(ntotal+1:))

         if ( iopt3 .eq. FUNCTYPE_BLOCK .or. iopt3 .eq. FUNCTYPE_LINEAR ) then

            ! read time-dependent items on breakpoints

            if ( gettoken( nobrk, ierr2 ) .ne. 0 ) goto 105
            write( lunut, 2030 ) nobrk

            ! read the scale factors

            allocate(scale(max(ndim3,ndim2)))
            do i = 1, ndim3
               if ( gettoken( scale(i), ierr2 ) .ne. 0 ) goto 105
            enddo

            ! read the data itself

            ntt = nvarnw*ndim2
            allocate(times(nobrk))
!jvb        allocate(values(ndim1,ndim2,nobrk))
            allocate(values(ndim2,nvarnw,nobrk))
            call fmread ( nvarnw , ipnt   , ndim2  , ndim3  , scale  ,
     &                    nobrk  , times  , values , dtflg  , dtflg3 ,
     &                    ifact  , iwidth , ioutpt , ierr   )
            deallocate(scale)

         elseif ( iopt3 .eq. FUNCTYPE_HARMONIC .or. iopt3 .eq. FUNCTYPE_FOURIER ) then

            ! Read items as functions

            write ( lunut , 2050 )

            ! nr harmonics

            if ( gettoken( nhar, ierr2 ) .ne. 0 ) goto 105

            nobrk = nhar + 1
            allocate(times(nobrk),phase(nobrk),values(nvarnw,ndim2,nobrk),stat=ierr_alloc)
            if ( ierr_alloc .ne. 0 ) then
               write(lunut,2130) ierr_alloc,nvarnw,ndim2,nhar
               call srstop(1)
            endif

            vrswrk = 5.0
            call rwfun2 ( iopt3  , nhar   , nvarnw , ndim2  , ipnt   ,
     &                    values , times  , phase  , ifact  , dtflg  ,
     &                    iwidth , lunuit , dtflg3 , vrsion , ioutpt ,
     &                    vrswrk , ierr   )

         endif

         dlwqdata%functype     = iopt3
         dlwqdata%subject      = subject
         dlwqdata%igrid        = 1
         dlwqdata%extern       = .false.
         dlwqdata%filetype     = FILE_NONE
         dlwqdata%filename     = ' '
         dlwqdata%lun          = 0
         dlwqdata%param_named  = .false.
         dlwqdata%loc_named    = .false.
         dlwqdata%loc_defaults = .false.
         dlwqdata%scaled       = .false.
         dlwqdata%param_scaled = .false.
         dlwqdata%loc_scaled   = .false.
         if ( subject .eq. SUBJECT_BOUNDARY     .or.
     &        subject .eq. SUBJECT_WASTE            ) then
            dlwqdata%iorder          = ORDER_PARAM_LOC
            dlwqdata%no_loc          = nvarnw
            dlwqdata%loc_pointered   = .true.
            dlwqdata%loc_pointers    => ipnt
            dlwqdata%no_param        = ndim2
            dlwqdata%param_pointered = .false.
         elseif ( subject .eq. SUBJECT_FUNCTION .or.
     &            subject .eq. SUBJECT_SEGFUNC      ) then
            dlwqdata%iorder          = ORDER_LOC_PARAM
            dlwqdata%no_param        = nvarnw
            dlwqdata%param_pointered = .true.
            dlwqdata%param_pointers  => ipnt
            dlwqdata%no_loc          = ndim2
            dlwqdata%loc_pointered   = .false.
         else
            write(lunut,*) ' ERROR: unexpected data item for new file processing'
            write(lunut,*) ' subject type:',dlwqdata%subject
            write(  *  ,*) ' ERROR: unexpected data item for new file processing'
            call srstop(1)
         endif

         dlwqdata%no_brk =  nobrk
         dlwqdata%times  => times ; nullify(times)
         dlwqdata%phase  => phase ; nullify(phase)
         dlwqdata%values => values; nullify(values)
         nullify(ipnt)

         ! add data block to collection

         idata = dlwqdataCollAdd( datacoll, dlwqdata )

         ! Check for completion of input

         ntotal = ntotal + nvarnw

      enddo

      if ( ntotal .gt. ndim1 ) then
           write ( lunut , 2060 ) ntotal , ndim1
           ierr = ierr + 1
      endif

      ! check complete pointer structure

      do i = 1,ndim1
         found = .false.
         do i2 = 1, ndim1
            if ( abs(ipnt_all(i2)) .eq. i ) then
               if ( found ) then
                  write ( lunut , 2070 ) i
                  ierr = ierr+1
               else
                  found = .true.
               endif
            endif
         enddo
         if ( .not. found ) then
            write ( lunut , 2080 ) i
            ierr = ierr+1
         endif
      enddo

      if (timon) call timstop( ithndl )
      return

  105 ierr = ierr+1
      if (timon) call timstop( ithndl )
      return

!       Output formats

 2000 format (    ' Time variable data.')
 2010 format (  /,' Option selected : ',I2 )
 2020 format (  /,' ERROR, option not implemented')
 2030 format (   ' Number of breakpoints:',I4 )
 2050 format (    ' Block with periodic functions.')
 2060 format (    ' ERROR, too many (',I5,') items, ',I5,' expected!')
 2070 format (    ' ERROR, duplicate item:',I5)
 2080 format (    ' ERROR, non-initialised item:',I5)
 2130 format(   ' ERROR: allocating work arrays for harmonics:',I10,
     &        /,' number of items    :',I10,
     &        /,' number of values   :',I10,
     &        /,' number of harmonics:',I10)

      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
SUBROUTINE ev_time
 !=======================================
 ! This subroutine is called in:
 ! main_MD (3 places)
 !=======================================
 USE PARACLUSTER
 USE CLUSTER
 USE ENFORCE
 USE DISTANCE
 USE META
 USE SUBSTRATE
 USE ENVIRONMENT
 USE module_sticky

 IMPLICIT NONE

 INTEGER :: i, np, l_cv
 REAL(8) :: secd, firstd, rvec
 REAL(8) :: delta
 REAL(8) :: partialx(nsiz), partialy(nsiz), partialz(nsiz)
 REAL(8), DIMENSION(nsiz, num_cv) :: fxmeta, fymeta, fzmeta !Used to check the forces when collvar_wanted is .true.

 ! Metadynamics/CVs and MgO substrate stuff use their own shared
 ! parallel soubroutine, called only once per time step, when these functions are switched on.
 ! However standard MD does not use it.

 !%%%%%%%%%%%%%%%%%%%%%%%%%
 !=========================
 ! Main MD loop (over npas)
 !=========================
 DO ipas=1, npas
    tempo=(ipas+(nd_proc-1)*npas)*tstep*1.d12   !!TIME in picoseconds

    partialx(1:natom)=0.d0
    partialy(1:natom)=0.d0
    partialz(1:natom)=0.d0

    IF(type_potential=='rgl') THEN
       !Verlet cage to save time in force_rgl
       firstd=0.d0
       secd=0.d0

       DO i=1,natom
          partialx(i)=partialx(i)+du(i)    !du dv dw of natom must be initialised in case of growth!
          partialy(i)=partialy(i)+dv(i)
          partialz(i)=partialz(i)+dw(i)
          delta=partialx(i)**2+partialy(i)**2+partialz(i)**2

          IF ((delta.GT.secd).AND.(delta.LE.firstd)) THEN
             secd=delta
          ENDIF
          IF (delta.GT.firstd) THEN
             secd=firstd
             firstd=delta
          ENDIF
       ENDDO
       rvec=DSQRT(firstd)+DSQRT(secd)

       IF (rvec.GT.rshell) THEN
          DO i=1,natom
             partialx(i)=0.d0
             partialy(i)=0.d0
             partialz(i)=0.d0
          ENDDO
          CALL bigvoi
       ENDIF
    ENDIF !! bigvoi is called only for RGL-type potential

    !IN TIME: call to voisin, force, and vel at time t, finally a call to therma for the new position at time t+dt
    !
    CALL voisin   !!calculating the number of first neighbours
    !
    CALL force_choice
    !
    !
    ! RGL force, meta-force and substrate force are then in eV/arete <-- force in eV/A are divided by aretebim
    !KEEP THE UNITS OF THE FORCES in eV/arete as the time step is done with respect that choice.
    !
    !-----------------------------------------------
    ! Start of Metadynamics/CVs related calculations
    !-----------------------------------------------
    IF ((collvar_wanted).OR.(metadyn.eq.'ya').OR.mgo_substrate .OR. impl_env) THEN
       ! Calling the subroutine that calculates pair distances
       ! (It needs to be done for all the implemented CVs)
       CALL pair_distances

       if ((collvar_wanted).or.(metadyn.eq.'ya')) then
          CALL collective

          if ((metadyn.eq.'ya') .and. (ipas.GE.metaperiod)) then
             CALL metaforce

             ! Adding the meta-forces to the "real" forces
             DO l_cv = 1, num_cv
                SELECT CASE (collvar_case(l_cv))
                   !=================================================
                   CASE (1) ! CN
                      ! Adding the force contribution of metadynamics
                      fx = fx + fgx
                      fy = fy + fgy
                      fz = fz + fgz
                      !=================================================
                   CASE (2) ! C2N
                      ! Adding the force contribution of metadynamics
                      fx = fx + fg2nx
                      fy = fy + fg2ny
                      fz = fz + fg2nz
                      !=================================================
                   CASE (3) ! SFN
                      ! Adding the force contribution of metadynamics
                      fx = fx + fgsfx
                      fy = fy + fgsfy
                      fz = fz + fgsfz
                      !=================================================
                   CASE (4) ! CNN
                      ! Adding the force contribution of metadynamics
                      fx = fx + fgcnnx
                      fy = fy + fgcnny
                      fz = fz + fgcnnz
                      !=================================================
                   CASE (5) ! CN_bim
                      ! Adding the force contribution of metadynamics
                      fx = fx + cn_fgx
                      fy = fy + cn_fgy
                      fz = fz + cn_fgz
                      !=================================================
                   CASE (6) ! d_com
                      ! Adding the force contribution of metadynamics
                      fx = fx + d_com_fgx
                      fy = fy + d_com_fgy
                      fz = fz + d_com_fgz
                      !=================================================
                END SELECT
             END DO ! on l_cv
          end if

          if (collvar_wanted) then
             !-------------------------------------------------------------------------------------------
             ! Writing the value of collective variables and forces at each step (AS A CHECK)
             !
             if (ipas.eq.1) then
                write (101, *) '# ipas, ', trim(collvar_name(1)), ', ', trim(collvar_name(2))
                IF (metadyn.eq.'ya') THEN
                   WRITE(102, *) '# ipas,  Average |force (real+MT)|,  Average |MT-force(1)|,  Average |MT-force(2)|'
                END IF
             end if
             IF ((metadyn.eq.'ya').AND.(ipas.GE.metaperiod)) THEN
                DO l_cv = 1, num_cv
                   SELECT CASE (collvar_case(l_cv))
                      !=================================================
                      CASE (1) ! CN
                         fxmeta(1:natom,l_cv) = fgx
                         fymeta(1:natom,l_cv) = fgy
                         fzmeta(1:natom,l_cv) = fgz
                      !=================================================
                      CASE (2) ! C2N
                         fxmeta(1:natom,l_cv) = fg2nx
                         fymeta(1:natom,l_cv) = fg2ny
                         fzmeta(1:natom,l_cv) = fg2nz
                      !=================================================
                      CASE (3) ! SFN
                         fxmeta(1:natom,l_cv) = fgsfx
                         fymeta(1:natom,l_cv) = fgsfy
                         fzmeta(1:natom,l_cv) = fgsfz
                      !=================================================
                      CASE (4) ! CNN
                         fxmeta(1:natom,l_cv) = fgcnnx
                         fymeta(1:natom,l_cv) = fgcnny
                         fzmeta(1:natom,l_cv) = fgcnnz
                      !=================================================
                      CASE (5) ! CN_bim
                         fxmeta(1:natom,l_cv) = cn_fgx
                         fymeta(1:natom,l_cv) = cn_fgy
                         fzmeta(1:natom,l_cv) = cn_fgz
                      !=================================================
                      CASE (6) ! d_com
                         fxmeta(1:natom,l_cv) = d_com_fgx
                         fymeta(1:natom,l_cv) = d_com_fgy
                         fzmeta(1:natom,l_cv) = d_com_fgz
                      !=================================================
                   END SELECT
                END DO ! on l_cv
             END IF
             if (num_cv.eq.1) then
                write (101, *) ipas, coll_var(1)
                IF ((metadyn.eq.'ya').AND.(ipas.GE.metaperiod)) THEN
                   WRITE(102, '(1I12, 2F16.8)') ipas, SUM(SQRT(fx*fx+fy*fy+fz*fz))/natom, &
                             &SUM(SQRT(fxmeta(1:natom,1)**2 +fymeta(1:natom,1)**2 +fzmeta(1:natom,1)**2))/natom
                END IF
             else
                write (101, *) ipas, coll_var(1), coll_var(2)
                IF ((metadyn.eq.'ya').AND.(ipas.GE.metaperiod)) THEN
                   WRITE(102, '(1I12, 3F16.8)') ipas, SUM(SQRT(fx*fx+fy*fy+fz*fz))/natom, &
                             &SUM(SQRT(fxmeta(1:natom,1)**2 +fymeta(1:natom,1)**2 +fzmeta(1:natom,1)**2))/natom, &
                             &SUM(SQRT(fxmeta(1:natom,2)**2 +fymeta(1:natom,2)**2 +fzmeta(1:natom,2)**2))/natom
                END IF
             end if
             !-------------------------------------------------------------------------------------------
          end if ! (collvar_wanted)

       end if ! on the check related to collvar_wanted and metadyn
    ENDIF
    !-----------------------------------------------
    ! End of Metadynamics/CVs related calculations
    !-----------------------------------------------

    !--------------------------------------------------------------
    ! Adding forces from the substrate MgO
    IF (mgo_substrate) THEN
       IF (metal_on_top) THEN
          CALL force_mgo_mot
       ELSE
          CALL force_mgo
       ENDIF

       fx = fx + mgo_fx
       fy = fy + mgo_fy
       fz = fz + mgo_fz

       ! Summing the energy of the sub to the pot energy of the system

       ener = ener + ener_sub
    ENDIF
    !-------------------------------------------------------------

    !--------------------------------------------------------------
    ! Adding forces for the sticky atoms
    IF (sticky_atoms_wanted) THEN
       CALL sticky_force

       fx = fx + sticky_fx
       fy = fy + sticky_fy
       fz = fz + sticky_fz
    ENDIF
    !-------------------------------------------------------------

    !--------------------------------------------------------------
    ! Adding forces for the atoms - environment interaction
    IF (impl_env) THEN

       CALL fenvironment

   !DO i = 1, natom
   !write(164,*) ipas, i, fx(i), fy(i), fz(i), env_fx(i), env_fy(i), env_fz(i)
   !END DO

       fx = fx + env_fx
       fy = fy + env_fy
       fz = fz + env_fz

   !DO i = 1, natom
   !write(224,*) ipas, i, fx(i), fy(i), fz(i)
   !END DO

       ener = ener + ener_env
    ENDIF
    !-------------------------------------------------------------

    CALL vel

    !write(70,*) ipas*tstep , etot
    !write(71,*) ipas*tstep , temp

    IF ((canonical == 'ya') .and. (vel_af) ) THEN
	do i=1,natom
		v_acf(ipas,i,1)=vx(i)
		v_acf(ipas,i,2)=vy(i)
		v_acf(ipas,i,3)=vz(i)

	end do
    !write(32,*) ipas,vx(1),vy(1),vz(1)
    !CALL vel_acf_calc
    END IF
    !
    ! WRITING PHYSICAL QUANTITIES at time t
    !
    IF(quenching=='ya') CALL scrittoq   !!output file for quenching
    !
    IF((caloric=='ya').OR.(canonical=='ya')) CALL scrittoc     !!output file at time t
    !
    IF(coalescence=='ya') CALL writing_coal     !!output file at time t
    !
    IF(deposizione=='ya') CALL scrittod !!output file for natom
    !
    if(metadyn.eq.'ya') call scrittometa !output file for metadynamics
    !
    !Calculate the new position time t+tstep
    !
    CALL therma

    !
 ENDDO !!end loop on ipas

 !if ((canonical == 'ya') .and. (vel_af)) then
!	 call dff_test
 !end if
 !%%%%%%%%%%%%%%%%%%%%%%%%%
 !=========================
 ! End of the main MD loop

END SUBROUTINE ev_time

max_stars_repo_name: kcl-tscm/LoDiS

max_issues_repo_name: kcl-tscm/LoDiS

max_forks_repo_name: kcl-tscm/LoDiS