C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
                        SUBROUTINE DIGFLT
C-----------------------------------------------------------------
      INCLUDE "EXCHM.h"
      INCLUDE "parmeta"
      INCLUDE "parm.tbl"
      INCLUDE "mpp.h"
      INCLUDE "mpif.h"
#include "sp.h"
C-----------------------------------------------------------------------
                             P A R A M E T E R
     & (PIQ=3.141592654,LP1=LM+1,JAM=6+2*(JM-10))
C-----------------------------------------------------------------------
C     NTIM IS A HALF-SPAN IN TIME STEP UNITS
C     ACTUAL TIME IS NTIM*DT
C
                             P A R A M E T E R
     & (NTIM=110)
C
                             P A R A M E T E R
     1 (CM1=2937.4,CM2=4.9283,CM3=23.5518,EPS=0.622)
C-----------------------------------------------------------------------
                             L O G I C A L
     & RUN,FIRST,RESTRT,SIGMA,NEST
C-----------------------------------------------------------------------
                             R E A L
     & PDQ(IDIM1:IDIM2,JDIM1:JDIM2)
     &,TQ(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,UQ(IDIM1:IDIM2,JDIM1:JDIM2,LM),VQ(IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,RELHUM(IDIM1:IDIM2,JDIM1:JDIM2),PDSL0(IDIM1:IDIM2,JDIM1:JDIM2)
C-----------------------------------------------------------------------
      INCLUDE "CTLBLK.comm"
C-----------------------------------------------------------------------
      INCLUDE "MASKS.comm"
C-----------------------------------------------------------------------
      INCLUDE "DYNAM.comm"
C-----------------------------------------------------------------------
      INCLUDE "CONTIN.comm"
C-----------------------------------------------------------------------
      INCLUDE "VRBLS.comm"
C-----------------------------------------------------------------------
      INCLUDE "PVRBLS.comm"
C-----------------------------------------------------------------------
      DATA KNT/0/,IUNRH/51/,IUNDF/52/
C-----------------------------------------------------------------------
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C-----------------------------------------------------------------------
      NBOCO=0
      REWIND IUNRH
      REWIND IUNDF
C***
C***  COMPUTE AND WRITE OUT THE RELATIVE HUMIDITY 
C***
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        PDSL0(I,J)=RES(I,J)*PD(I,J)
      ENDDO
      ENDDO
C
      DO L=1,LM
C
        CALL ZERO2(RELHUM)
C
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          IF(HTM(I,J,L).GT.0.5)THEN
            CLOGES=-CM1/T(I,J,L)-CM2*ALOG10(T(I,J,L))+CM3
            ESE=10.**(CLOGES+2.)
            PRS=AETA(L)*PDSL0(I,J)+PT
            QIJ=Q(I,J,L)
            E=PRS*QIJ/(EPS*(1.-QIJ)+QIJ)
            RELHUM(I,J)=AMIN1(E/ESE,0.98)
          ELSE
            RELHUM(I,J)=0.
          ENDIF
        ENDDO
        ENDDO
C
C***  SMOOTH THE INITIAL RH FIELD THEN SAVE IT
C
ccccccc CALL FILT25(RELHUM,HTM(IDIM1,JDIM1,L),5,l,0)
C
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          RELHUM(I,J)=AMIN1(RELHUM(I,J),0.98)
          RELHUM(I,J)=AMAX1(RELHUM(I,J),0.02)
        ENDDO
        ENDDO
C***
C***  ASSEMBLE GLOBAL RELHUM
C***
        CALL LOC2GLB(RELHUM,TEMP1)
C
        IF(MYPE.EQ.0)WRITE(IUNRH)TEMP1
C
C***  SMOOTH THE INITIAL TEMPERATURE FIELD BEFORE EXECUTING
C***  THE DIGITAL FILTER
C
ccccccc CALL FILT25(T(IDIM1,JDIM1,L),HTM(IDIM1,JDIM1,L),5,l,1)
C
      ENDDO
C
C***  SAVE CURRENT PROG VARIABLES IN WORK FILE FOR LATER USE
C
      CALL LOC2GLB(PD,TEMP1)
      IF(MYPE.EQ.0)WRITE(IUNDF)TEMP1
C
      DO L=1,LM
        CALL LOC2GLB(T(IDIM1,JDIM1,L),TEMP1)
        IF(MYPE.EQ.0)WRITE(IUNDF)TEMP1
C
        CALL LOC2GLB(U(IDIM1,JDIM1,L),TEMP1)
        IF(MYPE.EQ.0)WRITE(IUNDF)TEMP1
C
        CALL LOC2GLB(V(IDIM1,JDIM1,L),TEMP1)
        IF(MYPE.EQ.0)WRITE(IUNDF)TEMP1
      ENDDO
C
C***  SET UP ARRAYS TO HOLD SUMS FOR DIGITAL FILTER
C
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        PDQ(I,J)=PD(I,J)
      ENDDO
      ENDDO
C
      DO L=1,LM
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          TQ(I,J,L)=T(I,J,L)
          UQ(I,J,L)=U(I,J,L)
          VQ(I,J,L)=V(I,J,L)
        ENDDO
        ENDDO
      ENDDO
C
      HSPAN=FLOAT(NTIM)*DT/3600.
C
      IF(MYPE.EQ.0)WRITE(6,100) HSPAN
 100  FORMAT(' ','INITIALIZATION CALLED WITH HALF-SPAN',F5.1,' HOURS')
C
      CALL MPI_BARRIER(MPI_COMM_COMP,IRTN)
C-----------------------------------------------------------------------
C
C      RUN THE FORECAST MODEL FORWARD AND BACKWARD (DIGITAL FILTER)
C
C*****    FIRST, FORWARD INTEGRATION *************
C
C-----------------------------------------------------------------------
C      CALCULATE NORM NEEDED FOR LANCZOS WINDOW
C
C      'TETC' DEFINES CUT-OFF FREQUENCY FOR THE WINDOW
C   (FACTOR 2 APPEARS ONLY TO SHOW GENERAL FORMULA)
C   (IT SHOULD BE TETC=PIQ/FLOAT(NTIM))
C
C      'NTIM' IS A NUMBER OF TIME STEPS IN A HALF-SPAN
C
C      'TIME' CORRESPONDS TO NUMBER OF TIME STEPS OF A SPAN
C
      TIME=2.*FLOAT(NTIM)
      QNT=FLOAT(NTIM+1)
      TETC=2.*PIQ/TIME
      CSTT=TETC/PIQ
      SQ=CSTT
C
      DO NT=1,NTIM
        FNT=FLOAT(NT)
        AS1=PIQ*FNT/QNT
        AS2=FNT*TETC
        ASS1=SIN(AS1)/AS1
        ASS2=SIN(AS2)/AS2
        SQ=SQ+2.*CSTT*ASS1*ASS2
      ENDDO
C
C----------------------------------------------------------------
C
C***  NORMALIZATION OF THE WINDOW FUNCTION
C
      CSTT=CSTT/SQ
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        PDQ(I,J)=PDQ(I,J)*CSTT
      ENDDO
      ENDDO
C
      DO L=1,LM
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          TQ(I,J,L)=TQ(I,J,L)*CSTT        
          UQ(I,J,L)=UQ(I,J,L)*CSTT        
          VQ(I,J,L)=VQ(I,J,L)*CSTT        
        ENDDO
        ENDDO
      ENDDO
C
C----------------------------------------------------------------
C--------------------------------------------------------------------
C
      NTSD=0
C
C-----------------------------------------------------------------------
C********ENTRY INTO THE TIME LOOP***************************************
 2010 NTSD=NTSD+1
      KNT=KNT+1
      IF(MYPE.EQ.0)WRITE(6,2015)NTSD,NTIM
 2015 FORMAT(' NTSD=',I5,'  NTSTM=',I4)
C***********************************************************************
C
C***  DIVERGENCE AND HORIZONTAL PART OF THE OMEGA-ALPHA TERM
C
      IF(NTSD.GT.1)CALL EXCH(T,LM,U,LM,V,LM,2,2)    !Exchange T, U, and V
C
      CALL DIVHOA
C-----------------------------------------------------------------------
C
C***  PRESS. TEND.,ETA DOT & VERTICAL OMEGA-ALPHA
C
      CALL PDTEDT
C-----------------------------------------------------------------------
C
C***  VERTICAL ADVECTION
C
      IF(MOD(NTSD-1,IDTAD).EQ.0)THEN
        CALL EXCH(ETADT,LM-1,1,1)
C
        CALL VTADVF
C
        CALL EXCH(T,LM,U,LM,V,LM,Q2,LM,1,1)
      ENDIF
C-----------------------------------------------------------------------
C
C***  UPDATE SURFACE PRESSURE (MINUS PTOP)
C
      CALL PDNEW
C-----------------------------------------------------------------------
C
C***  UPDATE H BOUNDARY POINTS
C
      IF(MOD(NTSD,IDTAD).EQ.0)THEN
        CALL EXCH(T,LM,Q2,LM,1,1)
      ENDIF
      CALL EXCH(PD,1,1,1)
C
      CALL BOCOHF
C-----------------------------------------------------------------------
C
C***  PRESSURE GRADIENT AND CORIOLIS FORCE TERMS
C
      CALL EXCH(PD,1,T,LM,2,2)            !Exchange PD and T
C
      CALL PGCOR
C
      CALL EXCH(PDSL,1,5,5)
C-----------------------------------------------------------------------
C
C***  UPDATE V BOUNDARY POINTS
C
      CALL EXCH(U,LM,V,LM,1,1)           !Exchange U and V
C
      CALL BOCOV
C-----------------------------------------------------------------------
C
C***  HORIZONTAL ADVECTION
C
      IF(MOD(NTSD,IDTAD).EQ.0)THEN
        CALL EXCH(T,LM,U,LM,V,LM,4,4)         !Exchange T, U, and V
        CALL EXCH(Q2,LM,5,5)
C
        CALL HZADV
C
        CALL EXCH(U,LM,V,LM,2,2)         !Exchange U and V
      ENDIF
C-----------------------------------------------------------------------
C
C***  CALCULATE WINDOW PARAMETER
C
      BNT=FLOAT(NTSD)
      BS1=BNT*PIQ/QNT
      BS2=BNT*TETC
      BSS1= SIN(BS1)/BS1
      BSS2= SIN(BS2)/BS2
      HNTSD=CSTT*BSS1*BSS2
C+++     HNTSD=HNTSD*1.042
C
      DO L=1,LM
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          TQ(I,J,L)=TQ(I,J,L)+HNTSD*T(I,J,L)
          UQ(I,J,L)=UQ(I,J,L)+HNTSD*U(I,J,L)
          VQ(I,J,L)=VQ(I,J,L)+HNTSD*V(I,J,L)
        ENDDO
        ENDDO
      ENDDO
C
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        PDQ(I,J)=PDQ(I,J)+HNTSD*PD(I,J)
      ENDDO
      ENDDO
C
C-----------------------------------------------------------------------
C
      IF(NTSD.EQ.NTIM) GO TO 2013
C
C***********************************************************************
      GO TO 2010
C********EXIT FROM THE TIME LOOP****************************************
C
 2013 CONTINUE
C
      CALL MPI_BARRIER(MPI_COMM_COMP,IRTN)
C-----------------------------------------------------------------------
C
C      NOW BACKWARD INTEGRATION, STARTING FROM THE INITIAL TIME
C
C-----------------------------------------------------------------------
C
C***  CHANGE (SIGN ONLY OF) IMPORTANT TIME CONSTANTS
C
      DT   =-DT
      CPGFV=-CPGFV
      EN   =-EN
      ENT  =-ENT
      F4D  =-F4D
      F4Q  =-F4Q
      EF4T =-EF4T
C
      DO JK=1,JAM
        EM (JK)=-EM (JK)
        EMT(JK)=-EMT(JK)
      ENDDO
C
      DO L=1,LM
        F4Q2(L)=-F4Q2(L)
      ENDDO
C
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        WPDAR(I,J)=-WPDAR(I,J)
        CPGFU(I,J)=-CPGFU(I,J)
        CURV (I,J)=-CURV (I,J)
        FCP  (I,J)=-FCP  (I,J)
        FAD  (I,J)=-FAD  (I,J)
        F    (I,J)=-F    (I,J)
      ENDDO
      ENDDO
C
      CALL EXCH(WPDAR,1,CPGFU,1,2,2)
      CALL EXCH(CURV,1,FCP,1,2,2)
      CALL EXCH(FAD,1,F,1,2,2)
C
C------------------- END OF CHANGE -------------------------
C
C      DEFINE INITIAL DATA FOR BACKWARD INTEGRATION
C      (PDQ,TQ,UQ,VQ ARE SUMS NEEDED FOR DIGITAL FILTER)
C
      IF(MYPE.EQ.0)THEN
        REWIND IUNDF
        READ(IUNDF)TEMP1
      ENDIF
      CALL DSTRB(TEMP1,PD,1,1,1)
C
      DO L=1,LM
        IF(MYPE.EQ.0)THEN
          READ(IUNDF)TEMP1
          READ(IUNDF)TEMP2
          READ(IUNDF)TEMP3
        ENDIF
C
        CALL DSTRB(TEMP1,T,1,LM,L)
        CALL DSTRB(TEMP2,U,1,LM,L)
        CALL DSTRB(TEMP3,V,1,LM,L)
      ENDDO
C
      CALL EXCH(T,LM,U,LM,V,LM,5,5)    !Exchange T, U, and V
      CALL EXCH(PD,1,5,5)
C--------------------------------------------------------------------
      NTSD=0
      FIRST=.TRUE.
      TSPH=-3600./DT
      IF(MYPE.EQ.0)THEN
        IF(.NOT.NEST)THEN
          REWIND NBC
          READ(NBC)
          READ(NBC)BCHR
        ELSE
          READ(NBC,REC=2)BCHR
        ENDIF
      ENDIF
C
      CALL MPI_BCAST(BCHR,1,MPI_REAL,0,MPI_COMM_COMP,IRTN)
      NBOCO=INT(BCHR*TSPH+0.5)
C
C-----------------------------------------------------------------------
C********ENTRY INTO THE TIME LOOP***************************************
 2020 NTSD=NTSD+1
      KNT=KNT+1
      IF(MYPE.EQ.0)WRITE(6,2015)NTSD,NTIM
C***********************************************************************
C
C***  DIVERGENCE AND HORIZONTAL PART OF THE OMEGA-ALPHA TERM
C
      IF(NTSD.GT.1)CALL EXCH(T,LM,U,LM,V,LM,2,2)    !Exchange T, U, and  V
C
      CALL DIVHOA
C-----------------------------------------------------------------------
C
C***  PRESS. TEND.,ETA DOT & VERTICAL OMEGA-ALPHA
C
      CALL PDTEDT
C-----------------------------------------------------------------------
C
C***  VERTICAL ADVECTION
C
      IF(MOD(NTSD-1,IDTAD).EQ.0)THEN
        CALL EXCH(ETADT,LM-1,1,1)
C
        CALL VTADVF
C
        CALL EXCH(T,LM,U,LM,V,LM,Q2,LM,1,1)          !Exchange T, U, and V
      ENDIF
C-----------------------------------------------------------------------
C
C***  UPDATE SURFACE PRESSURE (MINUS PTOP)
C
      CALL PDNEW
C-----------------------------------------------------------------------
C
C***  UPDATE H BOUNDARY POINTS
C
      IF(MOD(NTSD,IDTAD).EQ.0)THEN
        CALL EXCH(T,LM,Q2,LM,1,1)
      ENDIF
      CALL EXCH(PD,1,1,1)
C
      CALL BOCOHF
C-----------------------------------------------------------------------
C
C***  PRESSURE GRADIENT AND CORIOLIS FORCE TERMS
C
      CALL EXCH(PD,1,T,LM,2,2)            !Exchange PD and T
C
      CALL PGCOR
C
      CALL EXCH(PDSL,1,5,5)
C-----------------------------------------------------------------------
C
C***  UPDATE V BOUNDARY POINTS
C
      CALL EXCH(U,LM,V,LM,1,1)           !Exchange U and V
C
      CALL BOCOV
C-----------------------------------------------------------------------
C
C***  HORIZONTAL ADVECTION
C
      IF(MOD(NTSD,IDTAD).EQ.0)THEN
C
        CALL EXCH(T,LM,U,LM,V,LM,4,4)         !Exchange T, U, and V
        CALL EXCH(Q2,LM,5,5)
C
        CALL HZADV
C
        CALL EXCH(U,LM,V,LM,2,2)         !Exchange U and V
      ENDIF
C-----------------------------------------------------------------------
C
C***  CALCULATE WINDOW PARAMETER
C
      BNT=FLOAT(NTSD)
      BS1=BNT*PIQ/QNT
      BS2=BNT*TETC
      BSS1= SIN(BS1)/BS1
      BSS2= SIN(BS2)/BS2
      HNTSD=CSTT*BSS1*BSS2
C+++     HNTSD=HNTSD*1.042
C
      DO L=1,LM
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          TQ(I,J,L)=TQ(I,J,L)+HNTSD*T(I,J,L)
          UQ(I,J,L)=UQ(I,J,L)+HNTSD*U(I,J,L)
          VQ(I,J,L)=VQ(I,J,L)+HNTSD*V(I,J,L)
        ENDDO
        ENDDO
      ENDDO
C
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        PDQ(I,J)=PDQ(I,J)+HNTSD*PD(I,J)
      ENDDO
      ENDDO
C
C-----------------------------------------------------------------------
      IF(NTSD.EQ.NTIM) GO TO 2022
C
C***********************************************************************
      GO TO 2020
C********EXIT FROM THE TIME LOOP****************************************
C
 2022 CONTINUE
C
C-----------------------------------------------------------------------
C
C***  CHANGE BACK (SIGN ONLY) IMPORTANT TIME CONSTANTS
C
      DT   =-DT
      CPGFV=-CPGFV
      EN   =-EN
      ENT  =-ENT
      F4D  =-F4D
      F4Q  =-F4Q
      EF4T =-EF4T
C
      DO JK=1,JAM
        EM (JK)=-EM (JK)
        EMT(JK)=-EMT(JK)
      ENDDO
C
      DO L=1,LM
        F4Q2(L)=-F4Q2(L)
      ENDDO
C
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        WPDAR(I,J)=-WPDAR(I,J)
        CPGFU(I,J)=-CPGFU(I,J)
        CURV (I,J)=-CURV (I,J)
        FCP  (I,J)=-FCP  (I,J)
        FAD  (I,J)=-FAD  (I,J)
        F    (I,J)=-F    (I,J)
      ENDDO
      ENDDO
C
      CALL EXCH(WPDAR,1,CPGFU,1,2,2)
      CALL EXCH(CURV,1,FCP,1,2,2)
      CALL EXCH(FAD,1,F,1,2,2)
C
C***  UPDATE INITIALIZED PROGNOSTIC VARIALBES
C
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        PD(I,J)=PDQ(I,J)
      ENDDO
      ENDDO
C
      DO L=1,LM
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          T(I,J,L)=TQ(I,J,L)
          U(I,J,L)=UQ(I,J,L)
          V(I,J,L)=VQ(I,J,L)
        ENDDO
        ENDDO
      ENDDO
C***
C***  RETRIEVE THE INITIAL RELATIVE HUMIDITY AND COMPUTE Q
C***  SO AS TO MAINTAIN THE RH GIVEN THE NEW TEMPERATURES
C***
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        PDSL0(I,J)=RES(I,J)*PD(I,J)
      ENDDO
      ENDDO
C
      IF(MYPE.EQ.0)REWIND IUNRH
C
      DO L=1,LM
        IF(MYPE.EQ.0)READ(IUNRH)TEMP1
        CALL DSTRB(TEMP1,RELHUM,1,1,1)
C
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          IF(HTM(I,J,L).GT.0.5)THEN
            CLOGES=-CM1/T(I,J,L)-CM2*ALOG10(T(I,J,L))+CM3
            ESE=10.**(CLOGES+2.)
            PRS=AETA(L)*PDSL0(I,J)+PT
            E=RELHUM(I,J)*ESE
            Q(I,J,L)=EPS*E/(PRS+E*(EPS-1.))
          ENDIF
        ENDDO
        ENDDO
      ENDDO
C
      CALL EXCH(T,LM,Q,LM,U,LM,V,LM,4,4)  !Exchange T, Q, U, and V
      CALL EXCH(PD,1,5,5)
C
C  RETURN BC FILE TO START FORECAST
C
      NTSD=0
      IF(MYPE.EQ.0)THEN
        IF(.NOT.NEST)THEN
          REWIND NBC
          READ(NBC)
          READ(NBC)BCHR
        ELSE
          READ(NBC,REC=2)BCHR
        ENDIF
      ENDIF
C
      CALL MPI_BCAST(BCHR,1,MPI_REAL,0,MPI_COMM_COMP,IRTN)
      NBOCO=INT(BCHR*TSPH+0.5)
C-----------------------------------------------------------
C------------------- END OF CHANGE -------------------------
C-----------------------------------------------------------
      RETURN
      END