      subroutine initpost(nh)

      USE mod_extra
      implicit none
      include 'param_o.h'
      include 'vrbls.comm'
!!JNT      include 'extra.comm'
      include 'mapot.comm'
      include 'const.h'
      include 'dynam.comm'
      include 'masks.comm'
      include 'pvrbls.comm'
     
      
      integer,parameter::npes=nsub*nsub*nm
      integer::nh,n,mype0,ix,jy,imaxV,imaxH,imin,jmaxV,jmaxH,jmin, &
               mype,i,j,l
      real, dimension(0:ilm+1,0:jlm+1):: dum1
      real, dimension(0:ilm+1,0:jlm+1,lm):: dum2
      real, dimension(0:ilm+1,0:jlm+1,nsoil):: dum3!!!!!!!!!!!!!!!!!!!!!!!!
      
    !  character (len=03):: c_mype, c_nh !!!!!!!!DRAGAN
      character (len=04):: c_mype, c_nh  !DRAGAN
      !character (len=02):: c_mype!DRAGAN
      
      
      character (len=80):: c_filename1!!!!!!!!!!!!!!!!!!!!!!!
     
     
      real::pbi
      integer::lmap1
      real,parameter::R=287.04

!      spl=(/50000./)
!      spl=(/100000./)
      spl=(/5000.,10000.,15000.,20000.,25000.,30000.,35000.,40000.,  &
            45000.,50000.,55000., &
            60000.,65000.,70000.,75000.,80000.,85000.,90000.,95000.,100000./)
      alsl=alog(spl)

      write(c_nh,1000) nh
 !1000 format(i3.3) !!!!!!!!!!!DRAGAN
 
 1000 format(i4.4)
 
!***************************************************************************   
        do mype = 0, npes-1
!***************************************************************************

          n = mype/nxy + 1
          mype0 = mype - (n-1)*nxy
          jy = (mype0/ixm) + 1
          ix = mod(mype0, ixm)+1
          list = 100 + mype + 1

!-----------------------------------------------------------------------   
          write(c_mype,1000) mype
!          c_filename1='../../RUN/data/data_out/GEFfcst_'//   &
           c_filename1='./test2/GEFfcst_'//   &
                  c_mype//"."//c_nh!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         open(unit=list, file=c_filename1, form='unformatted')
!-----------------------------------------------------------------------   
         imaxV = ix * ildom
         imaxH = imaxV + 1
         imin  = imaxH - ildom
         jmaxV = jy * jldom
         jmaxH = jmaxV + 1
         jmin  = jmaxH - jldom
!-----------------------------------------------------------------------   
         read(list) dum1 
         pd(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 read(list)dum2
         t(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
	 read(list)dum2
         u(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
 	 read(list)dum2
         v(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
 	 read(list)dum2
         q(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
!! 	   read(list)dum2
!!         q2(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
 	 read(list)dum2

        
         cwm(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
         read(list) dum1 
         ths(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
         read(list) dum1 
         u00(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)

          read(list)ul
       
!         read(list) dum1 
!         tshltr(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
          

!         read(list) dum1 
!         qshltr(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
 
!         read(list) dum1 
!         plm(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)

!!!!!!!!!!!!!!!!!!!!non-hydrostatic!!!!!!!!!!!!!!!

         w(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
 	 read(list)dum2
	 
	 dwdt(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
 	 read(list)dum2
	 
	 z(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
 	 read(list)dum2
!!!!!!!!!!!!!!!!!!!!!





      close(list)


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!         c_filename1='../../RUN/data/data_out/GEFfcstVEG_'//   &
!                  c_mype//"."//c_nh
         c_filename1='./test2/GEFfcstVEG_'//  &
	 c_mype//"."//c_nh
	 
	 open(unit=list, file=c_filename1, form='unformatted')
         
	 read(list) dum1 
         tshltr(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         qshltr(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 
	 read(list) dum1 
         plm(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 
	 
	 read(list) dum1 
         sfclhx(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	 read(list) dum1 
         sfcshx(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         qwbs(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         twbs(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         grnflx(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         vegfrc(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
         IVGTYP(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	read(list) dum1 
         ISLTYP(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
         glat(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	read(list) dum1 
         glon(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         sm(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
!	 read(list)dum2
!         htm(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
	 
!	read(list)dum2
!         vtm(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
	
	 read(list) dum1 
        sice(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         z0(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list)dum2
         omgalf(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
	
	 read(list) dum1 
         albedo(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	read(list) dum1 
        albase(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        mxsnal(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        czen(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	read(list) dum1 
        czmean(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	read(list) dum1 
        fis(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        q10(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        th10(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        qz0(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        thz0(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        u10(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        v10(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        uz0(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        vz0(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        ustar(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        ths(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm) 
	 
         read(list)dum2
         q2(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
	 
	 read(list)dum2
         tcucn(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
	 
	 read(list)dum2
         train(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
	 
	 read(list)dum2
         div(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
	 
	 read(list)dum2
         rtop(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)

      close(list)


!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


!         c_filename1='../../RUN/data/data_out/GEF_sw_rdtn'//   &
         c_filename1='./test2/GEF_sw_rdtn'//   &
                  c_mype//"."//c_nh
         open(unit=list, file=c_filename1, form='unformatted')
         read(list) dum1 
         RSWIN(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	 read(list) dum1 
         RSWOUT(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         RSWTOA(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         ASWIN(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	 read(list) dum1 
         ASWOUT(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	 read(list) dum1 
         ASWTOA(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         rswnet(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
         read(list)dum2
         RSWTT(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)

         close(list)
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


!         c_filename1='../../RUN/data/data_out/GEF_lw_rdtn'//   &
         c_filename1='./test2/GEF_lw_rdtn'//   &
                  c_mype//"."//c_nh
         open(unit=list, file=c_filename1, form='unformatted')
         read(list) dum1 
         RLWIN(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	 read(list) dum1 
         RLWOUT(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         RLWTOA(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         ALWIN(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	 read(list) dum1 
         ALWOUT(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	 read(list) dum1 
         ALWTOA(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         rlwnet(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
         read(list)dum2
         RLWTT(imin:imaxH, jmin:jmaxH,n,:) = dum2(1:ilm,1:jlm,:)
	 
        read(list) dum1 
         epsr(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)

        read(list) dum1 
        radot(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        sigt4(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)

         close(list)
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!        c_filename1='../../RUN/data/data_out/GEF_clouds'//   &
         c_filename1='./test2/GEF_clouds'//   &
                  c_mype//"."//c_nh
         open(unit=list, file=c_filename1, form='unformatted')
         read(list) dum1 
         CUPPT(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	 read(list) dum1 
         cfracl(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         cfracm(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         cfrach(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)

       read(list) dum1 
         acfrcv(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
       read(list) dum1 
         acfrst(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)

       read(list) dum1 
       cnvtop(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
       read(list) dum1 
       cnvbot(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
       read(list) dum1 
       htop(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
       read(list) dum1 
       hbot(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)


         close(list)
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!        c_filename1='../../RUN/data/data_out/GEF_precANDsoil'//   &
         c_filename1='./test2/GEF_precANDsoil'//   &
	          c_mype//"."//c_nh
         open(unit=list, file=c_filename1, form='unformatted')
         
	 read(list) dum1 
         acprec(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	 read(list) dum1 
         cuprec(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         aprec(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         acsnow(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)

      read(list) dum1 
      acsnom(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)

         read(list) dum1 
         prec(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
!- - - - - - - - - - - 
     
      read(list)dum3
      stc(imin:imaxH, jmin:jmaxH,n,:) = dum3(1:ilm,1:jlm,:)


      read(list)dum3
      smc(imin:imaxH, jmin:jmaxH,n,:) = dum3(1:ilm,1:jlm,:)

         read(list) dum1 
         SMSTAV(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
	 read(list) dum1 
         BGROFF(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)

         read(list) dum1 
         SSROFF(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
      read(list) dum1 
      cmc(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	 
      read(list) dum1 
      potevp (imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)

      read(list) dum1 
      potflx(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
         
	 
	 read(list) dum1 
        sfcevp(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm) 
	
	read(list) dum1 
        subshx(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        snopcx(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        si(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        sno(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        smstot(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        soiltb(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
!	read(list) dum1 
!        tg(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        sr(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        akhs(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm)
	
	read(list) dum1 
        akms(imin:imaxH, jmin:jmaxH,n) = dum1(1:ilm,1:jlm) 

   close(list)



!***********************************************************************
            end do
  
!***********************************************************************



!-----------------------------------------------------------------------   
       

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

!	open(1,file='../data_out/sfceta.dat',form='unformatted')
	open(1,file='./NCEP_GEF_initdata/sfceta.dat',form='unformatted')
	read(1)fis,res,htm,vtm,lmh,lmv
	close(1)


!	open(1,file='../data_out/dynam.dat',form='unformatted')
	open(1,file='./NCEP_GEF_initdata/dynam.dat',form='unformatted')
	read(1)fadv, fadv2,fadt, rd,  f4d,f4q2, ef4t, fkin, &
                    deta, rdeta, aeta, eta, dfl, daeta
	close(1)

!      print *,fadv,fadt,rd,f4d
      do l=1,lp1
!        print *,l,deta(l),dfl(l)
      enddo
!      stop

      do n=1,nm
      DO J=1,jm
      DO I=1,IM
        PDSL(I,J,n)    =RES(I,J,n)*PD(I,J,n)
        PBI            =PD(I,J,n)+PT
        PINT(I,J,n,LP1)  =PBI
        ALPINT(I,J,n,LP1)=ALOG(PBI)
	if(i.eq.56.and.j.eq.65.and.n.eq.6)then
!	  print *,i,j,n,pd(i,j,n),pdsl(i,j,n),res(i,j,n),pint(i,j,n,lp1)
        endif
      ENDDO
      ENDDO
      ENDDO

!
!     COMPUTE PRESSURE VALUES ABOVE THE SURFACE.
!
      DO L=LM,1,-1
      do n=1,nm
      DO J=1,jm
      DO I=1,IM
          PBI          =PDSL(I,J,n)*ETA(L)+PT
          PINT(I,J,n,L)  =PBI
          ALPINT(I,J,n,L)=ALOG(PBI)
	if(i.eq.1.and.j.eq.21.and.n.eq.1)then
!	  print *,i,j,n,l,pbi,pdsl(i,j,n),res(i,j,n),eta(l)
        endif
      ENDDO
      ENDDO
      ENDDO
      ENDDO

      do n=1,nm
      DO J=1,JM-1
      DO I=1,IM-1
        PDVP1(I,J,n)=0.25*(PDSL(I,J,n)+PDSL(I+1,J,n)  &
                       +PDSL(I,J+1,n)+PDSL(I+1,J+1,n))
      ENDDO
      ENDDO
      ENDDO

!C
!C       COMPUTE ALPINT AND PINT ONLY FOR THOSE POINTS THAT WE OWN
!C
!C       IF ( J .GE. JSTA .AND. J .LE. JEND ) THEN
!C
      do n=1,nm
      do j=1,jm
      do i=1,im
        LMAP1=LMH(I,J,n)+1
        DO L=LMAP1,LM
          ALPINT(I,J,n,L+1)=(DFL(L)-DFL(L+1))/(R*T(I,J,n,L))  &
                        +ALPINT(I,J,n,L)
          PINT(I,J,n,L+1)=EXP(ALPINT(I,J,n,L+1))
        ENDDO
      enddo
      enddo
      enddo

      do l=1,lm
!        print *,l,t(30,101,3,l),t(30,81,3,l),q(30,101,3,l)
      enddo

      end subroutine initpost
