        subroutine eta2p

       USE mod_vrbls
       USE mod_extra
       USE mod_masks
       USE mod_pvrbls
       USE mod_graph

	implicit none
       include 'param_o.h'
!!JNT       include 'extra.comm'
       include 'mapot.comm'
       include 'const.h'
       include 'dynam.comm'
       include 'params'

       integer,parameter::im_jm_nm=(im+2)*(jm+2)*nm
       real,parameter::g=9.80616,RD1=287.04,GAMMA=6.5E-3,   &
                       RGAMOG=RD1*GAMMA/G,PQ0=379.90516,A2=17.2693882  &
                       , A3=273.16,A4=35.86
       integer::nhold,lp,i,j,n,lsl,nn,l,lmb,il
!!JNT       real,dimension(0:im+1,0:jm+1,nm)::tsl,qsl,fsl,q2sl
!!JNT       real,dimension(0:im+1,0:jm+1,nm,lm)::iw
!!JNT       integer,dimension(0:im+1,0:jm+1,nm)::nl1x
!!JNT       integer,dimension(im_jm_nm)::ihold,jhold,nnhold

       real, allocatable, dimension(:,:,:)::tsl,qsl,fsl,q2sl
       real, allocatable, dimension(:,:,:,:)::iw
       integer, allocatable, dimension(:,:,:)::nl1x
       integer, allocatable, dimension(:)::ihold,jhold,nnhold
       real::pnl1,B,fac,pu,tu,tabv,pl,tl,tabo,ahf,tblo,petau,alpetu, &
             petal,alpetl,alpet2,fact,alpet1,trf,zu,qu,ai,bi,tmt0,tmt15, &
	     qw,qi,qint,qsat,utim,climit,lml,hh,tkl,qkl,cwmkl,pp,   &
             u00kl,fiq,iwu,qabv,rhu,bq,ahfq,ql,iwl,rhl,qblo,q2a,bq2,   &
             q2b,ahfq2
!!JNT       real,dimension(0:im+1,0:jm+1,nm)::alpetux,alpet2x,usl,vsl
       real, allocatable, dimension(:,:,:)::alpetux,alpet2x,usl,vsl

!        print *,'u,v,',u(10,10,6,6),v(10,10,6,6)
!        print *,'us,vs,',u(10,10,6,6)*u2us(10,10,6)+v2us(10,10,6)*v(10,10,6,6)
!        print *,'us,vs,',u(10,10,6,6)*u2vs(10,10,6)+v2vs(10,10,6)*v(10,10,6,6)
	
          UTIM=1.
          CLIMIT =1.0E-20

          ALLOCATE(tsl(0:im+1,0:jm+1,nm))
          ALLOCATE(qsl(0:im+1,0:jm+1,nm))
          ALLOCATE(fsl(0:im+1,0:jm+1,nm))
          ALLOCATE(q2sl(0:im+1,0:jm+1,nm))
          ALLOCATE(iw(0:im+1,0:jm+1,nm,lm))

          ALLOCATE(nl1x(0:im+1,0:jm+1,nm))

          ALLOCATE(ihold(im_jm_nm))
          ALLOCATE(jhold(im_jm_nm))
          ALLOCATE(nnhold(im_jm_nm))

          ALLOCATE(alpetux(0:im+1,0:jm+1,nm))
          ALLOCATE(alpet2x(0:im+1,0:jm+1,nm))
          ALLOCATE(usl(0:im+1,0:jm+1,nm))
          ALLOCATE(vsl(0:im+1,0:jm+1,nm))

              do n=1,nm
              DO J=1,jm
              DO I=1,IM
                IW(I,J,n,1)=0.
              ENDDO
              ENDDO
              ENDDO
!C
          DO L=2,LM
            do n=1,nm
            DO J=1,jm
            DO I=1,IM
              LML=LM-LMH(I,J,n)
              HH=HTM(I,J,n,L)
              TKL=T(I,J,n,L)
              QKL=Q(I,J,n,L)
              CWMKL=W(I,J,n,L)
              TMT0=(TKL-273.16)*HH
              TMT15=AMIN1(TMT0,-15.)*HH
              PP=PDSL(I,J,n)*AETA(L)+PT
              QW=HH*PQ0/PP*EXP(HH*A2*(TKL-A3)/(TKL-A4))
              QI=QW*(1.+0.01*AMIN1(TMT0,0.))
              U00KL=U00(I,J,n)+UL(L+LML)*(0.95-U00(I,J,n))*UTIM
!C
              IF(TMT0.LT.-15.0)THEN
                FIQ=QKL-U00KL*QI
                IF(FIQ.GT.0..OR.CWMKL.GT.CLIMIT) THEN
                  IW(I,J,n,L)=1.
                ELSE
                  IW(I,J,n,L)=0.
              ENDIF
!C
              IF(TMT0.GE.0.0)IW(I,J,n,L)=0.
              IF(TMT0.LT.0.0.AND.TMT0.GE.-15.0)THEN
                IW(I,J,n,L)=0.
                IF(IW(I,J,n,L-1).EQ.1.0.AND.CWMKL.GT.CLIMIT)IW(I,J,n,L)=1.
              ENDIF
            endif
           enddo
           enddo
           enddo
          enddo
       open(1,file='tmpfld.dat',form='unformatted')
        lsl=lsm
        DO 310 LP=1,LSL
        NHOLD=0

        TRF=2.*ALSL(LP)

        do 125 n=1,nm
        DO 125 J=1,jm
        DO 125 I=1,IM
!
        TSL(I,J,n)=-1.E6
        QSL(I,J,n)=-1.E6
        FSL(I,J,n)=-1.E6
!
!***  LOCATE VERTICAL INDEX OF MODEL INTERFACE JUST BELOW
!***  THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING.
!
        DO 115 L=2,LM
!	  if(i.eq.25.and.j.eq.22.and.n.eq.3)then
!	     print *,l,alpint(i,j,n,l),alsl(lp),alpint(i,j,n,lp1)
!           endif
        IF(ALPINT(I,J,n,L).GE.ALSL(LP))THEN
          NL1X(I,J,n)=L
          NHOLD=NHOLD+1
          IHOLD(NHOLD)=I
          JHOLD(NHOLD)=J
	  nnhold(nhold)=n
          GO TO 125
        ENDIF
  115   CONTINUE
        nl1x(i,j,n)=l
        NHOLD=NHOLD+1
        IHOLD(NHOLD)=I
        JHOLD(NHOLD)=J
        nnhold(nhold)=n
!
  125   CONTINUE

        DO 220 NN=1,NHOLD
        I=IHOLD(NN)
        J=JHOLD(NN)
        n=nnHOLD(NN)
        PNL1=PINT(I,J,n,NL1X(I,J,n))
        IF(NL1X(I,J,n).EQ.1)THEN
!---------------------------------------------------------------------
!***  EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
!---------------------------------------------------------------------
!
        PU=PINT(I,J,n,2)
        ZU=ZINT(I,J,n,2)
        TU=0.5*(T(I,J,n,1)+T(I,J,n,2))
        QU=0.50*(Q(I,J,n,1)+Q(I,J,n,2))
	IWU=0.5*(IW(I,J,n,1)+IW(I,J,n,2))

        TABV=TU*(SPL(LP)/PU)**RGAMOG

              TMT0=TU-273.16
              TMT15=AMIN1(TMT0,-15.)
              AI=0.008855
              BI=1.
              IF(TMT0.LT.-20.)THEN
                AI=0.007225
                BI=0.9674
              ENDIF
              QW=PQ0/PU*EXP(A2*(TU-A3)/(TU-A4))
              QI=QW*(BI+AI*AMIN1(TMT0,0.))
              QINT=QW*(1.-0.00032*TMT15*(TMT15+15.))
              IF(TMT0.LT.-15.)THEN
                  QSAT=QI
              ELSEIF(TMT0.GE.0.)THEN
                  QSAT=QINT
              ELSE
                IF(IWU.GT.0.0) THEN
                  QSAT=QI
                ELSE
                  QSAT=QINT
                ENDIF
              ENDIF
!CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT
!C             DELETE THIS LINE TO SWITCH BACK TO RH VS ICE
              QSAT=QW
!CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT
              RHU =QU/QSAT
!C
              IF(RHU.GT.1)THEN
                RHU=1
                QU =RHU*QSAT
              ENDIF

              IF(RHU.LT.0.01)THEN
                RHU=0.01
                QU =RHU*QSAT
              ENDIF

              TMT0=TABV-273.16
              TMT15=AMIN1(TMT0,-15.)
              AI=0.008855
              BI=1.
              IF(TMT0.LT.-20.)THEN
                AI=0.007225
                BI=0.9674
              ENDIF
              QW=PQ0/SPL(LP)*EXP(A2*(TABV-A3)/(TABV-A4))
              QI=QW*(BI+AI*AMIN1(TMT0,0.))
              QINT=QW*(1.-0.00032*TMT15*(TMT15+15.))
              IF(TMT0.LT.-15.)THEN
                  QSAT=QI
              ELSEIF(TMT0.GE.0.)THEN
                  QSAT=QINT
              ELSE
                IF(IWU.GT.0.0) THEN
                  QSAT=QI
                ELSE
                  QSAT=QINT
                ENDIF
              ENDIF

!CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT
!C             DELETE THIS LINE TO SWITCH BACK TO RH VS ICE
              QSAT=QW
!CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT
              QABV =RHU*QSAT
              QABV =AMAX1(1.e-12,QABV)

              B    =TABV
              BQ   =QABV
              FAC  =0.
              AHF  =0.
              AHFQ =0.
              Q2A  =0.50*(Q2(I,J,n,1)+Q2(I,J,n,2))
              BQ2  =Q2A
        ELSEIF(NL1X(I,J,n).EQ.LP1)THEN
!---------------------------------------------------------------------
!***  EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
!---------------------------------------------------------------------
!
        PL=PINT(I,J,n,LM-1)
        TL=0.5*(T(I,J,n,LM-2)+T(I,J,n,LM-1))
        QL=0.5*(Q(I,J,n,LM-2)+Q(I,J,n,LM-1))
        IWL=0.50*(IW(I,J,n,LM-2)+IW(I,J,n,LM-1))
        TBLO=TL*(SPL(LP)/PL)**RGAMOG

              TMT0=TL-273.16
              TMT15=AMIN1(TMT0,-15.)
              AI=0.008855
              BI=1.
              IF(TMT0.LT.-20.)THEN
                AI=0.007225
                BI=0.9674
              ENDIF
              QW=PQ0/PL*EXP(A2*(TL-A3)/(TL-A4))
              QI=QW*(BI+AI*AMIN1(TMT0,0.))
              QINT=QW*(1.-0.00032*TMT15*(TMT15+15.))
              IF(TMT0.LT.-15.)THEN
                  QSAT=QI
              ELSEIF(TMT0.GE.0.)THEN
                  QSAT=QINT
              ELSE
                IF(IWL.GT.0.0) THEN
                  QSAT=QI
                ELSE
                  QSAT=QINT
                ENDIF
              ENDIF
!CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT
!C             DELETE THIS LINE TO SWITCH BACK TO RH VS ICE
              QSAT=QW
!CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT
              RHL=QL/QSAT
              IF(RHL.GT.1)THEN
               RHL=1
               QL =RHL*QSAT
              ENDIF

              IF(RHL.LT..01)THEN
                RHL=.01
                QL =RHL*QSAT
              ENDIF
              TMT0=TBLO-273.16
              TMT15=AMIN1(TMT0,-15.)
              AI=0.008855
              BI=1.
              IF(TMT0.LT.-20.)THEN
                AI=0.007225
                BI=0.9674
              ENDIF
              QW=PQ0/SPL(L)  &
               *EXP(A2*(TBLO-A3)/(TBLO-A4))
              QI=QW*(BI+AI*AMIN1(TMT0,0.))
              QINT=QW*(1.-0.00032*TMT15*(TMT15+15.))
              IF(TMT0.LT.-15.)THEN
                  QSAT=QI
              ELSEIF(TMT0.GE.0.)THEN
                  QSAT=QINT
              ELSE
                IF(IWL.GT.0.0) THEN
                  QSAT=QI
                ELSE
                  QSAT=QINT
                ENDIF
              ENDIF
!CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT
!C             DELETE THIS LINE TO SWITCH BACK TO RH VS ICE
              QSAT=QW
!CMEB 12/22/98 SWITCH TO RH VS WATER NO MATTER WHAT
              QBLO =RHL*QSAT
              QBLO =AMAX1(1.e-12,QBLO)

              B    =TBLO
              BQ   =QBLO
              FAC  =0.
              AHF  =0.
              AHFQ  =0.
              Q2A  =0.50*(Q2(I,J,n,LMH(I,J,n)-1)+Q2(I,J,n,LMH(I,J,n)))
              BQ2  =Q2A
        ELSE
!---------------------------------------------------------------------
!***  INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
!---------------------------------------------------------------------
!
        B     =T(I,J,n,NL1X(I,J,n))
        BQ    =Q(I,J,n,NL1X(I,J,n))
        FAC  =2.*ALOG(PT+PDSL(I,J,n)*AETA(NL1X(I,J,n)))
        AHF  =(B-T(I,J,n,NL1X(I,J,n)-1))/     &
               (ALPINT(I,J,n,NL1X(I,J,n)+1)-ALPINT(I,J,n,NL1X(I,J,n)-1))
        AHFQ =(BQ-Q(I,J,n,NL1X(I,J,n)-1))/        &
                   (ALPINT(I,J,n,NL1X(I,J,n)+1)-ALPINT(I,J,n,NL1X(I,J,n)-1))

        Q2B   =0.50*(Q2(I,J,n,NL1X(I,J,n)-1)+Q2(I,J,n,NL1X(I,J,n)))
!C
       IF(NL1X(I,J,n).GT.2)THEN
           Q2A=0.50*(Q2(I,J,n,NL1X(I,J,n)-2)+Q2(I,J,n,NL1X(I,J,n)-1))
       ELSE
           Q2A=Q2B
       ENDIF
!C
        BQ2=Q2B*HTM(I,J,n,NL1X(I,J,n))
        AHFQ2=(BQ2-Q2A*HTM(I,J,n,NL1X(I,J,n)-1))/    &
             (ALPINT(I,J,n,NL1X(I,J,n)+1)-ALPINT(I,J,n,NL1X(I,J,n)-1))
        if(i.eq.25.and.j.eq.22.and.n.eq.3)then
!	   print *,'tsl,',lp,nl1x(i,j,n),ahf,ALPINT(I,J,n,NL1X(I,J,n)+1) &
!	          ,ALPINT(I,J,n,NL1X(I,J,n)-1)
	endif

        ENDIF

        TSL(I,J,n)=B+AHF*(TRF-FAC)
        QSL(I,J,n)=BQ+AHFQ*(TRF-FAC)
        QSL(I,J,n)=AMAX1(QSL(I,J,n),1.e-12)
        Q2SL(I,J,n)=BQ2+AHFQ2*(TRF-FAC)
        Q2SL(I,J,n)=AMAX1(Q2SL(I,J,n),0.)


     



        FSL(I,J,n)=(PNL1-SPL(LP))/(SPL(LP)+PNL1)   &
            *((ALSL(LP)+ALPINT(I,J,n,NL1X(I,J,n))-FAC)*AHF+B)*Rd*2.  &
            +ZINT(I,J,n,NL1X(I,J,n))*G
        if(i.eq.25.and.j.eq.22.and.n.eq.3)then
!	   print *,'tsl,',lp,tsl(i,j,n),b,ahf,trf,fac,nl1x(i,j,n)
	endif
  220   CONTINUE

          do 281 n=1,nm
          DO 281 J=1,jm
          DO 281 I=1,IM
!NOTE
!NOTE         29 JANUARY 1993, RUSS TREADON.
!NOTE          - AS FOR THE OTHER FIELDS WE INTERPOLATE ONLY
!NOTE            BETWEEN THE FAL AND THE MODEL TOP.  BELOW
!NOTE            SURFACE VALUES ARE FAL VALUES.
!
            LMB = LMV(I,J,n)
!
              PETAU=PT+PDVP1(I,J,n)*ETA(1)
              ALPETU=ALOG(PETAU)
            DO 280 IL=2,LMB
              PETAL=PT+PDVP1(I,J,n)*ETA(IL)
!             PETAU=PT+PDVP1(I,J)*ETA(IL-1)
              ALPETL=ALOG(PETAL)
!             ALPETU=ALOG(PETAU)
              ALPET2=SQRT(0.5E0*(ALPETL*ALPETL+ALPETU*ALPETU))
!
!          SEARCH FOR HIGHEST MID-LAYER ETA SURFACE (NOT SUBMERGED)
!          THAT IS BELOW THE GIVEN STANDARD PRESSURE LEVEL.
              IF(ALSL(Lp).LT.ALPET2)THEN
                NL1X(I,J,n)=IL-1
                ALPETUX(I,J,n)=ALPETU
                ALPET2X(I,J,n)=ALPET2
                GO TO 281
              ENDIF
!      If we arent on the last iterate of the 280 loop, reset  PETAU and ALPETU
            if ( il .eq. lmb ) goto 280
            PETAU=PETAL
            ALPETU=ALPETL
  280       CONTINUE
            NL1X(I,J,n)=LMB+1
            ALPETUX(I,J,n)=ALPETU
            ALPET2X(I,J,n)=ALPET2
 281     CONTINUE
!
!         BELOW GROUND USE FAL WINDS.
!
!$omp  parallel do
!$omp  private(alpet1,alpetl,alpetu,fact,petau)
          do 290 n=1,nm
          DO 290 J=1,jm
          DO 290 I=1,IM
            IF(NL1X(I,J,n).GT.LMV(I,J,n))THEN
              USL(I,J,n)=U(I,J,n,LMV(I,J,n))
              VSL(I,J,n)=V(I,J,n,LMV(I,J,n))
!
!          IF REQUESTED PRESSURE LEVEL IS NOT BELOW THE LOCAL GROUND
!          THEN WE HAVE TWO POSSIBILITIES.  IF THE REQUESTED PRESSURE
!          LEVEL IS BETWEEN THE LOCAL SURFACE PRESSURE AND TOP OF
!          MODEL PRESSURE, VERTICALLY INTERPOLATE BETWEEN NEAREST
!          BOUNDING ETA LEVELS TO GET THE WIND COMPONENTS.  IF THE
!          REQUESTED PRESSURE LEVEL IS ABOVE THE MODEL TOP, USE
!          CONSTANT EXTRAPOLATION OF TOP ETA LAYER (L=1) WINDS.
!
            ELSE
              IF(NL1X(I,J,n).GT.1)THEN
                ALPETL=ALPETUX(I,J,n)
                PETAU=PT+PDVP1(I,J,n)*ETA(NL1X(I,J,n)-1)
                ALPETU=ALOG(PETAU)
                ALPET1=SQRT(0.5*(ALPETL*ALPETL+ALPETU*ALPETU))
                FACT=(ALPET2X(I,J,n)-ALSL(Lp))/(ALPET2X(I,J,n)-ALPET1)
                USL(I,J,n)=U(I,J,n,NL1X(I,J,n))  &
                       +(U(I,J,n,NL1X(I,J,n)-1)-U(I,J,n,NL1X(I,J,n)))*FACT
                VSL(I,J,n)=V(I,J,n,NL1X(I,J,n))+(V(I,J,n,NL1X(I,J,n)-1)  &
                        -V(I,J,n,NL1X(I,J,n)))*FACT
              ELSE
                USL(I,J,n)=U(I,J,n,NL1X(I,J,n))
                VSL(I,J,n)=V(I,J,n,NL1X(I,J,n))
              ENDIF
!
!            ALPET2 IS MID-LAYER ETA SURFACE JUST BELOW STANDARD PRESSURE
!            LEVEL AND ALPET1 IS DASHED ETA SURFACE JUST ABOVE.
!            NOTE THAT IF THE STANDARD PRESSURE SURFACE IS SUBMERGED, THEN
!            ALPET2 AND ALPET1 ARE THE LOWEST AND 2ND LOWEST MID-LAYER
!            ETA SURFACES ABOVE THE TOPOGRAPHY (WITH OLDRD=.TRUE., ZJ).
!
            ENDIF
!	    print *,i,j,n,usl(i,j,n),vsl(i,j,n)
  290     CONTINUE

        do j=1,jm
	do i=1,im
!	  print *,i,j,q2sl(i,j,3)
        enddo
	enddo

!        print *,'usl,99,58,6,',usl(99,58,6),usl(99,59,6) 
!        print *,'vsl,99,58,6,',lp,vsl(99,58,6),vsl(99,59,6) 
        call bocoh(fsl,im,jm,nm)
	call bocov(usl,vsl,im,jm,nm)
!        print *,'usl,99,58,6,',usl(99,58,6),usl(99,59,6) 
!        print *,'vsl,99,58,6,',lp,vsl(99,58,6),vsl(99,59,6) 
 
	write(1)fsl
	write(1)usl
	write(1)vsl
        write(1)qsl  !DRAGAN 23.09.'10
        write(1)tsl
!!!avg2017        write(1)q2sl

!!!        write(1)w
!!!	write(1)dwdt   !non-hydrostatic 2016
!!!	write(1)z


!!      write(1)tshltr
!!      write(1)qshltr
!      write(1)plm


  310   continue

	close(1)

      DEALLOCATE(tsl,qsl,fsl,q2sl,iw)

      DEALLOCATE(nl1x,ihold,jhold,nnhold)

      DEALLOCATE(alpetux,alpet2x,usl,vsl)

	end subroutine eta2p
