C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
                             SUBROUTINE CNSTS(htmo,vtmo,lmho,lmvo)
C     ******************************************************************
C     *                                                                *
C     *  ROUTINE FOR INITIALIZATION OF CONSTANTS AND VARIABLES         *
C     *                                                                *
C     ******************************************************************
C 01-03-01 H CHUANG MODIFIED for MPI
C 01-05-31 H CHUANG MODIFIED TO INCLUDE NEW SOIL PARM 
C----------------------------------------------------------------------
c      use mskout
c
      INCLUDE "parmeta.resout"
      INCLUDE "parmeta.resin"
      INCLUDE "parmsoil"
      INCLUDE "mpif.h"
      INCLUDE "MPI_PARM.comm"
C----------------------------------------------------------------------
                             P A R A M E T E R
     & (D00=0.0,D5=5.E-1,D01=1.00E-2,H1=1.0,HM1=-1.0
     &, H90=90.0,H360=360.0,EPS=1.E-10)
                             P A R A M E T E R
     & (IM=IMO,JM=JMO,LM=LMO)
                             P A R A M E T E R
     & (ITB=076,JTB=134,ITBQ=152,JTBQ=440)
                             P A R A M E T E R
     & (IM2=IM-2,JM1=JM-1,JM2=JM-2,JM3=JM-3,JM4=JM-4,JM5=JM-5
     &, IMJM=IM*JM-JM/2,JAM=6+2*(JM-10)
     &, KHN=IM-1,KHS=-IM,KB=2*IMO+JMO-3
     &, KNE=IM    ,KNW=IM-1 ,KSW=-IM    ,KSE=-IM+1
     &, KHL00=1                    ,KHH00=IM*JM-JM/2
c     &, LP1=LM+1,nglats=190,nglons=384)
     &, LP1=LM+1,nglats=256,nglons=512)
C      In the statement above :
C      for t62 =====> nglats=94,   nglons=192
C      for t126 ====> nglats=190,  nglons=384
C-----------------------------------------------------------------------
                             L O G I C A L
     & SIGMA, NUSNOSST, LUSAF
C
cwas      REAL,ALLOCATABLE,DIMENSION(:,:,:)::HTMO
C      REAL,ALLOCATABLE,DIMENSION(:,:)::DX,WPDAR,CPGFU,CURV,FCP
C      REAL,ALLOCATABLE,DIMENSION(:,:)::FDIV,FAD,F,DDMPU,DDMPV
C      REAL,ALLOCATABLE,DIMENSION(:,:)::EPSR,GFFC,WFK,HDAC,HDACV
C-----------------------------------------------------------------------
                             C O M M O N  /SFCFLDS/
     & SST(IM,JM),ALBEDO(IM,JM),SI(IM,JM),SNO(IM,JM),SMC(IM,JM,NSOIL)
     &,STC(IM,JM,NSOIL),CMC(IM,JM)
C
                             C O M M O N /SFCTYPES/
     & IVGTPK(IM,JM),ISLTPK(IM,JM),ISPTPK(IM,JM)
C--------------------------------------------------------------------
                             C O M M O N /OUTETA/
     1 PDO(IM,JM),PDB(KB,2)
     5,TB (KB,LMO,2),QB (KB,LMO,2),UB (KB,LMO,2),VB (KB,LMO,2)
     6,PT,DETA(LMO),AETA(LMO),ETA(LMO+1),DFL(LMO+1),SM(IM,JM)
     7,RES(IM,JM),FIS(IM,JM)
C--------------------------------------------------------------------
                             C O M M O N /DAT/
     1 IDAT(3) 
Cmp
c	COMMON /MSKOUT/
c     1 HTMO(IM,JM,LMO),VTMO(IM,JM,LMO),
c     2 LMHO(IM,JM),LMVO(IM,JM)
C-----------------------------------------------------------------------
                             D I M E N S I O N
     & DXJ   (JM),WPDARJ(JM),CPGFUJ(JM),CURVJ (JM),FCPJ  (JM)
     &,FDIVJ (JM),EMJ   (JM),EMTJ  (JM),FADJ  (JM)
     &,DDMPUJ(JM),DDMPVJ(JM),HDACJ (JM)
       character *32 glabel

C-----------------------------------------------------------------------
                             D I M E N S I O N
     & KHLA  (JAM),KHHA  (JAM),KVLA  (JAM),KVHA  (JAM)
     &,KHL2  (JM),KHH2  (JM),KVL2  (JM),KVH2  (JM)
     &,KHL3  (JM),KHH3  (JM),KVL3  (JM),KVH3  (JM)
C-----------------------------------------------------------------------
                             C O M M O N /MASKS/
     & HBM2(IM,JM),VBM2(IM,JM),VBM3  (IM,JM),SICE(IM,JM)
C-----------------------------------------------------------------------
      COMMON /LATLONS/
     & HLATO(IMO,JMO),HLONO(IMO,JMO),VLATO(IMO,JMO),VLONO(IMO,JMO)
     & ,ihinx(imo,jmo),ihiny(imo,jmo),ivinx(imo,jmo),iviny(imo,jmo)
     & ,ih2hx(imo,jmo),ih2hy(imo,jmo),wh(imo,jmo,4),wv(imo,jmo,4)
C-----------------------------------------------------------------------
       COMMON /IMMSK/ MSKSCVH (1024,1024)
       COMMON /AFMSK/ MSKAF   (512,512)

      COMMON /SLAKE/KLAKE(21)
C-----------------------------------------------------------------------
      REAL MXSNAL(IM,JSTA:JEND),ALBASE(IM,JSTA:JEND)
                             D I M E N S I O N
     & RDETA (LM),F4Q2  (LM)
     &,EM    (JAM),EMT   (JAM)
C
     &,DX    (IM,JM),WPDAR(IM,JM),CPGFU(IM,JM),CURV(IM,JM),FCP(IM,JM)
     &,FDIV  (IM,JM),FAD(IM,JM)
     &,F(IM,JM),DDMPU(IM,JM),DDMPV(IM,JM)
     &,GLAT  (IM,JM),GLON  (IM,JM),albct(361,180,4),albc1(361,180)
     &,albc2 (361,180),albc3(361,180),albc4(361,180)
     &,ALBC  (361,180), i1d(360,180), ioneta(im,jm)
     &,AVLM  (361,180),AVAILM(IM,JM)
     &,DUM2D(IM,JM),IDUM2D(IM,JM)
     &,LMHO(IMO,JSTA:JEND),LMVO(IMO,JSTA:JEND)
     &,HTMO(IMO,JSTA:JEND,LMO),VTMO(IMO,JSTA:JEND,LMO)
CFEI &,fpar1(360,180),fpar2(360,180)
Cmp     &,fpar1(2500,1250),fpar2(2500,1250)
CYL
Cmp     &,kpds(25), kgds(22), jpds(25), jgds(22)
Cmp      logical*1 bitmap(2500,1250)
CYL

      COMMON /VEGGIE/ FPAR1(2500,1250),FPAR2(2500,1250)
	COMMON /TGRND/ TG(IM,JM)

C-----------------------------------------------------------------------
       equivalence (albc1,albct(1,1,1)),(albc2,albct(1,1,2))
     &,(albc3,albct(1,1,3)),(albc4,albct(1,1,4))
c
c      equivalence (fpar1,albct(1,1,1)),(fpar2,albct(1,1,2))
c
                             D I M E N S I O N
     & DETA2(LM),AETA2(LM),ETA2(LP1),DFRLG(LP1)
     &,QS0   (JTB),SQS   (JTB),THE0  (ITB),STHE  (ITB)
     &,                        THE0Q(ITBQ),STHEQ(ITBQ)
C
     &,EPSR  (IM,JM),GFFC   (IM,JM),WFK   (IM,JM)
     &,HDAC  (IM,JM),HDACV (IM,JM)
C
     &,PTBL  (ITB,JTB),TTBL  (JTB,ITB),TTBLQ (JTBQ,ITBQ)
C-----------------------------------------------------------------------
                             D I M E N S I O N
     & DETA1 (LM),AETA1 (LM),ETA1  (LP1), julm(13)
C-----------------------------------------------------------------------
                             D I M E N S I O N
     & SLDPTH(NSOIL),RTDPTH(NSOIL)
C--------------HORIZONTAL GRID CONSTANTS--------------------------------
                             N A M E L I S T /ECONDATA/
     & TLM0D,TPH0D,WBD,SBD,DLMD,DPHD
     &,DT,W,IDTAD,IDTCF,COAC,CODAMP,TDDAMP,DFC,DDFC,SIGMA,NUSNOSST,LUSAF
     &,SLDPTH,RTDPTH
C--------------STEREOGRAPHIC MAP CONSTANTS------------------------------
                             D A T A
     & julm/0,31,59,90,120,151,181,212,243,273,304,334,365/
                             D A T A
     & CMLD/-80.00/,DP30/022.86/,X1P/108./,Y1P/120./,IXM/145/,IYM/130/
C--------------UNIVERSAL CONSTANTS--------------------------------------
                             D A T A
     & A/6376000./,G/9.806/,TWOM/.00014584/,R/287.04/,CP/1004.6/
     &,PI/3.141592654/
C--------------SURFACE DATA---------------------------------------------
                             D A T A
     & ROS/1500./,CS/1339.2/,DS/.100/,AKS/.0000005  /,DZG/2.5000/
     &,TG0/258.16/,TGA/30./
     &,ROI/916.6/,CI/2060.0/,DI/.100/,AKI/.000001075/,DZI/2.0/
     &,TI0/271.16/,THL/210./,PLQ/70000./
C--------------LOGICAL FILE NAMES---------------------------------------
                             D A T A
     & NHIBU /54/,NFCST /53/,NBC   /55/,LIST  /06/
C***********************************************************************
      rewind 15
      READ(15,ECONDATA)
      if(mype.eq.0) write(6,ECONDATA)
C--------------DERIVED VERTICAL CONSTANTS-------------------------------
                             DO 110 L=1,LP1
                             DFRLG(L)=DFL(L)
 110                         DFL(L)=DFL(L)*G
C--------------DUMMY CONSTANTS------------------------------------------
      PT1=PT
      PT2=PT
      R1=R
      R2=R
                             DO 120 L=1,LP1
      ETA1(L)=ETA(L)
 120  ETA2(L)=ETA(L)
                             DO 130 L=1,LM
      DETA1(L)=DETA(L)
      AETA1(L)=AETA(L)
      DETA2(L)=DETA(L)
 130  AETA2(L)=AETA(L)
C--------------DERIVED GEOMETRICAL CONSTANTS----------------------------
      DTR=PI/180.
      TPH0=TPH0D*DTR
      WB=WBD*DTR
      SB=SBD*DTR
      DLM=DLMD*DTR
      DPH=DPHD*DTR
      TDLM=DLM+DLM
      TDPH=DPH+DPH
      RDLM=1./DLM
      RDPH=1./DPH
C
      WBI=WB+TDLM
      SBI=SB+TDPH
      EBI=WB+IM2*TDLM
      ANBI=SB+JM3*DPH
C
      STPH0=SIN(TPH0)
      CTPH0=COS(TPH0)
C---------------TIME STEPPING RELATED CONSTANTS-------------------------
      TSPH=3600./DT
      NDDAMP=TDDAMP*TSPH+.5
C
      DTAD=IDTAD
      DTCF=IDTCF
C--------------DERIVED HORIZONTAL GRID CONSTANTS------------------------
      DY=A*DPH
      CPGFV=-DT/(48.*DY)
      EN= DT/( 4.*DY)*DTAD
      ENT=DT/(16.*DY)*DTAD
C
      IF(MYPE.EQ.0)THEN     
              DO J=1,JMO
C
C              KHL0(J)=IM*(J-1)-(J-1)/2+1
C              KVL1(J)=IM*(J-1)-(J-1)/2+1
              KHL2(J)=IM*(J-1)-(J-1)/2+2
C              KVL3(J)=IM*(J-1)-(J-1)/2+2
C              KHL4(J)=IM*(J-1)-(J-1)/2+3
C              KVL5(J)=IM*(J-1)-(J-1)/2+3
C              KHL6(J)=IM*(J-1)-(J-1)/2+4
C
C              KVL0(J)=IM*(J-1)-J/2+1
C              KHL1(J)=IM*(J-1)-J/2+2
              KVL2(J)=IM*(J-1)-J/2+2
C              KHL3(J)=IM*(J-1)-J/2+3
C              KVL4(J)=IM*(J-1)-J/2+3
C              KHL5(J)=IM*(J-1)-J/2+4
C              KVL6(J)=IM*(J-1)-J/2+4
C
C              KHH0(J)=IM*J-J/2
              KHH2(J)=IM*J-J/2-1
C              KVH1(J)=IM*J-J/2-1
C              KHH4(J)=IM*J-J/2-2
C              KVH3(J)=IM*J-J/2-2
C              KHH6(J)=IM*J-J/2-3
C              KVH5(J)=IM*J-J/2-3
C
C              KHH1(J)=IM*J-(J+1)/2
C              KVH0(J)=IM*J-(J+1)/2
C              KHH3(J)=IM*J-(J+1)/2-1
              KVH2(J)=IM*J-(J+1)/2-1
C              KHH5(J)=IM*J-(J+1)/2-2
C              KVH4(J)=IM*J-(J+1)/2-2
C              KVH6(J)=IM*J-(J+1)/2-3
C
             END DO
      END IF
c             TPH=SB-DPH
      IF(MYPE.EQ.0)THEN
	     DO 140 J=1,JM
c              TPH=TPH+DPH
              TPH=SB+float(J-1)*DPH
              DXP=A*DLM*COS(TPH)
              DXJ(J)=DXP
CVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
C             WPDARJ(J)=-DT*W*100000./(32.*DXP*DY)
        WPDARJ(J)=-W*((A*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+DY**2)
     2             /(DT*32.*DXP*DY)*.88
C    2             /(DT*32.*DXP*DY)
CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
              CPGFUJ(J)=-DT/(48.*DXP)
              CURVJ(J)=.5*DT*TAN(TPH)/A
              FCPJ(J)=DT/(CP*192.*DXP*DY)
              FDIVJ(J)=1./(12.*DXP*DY)
              EMJ(J)= DT/( 4.*DXP)*DTAD
              EMTJ(J)=DT/(16.*DXP)*DTAD
              FADJ(J)=-DT/(48.*DXP*DY)*DTAD
CVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
C             ACDT=COAC*DT*DTCF
C             ACDT=COAC*DT
C    2            *SQRT((A*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+DY**2)
              ACDT=DT
     2            *SQRT((A*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+DY**2)
CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
C             CDDAMP=CODAMP*ACDT
              CDDAMP=0.04*CODAMP*ACDT
C             HDACJ(J)=ACDT/(4.*DXP*DY)
              HDACJ(J)=COAC*ACDT/(4.*DXP*DY)
              DDMPUJ(J)=CDDAMP/DXP
              DDMPVJ(J)=CDDAMP/DY
 140  CONTINUE
      ENDIF
C***
      CALL MPI_BCAST(DXJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
      CALL MPI_BCAST(WPDARJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
      CALL MPI_BCAST(CPGFUJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
      CALL MPI_BCAST(CURVJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
      CALL MPI_BCAST(FCPJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
      CALL MPI_BCAST(FDIVJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
      CALL MPI_BCAST(EMJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
      CALL MPI_BCAST(EMTJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
      CALL MPI_BCAST(FADJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
      CALL MPI_BCAST(HDACJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
      CALL MPI_BCAST(DDMPUJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
      CALL MPI_BCAST(DDMPVJ,JM,MPI_REAL,0,MPI_COMM_WORLD,IRTN)
C
C--------------SPREADING OF UPSTREAM HEIGHT-POINT ADVECTION FACTOR------
C
          JA=0
          DO 161 J=3,5
          JA=JA+1
          KHLA(JA)=2
          KHHA(JA)=IM-1-MOD(J+1,2)
          EMT(JA)=EMTJ(J)
 161      CONTINUE
          DO 162 J=JM4,JM2
          JA=JA+1
          KHLA(JA)=2
          KHHA(JA)=IM-1-MOD(J+1,2)
          EMT(JA)=EMTJ(J)
 162      CONTINUE
          DO 163 J=6,JM5
          JA=JA+1
          KHLA(JA)=2
          KHHA(JA)=2+MOD(J,2)
          EMT(JA)=EMTJ(J)
 163      CONTINUE
          DO 164 J=6,JM5
          JA=JA+1
          KHLA(JA)=IM-2
          KHHA(JA)=IM-1-MOD(J+1,2)
          EMT(JA)=EMTJ(J)
 164      CONTINUE
C--------------SPREADING OF UPSTREAM VELOCITY-POINT ADVECTION FACTOR----
      JA=0
              DO 171 J=3,5
          JA=JA+1
          KVLA(JA)=2
          KVHA(JA)=IM-1-MOD(J,2)
 171      EM(JA)=EMJ(J)
              DO 172 J=JM4,JM2
          JA=JA+1
          KVLA(JA)=2
          KVHA(JA)=IM-1-MOD(J,2)
 172      EM(JA)=EMJ(J)
              DO 173 J=6,JM5
          JA=JA+1
          KVLA(JA)=2
          KVHA(JA)=2+MOD(J+1,2)
 173      EM(JA)=EMJ(J)
              DO 174 J=6,JM5
          JA=JA+1
          KVLA(JA)=IM-2
          KVHA(JA)=IM-1-MOD(J,2)
 174      EM(JA)=EMJ(J)
C--------------CORIOLIS PARAMETER IN TLL SYSTEM & RELATED CONSTANTS-----
c          TPH=SB-DPH
          DO 180 J=JSTA,JEND
            TLM=WB-TDLM+MOD(J,2)*DLM
c            TPH=TPH+DPH
            TPH=SB+float(J-1)*DPH
            STPH=SIN(TPH)
            CTPH=COS(TPH)
C
            DO 180 I=1,IM
             TLM=TLM+TDLM
             FP=TWOM*(CTPH0*STPH+STPH0*CTPH*COS(TLM))
          F(I,J)=0.5*DT*FP
 180      CONTINUE
C--------------GEOGRAPHIC LAT AND LONG OF TLL GRID POINTS-----------
c          TPH=SB-DPH 
          DO 183 J=JSTA,JEND
C
            TLM=WB-TDLM+MOD(J+1,2)*DLM
c            TPH=TPH+DPH
            TPH=SB+float(J-1)*DPH
            STPH=SIN(TPH)
            CTPH=COS(TPH)
C
           DO 183 I=1,IM
            TLM=TLM+TDLM
            SINPHI=CTPH0*STPH+STPH0*CTPH*COS(TLM)
            GLAT(I,J)=ASIN(SINPHI)
            COSLAM=CTPH*COS(TLM)/(COS(GLAT(I,J))*CTPH0)
     1        -TAN(GLAT(I,J))*TAN(TPH0)
            IF(COSLAM.GT.H1)COSLAM=H1
            IF(COSLAM.LT.-1.)COSLAM=-1.
            FACT=H1
            IF(TLM.GT.D00)FACT=HM1
            GLON(I,J)=-TLM0D*DTR+FACT*ACOS(COSLAM)
  183     CONTINUE
C--------------READ ALBEDO AND MOISTURE AVAILABILTY FILES---------------
      READ(21)ALBC1
      READ(22)ALBC2
      READ(23)ALBC3
      READ(24)ALBC4
c Find Julian day of year to do the temporal interpolation
      juld = julm(idat(1)) + idat(2)
c      juld = julm(idat(2)) + idat(1)
      if(juld.le.32) then
        s1 = 32 - juld
        s2 = juld + 30
        wgt1 = s2/(s1+s2)
        wgt2 = s1/(s1+s2)
        do j = 1,180
        do i = 1,361
        albc(i,j) = wgt1 * albc1(i,j) + wgt2 * albc4(i,j)
        end do
        end do
      else if(juld.le.121.and.juld.gt.32) then
        s1 = 121 - juld
        s2 = juld - 32
        wgt1 = s2/(s1+s2)
        wgt2 = s1/(s1+s2)
        do j = 1,180
        do i = 1,361
        albc(i,j) = wgt1 * albc2(i,j) + wgt2 * albc1(i,j)
        end do
        end do
      else if(juld.le.213.and.juld.gt.121) then
        s1 = 213 - juld
        s2 = juld - 121
        wgt1 = s2/(s1+s2)
        wgt2 = s1/(s1+s2)
        do j = 1,180
        do i = 1,361
        albc(i,j) = wgt1 * albc3(i,j) + wgt2 * albc2(i,j)
        end do
        end do
      else if(juld.le.305.and.juld.gt.213) then
        s1 = 305 - juld
        s2 = juld - 213
        wgt1 = s2/(s1+s2)
        wgt2 = s1/(s1+s2)
        do j = 1,180
        do i = 1,361
        albc(i,j) = wgt1 * albc4(i,j) + wgt2 * albc3(i,j)
        end do
        end do
      else
        s1 = 365 - juld + 32
        s2 = juld - 305
        wgt1 = s2/(s1+s2)
        wgt2 = s1/(s1+s2)
        do j = 1,180
        do i = 1,361
        albc(i,j) = wgt1 * albc1(i,j) + wgt2 * albc4(i,j)
        end do
        end do
      end if
      READ(20)AVLM
C--------INTERPOLATE ALBEDO AND MOISTURE AVAILABILTY TO E GRID-------
      DO 185 J=JSTA,JEND
       DO 185 I=1,IM
        ALBASE(i,j)=11.0
        MXSNAL(i,j)=21.
        IF(SM(I,J) .LT. 0.9) THEN
         ELAT=H90+GLAT(I,J)/DTR
         ELON=H360-GLON(I,J)/DTR
         IF(ELON.GT.H360)ELON=ELON-H360
          ILON1=INT(ELON)
          DIF=ELON-ILON1
          IF(DIF.GT.D5)ILON1=ILON1+1
          IF(ILON1.EQ.D00)ILON1=360
          ILON2=ILON1+1
          ILAT1=INT(ELAT)
          DIF=ELAT-ILAT1
          IF(DIF.GT.D5)ILAT1=MIN(ILAT1+1,179)
          IF(ILAT1.EQ.180.OR.ILAT1.EQ.0)THEN
           WRITE(6,6788)I,J,GLAT(I,J),GLON(I,J),ELAT,ELON
 6788      FORMAT(' I J=',2I4,' GLAT=',E12.5,' GLON=',E12.5
     1     ,' ELAT=',E12.5,' ELON=',E12.5)
           STOP666
          ENDIF
          ILAT2=ILAT1+1
          W1=ELON-ILON1+D5
          IF(W1.LT.D00)W1=W1+H360
          W2=ELAT-ILAT1+D5
          AR1=W1*W2
          AR2=W1*(H1-W2)
          AR3=(H1-W1)*(H1-W2)
          AR4=(H1-W1)*W2
          IF ( (ALBC(ILON2,ILAT2) .GT. 7.) .AND.
     .           (ALBC(ILON2,ILAT1) .GT. 7.) .AND.
     .           (ALBC(ILON1,ILAT1) .GT. 7.) .AND.
     .           (ALBC(ILON1,ILAT2) .GT. 7.) ) THEN
C-----------------------------------------------------------------------
C ALL 4 SURROUNDING POINTS ARE LAND POINTS
              ALBASE(I,J)=AR1*ALBC(ILON2,ILAT2)+
     .                    AR2*ALBC(ILON2,ILAT1)+
     .                    AR3*ALBC(ILON1,ILAT1)+
     .                    AR4*ALBC(ILON1,ILAT2)
              MXSNAL(I,J)=AR1*AVLM(ILON2,ILAT2)+
     .                    AR2*AVLM(ILON2,ILAT1)+
     .                    AR3*AVLM(ILON1,ILAT1)+
     .                    AR4*AVLM(ILON1,ILAT2)
          ELSE
C-----------------------------------------------------------------------
C ONE OR MORE OF THE 4 SURROUNDING POINT ARE NOT LAND POINTS
C TAKE NEAREST NEIGHBOR LAND POINT IN THE FOLLOWING ORDER:
C (ILON2,ILAT2),(ILON2,ILAT1),(ILON1,ILAT1),(ILON1,ILAT2)
C
              ARX=-999.
              IF (ALBC(ILON2,ILAT2) .GT. 7.) THEN
                IF (AR1 .GT. ARX) THEN
                  ARX=AR1
                  ILONX=ILON2
                  ILATX=ILAT2
                ENDIF
              ENDIF
              IF (ALBC(ILON2,ILAT1) .GT. 7.) THEN
                IF (AR2 .GT. ARX) THEN
                  ARX=AR2
                  ILONX=ILON2
                  ILATX=ILAT1
                ENDIF
              ENDIF
              IF (ALBC(ILON1,ILAT1) .GT. 7.) THEN
                IF (AR3 .GT. ARX) THEN
                  ARX=AR3
                  ILONX=ILON1
                  ILATX=ILAT1
                ENDIF
              ENDIF
              IF (ALBC(ILON1,ILAT2) .GT. 7.) THEN
                IF (AR4 .GT. ARX) THEN
                  ARX=AR4
                  ILONX=ILON1
                  ILATX=ILAT2
                ENDIF
              ENDIF
C-----------------------------------------------------------------------
C Use nearest land neighbor:
              IF (ARX .GT. 0.) THEN
                ALBASE(I,J)=ALBC(ILONX,ILATX)
                MXSNAL(I,J)=AVLM(ILONX,ILATX)
              ELSE
C-----------------------------------------------------------------------
C NO SURROUNDING POINTS ARE LAND (ARX=-999):
C SET DEFAULT SNOWFREE ALBEDO=20 PERCENT, DEFAULT SNOWALB=55 PERCENT
                ALBASE(I,J)=20.
                MXSNAL(I,J)=55.
                PRINT *,'AT: LAT,LON',ELAT,ELON,MYPE
                PRINT *,'SNOWFREE ALBEDO SET TO DEFAULT VALUE OF 20%'
                PRINT *,'MAX SNOW ALBEDO SET TO DEFAULT VALUE OF 55%'
C-----------------------------------------------------------------------
C end of ALBASE,MXSNAL land points<4 block
              ENDIF
C-----------------------------------------------------------------------
C end of ALBASE,MXSNAL interpolation block
            ENDIF
C-----------------------------------------------------------------------
C end of land (SM=0) block
           ENDIF
C-----------------------------------------------------------------------
C CONVERT ALBEDO UNITS: PERCENT TO FRACTION
           ALBASE(I,J)=ALBASE(I,J)*D01
           MXSNAL(I,J)=MXSNAL(I,J)*D01
c           ALBEDO(I,J)=AR1*ALBC(ILON2,ILAT2)+AR2*ALBC(ILON2,ILAT1)+
c     1     AR3*ALBC(ILON1,ILAT1)+AR4*ALBC(ILON1,ILAT2)
c           ALBEDO(I,J)=ALBEDO(I,J)*D01
c           AVAILM(I,J)=AR1*AVLM(ILON2,ILAT2)+AR2*AVLM(ILON2,ILAT1)+
c     1     AR3*AVLM(ILON1,ILAT1)+AR4*AVLM(ILON1,ILAT2)
c           IF(AVAILM(I,J).LT.1.E-8)  AVAILM(I,J)=1.E-8
 185  CONTINUE
C--------------DERIVED VERTICAL GRID CONSTANTS--------------------------
      EF4T=.5*DT/CP
      F4Q =   -DT*DTAD
      F4D =-.5*DT*DTAD
      DO 195 L=1,LM
        RDETA(L)=1./DETA(L)
 195  F4Q2(L)=-.25*DT*DTAD/DETA(L)
C--------------BOUNDARY MASKS-------------------------------------------
C      IF(MYPE.EQ.0)THEN
      DO 210 J=JSTA,JEND
      DO 210 I=1,IM
      HBM2(I,J)=0.
      VBM2(I,J)=0.
 210  VBM3(I,J)=0.
C 
      JJS=JSTA
      JJE=JEND
      IF(MYPE.EQ.0)JJS=3
      IF(MYPE.EQ.(NUM_PROCS-1))JJE=JEND-2
      DO 220 J=JJS,JJE
      DO 220 I=2,IM-2+MOD(J,2)
 220  HBM2(I,J)=1.
C 
      DO 230 J=JJS,JJE
      DO 230 I=2,IM-1-MOD(J,2)
 230  VBM2(I,J)=1.
C
      IF(MYPE.EQ.0)JJS=4
      IF(MYPE.EQ.(NUM_PROCS-1))JJE=JEND-3
      DO 235 J=JJS,JJE
      DO 235 I=3-MOD(J,2),IM-2
235   VBM3(I,J)=1.
C      END IF
C--------------TOPOGRAPHY MASKS & MAXIMUM VERTICAL INDICES--------------
      IF(SIGMA)    GO TO 320
C-----------------------------------------------------------------------
Cmp          DO 260 K=KHL00,KHH00
C                             DO 270 L=1,LM
C                             IF(ETA(L+1).LT.RES(K))     GO TO 270
C                             LMH(K)=L
C                             GO TO 260
C 270                         CONTINUE
C 260  CONTINUE
C
C          DO 280 K=KHL00,KHH00
C                             DO 280 L=1,LM
C                             IF(ETA(L+1).LE.RES(K))     GO TO 280
C                             HTM(K,L)=0.
C 280                         CONTINUE
C
C              DO 290 J=2,JM1
C              KHL=KHL1(J)
C              KHH=KHH1(J)
C          DO 290 K=KHL,KHH
C                             DO 290 L=1,LM
C                             IF(ETA(L+1).LE.RES(K))     GO TO 290
C                             VTM(K+KHS,L)=0.
C                             VTM(K-1  ,L)=0.
C                             VTM(K    ,L)=0.
C                             VTM(K+KHN,L)=0.
C 290  CONTINUE
CC
C              DO 300 J=2,JM1
C              KVL=KVL1(J)
c              KVH=KVH1(J)
C          DO 300 K=KVL,KVH
C                             DO 310 L=1,LM
C                             IF(VTM(K,L).GT.0.1)   GO TO 310
C                             LMV(K)=L-1
C                             GO TO 300
C 310                         CONTINUE
C 300  CONTINUE
C--------------SPREADING OF LATITUDE DEPENDENT CONSTANTS----------------
C      allocate(DX(IM,JM))
C      allocate(WPDAR(IM,JM))
C      allocate(CPGFU(IM,JM))
C      allocate(CURV(IM,JM))
C      allocate(FCP(IM,JM))
C      ALLOCATE(FDIV(IM,JM))
C      ALLOCATE(FAD(IM,JM))
C      ALLOCATE(HDAC(IM,JM))
C      ALLOCATE(HDACV(IM,JM))          
 320  DO 200 J=JSTA,JEND
C       DO 200 I=1,IM-1+MOD(J,2)
       DO 200 I=1,IM
        DX(I,J)=DXJ(J)
        WPDAR(I,J)=WPDARJ(J)*HBM2(I,J)
        CPGFU(I,J)=CPGFUJ(J)
        CURV(I,J)=CURVJ(J)
        FCP(I,J)=FCPJ(J)
        FDIV(I,J)=FDIVJ(J)
        FAD(I,J)=FADJ(J)
        HDACV(I,J)=HDACJ(J)
        HDAC(I,J)=HDACJ(J)*1.25*HBM2(I,J)
 200  CONTINUE
C--------------INCREASING DIFFUSION ALONG THE BOUNDARIES----------------
C      DO 201 JA=1,JAM
C      KHL=KHLA(JA)
C      KHH=KHHA(JA)
C      DO 201 K=KHL,KHH
C 201  HDAC(I,J)=HDAC(I,J)* DFC
      JJS=JSTA
      JJE=JEND
      IF(MYPE.EQ.0)JJS=3
      IF(MYPE.EQ.(NUM_PROCS-1))JJE=JEND-2
      DO J=JJS,JJE
        IF (J.LE.5.OR.J.GE.JM4) THEN
          KHH=IM-2+MOD(J,2)
          DO I=2,KHH
           HDAC(I,J)=HDAC(I,J)* DFC
          ENDDO
        ELSE
          KHH=2+MOD(J,2)
          DO I=2,KHH
           HDAC(I,J)=HDAC(I,J)* DFC
          ENDDO
          KHH=IM-2+MOD(J,2)
          DO I=IM-2,KHH
           HDAC(I,J)=HDAC(I,J)* DFC
          ENDDO
        ENDIF
      ENDDO
C-----------------------------------------------------------------------
      DO 330 J=JSTA,JEND
      DO 330 I=1,IM
      DDMPU(I,J)=DDMPUJ(J)*VBM2(I,J)
      DDMPV(I,J)=DDMPVJ(J)*VBM2(I,J)
 330  HDACV(I,J)=HDACV(I,J)*VBM2(I,J)
C--------------INCREASING DIFFUSION ALONG THE BOUNDARIES----------------
      DO J=JJS,JJE
CMEB          KVL3(J)=IM*(J-1)-(J-1)/2+2 2 if j odd, 3 if j even
CMEB          KVL2(J)=IM*(J-1)-J/2+2     2 if j odd or even
CMEB          KVH2(J)=IM*J-(J+1)/2-1     IM-2 if j odd, IM-1 if j even
CMEB          KVH3(J)=IM*J-J/2-2         IM-2 if j odd or even
        IF (J.LE.5.OR.J.GE.JM4) THEN
          KVH=IM-1-MOD(J,2)
          DO I=2,KVH
            DDMPU(I,J)=DDMPU(I,J)*DDFC
            DDMPV(I,J)=DDMPV(I,J)*DDFC
            HDACV(I,J)=HDACV(I,J)* DFC
          ENDDO
        ELSE
          KVH=3-MOD(J,2)
          DO I=2,KVH
            DDMPU(I,J)=DDMPU(I,J)*DDFC
            DDMPV(I,J)=DDMPV(I,J)*DDFC
            HDACV(I,J)=HDACV(I,J)* DFC
          ENDDO
          KVH=IM-1-MOD(J,2)
          DO I=IM-2,KVH
            DDMPU(I,J)=DDMPU(I,J)*DDFC
            DDMPV(I,J)=DDMPV(I,J)*DDFC
            HDACV(I,J)=HDACV(I,J)* DFC
          ENDDO
        ENDIF
      ENDDO
C
CC
CC-MP
CC
CC insert from new GRDETA

C  READ NESDIS AND USAF SNOW MASKS
C
C   READ THE NESDIS-IMS LAND-WATER MASK
C                                          (SEA=0,LAND=1)
      DO J = 1, 1024
         READ(43,54321) (MSKSCVH(I,J),I=1,1024)
54321    FORMAT(80I1)
c        READ(43,'(80I1)') (MSKSCVH(I,J),I=1,1024)
      ENDDO
C
C   READ THE USAF AFGWC LAND/COAST/SEA MASK
C                           (SEA=1,LAND=2,COASTAL-LAND=4,OFFWORLD=9)
      READ(42)  MSKAF
C

C%%%%%%%%%%----SURFACE PARAMETERS----%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C  THE NEXT THREE SUBROUTINE CALLS INTERPOLATE TO THE ETA GRID:
C    A) - HI RES DAILY N.H. 47-KM USAF SNOW, WEEKLY NESDIS 190-KM ICE
C    B) - HI RES DAILY GLOBAL 1-DEG SST (WHICH IS OVERLAID WITH C BELOW)
C    C) - HI RES DAILY WESTERN N.H. 50-KM SST AND GREAT LAKES 14-KM SST
C
C  ROUTINE SNOHIRES WILL RETURN LUSAF=FALSE IF USAF SNOW I/O
C  PROBLEMS OR IF USAF SNOW ANAL TOO OLD.  WHEN SNOHIRES RETURNS
C  LUSAF=FALSE, THEN IT RETURNS SNOW COVER OF 0'S AND 1'S DERIVED FROM
C  WEEKLY 190-KM NESDIS/SAB ANAL.  WHEN SNOHIRES RETURNS LUSAF=TRUE THEN
C  DAILY USAF SNOWDEPTH IN THE RANGE OF 0 TO 10.9 METERS IS RETURNED.
C  OVER SEA, SNOHIRES RETURNS ICE FLAGS OF 0'S AND 1'S FROM NESDIS ANAL.
C
	if(mype.eq.0)write(6,*) 'enter snohires'
        CALL SNOHIRES (SI,SM,GLAT,GLON,LIST,LUSAF,DTR)
C
	if(mype.eq.0)write(6,*) 'enter ssthires'
        CALL SSTHIRES (SST,SM,GLAT,GLON,IDAT,LIST,DTR)
CC PHYSICAL CHECK FOR SST
      DO J=JSTA,JEND
       DO I=1,IM
        IF(SST(I,J).LT.150..OR.SST(I,J).GT.325.)PRINT*,
     +  'ENORMAL SST ',I,J,SST(I,J),MYPE
       END DO
      END DO  	
C
C   NEXT ROUTINE READS IN AND INTERPOLATES A WESTERN N.H. 14-KM SST FLD.
C   THIS 14-KM FLD SPANS A DOMAIN APPROX COVERING (10N-60N, 35W-160W).
C   NOTE: THIS 14-KM SST FLD WAS CREATED IN A PREVIOUS EXTERNAL JOBSTEP
C         AND IS IN REALITY MOSTLY A 50-KM (0.5 DEG) SST WITH 14-KM
C         SST OVERLAID ONLY OVER THE GREAT LAKES
C
      CALL SST14K   (SST,SM,GLAT,GLON)
C
c              TPH=SB-DPH
C
              DO 350 J=JSTA,JEND
C
              TLM=WB-TDLM+MOD(J+1,2)*DLM
c              TPH=TPH+DPH
              TPH=SB+float(J-1)*DPH
              STPH=SIN(TPH)
              CTPH=COS(TPH)
	      KHH=IM-MOD(J,2)
C
              DO 350 I=1,IM
              TLM=TLM+TDLM
              CTLM=COS(TLM)
              STLM=SIN(TLM)
C
              APH=ASIN(STPH0*CTPH*CTLM+CTPH0*STPH)
C
C  IF LUSAF=.FALSE., THEN SI ARE TRADITIONAL FLAGS OF 0'S AND 1'S
C  SO USE TRADITIONAL LOGIC TO SET SNOW DEPTH UP TO 3 MTRS
C  AS A FUNCTION OF LATITUDE AND TERRAIN HEIGHT
C
Cmp      IF ( .NOT. LUSAF)  SI(K)=SI(K)*SIN(APH)*3.
	IF (.NOT.LUSAF.AND.SM(I,J).LT.0.9)SI(I,J)=0.333*SI(I,J)*SIN(APH)
C
C........ IF(FIS(I,J).LT.0.0196.AND.APH/DTR.LT.66.) SI(I,J)=0.
C
C  NEXT SET LOWER BOUNDARY CONDITION OF SUBSOIL TEMPERATURE
C  AS AN ANNUAL FIXED FIELD BASED ON LATITUDE AND TERRAIN HEIGHT
C
C      TG(I,J)=TG0+TGA*COS(APH)+FIS(I,J)/3333.
C
  350 CONTINUE
C
C-----------------------------------------------------------------------
C      ALLOCATE(EPSR(IM,JM))
C      ALLOCATE(GFFC(IM,JM))
C      ALLOCATE(WFK(IM,JM))
      DO 340 J=JSTA,JEND
       DO 340 I=1,IM
        SICE(I,J)=0.
        SNO (I,J)=0.
        IF(ALBEDO(I,J).GT.0.99999999)ALBEDO(I,J)=0.99999999
        IF(SM(I,J).GT.0.9) THEN
C  SEA
        EPSR(I,J)=.97
        GFFC(I,J)=0.
c        WET(I,J)=1.
c        WFK(I,J)=1.
        ALBEDO(I,J)=.06
        ALBASE(I,J)=.06
        MXSNAL(I,J)=0.
        IF(SI (I,J).GT.0.    ) THEN
C  SEA-ICE    
	  SM(I,J)=0.
          SI(I,J)=0.
          SICE(I,J)=1.
C         TG(I,J)=TI0
          GFFC(I,J)=ROI*CI*AKI/DZI
          ALBEDO(I,J)=.60
          ALBASE(I,J)=.60
          MXSNAL(I,J)=0.
        END IF
      ELSE
C  LAND
        EPSR(I,J)=1.0
C  next two lines not needed but kept for impact studies
c        WFK(I,J)=.1125
c        WET(I,J)=AVAILM(I,J)*WFK(I,J)
C
        GFFC(I,J)=ROS*CS*AKS/DZG
        SICE(I,J)=0.
        SNO(I,J)=SI(I,J)*.20
C
cwas        IF(SI(I,J).GT.0.    ) THEN
C  SNOW      (change snowdepth in meters to water equiv in meters)
Cmp          SNO(K)=SI(K)*.10
cwas          SNO(I,J)=SI(I,J)*.20
cwas          WET(I,J)=WFK(I,J)
C.......ALBEDO(I,J)=.65 (now set in model using snowdepth, capped at .55)
cwas        END IF
      END IF
 340  CONTINUE
C--------------OPEN THE WATER POINTS OVER FALSELY FROZEN LAKES---------
c     DO KK=1,21
c       SM(KLAKE(KK))=1.0
c       SI(KLAKE(KK))=0.0
c       SICE(KLAKE(KK))=0.0
c       ALBEDO(KLAKE(KK))=.06
c       GFFC(KLAKE(KK))=0.
c     ENDDO
C
C-----------------------------------------------------------------------
C    MONTHLY GREEN VEG FRACTION, INTERPOLATE TO DAY OF YEAR.
C    THE MONTHLY FIELDS ASSUMED VALID AT 15TH OF MONTH.
C    VALUES RANGE FROM 0.001 TO 1.0 OVER LAND, 0.0 OVER WATER.
C    INTERPOLATE IN SPACE.
C
C    READ FROM A FILE THAT HAS 13 RECORDS, ONE PER MONTH WITH
C    JANUARY OCCURRING BOTH AS RECORD 1 AND AGAIN AS RECORD 13,
C    THE LATTER TO SIMPLIFY TIME INTERPOLATION FOR DAYS
C    BETWEEN DEC 16 AND JAN 15. WE TREAT JAN 1 TO JAN 15
C    AS JULIAN DAYS 366 TO 380 BELOW, I.E WRAP AROUND YEAR.
C *** THE FOLLOWING PART WAS REVISED BY F. CHEN 7/96 TO REFLECT
C        A NEW NESDIS VEGETATION FRACTION PRODUCT (FIVE-YEAR
C        CLIMATOLOGY WITH 0.144 DEGREE RESOLUTION
C        FROM 89.928S, 180W TO 89.928N, 180E)
C
       rewind 38
C ****  DO time interpolation ****
       juld=julm(idat(1))+idat(2)
       juld=juld+2
       if(juld.le.15) juld=juld+365
       mon2=idat(1)
       if(idat(2).gt.15) mon2=mon2+1
       if(mon2.eq.1) mon2=13
       mon1=mon2-1
C **** Assume data valid at 15th of month
       day2=julm(mon2)+15
       day1=julm(mon1)+15
       rday=juld
       wght1=(day2-rday)/(day2-day1)
       wght2=(rday-day1)/(day2-day1)

       CALL GETVEG(MON1,MON2,WGHT1,WGHT2)
       if(mype.eq.0)print*,'after call getveg'
C ** Space interpolation of fpar1 to E grid
       call putveg(TLM0D,TPH0D,DLMD,DPHD,hlato,hlono,
     &             fpar1,sm,sice,wfk)
       if(mype.eq.0)print *,'done with putveg, time = ',
     &    putveg_tim
c 
C    FIRST THE VEGETATION TYPES
      REWIND 30
      READ (30) I1D
c    SET DEFAULT FOR ETA LAND FAR FROM LAND IN GLOBAL
      DO J = JSTA, JEND
        DO I =1, IM
         IVGTPK(I,J)=7
        ENDDO
      ENDDO
      CALL PUTEM(HLATO,HLONO,I1D,SM,SICE,IVGTPK)
C    SECOND THE SOIL TYPES TYPES
       REWIND 31
       READ (31) I1D
C    SET DEFAULT FOR ETA LAND FAR FROM LAND IN GLOBAL
       DO J = JSTA, JEND
         DO I =1, IM
          ISLTPK(I,J)=2
         ENDDO
       ENDDO
       CALL PUTEM(HLATO,HLONO,I1D,SM,SICE,ISLTPK)
       DO J = JSTA, JEND
        DO I = 1,IM
         IF(ISLTPK(I,J).EQ.13) ISLTPK(I,J) = 9
        END DO
       END DO
c    Third the surface slope types
       REWIND 32
       READ (32) I1D
c       print*,'AFTER READING SURFACE SLOPE'
C        SET DEFAULT FOR ETA LAND FAR FROM LAND IN GLOBAL
       DO J = JSTA, JEND
         DO I =1, IM
          ISPTPK(I,J)=1
         ENDDO
       ENDDO
       CALL PUTEM(HLATO,HLONO,I1D,SM,SICE,ISPTPK)
c
C DETERMINE ALBEDO OVER LAND
      DO J=JSTA,JEND
        DO I=1,IM
          IF(SM(I,J).LT.0.9.AND.SICE(I,J).LT.0.9) THEN
C SNOWFREE ALBEDO
            IF ( (SNO(I,J) .EQ. 0.0) .OR.
     .           (ALBASE(I,J) .GE. MXSNAL(I,J) ) ) THEN
              ALBEDO(I,J) = ALBASE(I,J)
            ELSE
C MODIFY ALBEDO IF SNOWCOVER:
C BELOW SNOWDEPTH THRESHOLD...
              IF (SNO(I,J) .LT. SNUP) THEN
                RSNOW = SNO(I,J)/SNUP
                SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP))
C ABOVE SNOWDEPTH THRESHOLD...
              ELSE
                SNOFAC = 1.0
              ENDIF
C CALCULATE ALBEDO ACCOUNTING FOR SNOWDEPTH AND VGFRCK
              ALBEDO(I,J) = ALBASE(I,J)
     .          + (1.0-WFK(I,J))*SNOFAC*(MXSNAL(I,J)-ALBASE(I,J))
            ENDIF
          END IF
        ENDDO
      ENDDO
c
c     if(mype.eq.0) print*,'calling sfcgdas' 
c     CALL SFCGDAS
c     if(mype.eq.0) print*,'done calling sfcgdas'
C
      if(mype.eq.0) then
       if(int(jm/2).gt.jsta .and. int(jm/2).le.jend)then
        write(6,*) 'soil moisture pre SFCEDAS '
        write(6,*) 'SMC:', (SMC(I,JM/2,1),I=1,IM,25)
       end if
       write(6,*) 'calling SFCEDAS!!!!!'
      endif
      CALL SFCEDAS
      if(mype.eq.0) then
       if(int(jm/2).gt.jsta .and. int(jm/2).le.jend)then
        write(6,*) 'soil moisture post SFCEDAS '
        write(6,*) 'SMC:', (SMC(I,JM/2,1),I=1,IM,25)
       end if
      endif
      do j = jsta,jend
       do i = 1,im
        if(sice(i,j).eq.1.0) then
         tg(i,j) = ti0
        end if
       end do	
      end do
c
c  Set vegetation canopy water content equal to zero everywhere for now
      do j=jsta,jend
       do i=1,im
        cmc(i,j) = 0.0
       end do
      end do
C--------------FIND ALL OF THE ELEVATED WATER POINTS--------------------
C        NUMLAK=0
C        DO 366 I=1,IM
C	 NUMLAKI=NUMLAKI+1
C	 DO 366 J=1,JM
c        IF(SM(K).GT.0.9.AND.LMH(K).LT.LM)THEN
C          IF(SM(I,J).GT.0.9.AND.LMHO(I,J).LT.LM)THEN
C           NUMLAKJ=NUMLAKJ+1
C           LAKESI(NUMLAKI,NUMLAKJ)=I
C	   LAKESJ(NUMLAKI,NUMLAKJ)=J
C           WRITE(3,6743)NUMLAKI,LAKESI(NUMLAKI,NUMLAKJ)
C 6743      FORMAT(' LAKE POINT #',I4,' K INDEX IS',I6)
C          ENDIF
C 366    CONTINUE
C--------------MODIFICATION OF MAX WETNESS DUE TO SLOPES----------------
c             DO 341 J=2,JM1
c             KHL=KHL1(J)
c             KHH=KHH1(J)
c         DO 341 K=KHL,KHH
c     IF(SM(K).LT.0.9) THEN
c         SLOPX=ABS((FIS(K+KSE)+FIS(K+KNE)-FIS(K+KSW)-FIS(K+KNW))
c    2         /(4.E0*G*DX(K)))
c         SLOPY=ABS((FIS(K+KNW)+FIS(K+KNE)-FIS(K+KSW)-FIS(K+KSE))
c    2         /(4.E0*G*DY   ))
c         SLOP=AMAX1(SLOPX,SLOPY)
C         WFK(K)=.05S0+20.E0*SLOP
c         WFK(K)=WFK(K)*(1.E0+SLOP)
c         IF(SNO(K).GT.0.    )    WET(K)=WFK(K)
c     END IF
c341  CONTINUE
C-----------------------------------------------------------------------
      DO 360 J=JSTA,JEND
       DO 360 I=1,IM
        GFFC(I,J)=GFFC(I,J)*HBM2(I,J)
 360    EPSR(I,J)=EPSR(I,J)*HBM2(I,J)
C-----------------------------------------------------------------------
      CALL TABLE(PTBL,TTBL,PT
     1,          RDQ,RDTH,RDP,RDTHE,PL,THL,QS0,SQS,STHE,THE0)
      CALL TABLEQ(TTBLQ,RDPQ,RDTHEQ,PLQ,THL,STHEQ,THE0Q)
C-----------------------------------------------------------------------
C
      if(mype.eq.0) print *, 'begin write of nhb file '
      nfcsto = 13
      nbco = 16
      ii=imo/2
      jj=jmo/2
      if(mype.eq.0) then
      if(jj.gt.jsta.and.jj.le.jend)then
      print*,'writing sfc fields at ',ii,jj
      print*,'lmho=',lmho(ii,jj)
      print*,'lmvo=',lmvo(ii,jj)
      print*,'hbm2=',hbm2(ii,jj)
      print*,'vbm2=',vbm2(ii,jj)
      print*,'sm=',sm(ii,jj)
      print*,'sice=',sice(ii,jj)
C      do l=1,lm
C       print*,'htmo=',htmo(ii,jj,l)
C      end do 
C      do l=1,lm
C       print*,'vtmo=',vtmo(ii,jj,l)
C      end do
      print*,'tg=',tg(ii,jj)
      print*,'gffc=',gffc(ii,jj)
      print*,'glat=',glat(ii,jj)
      print*,'glon=',glon(ii,jj)
      print*,'wfk=',wfk(ii,jj)
      print*,'epsr=',epsr(ii,jj)
      print*,'albedo=',albedo(ii,jj)
      print*,'albase=',albase(ii,jj)
      print*,'mxsnal=',MXSNAL(ii,jj)
      print*,'sst=',sst(ii,jj)
      print*,'hdac=',hdac(ii,jj)
      print*,'hdacv=',hdacv(ii,jj)
      end if
      endif
      do j=jsta,jend 
       do i=1,imo
        if(albedo(i,j).le.0.)print*,'bad albedo',i,j,albedo(i,j)
        if(albase(i,j).le.0.)print*,'bad albase',i,j,albase(i,j)
       end do
      end do
      IF(MYPE.EQ.0)THEN
       WRITE(NHIBU)
     1 NFCSTO,NBCO,LIST
     2,DT,IDTAD,SIGMA
     3,KHLA,KHHA,KVLA,KVHA,KHL2,KHH2,KVL2,KVH2  
C
      END IF
c      ALLOCATE(LMHO(IMO,JSTA:JEND))
c      ALLOCATE(LMVO(IMO,JSTA:JEND))
c      ALLOCATE(HTMO(IMO,JSTA:JEND,LMO))
c      ALLOCATE(VTMO(IMO,JSTA:JEND,LMO))
      CALL ICOLLECT_LOC(LMHO,IDUM2D)
      IF(MYPE.EQ.0)WRITE(NHIBU)idum2d
      CALL ICOLLECT_LOC(LMVO,IDUM2D)
      IF(MYPE.EQ.0)WRITE(NHIBU)idum2d 
c       WRITE(NHIBU)ioneta
      CALL COLLECT_DUM(HBM2,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(VBM2,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(VBM3,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(SM,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(SICE,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
C 
      DO L=1,LM
        CALL COLLECT_LOC(HTMO(1,JSTA,L),dum2d)
        IF(MYPE.EQ.0)WRITE(NHIBU)((dum2d(I,J),I=1,IM),J=1,JM)
      END DO
C
      DO L=1,LM
        CALL COLLECT_LOC(VTMO(1,JSTA,L),dum2d)
        IF(MYPE.EQ.0)WRITE(NHIBU)((dum2d(I,J),I=1,IM),J=1,JM)
      END DO
c      DEALLOCATE(LMHO)
c      DEALLOCATE(LMVO)
c      DEALLOCATE(HTMO)
c      DEALLOCATE(VTMO)
C      
      IF(MYPE.EQ.0)WRITE(NHIBU)
     1 DY,CPGFV,EN,ENT,R,PT,TDDAMP
     2,F4D,F4Q,EF4T,DETA,RDETA,AETA,F4Q2,ETA,DFL
     3,EM,EMT
C
      CALL COLLECT_DUM(DX,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(WPDAR,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(CPGFU,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(CURV,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(FCP,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(FDIV,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(FAD,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(F,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(DDMPU,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(DDMPV,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(GLAT,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)PT2,dum2d
      CALL COLLECT_DUM(GLON,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      IF(MYPE.EQ.0)THEN
       WRITE(NHIBU)PLQ,RDPQ,RDTHEQ,STHEQ,THE0Q
       WRITE(NHIBU)
     1 ROS,CS,DS,ROI,CI,DI
     2,PL,THL,RDQ,RDTH,RDP,RDTHE
     3,DETA2,AETA2,DFRLG
     4,QS0,SQS,STHE,THE0
      END IF
c      CALL COLLECT_DUM(WFK,dum2d)
      CALL COLLECT_LOC(MXSNAL,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(EPSR,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(TG,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(GFFC,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
CC PHYSICAL CHECK FOR SST
      DO J=JSTA,JEND
       DO I=1,IM
        IF(SST(I,J).LT.150..OR.SST(I,J).GT.325.)PRINT*,
     +  'ENORMAL SST ',I,J,SST(I,J)
       END DO
      END DO  	 
      CALL COLLECT_DUM(SST,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
c      CALL COLLECT_DUM(ALBEDO,dum2d)
      CALL COLLECT_LOC(ALBASE,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(HDAC,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      CALL COLLECT_DUM(HDACV,dum2d)
      IF(MYPE.EQ.0)WRITE(NHIBU)dum2d
      IF(MYPE.EQ.0)WRITE(NHIBU)TTBLQ
      IF(MYPE.EQ.0)WRITE(NHIBU)
     1 PTBL,TTBL
     2,R1,PT1,TSPH
     3,WBD,SBD,TLM0D,TPH0D,DLMD,DPHD,CMLD,DP30,X1P,Y1P,IXM,IYM
     4,DETA1,AETA1,ETA1
C-----------------------------------------------------------------------
C    Time to write the soil model stuff to nhibu
c
C    FIRST THE VEGETATION TYPES
c       REWIND 30
c       READ (30) I1D
C        SET DEFAULT FOR ETA LAND FAR FROM LAND IN GLOBAL
c       DO J = JSTA, JEND
c         DO I =1, IM
c          IONETA(I,J)=7
c         ENDDO
c       ENDDO
c       CALL PUTEM(HLATO,HLONO,I1D,SM,SICE,IONETA)
c       CALL ICOLLECT_DUM(IONETA,IDUM2D)
c       IF(MYPE.EQ.0)WRITE(NHIBU) IDUM2D
       CALL ICOLLECT_DUM(IVGTPK,IDUM2D)
       IF(MYPE.EQ.0)WRITE(NHIBU) IDUM2D 
C
c       WRITE(6,*) 'VEG TYPE at',ii,jj,'=', IONETA(ii,jj)
C    SECOND THE SOIL TYPES TYPES
c       REWIND 31
c       READ (31) I1D
C        SET DEFAULT FOR ETA LAND FAR FROM LAND IN GLOBAL
c       DO J = JSTA, JEND
c         DO I =1, IM
c          IONETA(I,J)=2
c         ENDDO
c       ENDDO
c       CALL PUTEM(HLATO,HLONO,I1D,SM,SICE,IONETA)
c       DO J = JSTA, JEND 
c        DO I = 1,IM
c	 IF(IONETA(I,J).EQ.13) IONETA(I,J) = 9
c        END DO
c       END DO
c       CALL ICOLLECT_DUM(IONETA,IDUM2D)
c       IF(MYPE.EQ.0)WRITE(NHIBU) IDUM2D
       CALL ICOLLECT_DUM(ISLTPK,IDUM2D)
       IF(MYPE.EQ.0)WRITE(NHIBU) IDUM2D
c    Third the surface slope types
c       REWIND 32
c       READ (32) I1D
c       print*,'AFTER READING SURFACE SLOPE'
C        SET DEFAULT FOR ETA LAND FAR FROM LAND IN GLOBAL
c       DO J = JSTA, JEND 
c         DO I =1, IM
c          IONETA(I,J)=1
c         ENDDO
c       ENDDO
c       CALL PUTEM(HLATO,HLONO,I1D,SM,SICE,IONETA)
c       CALL ICOLLECT_DUM(IONETA,IDUM2D)
c       IF(MYPE.EQ.0)WRITE(NHIBU) IDUM2D
       CALL ICOLLECT_DUM(ISPTPK,IDUM2D)
       IF(MYPE.EQ.0)WRITE(NHIBU) IDUM2D
c
        CALL COLLECT_DUM(wfk,DUM2D)
         IF(MYPE.EQ.0)WRITE(NHIBU) DUM2D
C************* End of time and space interpolation of FPAR ***
c
c    Fifth the soil depth
c    Variable is SLDPTH
         IF(MYPE.EQ.0)write(nhibu) SLDPTH
c    Sixth the rooting density (not yet used in the FCST and POST)
c    Variable is RDTPTH
         IF(MYPE.EQ.0)write(nhibu) RTDPTH
C-----------------------------------------------------------------------
C      deallocate(DX)
C      deallocate(WPDAR)
C      deallocate(CPGFU)
C      deallocate(CURV)
C      deallocate(FCP)
C      DEALLOCATE(FDIV)
C      DEALLOCATE(FAD)
C      DEALLOCATE(F)
C      DEALLOCATE(DDMPU)
C      DEALLOCATE(DDMPV)
C      DEALLOCATE(EPSR)
C      DEALLOCATE(GFFC)
C      DEALLOCATE(WFK)
C      DEALLOCATE(HDAC)
C      DEALLOCATE(HDACV)
                             RETURN
                             END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      BLOCK DATA PATCH
C-----------------------------------------------------------------------
      PARAMETER (NP=21)
C-----------------------------------------------------------------------
      COMMON /SLAKE/KLAKE(NP)
C
      DATA KLAKE/24924,25074,25224,25225,25374,25375,25524,25525,25526
     1,          25674,25675,25676,25824,25825,25975,25976,26125,26126
     2,          26275,26276,26426/
C
      END
