CFPP$ NOCONCUR R
                             SUBROUTINE PROFS
C     ******************************************************************
C     *                                                                *
C     *  ROUTINE FOR VERTICAL TURBULENT MIXING OF LARGE SCALE          *
C     *  HISTORICAL VARIABLES                                          *
C     *                                                                *
C     *  SEE THE COMMENTS IN PROFQ2                                    *
C     *                                                                *
C     ******************************************************************
                             P A R A M E T E R
     & (H1=1.E0,H15=15.E0
     &, NCOL=256
CVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
     &, GHMX=.0240E0,GMMX=.3600E0
C    &, GHMX=.0320E0,GMMX=.4800E0
CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
     &, D00=0.E0,D25=.25E0,D50=.5E0
     &, EPSQ2=1.E-4,BETA=.00367E0,EPSDEN=1.E-12
     &, AKU=200.E0,CONST1=0.024/(9.8*BETA),CONST2=0.28/(9.8*BETA)
     &, A1=.92E0,B1=16.6E0,A2=.74E0,B2=10.1E0,C1=.08E0)
                             P A R A M E T E R
     & (BTAGM4=-4.E0*BETA*9.8E0
     &, DE1=A1*(1.E0-3.E0*C1),DE2=A2
     &, A1A16=A1*A1*6.E0,A1A112=A1*A1*12.E0
     &, A1A26=A1*A2*6.E0,A1A29=A1*A2*9.E0,A1A212=A1*A2*12.E0
     &, A2B23=A2*B2*3.E0
     &, A912M=-A1A112-A1A29,AB312M=-A2B23-A1A212,BACKK=0.1
     &, PAR1=A2*(1.-6.*A1/B1),PAR2=3.*A2*(B2+6.*A1),PAR3=9.*A1*A2
     &, SLOPE=135./(-3.19),SHRKM=1.E-6
     &, R=287.04E0
     &, D608=.608E0
     &, PAR4=A1*(1.-3.*C1-(6.*A1/B1)),PAR5=9.*A1*(2.*A1+A2))
C----------------------------------------------------------------------
      INCLUDE "parmeta"
      INCLUDE "parm.tbl"
C----------------------------------------------------------------------
                             P A R A M E T E R
     & (LA=13)
                             P A R A M E T E R
     & (IMJM=IM*JM-JM/2,JAM=6+2*(JM-10)
     &, KHN=IM-1,KHS=-IM,KVN=IM,KVS=-IM+1
     &, KHL01=IM+1,KHL02=2*IM
     &, KHH01=IM*(JM-1)-(JM-1)/2,KHH02=IM*(JM-2)-(JM-2)/2
     &, KVL02=2*IM
     &, KVH02=IM*(JM-2)-(JM-1)/2
     &, LM1=LM-1,LP1=LM+1,NTOT=KHH02-KHL02+1,NPASS=NTOT/NCOL
     &, LSCRCH=4*LM+1+LA +1
     &, L1=LA+1,L2=LA+LM+1,L3=LA+2*LM+1,L4=LA+3*LM+1)
C-----------------------------------------------------------------------
                             L O G I C A L
     & RUN,FIRST,RESTRT,SIGMA
C----------------------------------------------------------------------
      INCLUDE "CTLBLK"
C-----------------------------------------------------------------------
      INCLUDE "LOOPS"
C-----------------------------------------------------------------------
      INCLUDE "MASKS"
C-----------------------------------------------------------------------
      INCLUDE "PHYS"
C-----------------------------------------------------------------------
      INCLUDE "SCRCH"
C-----------------------------------------------------------------------
      INCLUDE "VRBLS"
C-----------------------------------------------------------------------
      INCLUDE "ACMSFC"
C-----------------------------------------------------------------------
      INCLUDE "PVRBLS"
      INCLUDE "DIANOS"
C-----------------------------------------------------------------------
                             D I M E N S I O N
     & UH    (IMJM),VH    (IMJM)
     &,ELH   (IMJM),RDZ   (IMJM),ELOQ2 (IMJM)
     &,UL    (IMJM),VL    (IMJM),GH    (IMJM),GM    (IMJM)
     &,C11   (IMJM),C12   (IMJM),C21   (IMJM),C22   (IMJM)
     &,DETF  (IMJM),AKMH  (IMJM)
     &,FFS   (IMJM),BKH   (IMJM),FTH   (IMJM),FQH   (IMJM)
     &,TF    (IMJM),BKL   (IMJM),CKL   (IMJM),CKLQ  (IMJM)
     &,FTHN  (IMJM),FTL   (IMJM),FQHN  (IMJM),FQL   (IMJM)
     &,ZINTV (IMJM)
     &,FUH   (IMJM),FVH   (IMJM),FUHN  (IMJM),FVHN  (IMJM)
     &,FUL   (IMJM),FVL   (IMJM)
C
     &,AKH   (IMJM,LM1),AKM   (IMJM,LM1)
     &,APE   (IMJM,LM ),THV   (IMJM,LM ),EL    (IMJM,LM )
C
     &,ZINT  (IMJM,LP1)
C-----------------------------------------------------------------------
                             D I M E N S I O N
     &  CL(NCOL,2:LM),CM(NCOL,LM),CU(NCOL,LM1),RV(NCOL,LM)
     &, SOLV(NCOL,LM),K0(NPASS+1)
C-----------------------------------------------------------------------
                             E Q U I V A L E N C E
     & (DUM   (1,01),CKLQ  (1)                                        )
     &,(DUM   (1,02),UH    (1)          ,FFS   (1),ZINTV (1)          )
     &,(DUM   (1,03),VH    (1)          ,BKH   (1)                    )
     &,(DUM   (1,04),ELH   (1)                                        )
     &,(DUM   (1,05)          ,C11   (1),FTH   (1),FUH   (1)          )
     &,(DUM   (1,06),RDZ   (1),C12   (1),FQH   (1),FVH   (1)          )
     &,(DUM   (1,07),ELOQ2 (1),C21   (1),TF    (1)                    )
     &,(DUM   (1,08),UL    (1),C22   (1),BKL   (1)                    )
     &,(DUM   (1,09),VL    (1),DETF  (1),CKL   (1)                    )
     &,(DUM   (1,10),GH    (1)          ,FTHN  (1),FUHN  (1)          )
     &,(DUM   (1,11),GM    (1)          ,FTL   (1),FUL   (1)          )
     &,(DUM   (1,12)          ,AKMH  (1),FQHN  (1),FVHN  (1)          )
     &,(DUM   (1,13)                    ,FQL   (1),FVL   (1)          )
                             E Q U I V A L E N C E
     & (DUM   (1,L1),THV   (1,1),AKM   (1,1)                          )
     &,(DUM   (1,L2),APE   (1,1)                                      )
     &,(DUM   (1,L3),EL    (1,1),AKH   (1,1)                          )
     &,(DUM   (1,L4),ZINT  (1,1)                                      )
C--------------GM,GH & ESM,ESH (EQS (33A)-(35), MY82)-----------------
CCMIC$ DO ALL VECTOR AUTOSCOPE
C         DO 210 K=KHL01,KHH01
C     UH(K)=U(K+KHS,1)+U(K-1  ,1)+U(K    ,1)+U(K+KHN,1)
C210  VH(K)=V(K+KHS,1)+V(K-1  ,1)+V(K    ,1)+V(K+KHN,1)
C-----------------------------------------------------------------------
C                   J. GERRITY    JUNE 11 1992
C                 PROPOSED MODIFICATIONS TO PROFS
C                 RELATED TO REVISED TKE EQUATION
C
C    THE BASIC IDEA OF THESE CHANGES IS TO FOLLOW THE OUTLINE
C    PROVIDED BY MARK HELFAND OF NASA GLA (1984) FOR TREATING THE
C    TURBULENT KINETIC ENERGY EQUATION AND THE EXPRESSIONS FOR THE
C    MIXING COEFFICIENTS...
C
C-----------------------------------------------------------------------
                             DO 200 L=1,LM1
CMIC$ DO ALL VECTOR AUTOSCOPE
          DO 220 K=KHL01,KHH01
        Q2KL=Q2(K,L)
C
        ELHK=EL(K,L)+EPSQ2
        ROOTQ2=SQRT(Q2(K,L))
        ELOQ2K=ELHK*ELHK/(Q2KL+EPSQ2)
C
        RDZK=1./(ZINT(K,L)-ZINT(K,L+2))
        RDZK=D50*RDZK
        GHK=THV(K,L)-THV(K,L+1)
        GHK=GHK*RDZK
        GHK=GHK*ELOQ2K*BTAGM4
C
C   INSERT GALPERIAN'S LIMITS ON GH -.28<=GH<=.0233 JPG
        IF (GHK .LT. -0.28) GHK=-0.28
        IF (GHK .GT.  0.0233) GHK=0.0233
C   INSERT GALPERIAN'S LIMITS ON GH -.28<=GH<=.0233 JPG
C 
C***
C***  SET ESH AND ESM ACCORDING TO HELFAND'S SUGGESTION
C***  FOR GROWING TURBULENCE
C***
            ESHE=PAR1/(1.-PAR2*GHK)
            ESME=(PAR5*GHK*ESHE+PAR4)/(1.-PAR3*GHK)
          AKMH(K) =ESME*ELHK*ROOTQ2+BACKK
          AKH(K,L)=ESHE*ELHK*ROOTQ2+BACKK
 220  CONTINUE
C--------------CALCULATION OF KM AT VELOCITY POINTS--------------------
C--------------REMOVING NEGATIVE EXCHANGE COEFFICIENTS-----------------
CMIC$ PARALLEL SHARED(AKMH,KVS,KVN,AKM,AKH,AKU,D00,VBM2,D25,
CMIC$1    KVL02,KVH02,KHL01,KHH01,L)  PRIVATE(K)
CMIC$ DO PARALLEL VECTOR
          DO 270 K=KVL02,KVH02
      AKM(K,L)=(AKMH(K+KVS)+AKMH(K    )+AKMH(K+1  )+AKMH(K+KVN))
     2        *VBM2(K)*D25
      AKM(K,L)=AMIN1(AKM(K,L),AKU)
      AKM(K,L)=AMAX1(AKM(K,L),D00)
 270  CONTINUE
C--------------REMOVING NEGATIVE EXCHANGE COEFFICIENTS-----------------
CMIC$ DO PARALLEL VECTOR
          DO 280 K=KHL01,KHH01
      AKH(K,L)=AMIN1(AKH(K,L),AKU)
      AKH(K,L)=AMAX1(AKH(K,L),D00)
 280  CONTINUE
CMIC$ END PARALLEL
C-----------------------------------------------------------------------
 200                         CONTINUE
C-----------------------------------------------------------------------
C ====================================================
      CALL TMPROF(ZINT,AKH,KHL01,KHH01)
C ====================================================
C
CDIA
C     START DIAGNOSTIC BLOCK.
C
C     EXCHANGE COEFFICIENTS.
      IF ( ((DIATKE).AND.(MOD(NTSD-1,NTKE).EQ.0)) .OR.
     X     ((DIAALL).AND.(MOD(NTSD-1,NALL).EQ.0)) ) THEN
         TITLE = 'KH HEAT EXCHANGE COEFFICIENT - PROFS'
         FTIME = (NTSD-1)*DT
         DO IPT = 1,NUMDIA
            KPT   = KDIAX(IPT)
            IDSTA = IDDIA(IPT)
            RLAT  = DIALAT(IPT)
            RLON  = DIALON(IPT)
            DO L = 1,LM
               DIAFLD(L) = AKH(KPT,L)
            END DO
            IF (DIATKE)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNTKE)
            IF (DIAALL)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNALL)
         END DO
C     
         TITLE = 'KM MOMENTUM EXCHANGE COEFFICIENT - PROFS'
         FTIME = (NTSD-1)*DT
         DO IPT = 1,NUMDIA
            KPT   = KDIAX(IPT)
            IDSTA = IDDIA(IPT)
            RLAT  = DIALAT(IPT)
            RLON  = DIALON(IPT)
            DO L = 1,LM
               DIAFLD(L) = AKM(KPT,L)
            END DO
            IF (DIATKE)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNTKE)
            IF (DIAALL)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNALL)
         END DO
      ENDIF
C
C     END DIAGNOSTIC BLOCK.
CDIA
C
C--------------SURFACE FLUX FACTOR--------------------------------------
CMIC$ DO ALL VECTOR AUTOSCOPE
          DO 350 K=KHL02,KHH02
      THZ0(K)=(H1-SM(K))*THS(K)+SM(K)*THZ0(K)
      QZ0 (K)=(H1-SM(K))*QS (K)+SM(K)*QZ0 (K)
      FFS (K)=AKHS(K)*D50
      CKLQ(K)=CKLQ(K)*FFS(K)
      IF (SNO(K).EQ.D00) CKLQ(K)=CKLQ(K)*WET(K)/WFK(K)
 350  CONTINUE
C
C-----------------------------------------------------------------------
C
CDIA
C     START DIAGNOSTIC BLOCK.
C
C     ARRAY THZ0.
      IF ( ((DIATKE).AND.(MOD(NTSD-1,NTKE).EQ.0)) .OR.
     X     ((DIAALL).AND.(MOD(NTSD-1,NALL).EQ.0)) ) THEN
         TITLE = 'ARRAY THZ0 - PROFS'
         FTIME = (NTSD-1)*DT
         DO IPT = 1,NUMDIA
            KPT   = KDIAX(IPT)
            IDSTA = IDDIA(IPT)
            RLAT  = DIALAT(IPT)
            RLON  = DIALON(IPT)
            DIAFLD(1) = THZ0(KPT)
            IF (DIATKE)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,1,DIAFLD,LUNTKE)
            IF (DIAALL)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,1,DIAFLD,LUNALL)
         END DO
c
         TITLE = 'ARRAY QZ0 - PROFS'
         FTIME = (NTSD-1)*DT
         DO IPT = 1,NUMDIA
            KPT   = KDIAX(IPT)
            IDSTA = IDDIA(IPT)
            RLAT  = DIALAT(IPT)
            RLON  = DIALON(IPT)
            DIAFLD(1) = QZ0(KPT)
            IF (DIATKE)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,1,DIAFLD,LUNTKE)
            IF (DIAALL)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,1,DIAFLD,LUNALL)
         END DO
C
         TITLE = 'ARRAY FFS - PROFS'
         FTIME = (NTSD-1)*DT
         DO IPT = 1,NUMDIA
            KPT   = KDIAX(IPT)
            IDSTA = IDDIA(IPT)
            RLAT  = DIALAT(IPT)
            RLON  = DIALON(IPT)
            DIAFLD(1) = FFS(KPT)
            IF (DIATKE)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,1,DIAFLD,LUNTKE)
            IF (DIAALL)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,1,DIAFLD,LUNALL)
         END DO
C
         TITLE = 'ARRAY CKLQ - PROFS'
         FTIME = (NTSD-1)*DT
         DO IPT = 1,NUMDIA
            KPT   = KDIAX(IPT)
            IDSTA = IDDIA(IPT)
            RLAT  = DIALAT(IPT)
            RLON  = DIALON(IPT)
            DIAFLD(1) = CKLQ(KPT)
            IF (DIATKE)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,1,DIAFLD,LUNTKE)
            IF (DIAALL)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,1,DIAFLD,LUNALL)
         END DO
C
      ENDIF
C
C     END DIAGNOSTIC BLOCK.
CDIA
C
C-----------------------------------------------------------------------
C***
C***  COMPUTE IMPLICIT VERTICAL DIFFUSION OF TEMPERATURE
C***  AND SPECIFIC HUMIDITY IN GROUPS OF NCOL COLUMNS
C***
      NSEND=NTOT/NCOL
      MODNC=MOD(NTOT,NCOL)
      KLAST=NCOL
      IF(MODNC.GT.0)THEN
        NSEND=NSEND+1
        KLAST=MODNC
      ENDIF
      DO  5000 NN=1,NSEND
      K0(NN) = KHL02 + (NN-1)*NCOL -1
 5000 CONTINUE
CMIC$ DO ALL SHARED (K0,ZINT,DTQ2,HTM,AKH,FFS,THZ0,NCOL,APE,T,CKLQ)
CMIC$1  PRIVATE(NN,L,K,KK,K2,CL,CM,CU,RV,SOLV) AUTOSCOPE
      DO 5100 NN=1,NSEND
      K2=NCOL
      IF(NN.EQ.NSEND)K2=KLAST
C***
C***  WE BEGIN BY COMPUTING THE COEFFICIENTS OF THE TRIDIAGONAL MATRIX
C***  FOR POTENTIAL TEMPERATURE AND WATER VAPOR DIFFUSION.
C***
C
C.....The loops below, up to 'ADD EFFECT OF SURFACE FLUX...' have been
C.....slightly re-arranged from the 'natural' order to take advantage
C.....of multi-tasking.
C
C
C***                    CL LOWER DIAGONAL
C***
C***                    CU UPPER DIAGONAL
C***
C***                    CM MAIN  DIAGONAL
C***
      DO 1000 L=2,LM1
CDIR$ IVDEP
      DO 1000 KK=1,K2
      K=K0(NN)+KK
      CL(KK,L)=HTM(K,L)*(-2.*AKH(K,L-1)*DTQ2)/
     1        ((ZINT(K,L)-ZINT(K,L+1))*(ZINT(K,L-1)-ZINT(K,L+1)))
      CU(KK,L)=HTM(K,L+1)*(-2.*AKH(K,L)*DTQ2)/
     1        ((ZINT(K,L)-ZINT(K,L+1))*(ZINT(K,L)-ZINT(K,L+2)))
      CM(KK,L)=1.+(2.*DTQ2/(ZINT(K,L)-ZINT(K,L+1)))*
     1        (HTM(K,L)*AKH(K,L-1)/(ZINT(K,L-1)-ZINT(K,L+1))
     2        +HTM(K,L+1)*AKH(K,L)/(ZINT(K,L)-ZINT(K,L+2)))
 1000 CONTINUE
C
C
C                    End cases; L=LM, L=1.
C
      DO 1020 KK=1,K2
      L=LM
      K=K0(NN)+KK
      CL(KK,L)=HTM(K,L)*(-2.*AKH(K,L-1)*DTQ2)/
     1        ((ZINT(K,L)-ZINT(K,L+1))*(ZINT(K,L-1)-ZINT(K,L+1)))
      CU(KK,1)=HTM(K,1+1)*(-2.*AKH(K,1)*DTQ2)/
     1        ((ZINT(K,1)-ZINT(K,1+1))*(ZINT(K,1)-ZINT(K,1+2)))
      CM(KK,1)=1.+(2.*AKH(K,1)*DTQ2)/
     1        ((ZINT(K,1)-ZINT(K,1+1))*(ZINT(K,1)-ZINT(K,1+2)))
      CM(KK,L)=1.+(2.*DTQ2/(ZINT(K,L)-ZINT(K,L+1)))*
     1        (HTM(K,L)*AKH(K,L-1)/(ZINT(K,L-1)-ZINT(K,L+1)))
 1020 CONTINUE
C***
C***           NOW SET UP RIGHT HAND SIDE VECTOR  RV
C***
      DO 1050 L=1,LM
CDIR$ IVDEP
      DO 1050 KK=1,K2
      K=K0(NN)+KK
      RV(KK,L)=T(K,L)*APE(K,L)
 1050 CONTINUE
C***
C***                LMH LAST LAYER ABOVE GROUND
C***           ADD EFFECT OF SURFACE FLUX TO CM FOR L=LMH
C***
      DO 1060 KK=1,K2
      K=K0(NN)+KK
      L=LMH(K)
      RV(KK,L)=RV(KK,L)+2.*DTQ2*FFS(K)*THZ0(K)/(ZINT(K,L)-ZINT(K,L+1))
      CM(KK,L)=CM(KK,L)+2.*DTQ2*FFS(K)/(ZINT(K,L)-ZINT(K,L+1))
 1060 CONTINUE
C***
C***        NOW CALL SOLUTION OF DIFFUSION EQUATION
C***
      CALL TRIDI(NCOL,1,K2,LM,CL,CM,CU,RV,SOLV)
C***
C***  CONVERT FROM POTENTIAL TEMPERATURE TO THERMODYNAMIC TEMPERATURE
C***
      DO 1070 L=1,LM
CDIR$ IVDEP
      DO 1070 KK=1,K2
      K=K0(NN)+KK
      T(K,L)=SOLV(KK,L)/APE(K,L)
 1070 CONTINUE
C***
C***  NOW SOLVE FOR DIFFUSION OF WATER VAPOR
C***
C***   SET UP RIGHT HAND SIDE VECTOR FOR SPECIFIC HUMIDITY
C***
      DO 1090 L=1,LM
CDIR$ IVDEP
      DO 1090 KK=1,K2
      K=K0(NN)+KK
      RV(KK,L)=Q(K,L)
 1090 CONTINUE
C***
C***                    LMH LAST LAYER ABOVE GROUND
C***  MODIFY ONLY CM(K,LMH) TO REPLACE HEAT FLUX WITH VAPOR FLUX
C***
      DO 1100 KK=1,K2
      K=K0(NN)+KK
      L=LMH(K)
      RV(KK,L)=RV(KK,L)+2.*DTQ2*CKLQ(K)*QZ0(K)/
     1                   (ZINT(K,L)-ZINT(K,L+1))
c     CM(KK,L)=CM(KK,L)+2*DT*FFS(K)/(ZINT(K,L)-ZINT(K,L+1))
      CM(KK,L)=CM(KK,L)+2.*DTQ2*(CKLQ(K)-FFS(K))/
     1        (ZINT(K,L)-ZINT(K,L+1))
 1100 CONTINUE
C***
C***        NOW CALL SOLUTION OF DIFFUSION EQUATION
C***
      CALL TRIDI(NCOL,1,K2,LM,CL,CM,CU,RV,SOLV)
C***
C***  SAVE THE SOLUTION IN Q(K,L) SPECIFIC HUMIDITY SOLUTION
C***
      DO 1110 L=1,LM
CDIR$ IVDEP
      DO 1110 KK=1,K2
      K=K0(NN)+KK
      Q(K,L)=SOLV(KK,L)
 1110 CONTINUE
 5100 CONTINUE
C-----------------------------------------------------------------------
C
CDIA
C     START DIAGNOSTIC BLOCK.
C
C     TEMPERATURE AND SPECIFIC HUMIDITY AFTER VERTICAL DIFFUSION
C
      IF ( ((DIATMP).AND.(MOD(NTSD-1,NTMP).EQ.0)) .OR.
     X     ((DIAALL).AND.(MOD(NTSD-1,NALL).EQ.0)) ) THEN
         TITLE = 'T AFTER VERTICAL DIFFUSION - PROFS'
         FTIME = (NTSD-1)*DT
         DO IPT = 1,NUMDIA
            KPT   = KDIAX(IPT)
            IDSTA = IDDIA(IPT)
            RLAT  = DIALAT(IPT)
            RLON  = DIALON(IPT)
            DO L = 1,LM
               DIAFLD(L) = T(KPT,L)
            END DO
            IF (DIATMP)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNTMP)
            IF (DIAALL)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNALL)
         END DO
      ENDIF
C
      IF ( (DIAMST).AND.(MOD(NTSD-1,NMST).EQ.0) ) THEN
         TITLE = 'Q AFTER VERTICAL DIFFUSION - PROFS'
         FTIME = (NTSD-1)*DT
         DO IPT = 1,NUMDIA
            KPT   = KDIAX(IPT)
            IDSTA = IDDIA(IPT)
            RLAT  = DIALAT(IPT)
            RLON  = DIALON(IPT)
            DO L = 1,LM
               DIAFLD(L) = Q(KPT,L)
            END DO
            IF (DIAMST)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNMST)
            IF (DIAALL)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNALL)
         END DO
      ENDIF
C
C     END DIAGNOSTIC BLOCK.
CDIA
C
C
C--------------INTERFACE HEIGHT AT V POINTS & SURF. FLUX FACTOR---------
C
CMIC$ PARALLEL SHARED(ZINT,KVL02,KVH02,LP1,KVS,KVN,D25,AKMS,D50,FFS)
CMIC$1   PRIVATE(L,K,ZINTV)
CMIC$ DO PARALLEL
                             DO 560 L=1,LP1
CDIR$ IVDEP
          DO 570 K=KVL02,KVH02
 570  ZINTV(K)=ZINT(K+KVS,L)+ZINT(K    ,L)+ZINT(K+1  ,L)+ZINT(K+KVN,L)
CDIR$ IVDEP
          DO 560 K=KVL02,KVH02
 560  ZINT(K,L)=ZINTV(K)*D25
C
CMIC$ DO PARALLEL VECTOR
          DO 580 K=KVL02,KVH02
 580  FFS(K)=AKMS(K)*D50
CMIC$ END PARALLEL
C-----------------------------------------------------------------------
C***
C***  SET UP COEFFICIENTS CL,CM AND CU FOR MOMENTUM
      NSEND=NTOT/NCOL
      MODNC=MOD(NTOT,NCOL)
      KLAST=NCOL
      IF(MODNC.GT.0)THEN
        NSEND=NSEND+1
        KLAST=MODNC
      ENDIF
C***
C***  FIRST FOR U COMPONENT
C***
CMIC$ DO ALL SHARED (K0,ZINT,DTQ2,HTM,AKM,FFS,NCOL,VTM,LMV,UZ0,VZ0,
CMIC$1  D25,SFCUVX,U,V)
CMIC$1  PRIVATE(NN,L,K,KK,K2,CL,CM,CU,RV,SOLV) AUTOSCOPE
      DO 5200 NN=1,NSEND
      K2=NCOL
      IF(NN.EQ.NSEND)K2=KLAST
C***
C***               CL LOWER DIAGONAL
C***
C***
C***               CU UPPER DIAGONAL
C***
C***
C***               CM MAIN  DIAGONAL
C***
      DO 2000 L=2,LM1
CDIR$ IVDEP
      DO 2000 KK=1,K2
      K=K0(NN)+KK
      CL(KK,L)=VTM(K,L)*(-2.*AKM(K,L-1)*DTQ2)/
     1        ((ZINT(K,L)-ZINT(K,L+1))*(ZINT(K,L-1)-ZINT(K,L+1)))
      CU(KK,L)=VTM(K,L+1)*(-2.*AKM(K,L)*DTQ2)/
     1        ((ZINT(K,L)-ZINT(K,L+1))*(ZINT(K,L)-ZINT(K,L+2)))
      CM(KK,L)=1.+(2.*DTQ2/(ZINT(K,L)-ZINT(K,L+1)))*
     1        (VTM(K,L)*AKM(K,L-1)/(ZINT(K,L-1)-ZINT(K,L+1))
     2        +VTM(K,L+1)*AKM(K,L)/(ZINT(K,L)-ZINT(K,L+2)))
 2000 CONTINUE
C***
C***               End cases.
C***
      DO 2020 KK=1,K2
      L=LM
      K=K0(NN)+KK
      CL(KK,L)=VTM(K,L)*(-2.*AKM(K,L-1)*DTQ2)/
     1        ((ZINT(K,L)-ZINT(K,L+1))*(ZINT(K,L-1)-ZINT(K,L+1)))
      CU(KK,1)=VTM(K,1+1)*(-2.*AKM(K,1)*DTQ2)/
     1        ((ZINT(K,1)-ZINT(K,1+1))*(ZINT(K,1)-ZINT(K,1+2)))
      CM(KK,1)=1.+(2.*AKM(K,1)*DTQ2)/
     1        ((ZINT(K,1)-ZINT(K,1+1))*(ZINT(K,1)-ZINT(K,1+2)))
      CM(KK,L)=1.+(2.*DTQ2/(ZINT(K,L)-ZINT(K,L+1)))*
     1        (VTM(K,L)*AKM(K,L-1)/(ZINT(K,L-1)-ZINT(K,L+1)))
 2020 CONTINUE
C***
C***              NOW SET UP RIGHT HAND SIDE VECTOR  RV
C***
      DO 2050 L=1,LM
CDIR$ IVDEP
      DO 2050 KK=1,K2
      K=K0(NN)+KK
      RV(KK,L)=U(K,L)
 2050 CONTINUE
C***
C***              ADD EFFECT OF SURFACE FLUX TO CM FOR L=LMV
C***              LMV LAST LAYER ABOVE GROUND
C***
      DO 2060 KK=1,K2
      K=K0(NN)+KK
      L=LMV(K)
      CM(KK,L)=CM(KK,L)+2.*DTQ2*FFS(K)/(ZINT(K,L)-ZINT(K,L+1))
      RV(KK,L)=RV(KK,L)+2.*DTQ2*FFS(K)*UZ0(K)/
     1            (ZINT(K,L)-ZINT(K,L+1))
 2060 CONTINUE
C***
C***          NOW CALL SOLUTION OF DIFFUSION EQUATION
C***
      CALL TRIDI(NCOL,1,K2,LM,CL,CM,CU,RV,SOLV)
C***
C***  TRANSFER SOLUTION TO U(K,L)
C***  NOW SOLVE FOR DIFFUSION OF V COMPONENT
C***  NEED NOT MODIFY CL,CM NOR CU
C***
C***  SET UP RIGHT HAND SIDE VECTOR FOR V COMPONENT
C***
      DO 2090 L=1,LM
CDIR$ IVDEP
      DO 2090 KK=1,K2
      K=K0(NN)+KK
      U(K,L)=SOLV(KK,L)
      RV(KK,L)=V(K,L)
 2090 CONTINUE
C***
C***                LMV LAST LAYER ABOVE GROUND
C***
      DO 2100 KK=1,K2
      K=K0(NN)+KK
      L=LMV(K)
C     RV(KK,L)=RV(KK,L)+2.*DT*FFS(K)*UZ0(K)/(ZINT(K,L)-ZINT(K,L+1))
      RV(KK,L)=RV(KK,L)+2.*DTQ2*FFS(K)*VZ0(K)/
     1                  (ZINT(K,L)-ZINT(K,L+1))
 2100 CONTINUE
C***
C***        NOW CALL SOLUTION OF DIFFUSION EQUATION
C***
      CALL TRIDI(NCOL,1,K2,LM,CL,CM,CU,RV,SOLV)
C***
C***  SAVE THE SOLUTION IN V(K,L)
C***
      DO 2110 L=1,LM
CDIR$ IVDEP
      DO 2110 KK=1,K2
      K=K0(NN)+KK
      V(K,L)=SOLV(KK,L)
 2110 CONTINUE
 5200 CONTINUE
C
C   COMPUTE SURFACE MOMENTUM FLUX FOR GCIP.
C
      DO K=KVL02,KVH02
        L=LMV(K)
        PLMV=AETA(L)*(RES(K)*PD(K)+RES(K+1)*PD(K+1)+
     &    RES(K+KVS)*PD(K+KVS)+RES(K+KVN)*PD(K+KVN))*D25+PT
        TLMV=(T(K,L)+T(K+1,L)+T(K+KVS,L)+T(K+KVN,L))*D25
        QLMV=(Q(K,L)+Q(K+1,L)+Q(K+KVS,L)+Q(K+KVN,L))*D25
        DENS=PLMV/(R*TLMV*(QLMV*D608+H1))
        DELU=U(K,L)-UZ0(K)
        DELV=V(K,L)-VZ0(K)
        SFCUVX(K) = ((DELU*DELU)+(DELV*DELV)) * AKMS(K) * DENS +
     &          SFCUVX(K)
      ENDDO
C-----------------------------------------------------------------------
C
CDIA
C     START DIAGNOSTIC BLOCK.
C
C     U AND V AFTER VERTICAL DIFFUSION.
C
      IF ( ((DIAUNV).AND.(MOD(NTSD-1,NUNV).EQ.0)) .OR.
     X     ((DIAALL).AND.(MOD(NTSD-1,NALL).EQ.0)) ) THEN
         TITLE = 'U AFTER VERTICAL DIFFUSION - PROFS'
         FTIME = (NTSD-1)*DT
         DO IPT = 1,NUMDIA
            KPT   = KDIAX(IPT)
            IDSTA = IDDIA(IPT)
            RLAT  = DIALAT(IPT)
            RLON  = DIALON(IPT)
            DO L = 1,LM
               DIAFLD(L) = U(KPT,L)
            END DO
            IF (DIAUNV)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNUNV)
            IF (DIAALL)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNALL)
         END DO
C
         TITLE = 'V AFTER VERTICAL DIFFUSION - PROFS'
         FTIME = (NTSD-1)*DT
         DO IPT = 1,NUMDIA
            KPT   = KDIAX(IPT)
            IDSTA = IDDIA(IPT)
            RLAT  = DIALAT(IPT)
            RLON  = DIALON(IPT)
            DO L = 1,LM
               DIAFLD(L) = V(KPT,L)
            END DO
            IF (DIAUNV)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNUNV)
            IF (DIAALL)
     X       CALL DWRITE(TITLE,IDSTA,RLAT,RLON,FTIME,LM,DIAFLD,LUNALL)
         END DO
      ENDIF
C
C     END DIAGNOSTIC BLOCK.
CDIA
C
C-----------------------------------------------------------------------
                             RETURN
                             END
