      pARAMETER (NN=3000000,KINLEN=3000000,ILENF=3000000,KWORD=32)
      PARAMETER (NX=201,NY=101,NLEV=16,NSFC=12) !20-31Mar2011
!      PARAMETER (NX=201,NY=101,NLEV=15,NSFC=16) !for the rest
!      PARAMETER (NX=201,NY=101,NLEV=11,NSFC=16)
      PARAMETER (NCYCLE=NLEV*4+NSFC) !20-31Mar2011
!      PARAMETER (NCYCLE=NLEV*5+12) !for the rest
!      PARAMETER (NCYCLE=NLEV*5+NSFC+1)
C
      DIMENSION IGRIB(NN)
      DIMENSION FPDATA(ILENF),FPDATA2D(NX,NY)
      DIMENSION KSEC0(2), KSEC1(1024), KSEC2(1024), PSEC2(128)
     1         ,KSEC3(2),PSEC3(2),KSEC4(512),KGDS(200)
      dimension PR(NLEV),amin(NLEV),amax(NLEV)
      INTEGER ILEN,KEOF
C
      CHARACTER CARG1*6,CARG2*2, CARG3*3
      CHARACTER  CED*4
      CHARACTER  FFNAME*128,FTMP*64,FFNAME2*128
      character fn_grib*128,outdir*128
      character atime*11
      

      DIMENSION IDAT(3)
      dimension month(12)
      dimension HD(NX,NY,NLEV),QD(NX,NY,NLEV),TD(NX,NY,NLEV)
      dimension UD(NX,NY,NLEV),VD(NX,NY,NLEV),SLPD(NX,NY,NSFC)

      real esat,es
C--------------cspir----------
        data month/31,28,31,30,31,30,31,31,30,31,30,31/
C--------------cspir----------
        data PR/1000.,925.,850.,700.,500.,400.,300.,250.
     $          ,200.,150.,100.,70.,50.,30.,20.,10./
!        data PR/1000.,925.,850.,700.,500.,400.,300.,250.
!     $          ,200.,150.,100.,50.,30.,20.,10./
!        data PR/1000.,925.,850.,700.,500.,400.,300.,250.
!     $          ,200.,150.,100./
c
      common /estab/esat(15000:45000),es(15000:45000)
c
C
C --Check the command input
C
cpk      read(5,'(a)')   FN_GRIB
cpk      read(5,'(a)')   outdir
cpk      read(5,'(3I2)') idat
cpk      read(5,'(I2)')  ihrst      
cpk      read(5,'(I3)')  ihr     


      CALL GETARG (1,fn_grib)
      CALL GETARG (2,outdir)
      CALL GETARG (3,CARG1)
      READ(CARG1,'(3I2)')IDAT
      CALL GETARG (4,CARG2)
      READ(CARG2,'(I2)')IHRST
      CALL GETARG (5,CARG3)
      READ(CARG3,'(I3)')IHR

      print*,fn_grib
      print*,outdir
      print*,idat
      print*,ihrst
      print*,ihr

      call es_ini

C
C --End command input
C
C
      print*,'Decoding file: ',FN_GRIB(1:INDEX(FN_GRIB,' '))

C
      do JJ=1,NY
      do II=1,NX
      do L = 1, NLEV
         HD(ii,jj,L)=-999.9
         UD(ii,jj,L)=-999.9
         VD(ii,jj,L)=-999.9
         QD(ii,jj,L)=-999.9
         TD(ii,jj,L)=-999.9
      enddo
      do L = 1, NSFC
         SLPD(ii,jj,L)=-999.9
      enddo
      enddo
      enddo
C
      I=0 

      CALL ROPEN(FN_GRIB,IRET) 
      
      IF (IRET.EQ.1) THEN
      WRITE(*,*)'ERROR ON OPEN '//FN_GRIB
      GOTO 200 
      END IF

 100  CALL READNEWG(IGRIB,ILEN,KEOF)
      CALL SWAP32(IGRIB,NN)

      IF (KEOF.EQ.1)THEN
      WRITE(*,*)'END OF FILE '//FN_GRIB
       GOTO 200
      END IF
      IF (KEOF.EQ.2) GOTO 100
      KRET=0
      JLENF=ILENF
      CALL GRIBEX (KSEC0,KSEC1,KSEC2,PSEC2,KSEC3,PSEC3,KSEC4,
     &             FPDATA,JLENF,IGRIB,KINLEN,KWORD,'I',KRET)
C
C     *** Avoid some strange records WDIR(34), VGRD(31) ***
C
      if (KSEC1(6).eq.31 .or. KSEC1(6).eq.34) GOTO 100
C---------------cspir-------------------
          ihrecm=ihr
          iday=idat(1)
          imonth=idat(2)
          iyear=idat(3)
7788   if(ihrecm.GE.24) then
          ihrecm=ihrecm-24
          iday=iday+1
          goto 7788
         endif
        if (iday .gt. month(idat(2))) then
          iday = iday - month(idat(2))
          imonth = imonth + 1
        endif
        if (idat(2) .eq. 13) then
          imonth=1
          iyear = iyear + 1
          if (iyear.ge.10) iyear = iyear - 10
        endif
C---------------cspir-------------------
      IF (KSEC1(10).eq.iyear .and.
     &    KSEC1(11).eq.imonth .and.
     &    KSEC1(12).eq.iday .and.
     &    KSEC1(16).eq.ihrst .and.        !!CAREFULL
     &    KSEC1(13).eq.ihrecm ) then      !!CAREFULL

          print*,'Valid date of data: ',KSEC1(10),'/',KSEC1(11),'/',
     &            KSEC1(12),' on ',KSEC1(13),' at ',KSEC1(16)

      ELSE
Ctas          print*,'Check the input dates'
          GOTO 100
      ENDIF
      
 101  CONTINUE
c  *** KSEC1(7) 100 is isobaric, 1 is sfc, 112 is soil ***
      IF((KSEC1(7).NE.100).AND.(KSEC1(7).NE.1).AND.
     $   (KSEC1(7).NE.112)) GO TO 100

      L=9999
      DO LV=1,NLEV
         IF((KSEC1(7).EQ.100).AND.(KSEC1(8).EQ.PR(LV))) L=LV
      ENDDO
      IF ((KSEC1(7).EQ.1).AND.(KSEC1(8).EQ.0)) L=0
      IF ((KSEC1(7).EQ.112).AND.(KSEC1(8).EQ.0))  L=1
      IF ((KSEC1(7).EQ.112).AND.(KSEC1(8).EQ.7))  L=2
      IF ((KSEC1(7).EQ.112).AND.(KSEC1(8).EQ.28)) L=3
      IF ((KSEC1(7).EQ.112).AND.(KSEC1(8).EQ.100))L=4
C
      IF (L.EQ.9999) GOTO 100
C
      KRET=0
      JLENF=KSEC2(2)*KSEC2(3)  ! nx*ny
      CALL GRIBEX (KSEC0,KSEC1,KSEC2,PSEC2,KSEC3,PSEC3,KSEC4,
     &             FPDATA,JLENF,IGRIB,KINLEN,KWORD,'D',KRET)
C
C     Convert 1D to 2D
C
      indx=0
      do jj=KSEC2(3),1,-1
      do ii=1,KSEC2(2)
         indx=indx+1
         FPDATA2D(ii,jj)=FPDATA(indx)
      enddo
      enddo
      print*,'***************************************************'
C
C 129 is Geopotential (m**2 s**-2)
      if((KSEC1(6).eq.129).AND.(KSEC1(7).eq.100)) then
      I=I+1
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
          HD(ii,jj,l)=fpdata2d(ii,jj)/9.81
        end do
        end do
        print*,I,' Geopotential at',KSEC1(8),L,'level'
        print*,'NX=',KSEC2(2),'  NY=',KSEC2(3)
      endif
C
C 156 is Height (geopotential) (m)
      if(KSEC1(6).eq.156) then
      print*,'Height',KSEC1(8),L
      I=I+1
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
          HD(ii,jj,l)=fpdata2d(ii,jj)
        end do
        end do
      end if
C
C 131 is U-velocity (m s**-1)
      if(KSEC1(6).eq.131) then
      I=I+1
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
          UD(ii,jj,l)=fpdata2d(ii,jj)
        end do
        end do
        print*,I,' U-velocity at',KSEC1(8),L,'level'
        print*,'NX=',KSEC2(2),'  NY=',KSEC2(3)
      end if
C
C 132 is V-velocity (m s**-1)
      if(KSEC1(6).eq.132) then
      I=I+1
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
          VD(ii,jj,l)=fpdata2d(ii,jj)
        end do
        end do
        print*,I,' V-velocity at',KSEC1(8),L,'level'
        print*,'NX=',KSEC2(2),'  NY=',KSEC2(3)
      end if
C
C 133 is Specific humidity (Kg/Kg)
      if(KSEC1(6).eq.133) then
      print*,'Specific Humidity',KSEC1(8),L
      I=I+1
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
          QD(ii,jj,l)=fpdata2d(ii,jj)
        end do
        end do
        print*,I,' Specific Humidity at',KSEC1(8),L,'level'
        print*,'NX=',KSEC2(2),'  NY=',KSEC2(3)
      end if
C
C130  is Temperature (K)
      if(KSEC1(6).eq.130) then
      I=I+1
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
          TD(ii,jj,l)=fpdata2d(ii,jj)
        end do
        end do
        print*,I,' Temperature at',KSEC1(8),L,'level'
        print*,'NX=',KSEC2(2),'  NY=',KSEC2(3)
      end if
C
C 157 is Relative humidity (%)
cpk      if(KSEC1(6).eq.157) then
cpk      print*,'Relative humidity',KSEC1(8),L
cpk      I=I+1
cpk        do jj=1,KSEC2(3)
cpk        do ii=1,KSEC2(2)
cpk          if (fpdata2d(ii,jj).gt.100.0)THEN
cpk              print*,'!!!! RH > 100% in ',ii,jj,fpdata2d(ii,jj)
cpk              fpdata2d(ii,jj)=100.0
cpk          endif
cpk          if ((fpdata2d(ii,jj).lt.0.0).and.(IHR.eq.0)) THEN
cpk              print*,'!!!! RH < 0% in ',ii,jj,fpdata2d(ii,jj)
cpk              fpdata2d(ii,jj)=0.0
cpk          endif
cpk          QD(ii,jj,l)=fpdata2d(ii,jj)
cpk        end do
cpk        end do
cpk      end if
C
C 172 is Land/Sea Mask (0/1)
      if(KSEC1(6).eq.172) then
      I=I+1
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
          SLPD(ii,jj,1)=fpdata2d(ii,jj)
        end do
        end do
        print*,I,' Land/Sea Mask (0/1) at',KSEC1(8),L,'level'
        print*,'NX=',KSEC2(2),'  NY=',KSEC2(3)
      end if
C
C 151 is Mean Sea-Level Pressure (Pa)
      if(KSEC1(6).eq.151) then
      I=I+1
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
          SLPD(ii,jj,2)=fpdata2d(ii,jj)
        end do
        end do
        print*,I,'  Mean Sea-Level Pressure at',KSEC1(8),L,'level'
        print*,'NX=',KSEC2(2),'  NY=',KSEC2(3)
      end if
C
C 134 is Surface Pressure (Pa) !!dummy
      if(KSEC1(6).eq.151) then
      I=I+1
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
cs          SLPD(ii,jj,3)=fpdata2d(ii,jj)
          SLPD(ii,jj,3)=-999.9
        end do
        end do
        print*,I,' Surface Pressure at',KSEC1(8),L,'level'
        print*,'NX=',KSEC2(2),'  NY=',KSEC2(3)
      end if
C
C 129 is Geopotential (m2/s2)
      if((KSEC1(6).eq.129).AND.(KSEC1(7).eq.1)) then
      print*,'SFC.Geopotential ',KSEC1(8),L
      I=I+1
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
          SLPD(ii,jj,4)=fpdata2d(ii,jj)/9.81
        end do
        end do
        print*,I,' Surface Geopotential at',KSEC1(8),L,'level'
        print*,'NX=',KSEC2(2),'  NY=',KSEC2(3)
      end if
C 139,170,183,236 is Soil Temperature (K) !!dummy!!
      if(KSEC1(7).eq.112) then
      if ((KSEC1(6).eq.139).OR.(KSEC1(6).eq.170).OR.
     $    (KSEC1(6).eq.183).OR.(KSEC1(6).eq.236)) then
      I=I+1
!      I=I+4
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
          SLPD(ii,jj,3+2*L)=fpdata2d(ii,jj)
cpk          SLPD(ii,jj,3+2*L+1)=0.14  !SMC
cpk          SLPD(ii,jj,3+2*L+2)=fpdata2d(ii,jj)
cpk          SLPD(ii,jj,3+2*L+3)=0.13  !SMC
        end do
        end do
        print*,I,'Soil Temp.',KSEC1(8),L,'level'
        print*,'NX=',KSEC2(2),'  NY=',KSEC2(3)
        print*,'Index ',3+2*L
      end if
      end if

C 140,171,184,237 is Soil Moisture(39,40,41,42)  !!dummy!!
      if(KSEC1(7).eq.112) then
      if ((KSEC1(6).eq.039).OR.(KSEC1(6).eq.040).OR.
     $    (KSEC1(6).eq.041).OR.(KSEC1(6).eq.042)) then
      I=I+1
!      I=I+4
        do jj=1,KSEC2(3)
        do ii=1,KSEC2(2)
          SLPD(ii,jj,4+2*L)=fpdata2d(ii,jj)
          if (SLPD(ii,jj,4+2*L).lt.0.) SLPD(ii,jj,4+2*L)=0.
        end do
        end do
        print*,I,'Soil Water ',KSEC1(8),L,'level'
        print*,'NX=',KSEC2(2),'  NY=',KSEC2(3)
        print*,'Index ',4+2*L
      end if
      end if

C     WRITE(*,*) ' *********** Count: ', I,NCYCLE,' *****************'
      print*,'Count ',I,NCYCLE

C
      IF(MOD(I,NCYCLE).EQ.0) THEN

c ------- convert rel.hum. to mixing ratio ----
cpk      do kk=1,nlev
cpk      do jj=1,KSEC2(3)
cpk      do ii=1,KSEC2(2)
cpk       it=TD(ii,jj,kk)*100
cpk       it=min(45000,max(15000,it))
cpk       xe=esat(it)
cpk       xmrsat=0.00622*xe/(pr(kk)-xe)
cpk       QD(ii,jj,kk)=QD(ii,jj,kk)*xmrsat
cpk      enddo
cpk      enddo
cpk      enddo

      if (IDAT(3) .ge. 100) IDAT(3) = IDAT(3) - 100
      print*,'IDAT= ',IDAT,' IHRST= ',IHRST,' IHR= ',IHR

      n=index(outdir,' ')-1
      write(atime,'(i2.2,i2.2,i2.2,i2.2,i3.3)')
     &idat(3),idat(2),idat(1),ihrst,ihr
      ffname=outdir(1:n)//'/'//atime//'.ETA_ecm'
      ffname2=outdir(1:n)//'/'//atime//'.ecm_wam'
      print*,'Output ',ffname(1:64)

C
      WRITE(CED,'(I4.4)')IHR
      FTMP=outdir(1:n)//'/minmax-ecm.'
      ISTR2=INDEX(FTMP,' ')-1
      FTMP(ISTR2+1:ISTR2+4)=CED(1:4)
      FTMP(ISTR2+5:ISTR2+8)='.dat'
      OPEN (12,file=FTMP(1:ISTR2+8),form='formatted')
 
      write(12,*)'IDAT: ',IDAT,' IHRST: ',IHRST,' IHR: ',IHR
      DO L=1,NLEV
         write(12,*)'======================================='
         write(12,*)'Level: ',L,' --> ',PR(L),' hPa'
         amin(L) = 100000.
         amax(L) = -100000.
         do jj=1,KSEC2(3)
         do ii=1,KSEC2(2)
           if (HD(ii,jj,l) .ge. amax(l)) amax(l)=HD(ii,jj,l)
           if (HD(ii,jj,l) .le. amin(l)) amin(l)=HD(ii,jj,l)
         enddo
         enddo
         write(12,*)'G.H : amin,amax',amin(l),amax(l)
C
         amin(l) = 100000.
         amax(l) = -100000.
         do jj=1,KSEC2(3)
         do ii=1,KSEC2(2)
           if (UD(ii,jj,l) .ge. amax(l)) amax(l)=UD(ii,jj,l)
           if (UD(ii,jj,l) .le. amin(l)) amin(l)=UD(ii,jj,l)
         enddo
         enddo
         write(12,*)'U : amin,amax',amin(l),amax(l)
C
         amin(l) = 100000.
         amax(l) = -100000.
         do jj=1,KSEC2(3)
         do ii=1,KSEC2(2)
           if (VD(ii,jj,l) .ge. amax(l)) amax(l)=VD(ii,jj,l)
           if (VD(ii,jj,l) .le. amin(l)) amin(l)=VD(ii,jj,l)
         enddo
         enddo
         write(12,*)'V : amin,amax',amin(l),amax(l)
C
         amin(l) = 100000.
         amax(l) = -100000.
         do jj=1,KSEC2(3)
         do ii=1,KSEC2(2)
           if (QD(ii,jj,l) .ge. amax(l)) amax(l)=QD(ii,jj,l)
           if (QD(ii,jj,l) .le. amin(l)) amin(l)=QD(ii,jj,l)
         enddo
         enddo
         write(12,*)'Q : amin,amax',amin(l),amax(l)
C
         amin(l) = 100000.
         amax(l) = -100000.
         do jj=1,KSEC2(3)
         do ii=1,KSEC2(2)
           if (TD(ii,jj,l) .ge. amax(l)) amax(l)=TD(ii,jj,l)
           if (TD(ii,jj,l) .le. amin(l)) amin(l)=TD(ii,jj,l)
         enddo
         enddo
         write(12,*)'T : amin,amax',amin(l),amax(l)

         ENDDO

         write(12,*)
         write(12,*)'SURFACE FIELDS'
         DO L=1,NSFC
            amin(l) = 1000000.
            amax(l) = -1000000.
            do jj=1,KSEC2(3)
            do ii=1,KSEC2(2)
              if (SLPD(ii,jj,l).ge.amax(l))amax(l)=SLPD(ii,jj,l)
              if (SLPD(ii,jj,l).le.amin(l))amin(l)=SLPD(ii,jj,l)
            enddo
            enddo
            write(12,*)L,'SFC FIELD : amin,amax',amin(l),amax(l)
         ENDDO

      close(12)
C
      OPEN (UNIT=16,FILE=FFNAME,STATUS='UNKNOWN',
     &      FORM='UNFORMATTED')
      WRITE(16)NLEV
      WRITE(16)UD
      WRITE(16)VD
      WRITE(16)HD
      WRITE(16)QD
      WRITE(16)PR
      WRITE(16)SLPD
      CLOSE (16)
      print*,'Write to file: ',FFNAME  !!(1:n+20)
C  --wave coupling--
      IF (IHR.EQ.0) THEN
         OPEN (UNIT=18,FILE=FFNAME2,STATUS='UNKNOWN',
     &         FORM='UNFORMATTED')
         WRITE(18) ((UD(ii,jj,2), ii=1,KSEC2(2)), jj=1,KSEC2(3))
         WRITE(18) ((VD(ii,jj,2), ii=1,KSEC2(2)), jj=1,KSEC2(3))
         WRITE(18) ((SLPD(ii,jj,1), ii=1,KSEC2(2)), jj=1,KSEC2(3))
         CLOSE (18)
         print*,'Write to file: ',FFNAME2
      ENDIF
c ------ saving geographical info -------
      do ii=1,200
        KGDS(ii)=KSEC2(ii)
      enddo
      itmp=KGDS(4)
      KGDS(4)=KGDS(7)
      KGDS(7)=itmp    

      print*,'writing kgds',(kgds(ii),ii=1,14)

      OPEN (UNIT=17,FILE="../../data/prep/gdsinfo.ETA_ecm",
     $form='unformatted',access='sequential')
      write(17)KGDS      
      CLOSE(17)

      END IF
C
      GOTO 100 
     
 200  CONTINUE
      CALL RCLOSE() 
C
 1001 CONTINUE
      print*,'Number of messages: ',I
C
      STOP
      END
C
C###################################################################
       SUBROUTINE SWAP32(A,N)
C
C      REVERSE ORDER OF BYTES IN INTEGER*4 WORD, or REAL*4
C
       INTEGER*4   A(N)
C
       CHARACTER*1 JTEMP(4)
       CHARACTER*1 KTEMP
C
       EQUIVALENCE (JTEMP(1),ITEMP)
C
       SAVE
C
       DO I = 3,N
         ITEMP    = A(I)
         KTEMP    = JTEMP(4)
         JTEMP(4) = JTEMP(1)
         JTEMP(1) = KTEMP
         KTEMP    = JTEMP(3)
         JTEMP(3) = JTEMP(2)
         JTEMP(2) = KTEMP
         A(I)     = ITEMP
      ENDDO
      RETURN
      END
C######################################################
      SUBROUTINE HIBUHR(KYEAR,KMONTH,KDAY,KHOUR,IHR)
C     ******************************************************************
C     *                                                                *
C     *  OVO JE RUTINA ZA IZRACUNAVANJE PROTEKLOG VREMENA U CASOVIMA   *
C     *  POCEVSI OD 1.1.1978                                           *
C     *                                                                *
C     *  PROGRAMER: Z.JANJIC,  SHMZ,  BEOGRAD,  1981                   *
C     *                                                                *
C     ******************************************************************
      DIMENSION MONTH(12)
      DATA MONTH/31,28,31,30,31,30,31,31,30,31,30,31/
C
      IDAY=0
      KYEAR1=KYEAR
C
      if (KYEAR1.lt.30) KYEAR1=KYEAR1+100  ! ioannis
C
      DO 100 K=78,KYEAR1-1
      IF(MOD(K,4).EQ.0) IDAY=IDAY+1
 100  IDAY=IDAY+365
C
      IF(MOD(KYEAR,4).EQ.0) MONTH(2)=29
C
      IF(KMONTH.EQ.1) GO TO 120
C
      DO 110 K=1,KMONTH-1
 110  IDAY=IDAY+MONTH(K)
C
 120  IHR=(IDAY+KDAY)*24+KHOUR
C
      RETURN
      END
C############################################################
      SUBROUTINE PROSELECT(KLON,KLAT,FIELD2D,FIELD1D)
      real FIELD2D(KLON,KLAT),IBITMAP(KLON,KLAT),FIELD1D(900000)
c      real FIELD2D(321,141),IBITMAP(321,141),FIELD1D(900000)
c      KLON=321
c      KLAT=141
      ZR=(1._8-SQRT(2._8/3._8))/2._8
      Z=REAL(KLON,KIND=8)*ZR
      ILON=NINT(Z)
      Z=REAL(KLAT,KIND=8)*ZR
      ILAT=NINT(Z)
      print*,ILON,ILAT

      DO JLAT=1,KLAT
       DO JLON=1,KLON
         FIELD2D(JLON,JLAT)=-999.9
         IBITMAP(JLON,JLAT)=0
       ENDDO
      ENDDO

      indx=0
C ***** South Boundary
      DO JLAT=1,ILAT
       DO JLON=1,KLON
        indx=indx+1
        FIELD2D(JLON,JLAT)=FIELD1D(indx)
        IBITMAP(JLON,JLAT)=1
       ENDDO
      ENDDO

C ***** Lateral Boundaries
      DO JLAT=ILAT+1,KLAT-ILAT
       DO JLON=1,ILON
        indx=indx+1
        FIELD2D(JLON,JLAT)=FIELD1D(indx)
        IBITMAP(JLON,JLAT)=1
       ENDDO

       DO JLON=KLON-ILON+1,KLON
        indx=indx+1
        FIELD2D(JLON,JLAT)=FIELD1D(indx)
        IBITMAP(JLON,JLAT)=1
       ENDDO
      ENDDO

C ***** North Boundary
      DO JLAT= KLAT-ILAT+1,KLAT
       DO JLON=1,KLON
        indx=indx+1
        FIELD2D(JLON,JLAT)=FIELD1D(indx)
        IBITMAP(JLON,JLAT)=1
        ENDDO
      ENDDO

      RETURN
      END

C############################################################
      subroutine es_ini
c
      common /estab/esat(15000:45000),es(15000:45000)
c
c *** Create tables of the saturation vapour pressure with up to
c        two decimal figures of accuraccy:
c
      do it=15000,45000
         t=it*0.01
         p1 = 11.344-0.0303998*t
         p2 = 3.49149-1302.8844/t
         c1 = 23.832241-5.02808*alog10(t)
         esat(it) = 10.**(c1-1.3816E-7*10.**p1+
     .               8.1328E-3*10.**p2-2949.076/t)
         es(it) = 610.78*exp(17.269*(t-273.16)/(t-35.86))
      enddo
c
      return
      end
c
