      SUBROUTINE TURBL
C     ******************************************************************
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    TURBL       VERTICAL TURBULENT EXCHANGE
C   PRGRMMR: JANJIC          ORG: W/NP2      DATE: 95-03-20       
C     
C ABSTRACT:
C     TURBL UPDATES THE TURBULENT KINETIC ENERGY WITH THE PROD-
C     UCTION/DISSIPATION TERM AND THE VERTICAL DIFFUSION TERM
C     DIFFUSION TERM (USING AN IMPLICIT FORMULATION).  EXCHANGE
C     COEFFICIENTS FOR THE SURFACE AND FOR ALL LAYER INTERFACES
C     ARE THEN COMPUTED AND THE EXCHANGE IS EXECUTED.
C     
C PROGRAM HISTORY LOG:
C   95-03-15  JANJIC     - ORIGINATOR
C   95-03-28  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
C   96-03-29  BLACK      - ADDED EXTERNAL EDGE; REMOVED SCRCH COMMON
C   96-07-19  MESINGER   - ADDED Z0 EFFECTIVE
C   98-??-??  TUCCILLO   - MODIFIED FOR CLASS VIII PARALLELISM
C   98-10-27  BLACK      - PARALLEL CHANGES INTO MOST RECENT CODE
C   14-07-XX  KATSAFADOS - COUPLING WITH OCEAN WAVE
C     
C USAGE: CALL TURBL FROM MAIN PROGRAM EBU
C   INPUT ARGUMENT LIST:
C       NONE     
C  
C   OUTPUT ARGUMENT LIST: 
C     NONE
C     
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C     UNIQUE: MIXLEN
C             PRODQ2
C             DIFCOF
C             SFCDIF
C             VDIFH
C             VDIFQ
C             VDIFV
C  
C     LIBRARY: NONE
C  
C   COMMON BLOCKS: CTLBLK
C                  LOOPS
C                  MASKS
C                  DYNAM
C                  PHYS2
C                  VRBLS
C                  PVRBLS
C                  INDX
C                  Z0EFFT
C   
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE : IBM SP
C$$$  
C***********************************************************************
C-----------------------------------------------------------------------
C
      INCLUDE "EXCHM.h"
      INCLUDE "parmeta"
      INCLUDE "mpif.h"   !added for MPI_COMM_WORLD pkatsaf
      INCLUDE "mpp.h"
#include "sp.h"
      INCLUDE "wave.inc"
C-----------------------------------------------------------------------
                             P A R A M E T E R
     & (KTMQ2=1,CAPA=0.28589641,G=9.8,RG=1./G,ROG=287.04/G
     &, EPSZ=1.E-4,EPSQ2=0.2
     &, IMJM=IM*JM-JM/2,LM1=LM-1,LP1=LM+1,JAM=6+2*(JM-10)
     &, ITB=76,JTB=134,ITBQ=152,JTBQ=440
     &, NHRZ=(IDIM2-IDIM1+1)*(JDIM2-JDIM1+1))
C-----------------------------------------------------------------------
      INCLUDE "CTLBLK.comm"
C-----------------------------------------------------------------------
      INCLUDE "LOOPS.comm"
C-----------------------------------------------------------------------
      INCLUDE "MASKS.comm"
C-----------------------------------------------------------------------
      INCLUDE "DYNAM.comm"
C-----------------------------------------------------------------------
      INCLUDE "PHYS2.comm"
C-----------------------------------------------------------------------
      INCLUDE "VRBLS.comm"
C-----------------------------------------------------------------------
      INCLUDE "PVRBLS.comm"
C-----------------------------------------------------------------------
      INCLUDE "INDX.comm"
C-----------------------------------------------------------------------
      INCLUDE "Z0EFFT.comm"
C-----------------------------------------------------------------------
      INCLUDE "WAVES.comm"
C-----------------------------------------------------------------------
                             L O G I C A L
     & RUN,FIRST,RESTRT,SIGMA
C-----------------------------------------------------------------------
                             R E A L
     & CKLQ(IDIM1:IDIM2,JDIM1:JDIM2)
     &,CT  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,APE (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,AKH (IDIM1:IDIM2,JDIM1:JDIM2,LM1)
     &,AKM (IDIM1:IDIM2,JDIM1:JDIM2,LM1)
     &,ZINT(IDIM1:IDIM2,JDIM1:JDIM2,LP1)
     &,UZ0H(IDIM1:IDIM2,JDIM1:JDIM2)
     &,VZ0H(IDIM1:IDIM2,JDIM1:JDIM2)
! coupling
     &,ACPRECpast(IDIM1:IDIM2,JDIM1:JDIM2)
! coupling
C
                             R E A L
     & AKMCOL(IDIM1:IDIM2,JDIM1:JDIM2,LM1) 
     &,AKHCOL(IDIM1:IDIM2,JDIM1:JDIM2,LM1) 
     &,AKMSV (IDIM1:IDIM2,JDIM1:JDIM2) 
     &,ZCOL  (IDIM1:IDIM2,JDIM1:JDIM2,LP1) 
     &,UCOL  (IDIM1:IDIM2,JDIM1:JDIM2,LM) 
     &,VCOL  (IDIM1:IDIM2,JDIM1:JDIM2,LM) 
C
                             R E A L
     & AKH_T   (LM1,IDIM1:IDIM2,JDIM1:JDIM2) 
     &,AKM_T   (LM1,IDIM1:IDIM2,JDIM1:JDIM2) 
     &,APECOL_T(LM,IDIM1:IDIM2,JDIM1:JDIM2) 
     &,ZCOL_T  (LP1,IDIM1:IDIM2,JDIM1:JDIM2) 
     &,ZCOL_T2 (LP1,IDIM1:IDIM2,JDIM1:JDIM2) 
     &,UCOL_T  (LM,IDIM1:IDIM2,JDIM1:JDIM2) 
     &,VCOL_T  (LM,IDIM1:IDIM2,JDIM1:JDIM2) 
     &,TCOL_T  (LM,IDIM1:IDIM2,JDIM1:JDIM2) 
     &,QCOL_T  (LM,IDIM1:IDIM2,JDIM1:JDIM2) 
     &,Q2COL_T (LM,IDIM1:IDIM2,JDIM1:JDIM2)
C
                             R E A L 
     & GM(LM1),GH(LM1),EL(LM1),ZEFF(4)
C-----------------------------------------------------------------------
C***
C***  THE FOLLOWING ARE USED FOR TIMIMG PURPOSES ONLY
C***
      real*8 timef
      real nhb_tim,mpp_tim,init_tim
      common/timing/surfce_tim,nhb_tim,res_tim,exch_tim
C***********************************************************************
      integer :: npest,npestm1
cpk   Coupling
      real WSpeed10, DIRW10
      real RRT, VVD, DDR, RRpast, RRnow, RRtot
C-----------------------------------------------------------------------
      IF(MYPE.EQ.0) THEN
         IF (NTSD.EQ.1) THEN
            write(6,*)'TURBL ZETA=DFL/G ',DFL/9.80
            RRpast=0.
            RRnow=0.
            RRtot=0.
            CALL ZERO2(ACPRECpast)
         ENDIF
      ENDIF
C-----------------------------------------------------------------------
      IF(COUPLE_WAVE .AND. MOD(NTSD,INT(DT_WAVE_RATIO)).EQ.1) THEN
      call mpi_comm_size(MPI_COMM_WORLD,npest,ierr)
      npestm1=npest-1
      IF(MYPE.EQ.4) THEN  !4for2x3,4x4, 2for2x2
         call mpi_recv
     & (idtpro,12,MPI_CHARACTER,npestm1,tagi,MPI_COMM_WORLD,status,ierr)
         call mpi_recv
     &  (betag,ngx*ngy,MPI_REAL,npestm1,tagw,MPI_COMM_WORLD,status,ierr)
         call mpi_recv
     &  (z0g,ngx*ngy,MPI_REAL,npestm1,tagz,MPI_COMM_WORLD,status,ierr)
         call mpi_recv
     &  (swhg,ngx*ngy,MPI_REAL,npestm1,tagh,MPI_COMM_WORLD,status,ierr)
         call mpi_recv
     &  (wdir,ngx*ngy,MPI_REAL,npestm1,tagd,MPI_COMM_WORLD,status,ierr)
         call mpi_recv
     &  (wmask,ngx*ngy,MPI_REAL,npestm1,tagm,MPI_COMM_WORLD,status,ierr)
         call mpi_recv
     &  (ageg,ngx*ngy,MPI_REAL,npestm1,taga,MPI_COMM_WORLD,status,ierr)
         call mpi_recv
     &  (ustrg,ngx*ngy,MPI_REAL,npestm1,tags,MPI_COMM_WORLD,status,ierr)
!         CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
!      ENDIF
!      IF(MYPE.EQ.0)
!      write(6,*)'====================================================='
!      write(6,*)'SUBROUTINE TURBL recv'
!      write(6,*)'ntsd=',ntsd,' wam date/time=',idtpro
!      write(6,*)'betag(ngx/2,ngy/2)=',betag(ngx/2,ngy/2)
!      write(6,*)'z0g(ngx/2,ngy/2)=',z0g(ngx/2,ngy/2)
!      write(6,*)'swhg(ngx/2,ngy/2)=',swhg(ngx/2,ngy/2)
!      write(6,*)'====================================================='
!      write(6,*)'TURBL recv TIMESTEP betag,z0g,swhg(ngx/2,ngy/2),
!     &           wam date/time',
!     &           ntsd,betag(ngx/2,ngy/2),z0g(ngx/2,ngy/2),
!     &           swhg(ngx/2,ngy/2),idtpro
      aminu=100000.
      amaxu=-100000.
      aminz=100000.
      amaxz=-100000.
      aminh=100000.
      amaxh=-100000.
      amind=100000.
      amaxd=-100000.
      amins=100000.
      amaxs=-100000.
      do i=1,ngx
      do j=1,ngy
         if (betag(i,j).lt.1.E+3) then
            if (betag(i,j).lt.aminu) aminu=betag(i,j)
            if (betag(i,j).gt.amaxu) amaxu=betag(i,j)
         endif
         if (z0g(i,j).lt.1.E+3) then
            if (z0g(i,j).lt.aminz) aminz=z0g(i,j)
            if (z0g(i,j).gt.amaxz) amaxz=z0g(i,j)
         endif
         if (swhg(i,j).ge.1.E-5 .and. swhg(i,j).lt.1.E+3) then
            if (swhg(i,j).lt.aminh) aminh=swhg(i,j)
            if (swhg(i,j).gt.amaxh) amaxh=swhg(i,j)
         endif
         if (wdir(i,j).lt.amind) amind=wdir(i,j)
         if (wdir(i,j).gt.amaxd) amaxd=wdir(i,j)
         if (ustrg(i,j).lt.amins) amins=ustrg(i,j)
         if (ustrg(i,j).gt.amaxs) amaxs=ustrg(i,j)
      enddo
      enddo
!      write(6,*)'TURBL recv betag min-max ',aminu,amaxu
!      write(6,*)'TURBL recv z0g min-max ',aminz,amaxz
!      write(6,*)'TURBL recv swhg min-max ',aminh,amaxh
      write(6,*)'TURBL recv ustrg min-max ',amins,amaxs
!      ENDIF
C
      CALL ll2imjm(betag,wavepar_global,1.E-5,1.E+3,0.018)
      CALL ll2imjm(z0g,z0g_global,1.E-5,1.E+3,1.E-5)
      CALL ll2imjm(swhg,swhg_global,1.E-5,1.E+3,1.E-3)
      CALL ll2imjm(wdir,wdir_global,1.E-5,1.E+3,1.E-5)
      CALL ll2imjm(wmask,wmask_global,1.E-5,1.E+3,1.)
      CALL ll2imjm(ageg,ageg_global,1.E-5,1.E+3,1.E-5)
      CALL ll2imjm(ustrg,ustrg_global,1.E-5,1.E+3,1.E-5)
C
!      IF(MYPE.EQ.0) THEN
      aminb=100000.
      amaxb=-100000.
      aminc=100000.
      amaxc=-100000.
      amind=100000.
      amaxd=-100000.
      aminr=100000.
      amaxr=-100000.
      amins=100000.
      amaxs=-100000.
      do i=1,im
      do j=1,jm
         if (wavepar_global(i,j).lt.1.E+3) then
            if (wavepar_global(i,j).lt.aminb) aminb=wavepar_global(i,j)
            if (wavepar_global(i,j).gt.amaxb) amaxb=wavepar_global(i,j)
         endif
         if (z0g_global(i,j).lt.1.E+3) then
            if (z0g_global(i,j).lt.aminc) aminc=z0g_global(i,j)
            if (z0g_global(i,j).gt.amaxc) amaxc=z0g_global(i,j)
         endif
         if (swhg_global(i,j).gt.1.E-5 .and. swhg_global(i,j).lt.1.E+3)
     &   then
            if (swhg_global(i,j).lt.amind) amind=swhg_global(i,j)
            if (swhg_global(i,j).gt.amaxd) amaxd=swhg_global(i,j)
         endif
         if (wdir_global(i,j).lt.aminr) aminr=wdir_global(i,j)
         if (wdir_global(i,j).gt.amaxr) amaxr=wdir_global(i,j)
         if (wmask_global(i,j).ge.0.5) wmask_global(i,j)=1.   !WAM sea
         if (wmask_global(i,j).lt.0.5) wmask_global(i,j)=0.
         if (ustrg_global(i,j).lt.amins) amins=ustrg_global(i,j)
         if (ustrg_global(i,j).gt.amaxs) amaxs=ustrg_global(i,j)
      enddo
      enddo
!      write(6,*)'TURBL wavepar_global ntsd,min-max ',ntsd,aminb,amaxb
      write(6,*)'====================================================='
      write(6,*)'SUBROUTINE TURBL recv'
      write(6,*)'ntsd=',ntsd,' wam date/time=',idtpro
      write(6,*)'betag min-(ngx/2,ngy/2)-max=',
     &           aminu,betag(ngx/2,ngy/2),amaxu
      write(6,*)'betag(wavepar_global) min-(im/2,jm/2)-max=', 
     &           aminb,wavepar_global(im/2,jm/2),amaxb
      write(6,*)
      write(6,*)'z0g min-(ngx/2,ngy/2)-max=',
     &           aminz,z0g(ngx/2,ngy/2),amaxz
      write(6,*)'z0g(z0g_global) min-(im/2,jm/2)-max=', 
     &           aminc,z0g_global(im/2,jm/2),amaxc
      write(6,*)
      write(6,*)'swhg min-(ngx/2,ngy/2)-max=',
     &           aminh,swhg(ngx/2,ngy/2),amaxh
      write(6,*)'swhg(swhg_global) min-(im/2,jm/2)-max=', 
     &           amind,swhg_global(im/2,jm/2),amaxd
      write(6,*)'wdir(wdir_global) min-(im/2,jm/2)-max=', 
     &           aminr,wdir_global(im/2,jm/2),amaxr
      write(6,*)'ustrg(ustrg_global) min-(im/2,jm/2)-max=', 
     &           amins,ustrg_global(im/2,jm/2),amaxs
      write(6,*)'====================================================='
      TEMP1=wavepar_global
      TEMP2=z0g_global
      TEMP3=swhg_global
      TEMP4=wdir_global
      TEMP5=wmask_global
      TEMP6=ageg_global
      TEMP7=ustrg_global
      ENDIF  !MYPE=2
!
      CALL DSTRBW(TEMP1,wavepar,1,1,1)
      CALL DSTRBW(TEMP2,z0gpar,1,1,1)
      CALL DSTRBW(TEMP3,swhgpar,1,1,1)
      CALL DSTRBW(TEMP4,wdirpar,1,1,1)
      CALL DSTRBW(TEMP5,wmaskpar,1,1,1)
      CALL DSTRBW(TEMP6,agegpar,1,1,1)
      CALL DSTRBW(TEMP7,ustrgpar,1,1,1)
!
      aminw1=100000.
      amaxw1=-100000.
      aminw2=100000.
      amaxw2=-100000.
      aminw3=100000.
      amaxw3=-100000.
      aminw4=100000.
      amaxw4=-100000.
!      ENDIF
      DO J=MYJS_P1,MYJE_P1          ! This line is correct
      DO I=MYIS_P1,MYIE_P1
         if (wavepar(i,j).gt.1.E-5 .and. wavepar(i,j).lt.1.E+3) then
            if (wavepar(i,j).lt.aminw1) aminw1=wavepar(i,j)
            if (wavepar(i,j).gt.amaxw1) amaxw1=wavepar(i,j)
         endif
         if (z0gpar(i,j).gt.1.E-5 .and. z0gpar(i,j).lt.1.E+3) then
            if (z0gpar(i,j).lt.aminw2) aminw2=z0gpar(i,j)
            if (z0gpar(i,j).gt.amaxw2) amaxw2=z0gpar(i,j)
         endif
         if (swhgpar(i,j).gt.1.E-5 .and. swhgpar(i,j).lt.1.E+3) then
            if (swhgpar(i,j).lt.aminw3) aminw3=swhgpar(i,j)
            if (swhgpar(i,j).gt.amaxw3) amaxw3=swhgpar(i,j)
         endif
         if (agegpar(i,j).gt.1.E-5 .and. agegpar(i,j).lt.1.E+3) then
            if (agegpar(i,j).lt.aminw4) aminw4=agegpar(i,j)
            if (agegpar(i,j).gt.amaxw4) amaxw4=agegpar(i,j)
         endif
      ENDDO
      ENDDO
!      write(6,*)'TURBL mype,DIMs ',mype,MYIS_P1,MYIE_P1,MYJS_P1,MYJE_P1
      write(6,*)'SUBROUTINE TURBL recv MYPE:', mype
      write(6,*)'betag(wavepar) MPI distribute min-max=', aminw1,amaxw1
      write(6,*)'z0g(z0gpar) MPI distribute min-max=', aminw2,amaxw2
      write(6,*)'swhg(swhgpar) MPI distribute min-max=', aminw3,amaxw3
      write(6,*)'wage(agegpar) MPI distribute min-max=', aminw4,amaxw4
      write(6,*)'====================================================='
!      write(6,*)'TURBL wavepar mype,ntsd,min-max ',mype,ntsd,aminw,amaxw
!      ENDIF
      ENDIF   !COUPLE_WAVE
C-----------------------------------------------------------------------
C - WAVES 
      CALL ZERO3(AKM,LM1)
      CALL ZERO3(ZINT,LP1)
      CALL ZERO3_T(AKH_T,LM1)
      CALL ZERO3_T(AKM_T,LM1)
      CALL ZERO2(UZ0H)
      CALL ZERO2(VZ0H)
C-----------------------------------------------------------------------
C***
C***  COMPUTE THE HEIGHTS OF THE LAYER INTERFACES AND THE EXNER FUNCTION
C***
!$omp parallel do 
      DO J=MYJS_P1,MYJE_P1          ! This line is correct
c     DO J=MYJS2_P1,MYJE2_P1        ! This line matches operations
      DO I=MYIS_P1,MYIE_P1
        ZINT(I,J,LP1)=EPSZ 
        IF(SIGMA)ZINT(I,J,LP1)=RG*FIS(I,J)
      ENDDO
      ENDDO
C
      DO 10 L=LM,1,-1
!$omp parallel do private (apests,pdsl)
      DO J=MYJS1_P1,MYJE1_P1          ! This line is correct
c     DO J=MYJS2_P1,MYJE2_P1        ! This line matches operations
      DO I=MYIS_P1,MYIE_P1
        PDSL=PD(I,J)*RES(I,J)
        APESTS=PDSL*AETA(L)+PT
C
        ZINT(I,J,L)=ZINT(I,J,L+1)+T(I,J,L)/APESTS
     1             *PDSL*DETA(L)*ROG*(Q(I,J,L)*0.608+1.)
        ZINT(I,J,L)=(ZINT(I,J,L)-DFRLG(L))*HTM(I,J,L)+DFRLG(L)
C
        APE(I,J,L)=(1.E5/APESTS)**CAPA
      ENDDO
      ENDDO
   10 CONTINUE
C-----------------------------------------------------------------------
C***
C***  REMOVE NEGATIVE Q2
C***
!$omp parallel do 
      DO 40 L=1,LM
      DO J=MYJS_P1,MYJE_P1
      DO I=MYIS_P1,MYIE_P1
        Q2(I,J,L)=AMAX1(Q2(I,J,L)*HBM2(I,J),EPSQ2)
      ENDDO
      ENDDO
   40 CONTINUE
C-----------------------------------------------------------------------
!$omp parallel do 
      DO J=MYJS2_P1,MYJE2_P1
      DO I=MYIS_P1,MYIE_P1
        UZ0H(I,J)=(UZ0(I+IHE(J),J)+UZ0(I+IHW(J),J)
     1            +UZ0(I,J+1)+UZ0(I,J-1))*HBM2(I,J)*0.25
        VZ0H(I,J)=(VZ0(I+IHE(J),J)+VZ0(I+IHW(J),J)
     1            +VZ0(I,J+1)+VZ0(I,J-1))*HBM2(I,J)*0.25
      ENDDO
      ENDDO
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C***               PREPARE THE EXCHANGE COEFFICIENTS
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C***
C***  COMPUTE VELOCITY COMPONENTS AT H POINTS
C***
!$omp parallel do private(rwmsk,wmsk)
      DO 60 L=1,LM
C
      DO J=MYJS2_P1,MYJE2_P1
      DO I=MYIS_P1,MYIE_P1
        WMSK=VTM(I+IHE(J),J,L)+VTM(I+IHW(J),J,L)
     1                        +VTM(I,J+1,L)+VTM(I,J-1,L)
        IF(WMSK.GT.0.)THEN
          RWMSK=1./WMSK
          UCOL(I,J,L)=(U(I+IHE(J),J,L)*VTM(I+IHE(J),J,L)
     1                +U(I+IHW(J),J,L)*VTM(I+IHW(J),J,L)
     2                +U(I,J+1,L)*VTM(I,J+1,L)+U(I,J-1,L)*VTM(I,J-1,L))
     3                *RWMSK
          VCOL(I,J,L)=(V(I+IHE(J),J,L)*VTM(I+IHE(J),J,L)
     1                +V(I+IHW(J),J,L)*VTM(I+IHW(J),J,L)
     2                +V(I,J+1,L)*VTM(I,J+1,L)+V(I,J-1,L)*VTM(I,J-1,L))
     3                *RWMSK
        ELSE
          UCOL(I,J,L)=0.
          VCOL(I,J,L)=0.
        ENDIF
      ENDDO
      ENDDO
   60 CONTINUE
C***
C***  FILL TRANSPOSED ARRAYS
C***
!$omp parallel sections
!$omp section
      CALL SGETMO(T,NHRZ,NHRZ,LM,TCOL_T,LM)
!$omp section
      CALL SGETMO(Q,NHRZ,NHRZ,LM,QCOL_T,LM)
!$omp section
      CALL SGETMO(APE,NHRZ,NHRZ,LM,APECOL_T,LM)
!$omp section
      CALL SGETMO(Q2,NHRZ,NHRZ,LM,Q2COL_T,LM)
!$omp section
      CALL SGETMO(ZINT,NHRZ,NHRZ,LP1,ZCOL_T,LP1)
!$omp section
      CALL SGETMO(UCOL,NHRZ,NHRZ,LM,UCOL_T,LM)
!$omp section
      CALL SGETMO(VCOL,NHRZ,NHRZ,LM,VCOL_T,LM)
!$omp end parallel sections
C----------------------------------------------------------------------
C***
C***  FIND THE MIXING LENGTH
C***
!$omp parallel do private(el,gh,gm,hpbl,lmhk,lmhm,lmhp,lpbl)
!$omp& private(ulm,vlm,wstar,zeff)
      amaxus=-100000.
      aminus=100000.
      
      DO 100 J=MYJS2_P1,MYJE2_P1
      DO 100 I=MYIS_P1,MYIE1_P1

      LMHK=LMH(I,J)
      LMHP=LMHK+1
      LMHM=LMHK-1
!      if (SM(I,J).gt.0.5) then
!         print*,"LMH ",i,j,LMH(I,J), LMHM, LMHP
!      endif
C
      CALL MIXLEN(LMHK,LPBL,HPBL,UCOL_T(1,I,J),VCOL_T(1,I,J)
     1,           TCOL_T(1,I,J),QCOL_T(1,I,J),Q2COL_T(1,I,J)
     2,           APECOL_T(1,I,J),ZCOL_T(1,I,J),GM,GH,EL)
!      CALL MIXLEN(LMHM,LPBL,HPBL,UCOL_T(1,I,J),VCOL_T(1,I,J) !coupling
!     1,           TCOL_T(1,I,J),QCOL_T(1,I,J),Q2COL_T(1,I,J)
!     2,           APECOL_T(1,I,J),ZCOL_T(1,I,J),GM,GH,EL)
       LPBL2(I,J)=LPBL  !coupling
       HPBL2(I,J)=HPBL  !coupling
C
C-----------------------------------------------------------------------
C***
C***  SOLVE FOR THE PRODUCTION/DISSIPATION OF 
C***  THE TURBULENT KINETIC ENERGY
C***
C
!	write(6,*) 'working ', I,J
      CALL PRODQ2(LMHK,DTQ2,USTAR(I,J),GM,GH,EL,Q2COL_T(1,I,J))
!      CALL PRODQ2(LMHM,DTQ2,USTAR(I,J),GM,GH,EL,Q2COL_T(1,I,J)) !coupling
C
C-----------------------------------------------------------------------
C***
C***  FIND THE EXCHANGE COEFFICIENTS IN THE FREE ATMOSPHERE
C***
      CALL DIFCOF(LMHK,GM,GH,EL,Q2COL_T(1,I,J)
     1,           ZCOL_T(1,I,J),AKM_T(1,I,J),AKH_T(1,I,J))
!      CALL DIFCOF(LMHM,GM,GH,EL,Q2COL_T(1,I,J)     !coupling
!     1,           ZCOL_T(1,I,J),AKM_T(1,I,J),AKH_T(1,I,J))

      do L=1,LM1
      if (abs(AKM_T(L,I,J)) .le. 2.e10) then
      else
      write(6,*) 'bad akm_t from DIFCOF: ',MYPE,
     +I,J,L,AKM_T(L,I,J)
      endif
      enddo
C-----------------------------------------------------------------------
C***
C***  CARRY OUT THE VERTICAL DIFFUSION OF
C***  TURBULENT KINETIC ENERGY
C***
C
      CALL VDIFQ(LMHK,KTMQ2,DTQ2,Q2COL_T(1,I,J),EL,ZCOL_T(1,I,J))
!      CALL VDIFQ(LMHM,KTMQ2,DTQ2,Q2COL_T(1,I,J),EL,ZCOL_T(1,I,J)) !coupling
C-----------------------------------------------------------------------
C***
C***  FIND THE Z0 EFFECTIVE
C***
      ZEFF(1)=ZEFFIJ(I,J,1)
      ZEFF(2)=ZEFFIJ(I,J,2)
      ZEFF(3)=ZEFFIJ(I,J,3)
      ZEFF(4)=ZEFFIJ(I,J,4)
C-----------------------------------------------------------------------
C***
C***  FIND THE SURFACE EXCHANGE COEFFICIENTS
C***
      ULM=UCOL(I,J,LMHK)
      VLM=VCOL(I,J,LMHK)
!      ULM=UCOL(I,J,LMHM)  !ULM-1 for coupling
!      VLM=VCOL(I,J,LMHM)  !ULM-1 for coupling

c --- ioannis ------------------
c --- store the necessary fields 
c --- in TEMPORARY variables 
c
      TMPUZ0H=UZ0H(I,J)
      TMPVZ0H=VZ0H(I,J)
      TMPTHZ0=THZ0(I,J)
      TMPQZ0=QZ0(I,J)
      TMPUSTAR=USTAR(I,J)
      TMPWSTAR=WSTAR
      TMPZ0=Z0(I,J)

      USTR=USTAR(i,j)*SM(i,j)
      if (USTR.lt.aminus) aminus=USTR
      if (USTR.gt.amaxus) amaxus=USTR
!      if (SM(I,J).EQ.1. .and. USTAR(I,J).GT.0.6)
!     &    print*,'TURBL USTAR US_WAM ',USTAR(I,J), US_WAM
C - WAVES 
      IF (COUPLE_WAVE) THEN
         WAVE_PAR = wavepar(I,J)
         Z0_WAM = z0gpar(I,J)
         IF (SM(i,j).gt.0.5 .and. wmaskpar(i,j).gt.0.5) THEN
            US_WAM = ustrgpar(I,J)
            WAGE = amax1(agegpar(I,J), 1.E-3)
!            print*,' TURBL US_WAM, USTR, WAGE ',US_WAM, USTR, WAGE
         ELSE
            US_WAM = USTR
         END IF
         IF (CTRL) THEN
            WAVE_PAR = CHARNOCK
            IF(MYPE.EQ.0 .AND. I.EQ.MYIS_P1 .AND. J.EQ.MYJS2_P1 )
     &         print*,'OFFLINE-Constant Charnock ',CTRL,CHARNOCK
         ENDIF
! get precipitation parameters for the revised roughness length
         IF (RAINEFF) THEN
!            RRT = RAINR(I,J) !from the GSMCOLUMN
            VVD = TERMV(I,J)
            DDR = DRODIA(I,J)
            RRnow=ACPREC(I,J)
            RRpast=ACPRECpast(I,J)
            RRtot=(RRnow-RRpast)/DTQ2
            RRT = RRtot
            IF (RRT.GT.5.E-3)print*,'RRT, RRtot ',RRT,RRtot,RRnow,RRpast
!            IF (RRT.GT.1.E-3) print*,'TURBL RR,VRAIN,IDR ', RRT, VVD,DDR
         ENDIF
      ELSE
         WAVE_PAR = CHARNOCK
            IF(MYPE.EQ.0 .AND. I.EQ.MYIS_P1 .AND. J.EQ.MYJS2_P1 )
     &         print*,'Constant Charnock ',CTRL,CHARNOCK
      ENDIF
! get precipitation parameters for the revised roughness length
      RRT = RAINR(I,J)
      VVD = TERMV(I,J)
      DDR = DRODIA(I,J)
      RRnow=ACPREC(I,J)
      RRpast=ACPRECpast(I,J)
      RRtot=(RRnow-RRpast)/DTQ2
      RRT = RRtot
      IF (RRT.GT.5.E-3)print*,'RRT, RRtot ', RRT, RRtot, RRnow, RRpast
!      IF (RRT.GT.1.E-3) print*,'TURBL RR,VRAIN,IDR ', RRT, VVD,DDR
C - WAVES 
      TMPAKMS=AKMS(I,J)
      TMPAKHS=AKHS(I,J)
      TMPCT=CT(I,J)
      TMPTH2M=TSHLTR(I,J)
      TMPTH10=TH10(I,J)
      TMPQ2M=QSHLTR(I,J)
      TMPQ10=Q10(I,J)
c -------------------------------

c ----ioannis------
c      if (MYIE_P1.ne.MYIE1_P1) then
c      write(6,*)'ioannis-new10mwind NONEQUAL',MYPE,MYIE_P1,MYIE1_P1
c      endif

c      if ((mype.eq.13).and.(i.eq.9).and.(j.eq.47)) then
c         write(6,*)'ioannis-10mwind=',MYIS_P1,MYIE_P1,MYIE1_P1,
c     $        MYJS2_P1,MYJE2_P1

c      write(6,*)'ioannis...10mwind=',NTSD,MYPE,I,J,FIS(i,j),Z0(i,j)
c     $,ZEFF(1),ZEFF(2),ZEFF(3),ZEFF(4),ULM,VLM,LMHK
c     $,U10(i,j),V10(i,j),I+IHE(J),I+IHW(J)
c     $,VTM(I+IHE(J),J,LMHK),VTM(I+IHW(J),J,LMHK),VTM(I,J+1,LMHK)
c     $,VTM(I,J-1,LMHK),' ...ENDioannis'
c      endif
c -----------------

C
      CALL SFCDIF(LMHK,SM(I,J),THS(I,J),QS(I,J)
     1,           UZ0H(I,J),VZ0H(I,J),THZ0(I,J),QZ0(I,J)
     2,           USTAR(I,J),WSTAR
     3,           Z0(I,J),ZEFF,AKMS(I,J),AKHS(I,J),HPBL,CT(I,J)
     4,           U10(I,J),V10(I,J),TSHLTR(I,J),TH10(I,J)
     5,           QSHLTR(I,J),Q10(I,J)
     6,           ULM,VLM,TCOL_T(1,I,J),QCOL_T(1,I,J)
     7,           APECOL_T(1,I,J),ZCOL_T(1,I,J),PD(I,J),PT
     8,           T(I,J,LMHK),WAVE_PAR,Z0_WAM,US_WAM,WAGE
     9,           RRT, VVD, DDR, RAINEFF)
!      CALL SFCDIF(LMHM,SM(I,J),THS(I,J),QS(I,J)   !coupling
!     1,           UZ0H(I,J),VZ0H(I,J),THZ0(I,J),QZ0(I,J)
!     2,           USTAR(I,J),WSTAR
!     3,           Z0(I,J),ZEFF,AKMS(I,J),AKHS(I,J),HPBL,CT(I,J)
!     4,           U10(I,J),V10(I,J),TSHLTR(I,J),TH10(I,J)
!     5,           QSHLTR(I,J),Q10(I,J)
!     6,           ULM,VLM,TCOL_T(1,I,J),QCOL_T(1,I,J)
!     7,           APECOL_T(1,I,J),ZCOL_T(1,I,J),PD(I,J),PT
!     8,           T(I,J,LMHM),WAVE_PAR,Z0_WAM)
C
c --- ioannis ------
c --- we avoid checking whether ULM or VLM are equal to zero
c --- because in some cases the wind speed may be equal to 
c --- zero not because of the interpolation method but due
c --- to extremely light winds
c
      ULM=UCOL(I,J,LMHK)
      VLM=VCOL(I,J,LMHK)
!      ULM=UCOL(I,J,LMHM)  !ULM-1 for coupling
!      VLM=VCOL(I,J,LMHM)  !ULM-1 for coupling
c     off for coupling
      DO L=LMHK,LMHK-8,-1
        IF((DFL(L)-DFL(LMHK)).lt.7357.5) THEN  !  7357.5 is 9.81*750m
        IF ((VTM(I+IHE(J),J,L)+VTM(I+IHW(J),J,L)
     $      +VTM(I,J+1,L)+VTM(I,J-1,L)).ne.0.0) THEN         
              ULM=UCOL(I,J,L)
              VLM=VCOL(I,J,L)
              goto 101          
        ENDIF
        ENDIF
      ENDDO
  101  continue
c ------------------

      CALL SFCDIF(LMHK,SM(I,J),THS(I,J),QS(I,J)
     1,           TMPUZ0H,TMPVZ0H,TMPTHZ0,TMPQZ0
     2,           TMPUSTAR,TMPWSTAR
     3,           TMPZ0,ZEFF,TMPAKMS,TMPAKHS,HPBL,TMPCT
     4,           U10(I,J),V10(I,J),TMPTH2M,TMPTH10
     5,           TMPQ2M,TMPQ10
     6,           ULM,VLM,TCOL_T(1,I,J),QCOL_T(1,I,J)
     7,           APECOL_T(1,I,J),ZCOL_T(1,I,J),PD(I,J),PT
     8,           T(I,J,LMHK),WAVE_PAR,Z0_WAM,US_WAM,WAGE
     9,           RRT, VVD, DDR, RAINEFF)
!      CALL SFCDIF(LMHM,SM(I,J),THS(I,J),QS(I,J)   !coupling
!     1,           UZ0H(I,J),VZ0H(I,J),THZ0(I,J),QZ0(I,J)
!     2,           USTAR(I,J),WSTAR
!     3,           Z0(I,J),ZEFF,AKMS(I,J),AKHS(I,J),HPBL,CT(I,J)
!     4,           U10(I,J),V10(I,J),TSHLTR(I,J),TH10(I,J)
!     5,           QSHLTR(I,J),Q10(I,J)
!     6,           ULM,VLM,TCOL_T(1,I,J),QCOL_T(1,I,J)
!     7,           APECOL_T(1,I,J),ZCOL_T(1,I,J),PD(I,J),PT
!     8,           T(I,J,LMHM),WAVE_PAR,Z0_WAM)

  100 CONTINUE
      ACPRECpast=ACPREC
C------------------------------------------------------------------------
C***
C***  FILL STANDARD ARRAYS FROM TRANSPOSED ARRAYS
C***
!$omp parallel sections
!$omp section
       CALL SGETMO(Q2COL_T,LM,LM,NHRZ,Q2,NHRZ)
!$omp section
       CALL SGETMO(AKH_T,LM1,LM1,NHRZ,AKH,NHRZ)
!$omp section
       CALL SGETMO(AKM_T,LM1,LM1,NHRZ,AKM,NHRZ)
!$omp end parallel sections
C-----------------------------------------------------------------------
C***
C***  UNCOMPUTED LOCATIONS MUST BE FILLED IN FOR THE POST-PROCESSOR
C***
      IIM=IM-MY_IS_GLB+1
      JJM=JM-MY_JS_GLB+1
C
C***  EASTERN GLOBAL BOUNDARY
C
      IF(MY_IE_GLB.EQ.IM)THEN
        DO J=1,JM
          IF(J.GE.MY_JS_GLB.AND.J.LE.MY_JE_GLB)THEN
            JJ=J-MY_JS_GLB+1
            TH10(IIM,JJ)=TH10(IIM-1,JJ)
            Q10(IIM,JJ)=Q10(IIM-1,JJ)
            U10(IIM,JJ)=U10(IIM-1,JJ)
            V10(IIM,JJ)=V10(IIM-1,JJ)
            TSHLTR(IIM,JJ)=TSHLTR(IIM-1,JJ)
            QSHLTR(IIM,JJ)=QSHLTR(IIM-1,JJ)
          ENDIF
        ENDDO
      ENDIF
C
C***  SOUTHERN GLOBAL BOUNDARY
C
      IF(MY_JS_GLB.EQ.1)THEN
        DO J=1,2
        DO I=1,IM
          IF(I.GE.MY_IS_GLB.AND.I.LE.MY_IE_GLB)THEN
            II=I-MY_IS_GLB+1
            TH10(II,J)=TH10(II,3)
            Q10(II,J)=Q10(II,3)
            U10(II,J)=U10(II,3)
            V10(II,J)=V10(II,3)
            TSHLTR(II,J)=TSHLTR(II,3)
            QSHLTR(II,J)=QSHLTR(II,3)
          ENDIF
        ENDDO
        ENDDO
      ENDIF
C
C***  NORTHERN GLOBAL BOUNDARY
C
      IF(MY_JE_GLB.EQ.JM)THEN
        DO J=JM-1,JM
        JJ=J-MY_JS_GLB+1
        DO I=1,IM
          IF(I.GE.MY_IS_GLB.AND.I.LE.MY_IE_GLB)THEN
            II=I-MY_IS_GLB+1
            TH10(II,JJ)=TH10(II,JJM-2)
            Q10(II,JJ)=Q10(II,JJM-2)
            U10(II,JJ)=U10(II,JJM-2)
            V10(II,JJ)=V10(II,JJM-2)
            TSHLTR(II,JJ)=TSHLTR(II,JJM-2)
            QSHLTR(II,JJ)=QSHLTR(II,JJM-2)
          ENDIF
        ENDDO
        ENDDO
      ENDIF
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      btim=timef()
      CALL EXCH(UZ0H,1,VZ0H,1,1,1)
      exch_tim=exch_tim+timef()-btim
C***
C***  AVERAGE UZ0 AND VZ0 BACK TO V POINTS
C***
!$omp parallel do 
      DO 125 J=MYJS2,MYJE2
      DO 125 I=MYIS,MYIE
      UZ0(I,J)=(UZ0H(I+IVE(J),J)*HBM2(I+IVE(J),J)
     1         +UZ0H(I+IVW(J),J)*HBM2(I+IVW(J),J)
     2         +UZ0H(I,J+1)*HBM2(I,J+1)+UZ0H(I,J-1)*HBM2(I,J-1))*0.25
      VZ0(I,J)=(VZ0H(I+IVE(J),J)*HBM2(I+IVE(J),J)
     1         +VZ0H(I+IVW(J),J)*HBM2(I+IVW(J),J)
     2         +VZ0H(I,J+1)*HBM2(I,J+1)+VZ0H(I,J-1)*HBM2(I,J-1))*0.25
  125 CONTINUE
C-----------------------------------------------------------------------
C***
C***  EXECUTE THE GROUND PROCESSES
C***
      CALL SURFCE(APE(IDIM1,JDIM1,1),ZINT(IDIM1,JDIM1,1)
     1,           CKLQ(IDIM1,JDIM1))
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C***                 EXECUTE THE VERTICAL EXCHANGE
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      btim=timef()

      CALL EXCH(AKMS,1,AKM,LM1,ZINT,LP1,1,1)
      exch_tim=exch_tim+timef()-btim
C
!$omp parallel do
      DO L=1,LM1
        DO J=MYJS2,MYJE2
        DO I=MYIS,MYIE
          AKMCOL(I,J,L)=(AKM(I+IVE(J),J,L)*HBM2(I+IVE(J),J)
     1        +AKM(I+IVW(J),J,L)*HBM2(I+IVW(J),J)
     2        +AKM(I,J+1,L)*HBM2(I,J+1)+AKM(I,J-1,L)*HBM2(I,J-1))
     3        *VTM(I,J,L)*VBM2(I,J)*0.25
C
      if (abs(AKMCOL(I,J,L)) .le. 2.e10) then
      else
      write(6,*) 'bad AKMCOL: ', I,J,L,AKMCOL(I,J,L)
      write(6,*) 'Js needed for AKM: ', MYJS2-1,MYJE2+1
      write(6,*) 'AKM values: '
      write(6,*) AKM(I+IVE(J),J,L),AKM(I+IVW(J),J,L),
     +AKM(I,J+1,L),AKM(I,J-1,L)
      write(6,*) 'HBM2 values: '
      write(6,*) HBM2(I+IVE(J),J),HBM2(I+IVW(J),J),
     +HBM2(I,J+1),HBM2(I,J-1)
      write(6,*) 'VTM,VBM2: ', VTM(I,J,L),VBM2(I,J)
      endif
          AKHCOL(I,J,L)=AKH(I,J,L)*HTM(I,J,L)*HBM2(I,J)
        ENDDO
        ENDDO
      ENDDO
C
!$omp parallel do
      DO J=MYJS2,MYJE2
      DO I=MYIS,MYIE
        THZ0(I,J)=(1.-SM(I,J))*THS(I,J)+SM(I,J)*THZ0(I,J)
        QZ0 (I,J)=(1.-SM(I,J))*QS (I,J)+SM(I,J)*QZ0 (I,J)
        AKMSV(I,J)=(AKMS(I+IVE(J),J)*HBM2(I+IVE(J),J)
     1             +AKMS(I+IVW(J),J)*HBM2(I+IVW(J),J)
     2             +AKMS(I,J+1)*HBM2(I,J+1)+AKMS(I,J-1)*HBM2(I,J-1))
     3             *VBM2(I,J)*0.25
      ENDDO
      ENDDO
C
!$omp parallel do
      DO L=1,LP1
        DO J=MYJS2,MYJE2
        DO I=MYIS,MYIE
          ZCOL(I,J,L)=0.25*(ZINT(I+IVE(J),J,L)+ZINT(I+IVW(J),J,L)
     1                     +ZINT(I,J+1,L)+ZINT(I,J-1,L))
        ENDDO
        ENDDO
      ENDDO
C***
C***  TRANSPOSE ARRAYS
C***
!$omp parallel sections
!$omp section
      CALL SGETMO(ZCOL,NHRZ,NHRZ,LP1,ZCOL_T2,LP1)
!$omp section
      CALL SGETMO(U,NHRZ,NHRZ,LM,UCOL_T,LM)
!$omp section
      CALL SGETMO(V,NHRZ,NHRZ,LM,VCOL_T,LM)
!$omp section
      CALL SGETMO(AKHCOL,NHRZ,NHRZ,LM1,AKH_T,LM1)
!$omp section
      CALL SGETMO(AKMCOL,NHRZ,NHRZ,LM1,AKM_T,LM1)
!$omp end parallel sections
C-----------------------------------------------------------------------
!$omp parallel do private(lmhk,lmvk)
      DO 200 J=MYJS2,MYJE2
      DO 200 I=MYIS,MYIE1
C
      LMHK=LMH(I,J)
      LMVK=LMV(I,J)
!      LMHM=LMHK-1   !coupling
!      LMVM=LMVK-1   !coupling
C***
C***  CARRY OUT THE VERTICAL DIFFUSION OF
C***  TEMPERATURE AND WATER VAPOR
C***
      CALL VDIFH(LMHK,KTMQ2,DTQ2,THZ0(I,J),QZ0(I,J),AKHS(I,J)
     1,          CT(I,J),CKLQ(I,J)
     2,          TCOL_T(1,I,J),QCOL_T(1,I,J),AKH_T(1,I,J)
     3,          APECOL_T(1,I,J),ZCOL_T(1,I,J))
!      CALL VDIFH(LMHM,KTMQ2,DTQ2,THZ0(I,J),QZ0(I,J),AKHS(I,J) !coupling
!     1,          CT(I,J),CKLQ(I,J)
!     2,          TCOL_T(1,I,J),QCOL_T(1,I,J),AKH_T(1,I,J)
!     3,          APECOL_T(1,I,J),ZCOL_T(1,I,J))
C     
C-----------------------------------------------------------------------
C***
C***  CARRY OUT THE VERTICAL DIFFUSION OF
C***  VELOCITY COMPONENTS
C***
       CALL VDIFV(LMVK,KTMQ2,DTQ2,UZ0(I,J),VZ0(I,J)
     1,          AKMSV(I,J),UCOL_T(1,I,J),VCOL_T(1,I,J)
     2,          AKM_T(1,I,J),ZCOL_T2(1,I,J))
!       CALL VDIFV(LMVM,KTMQ2,DTQ2,UZ0(I,J),VZ0(I,J)   !coupling
!     1,          AKMSV(I,J),UCOL_T(1,I,J),VCOL_T(1,I,J)
!     2,          AKM_T(1,I,J),ZCOL_T2(1,I,J))
C
C-----------------------------------------------------------------------
  200 CONTINUE
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C***
C***  FILL STANDARD ARRAYS FROM TRANSPOSED ARRAYS
C***
!$omp parallel sections
!$omp section
      CALL SGETMO(QCOL_t,LM,LM,NHRZ,Q,NHRZ)
!$omp section
      CALL SGETMO(TCOL_t,LM,LM,NHRZ,T,NHRZ)
!$omp section
      CALL SGETMO(UCOL_t,LM,LM,NHRZ,U,NHRZ)
!$omp section
      CALL SGETMO(VCOL_t,LM,LM,NHRZ,V,NHRZ)
!$omp end parallel sections
C
C***
C***  REMOVE NEGATIVE Q2
C***
!$omp parallel do 
      DO L=1,LM
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          Q2(I,J,L)=AMAX1(Q2(I,J,L)*HBM2(I,J),EPSQ2)
        ENDDO
        ENDDO
      ENDDO
C
!      print*,'TURBL min-(im/2,jm/2),max USTAR: ',
!     &        aminus,USTAR(im/2,jm/2),amaxus
C----------------------------------------------------------------------
!$omp parallel do
!cpk Neutral winds
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
!         DIRW10=ATAN2(V10(I,J),U10(I,J))
         DIRW10=57.2957795 * ATAN2(U10(I,J),V10(I,J)) + 180.
         WSpeed10 = 2.5*USTAR(I,J)*ALOG(10./Z0(I,J)) ! 10 m
         Uneu(I,J) = -WSpeed10 * SIN(DIRW10*.01745239)
         Vneu(I,J) = -WSpeed10 * COS(DIRW10*.01745239)
!         Uneu(I,J) = WSpeed10*COS(DIRW10)
!         Vneu(I,J) = WSpeed10*SIN(DIRW10)
         Uneu(I,J) = UCOL(I,J,LM)
         Vneu(I,J) = VCOL(I,J,LM)
      ENDDO
      ENDDO
C----------------------------------------------------------------------
                                 RETURN
                                 END

