      SUBROUTINE DYN_SPS(DTD)
C
C     ADVECTION AND DIFFUSION
C
      INCLUDE 'PRMT.INC'
      INCLUDE 'basp.h'
      INCLUDE 'garr.h'
      INCLUDE 'din.h'
      INCLUDE 'fam.h'                                  ! MS
      INCLUDE 'xx.h'
      INCLUDE 'darr.h'
      real*8 tmr1, tmr2
      real extime,cputime(2)
      COMMON /NCALDL/ NCAL,NDAL,NCDL
C
C     Velocities, temperature, humidity and surface pressure from GCM
C
      COMMON/MEMORY/ XDNM(NLONP2,NLEV,5,NLAT)
      COMMON /NEWVAL/UDNM(NLONP2,NLAT,KLM),VDNM(NLONP2,NLAT,KLM),
     $TTDNM(NLONP2,NLAT,KLM),TQDNM(NLONP2,NLAT,KLM),
     $TAUDNM(NLONP2,NLAT,KLM)
CSPS      DIMENSION Y0(L,M),Y1(L,M)
      real AMFIX (L),Y0D(NC),SERW(NY),USPS(L,M,NLO1),VSPS(L,MMA,NLO),
     $WSPS(LL,M,NLO),TSPS(L,0:MMA,0:NLO1),SIGK(L),
     $PSPS(L,0:MMA,0:NLO1)
      real, dimension (NLO,M,L,NY) :: YTST,PTST,DTST
      COMMON/DYNPREP/SERD(NY,M,NLO),WDD(NY,M,NLO),AMC(L,0:MMA,0:NLO1)
      DATA SIGK /3.573e-008,7.267e-008,
     $     1.478e-007, 3.007e-007, 6.115e-007, 1.244e-006,
     *     0.00000253, 0.00000445, 0.00000780, 0.00001370,
     *     0.00002400, 0.00004210, 0.00007390, 0.00012400,
     *     0.00020300, 0.00031800, 0.00048600, 0.00073000,
     *     0.00110000, 0.00165000, 0.00248000, 0.00372000,
     *     0.00560000, 0.00842000, 0.01270000, 0.01900000,
     *     0.02860000, 0.04290000, 0.06460000, 0.09000000,
     *     0.12000000, 0.17000000, 0.23000000, 0.30000000,
     *     0.38000000, 0.46000000, 0.54000000, 0.62000000,
     *     0.70000000, 0.78000000, 0.85000000, 0.91000000,
     *     0.95000000, 0.97500000, 0.99300000 /
       DATA Avog_N/6.0225E23/, Grava/981./, R_Ud/287.E4/
!      print *,' Dyn_sps, DT=',DTD
!      call minmax1(XXALLY(2:NLONP1,2:NLAT,:,1),
!     $NLON*NLATM1*NLEV,vmin,vmax,1)
!      call minmax1(XXALLY(:,:,:,1),NLONP2*NLAT*NLEV,vmin,vmax,1)
      extime = etime(cputime)
      TIMCR0=cputime(1)
      goto 15
C
C     TEMPERATURE
C
      DO 3 IR=1,L
         I=L-IR+1
         DO 9 J=0,NLO1
            DO 9 K=1,M
 9             TSPS(I,K,J)=XDNM(J+1,IR,3,K+1)
         DO 3 J=0,NLO1
            !print *,J,JNLOH(J)
            TSPS(I,MMA,J)=TSPS(I,M,JNLOH(J))
 3          TSPS(I,0,J)=TSPS(I,1,JNLOH(J))
      !pause 1
C
C     Pressure
C
      DO 7 J=0,NLO1
         DO 7 K=1,M
            PSPS(1,K,J)=XDNM(J+1,L,5,K+1)*1.E3
            DO 7 I=2,L
  7            PSPS(I,K,J)=PSPS(1,K,J)*SIGK(L-I+1)
      DO 10 J=0,NLO1
         !print *,J,JNLOH(J)
         PSPS(:,0,J)=PSPS(:,1,JNLOH(J))
 10      PSPS(:,MMA,J)=PSPS(:,M,JNLOH(J))
      !pause 2
C
C     AIR NUMBER DENSITY
C
  15  DO 1 I=1, L
   1     AMFIX(I) = 2.55E19*EXP(-ZKM(I)/7.) !2.5e16*PA(i)
C
C     TIME STEP CALCULATION
C
      DTA=DTD              ! ******* SPS ***** 2012
      N=1                  ! ******* SPS ***** 2012
      GOTO77               ! ******* SPS ***** 2012
      TAU=DTD
      DTAC=0.5*TAU
      DTA=0.5*DTD
      N=2
      DO WHILE(DTA.GT.DTAC)
         N=N+1
         DTA=DTD/REAL(N)
      ENDDO
      TNN=1./DTA
      TN=TNN
      DO WHILE(TNN.LE.TN)
         TNN=2.*TNN
         N=2*N
         DTA=0.5*DTA
      ENDDO
  77  TN=1./DTA
!      print *,' before KDI3'
!      CALL KDI3(AMC)
!      extime = etime(cputime)
!      TIMCR1=cputime(1)
!      PRINT *,' DTIM KDI3 =',TIMCR1-TIMCR0
!      print *,' after KDI3',C(4,15,15,15)
!      print *,L,M,NLO,LL,MMA,NLO1,NY,NDAL
!      PRINT *,' N=',N,' DT=',DTD,' DTA=',DTA !NLONP1-1,NLAT-1,NLEV
      YTST(:,:,:,:)=XXALLY(2:NLONP1,2:NLAT,1:NLEV,1:NY)
      PTST(:,:,:,:)=XXALLP(2:NLONP1,2:NLAT,1:NLEV,1:NY)
      DTST(:,:,:,:)=XXALLD(2:NLONP1,2:NLAT,1:NLEV,1:NY)
!      print *,' DYN_SPS B',XXALLY(36,23,1:3,8)/XXALLY(36,23,1:3,NC)
      DO 2 IDA=1,N
C
C        IMPLICIT DIFFUSION SCHEME
C
         CALL ZK3D(YTST,C,L,M,NLO,LL,MMA,NLO1,NY,NBN,KEN,
     $   SERD,WDD,NDA,NDAL,AMC,1,DYM,PTST,DTST,TN,JNLOH,NYPR,ALON)
!      extime = etime(cputime)
!      TIMCR2=cputime(1)
!      PRINT *,' DTIM ZK3D =',TIMCR2-TIMCR1
C
C        FAMILY CORRECTION
C
         DO 5 I=1,NLO
            DO 5 J=1,M
               DO 5 K=1,L
                  CALL SEM_UNID(YTST(I,J,K,:),NY,MS,NBS,NFS,NUS)
   5              CONTINUE
!      extime = etime(cputime)
!      TIMCR3=cputime(1)
!
!      PRINT *,' DTIM SEM =',TIMCR3-TIMCR2
   2     CONTINUE
      XXALLY(2:NLONP1,2:NLAT,1:NLEV,1:NY)=YTST(:,:,:,:)
!      print *,' DYN_SPS A',XXALLY(36,23,1:3,8)/XXALLY(36,23,1:3,NC)
!      call minmax1(XXALLY(:,:,:,1),NLONP2*NLAT*NLEV,vmin,vmax,1)
!      write(50,100)XXALLY(37:72,24:46,:,1:1)
!      write(50,100)XXALLP(:,:,:,1:4)
!      write(50,100)XXALLD(:,:,:,1:4)
 100  FORMAT(8(1PE10.3))
      RETURN
      END
C
C **********************************************************************
C
      SUBROUTINE PRAT(VB,WB,LS,ZS,DT,FST,ZST,RHOY,RHOX,RHO0,SM,
     *X,XX,XXX,XY,XXY,XYY)
      INTEGER LS,ZS
      PARAMETER(N=1000)
      DIMENSION VB(LS,ZS),WB(LS,ZS),RHOY(LS,ZS),RHOX(LS,ZS),SM(LS,ZS),
     *X(LS,ZS),XX(LS,ZS),XXX(LS,ZS),XY(LS,ZS),XXY(LS,ZS),XYY(LS,ZS),
     *F0(N),FM(N),FX(N),FXX(N),FY(N),FXY(N),FYY(N),RHO0(LS,ZS),
     *ZST(*),FST(*)
      CALL PRATU(WB,ZS,LS,DT,FST,RHOY,RHO0,X,XY,XYY,XX,XXY,XXX,
     *F0,FY,FYY,FX,FXY,FXX,FM,N,1,SM)
      CALL PRATU(VB,LS,ZS,DT,ZST,RHOX,RHO0,X,XX,XXX,XY,XXY,XYY,
     *F0,FX,FXX,FY,FXY,FYY,FM,N,0,SM)
      RETURN
      END
C
C **********************************************************************
C
      SUBROUTINE PRATU(VB,LZ,ZL,DT,ZST,RHO,RHO0,X,XX,XXX,XY,XXY,XYY,
     *                 F0,FX,FXX,FY,FXY,FYY,FM,N,IP,SM)
C234567890123456789012345678901234567890123456789012345678901234567890
      INTEGER LZ,ZL
      DIMENSION F0(N),FX(N),FY(N),FXX(N),FYY(N),FXY(N),FM(N),VB(*),
     *X(*),XX(*),XXX(*),XY(*),XXY(*),XYY(*),SM(*),RHO(*),RHO0(*),ZST(*)
      DO 40 IK=1,ZL
         DO 12 IJ=1,LZ
            SLPMAX = 0.0
            IF(IP.EQ.1)THEN
               J=ZL*(IJ-1)+IK
            ELSE
               J=LZ*(IK-1)+IJ
            ENDIF 
            IF(X(J).GT.0.00) SLPMAX = X(J)  
            S1MAX = 1.50*SLPMAX 
            S1NEW = AMIN1(+S1MAX,AMAX1(-S1MAX,XX(J))) 
            S2NEW = AMIN1((SLPMAX+SLPMAX-0.3334*ABS(S1NEW)), 
     2		    AMAX1(ABS(S1NEW)-SLPMAX,XXX(J))) 
	      XXY(J) = AMIN1(+SLPMAX,AMAX1(-SLPMAX,XXY(J))) 
	      XX(J) = S1NEW 
	      XXX(J) = S2NEW  
   12	   CONTINUE
         DO 15 IP1=1,N
            fm(ip1)=0.
            F0(IP1)=0.
            FX(IP1)=0.
            FY(IP1)=0.
            FXX(IP1)=0.
            FYY(IP1)=0.
            FXY(IP1)=0.
   15	   CONTINUE
         DO 30 IJ=1,LZ-1
            IF(IP.EQ.1)THEN
               J=ZL*(IJ-1)+IK
               JP=J+ZL
            ELSE
               J=LZ*(IK-1)+IJ
               JP=J+1
            ENDIF 
            IF(VB(J).GE.0.)THEN
               AMN=1
               JPP=J
            ELSE
               AMN=-1
               JPP=JP
            ENDIF
            FM(IJ)=AMN*VB(J)*DT*RHO(J)*(ZST(IK+1)-ZST(IK))
            EPS = ABS(FM(IJ))/SM(JPP)
            EPSQ = EPS*EPS 
            EPS1 = 1.0 - EPS 
            EPS1Q = EPS1*EPS1 
C234567890123456789012345678901234567890123456789012345678901234567890
            CALL SUB(EPS,EPSQ,EPS1,EPS1Q,FM(IJ),SM(JPP),X(JPP),
     *      XX(JPP),XXX(JPP),XY(JPP),XXY(JPP),XYY(JPP),F0(IJ),FX(IJ),
     *      FXX(IJ),FY(IJ),FXY(IJ),FYY(IJ),AMN)
  30	   CONTINUE
         DO 35 IJ=1,LZ-1
            ALF=1.
            IF(IP.EQ.1)THEN
               J=ZL*(IJ-1)+IK
               JP=J+ZL
C               IF(IJ.EQ.1)ALF=2.
            ELSE
               J=LZ*(IK-1)+IJ
               JP=J+1
            ENDIF 
            IF(VB(J).GE.0.)THEN
               AMN=1
               I=JP
            ELSE
               AMN=-1
               I=J
            ENDIF 
            SM(I)=SM(I)+FM(IJ)
            EPS=FM(IJ)/SM(I)*ALF
            EPS1=1.-EPS
C234567890123456789012345678901234567890123456789012345678901234567890
            CALL RECON(EPS,EPS1,X(I),XX(I),XXX(I),XY(I),XXY(I),XYY(I),
     *      F0(IJ),FX(IJ),FXX(IJ),FY(IJ),FXY(IJ),FYY(IJ),AMN)
   35	   CONTINUE
   40 CONTINUE
      RETURN  
      END
C
C **********************************************************************
C
      SUBROUTINE SUB(EPS,EPSQ,EPS1,EPS1Q,FM,SM,X,XX,XXX,XY,XXY,XYY,
     *F0,FX,FXX,FY,FXY,FYY,AMN)
      SM=SM-FM
      F0S = EPS*(X+AMN*EPS1*(XX+AMN*(EPS1-EPS)*XXX)) 
      FXS = EPSQ*(XX+3.*AMN*EPS1*XXX) 
      FXXS = EPS*EPSQ*XXX 
      FYS = EPS*(XY+AMN*EPS1*XXY) 
      FXYS = EPSQ*XXY 
      FYYS = EPS*XYY 
      X = X - F0S 
      XX = EPS1Q*(XX - 3.*AMN*EPS*XXX) 
      XXX = EPS1*EPS1Q*XXX 
      XY = XY - FYS 
      XYY = XYY - FYYS 
      XXY = EPS1Q*XXY 
      F0 = F0+F0S
      FX = FX+FXS
      FXX = FXX+FXXS
      FY = FY+FYS
      FXY = FXY+FXYS
      FYY = FYY+FYYS
      RETURN
      END
C
C **********************************************************************
C
      SUBROUTINE RECON(EPS,EPS1,X,XX,XXX,XY,XXY,XYY,
     *F0,FX,FXX,FY,FXY,FYY,AMN)
      TEMPTM = AMN*(EPS*X - EPS1*F0) 
      X = X + F0 
      XXX = EPS*EPS*FXX+EPS1*EPS1*XXX+5.*(EPS*EPS1*AMN*
     *(XX-FX)-AMN*(EPS1-EPS)*TEMPTM) 
      XX = EPS*FX+EPS1*XX + 3.0*TEMPTM 
      XXY = EPS*FXY+EPS1*XXY+3.*AMN*(-EPS1*FY+EPS*XY) 
      XY = XY + FY 
      XYY = XYY + FYY 
      RETURN
      END
C
C **********************************************************************
C
      SUBROUTINE PROGG(A,B,C,D,X,N,K)
C
C     System of N equations A(I)*X(I-1)-B(I)*X(I)+C(I)*X(I+1)=D(I)
C
C
      DIMENSION A(N),B(N),C(N),D(N),X(N),S(N),U(N)
      M=K+1
      BB=1./B(K)
      S(K)=C(K)*BB
      U(K)=-D(K)*BB
      DO 1 I=M,N
         BA=1./(B(I)-A(I)*S(I-1))
         S(I)=C(I)*BA
   1     U(I)=(A(I)*U(I-1)-D(I))*BA
      MN=N-K
      X(N)=U(N)
      DO 2 II=1,MN
         I=N-II
   2     X(I)=S(I)*X(I+1)+U(I)
      RETURN
      END
C
C *******************************************************************
C
      SUBROUTINE SEM_UNID(Y,NC,NS,NB,NF,NU)
C
C     REDISTRIBUTION OF FAMILY GASES FOR GRIDBOX
C
      DIMENSION Y(NC),NB(NS),NF(NC,NS),NU(NC,NS)
C  NS is number of families, NB is number of species belonging to
C  each family, NF indicates which species belongs to which family, and
C  NU takes into account how many molecules of the total family value a
C  specific species makes up (e.g., N2O5=2, NO=1). GM3/01
C
      DO 1 J=1,NS
         YY=0.
         DO 2 JI=1,NB(J)
 2          YY=YY+Y(NF(JI,J))*NU(JI,J)
         IF(YY.NE.0.)THEN
            YYY=Y(J)/YY
         ELSE
            YYY=1.
         ENDIF
         DO 4 JI=1,NB(J)
 4          Y(NF(JI,J))=YYY*Y(NF(JI,J))
 1          CONTINUE
      RETURN
      END

      SUBROUTINE ZK3D(Y3,C,L,M,N,LP,MP,NP,ND,KG,KE,S,W,NDA,NDL,AM,
     $ITM,DYM,P,D,TN,JNLOH,NYPR,alon)
      DIMENSION Y3(N,M,L,ND),Y2(0:LP,0:MP,0:NP),Y5(0:LP,0:MP,0:NP),
     $C(7,L,M,N),P(N,M,L,ND),D(N,M,L,ND),YY(N),AM(L,0:MP,0:NP),
     $KG(ND),KE(ND),S(ND,M,N),W(ND,M,N),NDA(NDL),A(N),B(N),CY(N),
     $DD(N),JNLOH(0:NP),NYPR(N),alon(0:NP)
      REAL, DIMENSION (L,M,N) :: Y1,B0,D0,AMMR
      real*8 tmr1, tmr2
      real extime,cputime(2)
      Y2=0.
      Y5=0.
      TN0=TN
      DO 75 I=1,L
         DO 75 K=1,M
            DO 75 J=1,N
  75           AMMR(I,K,J)=1./AM(I,K,J)
!      extime = etime(cputime)
!      TIMCR0=cputime(1)
      DO 9 JJ=1,NDL                          ! *** Species cycle ***
         NN=NDA(JJ)                          ! *** Only calculated
         KGG=KG(NN)                          ! *** Low boundary level
         TN=TN0
         IF(NYPR(NN).eq.2)TN=TN0*2.
         ITER=1
!      extime = etime(cputime)
!      TIMCR0=cputime(1)
!      PRINT *,' Gas =',NN,NYPR(NN),1./TN !TIMCR0-TIMCR4
!      pause
C
C        Preparation for different direction mixing
C
!      print *,' NN=',NN,' DT=',1./TN
         DO 10 J=1,N
            DO 10 K=1,M
               DO 11 I=1,L
                  B0(I,K,J)=-C(4,I,K,J)+D(J,K,I,NN)+TN
                  Y1(I,K,J)=amax1(1.E-20,Y3(J,K,I,NN)) ! *** Initial value
                  D0(I,K,J)=-(P(J,K,I,NN)+TN*Y1(I,K,J))*AMMR(I,K,J)
!      print *,I,P(J,K,I,NN),TN*Y1(I,K,J),D(J,K,I,NN)*Y1(I,K,J)
  11              Y2(I,K,J)=Y1(I,K,J)*AMMR(I,K,J)      ! *** Mixing ratio
!      pause
C
C              Low boundary conditions
C
C              Surface deposition for low level condition (KGG=1)
C
               DO I=KGG,1
                  B0(I,K,J)=B0(I,K,J)+W(NN,K,J)*C(3,I,K,J)
               ENDDO
C
C              Surface flux for low level condition (KGG=1),
C              first level mixing ratio for upper level condition
C
  10           Y2(KGG-1,K,J)=S(NN,K,J) ! *** Low boundary cond.
C
C        Poles
C
         DO 12 J=1,N
            !print *,j,alon(j),JNLOH(J),alon(JNLOH(J))
            Y2(:,0,J)=Y2(:,1,JNLOH(J))
  12        Y2(:,MP,J)=Y2(:,M,JNLOH(J))
         !pause 3
C
C        Longitude boundary
C
         Y2(:,:,0)=Y2(:,:,N)
         Y2(:,:,NP)=Y2(:,:,1)
C
C        Second MR array
C
         Y5=Y2
!      if(NN.eq.8)
!     $print *,' before iterations',Y5(1,23,36),Y2(1,23,36)
!      extime = etime(cputime)
!      TIMCR1=cputime(1)
!      PRINT *,' Prepar =',TIMCR1-TIMCR0
C
C        Start iterations
C
C        First - longitude mixing
C
   7     CONTINUE
         GOTO 81
         DO 1 I=KGG,L                      ! *** Altitude loop
            DO 1 K=1,M                     ! *** Latitude loop
               DO 2 J=1,N                  ! *** Preparation for factorization 
                  A(J)= C(1,I,K,J)
                  CY(J)=C(7,I,K,J)
                  B(J)=B0(I,K,J)
                  DD(J)=D0(I,K,J)-
     $            C(2,I,K,J)*Y2(I,K-1,J)-C(6,I,K,J)*Y2(I,K+1,J)-
     $            C(3,I,K,J)*Y2(I-1,K,J)-C(5,I,K,J)*Y2(I+1,K,J)
   2              CONTINUE
               DD(1)=DD(1)-A(1)*Y2(I,K,0)
               DD(N)=DD(N)-CY(N)*Y2(I,K,NP)
!               CALL PROGG(A,B,CY,DD,Y5(I,K,1:N),N,1)
               CALL PROGG(A,B,CY,DD,YY,N,1)
               Y5(I,K,1:N)=YY(1:N)
 1             continue
C
C        Poles
C
         DO J=1,N
            Y5(:,0,J)=Y5(:,1,JNLOH(J))
            Y5(:,MP,J)=Y5(:,M,JNLOH(J))
         ENDDO
C
C        Longitude boundary
C
         Y5(:,:,0)=Y5(:,:,N)
         Y5(:,:,NP)=Y5(:,:,1)
!      if(NN.eq.8)
!     $print *,' after longitude',Y5(1,23,36),Y2(1,23,36)
!      extime = etime(cputime)
!      TIMCR2=cputime(1)
!      PRINT *,' Longitude =',TIMCR2-TIMCR1
C
C        Second - latitude mixing
C
  81     DO 3 I=KGG,L           ! *** Altitude loop       
            DO 3 J=1,N                     ! *** Longitude loop
               DO 4 K=1,M
                  A(K)= C(2,I,K,J)
                  CY(K)=C(6,I,K,J)
                  B(K)=B0(I,K,J)
                  DD(K)=D0(I,K,J)-
     $            C(1,I,K,J)*Y5(I,K,J-1)-C(7,I,K,J)*Y5(I,K,J+1)-
     $            C(3,I,K,J)*Y5(I-1,K,J)-C(5,I,K,J)*Y5(I+1,K,J)
   4              CONTINUE
               DD(1)=DD(1)-A(1)*Y5(I,0,J)
               DD(M)=DD(M)-CY(M)*Y5(I,MP,J)
!               CALL PROGG(A,B,CY,DD,Y2(I,1:M,J),M,1)
               CALL PROGG(A,B,CY,DD,YY,M,1)
               Y2(I,1:M,J)=YY(1:M)
 3             continue
C
C        Poles
C
         DO J=1,N
            Y2(:,0,J)=Y2(:,1,JNLOH(J))
            Y2(:,MP,J)=Y2(:,M,JNLOH(J))
         ENDDO
C
C        Longitude boundary
C
         Y2(:,:,0)=Y2(:,:,N)
         Y2(:,:,NP)=Y2(:,:,1)
!      if(NN.eq.8)
!     $print *,' after latitude',Y5(1,23,36),Y2(1,23,36)
!      extime = etime(cputime)
!      TIMCR3=cputime(1)
!      PRINT *,' Latitude =',TIMCR3-TIMCR2
C
C        Third - altitude mixing
C
         DO 5 K=1,M                        ! *** Latitude loop        
            DO 5 J=1,N                     ! *** Longitude loop
               DO 6 I=KGG,L
                  A(I)= C(3,I,K,J)
                  CY(I)=C(5,I,K,J)
                  B(I)=B0(I,K,J)
                  DD(I)=D0(I,K,J)-
     $            C(1,I,K,J)*Y2(I,K,J-1)-C(7,I,K,J)*Y2(I,K,J+1)-
     $            C(2,I,K,J)*Y2(I,K-1,J)-C(6,I,K,J)*Y2(I,K+1,J)
   6              CONTINUE
               DD(KGG)=DD(KGG)-A(KGG)*Y2(KGG-1,K,J)
               DD(L)=DD(L)-A(L)*Y2(LP,K,J)
!               DO I=KGG,2,-1
!                  DD(I)=DD(I)-C(3,I,K,J)*Y5(I-1,K,J)
!               ENDDO
!               CALL PROGG(A,B,CY,DD,Y5(1:L,K,J),L,KGG)
               CALL PROGG(A,B,CY,DD,YY,L,KGG)
!               Y5(1:L,K,J)=YY(1:L)
               Y5(KGG:L,K,J)=YY(KGG:L)
 5       CONTINUE
C
C        Poles
C
         DO J=1,N
            Y5(:,0,J)=Y5(:,1,JNLOH(J))
            Y5(:,MP,J)=Y5(:,M,JNLOH(J))
         ENDDO
C
C        Longitude boundary
C
         Y5(:,:,0)=Y5(:,:,N)
         Y5(:,:,NP)=Y5(:,:,1)
!      if(NN.eq.8)
!     $print *,' after altitude',Y5(1,23,36),Y2(1,23,36)
!      extime = etime(cputime)
!      TIMCR4=cputime(1)
!      PRINT *,' Altitude =',TIMCR4-TIMCR0
!      pause
         IF(ITER.GE.ITM)GOTO8
         Y2=Y5
         DY=0.
!         DO 8 K=1,M
!            DO 8 J=1,N
!               DO 8 I=1,L                     !KE(J)
!                  YNEW=Y5(I,K,J)*AM(I,K,J)
!                  DY=AMAX1(DY,ABS(YNEW-Y3(J,K,I,NN))/Y3(J,K,I,NN))
!   8              Y3(J,K,I,NN)=YNEW
         IF(DY.LE.DYM)GOTO8
         ITER=ITER+1
         GOTO7
   8     DO 9 K=1,M
            DO 9 J=1,N
               DO 9 I=1,L                     !KE(J)
                  Y3(J,K,I,NN)=Y5(I,K,J)*AM(I,K,J)
!         DO 71 I=1,L
!            DO 71 K=1,M
!               DO 71 J=1,N
! 71                PRINT 100,NN,I,K,J,Y1(I,K,J),Y2(I,K,J),Y5(I,K,J),
!     $             AM(I,K,J),C(4,I,K,J),C(1,I,K,J)
   9              CONTINUE
       TN=TN0
!      extime = etime(cputime)
!      TIMCR4=cputime(1)
!      PRINT *,' Altitude =',TIMCR4-TIMCR0
!      pause
 100  FORMAT(4I4,7(1PE10.2))
 101  FORMAT(3I3,7(1PE10.2))
 102  FORMAT(8(1PE10.2))
 103  FORMAT(a8,3I4,7(1PE10.2))
      RETURN
      END 

      SUBROUTINE MINMAX1(Y,N,VAL_MIN,VAL_MAX,KPR)
C
C     Minimum and maximum values for 1-D array
C     KPR=0 - no print
C     KPR=1 - print
C
      DIMENSION Y(N)
      VAL_MIN=1.E20
      VAL_MAX=-1.E20
      DO 1 J=1,N
         IF(Y(J).LT.VAL_MIN)VAL_MIN=Y(J)
         IF(Y(J).GT.VAL_MAX)VAL_MAX=Y(J)
 1       CONTINUE
      IF(KPR.NE.0)
     $PRINT *,'MINMAX1,  VAL_MIN=',VAL_MIN,' VAL_MAX=',VAL_MAX
      RETURN
      END

      SUBROUTINE MINMAX3(Y,N,M,L,VAL_MIN,VAL_MAX)
      DIMENSION Y(N,M,L)
      VAL_MIN=1.E20
      VAL_MAX=-1.E20
!      print *,' N=',N,' M=',M,' L=',L
      DO 1 J=1,N
         DO 1 K=1,M
            DO 1 I=1,L
!      print *,i,k,j,Y(J,K,I),VAL_MIN,VAL_MAX
!      pause
               IF(Y(J,K,I).LT.VAL_MIN)VAL_MIN=Y(J,K,I)
               IF(Y(J,K,I).GT.VAL_MAX)VAL_MAX=Y(J,K,I)
   1           CONTINUE
      PRINT *,' VAL_MIN=',VAL_MIN,' VAL_MAX=',VAL_MAX
      RETURN
      END

      SUBROUTINE DYN_PREP
      INCLUDE 'PRMT.INC'
      INCLUDE 'basp.h'
      INCLUDE 'garr.h'
      INCLUDE 'din.h'
      INCLUDE 'darr.h'
      COMMON/DYNPREP/SERD(NY,M,NLO),WDD(NY,M,NLO),AMC(L,0:MMA,0:NLO1)
      COMMON /HISDAY/ TXYZ(NLONP2,NLAT,NLEV),T_DAY(NLONP2,NLAT,NLEV)
      COMMON /HISPS/ PSDAY4(NLONP2,NLAT),PSDAY(NLONP2,NLAT)
      REAL SIGK(NLEV)
      DATA SIGK /3.573e-008,7.267e-008,
     $     1.478e-007, 3.007e-007, 6.115e-007, 1.244e-006,
     *     0.00000253, 0.00000445, 0.00000780, 0.00001370,
     *     0.00002400, 0.00004210, 0.00007390, 0.00012400,
     *     0.00020300, 0.00031800, 0.00048600, 0.00073000,
     *     0.00110000, 0.00165000, 0.00248000, 0.00372000,
     *     0.00560000, 0.00842000, 0.01270000, 0.01900000,
     *     0.02860000, 0.04290000, 0.06460000, 0.09000000,
     *     0.12000000, 0.17000000, 0.23000000, 0.30000000,
     *     0.38000000, 0.46000000, 0.54000000, 0.62000000,
     *     0.70000000, 0.78000000, 0.85000000, 0.91000000,
     *     0.95000000, 0.97500000, 0.99300000 /
C
C     LOW BONDARY CONDITIONS DEPENDING ON TIME
C
      CALL BCR(NLBC05,FLVBC,NY,ITYEAR,SER,AY,TI,NBN)
C      print *,ITYEAR,TI
C      print *,SER
C      print *,NBN
C      pause
C
C     IF SUPPOSE UNIFORM LOW BOUNDARY CONDITION
C
      DO 1 J=1,NLO
         DO 1 K=1,M
            DO 4 JG=1,NY
               WDD(JG,K,J)=WD(JG)
   4           SERD(JG,K,J)=SER(JG)
            SERD(NNOY,K,J)=BONL(K,1)*1.E-12*5.
            SERD(NCO,K,J)=BONL(K,2)*1.E-9
   1        CONTINUE
C
C     PARAMETERS FOR MIXING
C
      AZZ=EKZ
      AYY=1.E9
      AXX=1.E8
      DO 5 K=1,M
         DO 5 J=1,NLO
   5        CALL KZZ(EKZ,ZKM,L,FI(K),AZZ(1:L,K,J))
      AZZ(0,:,:)=AZZ(1,:,:)
      AZZ(LL,:,:)=AZZ(L,:,:)
!      CALL PRNS(AZZ(1:L,1:M,36),1,L,M,ZKM,FI,2,1,6,' Kzz    ')
C
C     AIR NUMBER DENSITY
C
      DO 3 I=1,L
         II=L-I+1
         DO 3 J=0,NLO1
            JJ=J+1
            DO 2 K=0,M
               KK=K+1
   2           AMC(I,K,J)=SIGK(II)*PSDAY4(JJ,KK)*1000./
     $         (TXYZ(JJ,KK,II)*1.3806e-19)
            !print *,J,JNLOH(J)+1
            AMC(I,MMA,J)=SIGK(II)*PSDAY4(JNLOH(J)+1,MMA)*1000./
     $      (TXYZ(JNLOH(J)+1,MMA,II)*1.3806e-19)
   3        CONTINUE
      CALL KDI3(AMC)
      !pause 4
      RETURN
      END

      SUBROUTINE KYY(EKYYI,ZKM,L,FI,AK)
      DIMENSION AK(L),ZKM(L),FYY(9),ZFZ(9)
      DATA FYY/1.0E1,1.0E1,1.0E1,1.0E1,6.0E0,4.0E0,2.5E0,1.6E0,1./
      DATA ZFZ/0.,2.,4.,6.,8.,10.,12.,14.,16./, DTR/0.0174533/
      DO 1 I=1,L
   1     AK(I)=EKYYI
      ZUP=14.-6.*ABS(FI)/90.
      NI=INT((ZUP+1.1)/2.)+2
      II=10-NI
      I=1
      DO WHILE(ZKM(I).LE.ZFZ(NI)+0.1)
         AMN=A18KWI(ZKM(I),ZFZ,FYY(II),1,NI)
         AK(I)=EKYYI*AMN
         I=I+1
      ENDDO
   3  RETURN
      END

      SUBROUTINE KZZ(EKZZI,ZKM,L,FI,AK)
      DIMENSION AK(L),ZKM(L),FZZ(9),ZFZ(9)
      DATA FZZ/3*5.0E1,4.0E1,3.0E1,2.0E1,1.0E1,2.5,1./
      DATA ZFZ/0.,2.,4.,6.,8.,10.,12.,14.,16./, DTR/0.0174533/
      DO 1 I=1,L
   1     AK(I)=EKZZI 
      ZUP=14.-6.*ABS(FI)/90.
      NI=INT((ZUP+1.1)/2.)+2
      II=10-NI
      I=1
      DO WHILE(ZKM(I).LE.ZFZ(NI)+0.1)
         AMN=A18KWI(ZKM(I),ZFZ,FZZ(II),1,NI)
         AK(I)=2.E3*AMN 
         I=I+1
      ENDDO
C      GOTO3
      AFI=17.*COS(FI*DTR)*COS(FI*DTR)
      ZFI1=30.-AFI
      ZFI2=55.-AFI
      DO 2 I=1,L
         IF(ZKM(I).LT.30.)GOTO2
         IF(ZKM(I).GT.ZFI1.AND.ZKM(I).LE.ZFI2)THEN
            AK(I)=(ZKM(I)-ZFI1)*720.+EKZZI
         ENDIF
         IF(ZKM(I).GT.ZFI2)THEN
            AK(I)=(ZKM(I)-ZFI2)*1000.+18000.+EKZZI
         ENDIF
    2    CONTINUE
    3 RETURN
      END



