SUBROUTINE getcoamps(nx_ext, ny_ext, nz_ext,                            & 1,8
           dir_extd,extdinit,extdfcst,                                  &
           iproj_ext,scale_ext,                                         &
           trlon_ext,latnot_ext,x0_ext,y0_ext,                          &
           lat_ext,lon_ext,                                             &
           p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,                      &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                          &
           tsfc_ext,tsoil_ext,                                          &
           wetsfc_ext,wetdp_ext,wetcanp_ext,istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
! This programs reads coamps2.0 sigma level output
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    dir_extd      directory where the coamps files are located
!    extdinit      Initialized time in mm-dd-yyyy:hh:mm:ss format
!    extdfcst      Forecast hour in HHH:MM:SS format
!
!  OUTPUT:
!
!    iproj_ext     Map projection number of external data
!    scale_ext     Scale factor of external data
!    trlon_ext     True longitude of external data (degrees E)
!    latnot_ext(2) True latitude(s) of external data (degrees N)
!    x0_ext        x coordinate of origin of external data
!    y0_ext        y coordinate of origin of external data
!    lat_ext       latitude of external data points (degrees N)
!    lon_ext       longitude of external data points (degrees E)
!
!    p_ext         pressure (Pascal)
!    hgt_ext       height (m)
!    t_ext         temperature (K)
!    tsfc_ext      ground/sea-surface temperature (K)
!    tsoil_ext     Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc_ext    Surface soil moisture in the top 1 cm layer
!                                              (fraction,0--1.0))
!    wetdp_ext     Deep soil moisture in the deep 1 m layer
!    wetcanp_ext   Canopy water amount
!    qv_ext        specific humidity(mixing ratio used) (kg/kg)
!    qc_ext        Cloud water mixing ratio (kg/kg)
!    qr_ext        Rain  water mixing ratio (kg/kg)
!    qi_ext        Ice         mixing ratio (kg/kg)
!    qs_ext        Snow        mixing ratio (kg/kg)
!    qh_ext        Hail        mixing ratio (kg/kg)
!    u_ext         u wind component (m/s)
!    v_ext         v wind component (m/s)
!    istatus       status indicator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!

  INTEGER :: nx_ext, ny_ext, nz_ext
!
  CHARACTER (LEN=60) :: dir_extd
  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=200) :: FILE,cfile
!
!  COAMPS input fields
!
  REAL :: DATA(500)
  REAL :: fp(nx_ext,ny_ext,nz_ext),ft(nx_ext,ny_ext,nz_ext),            &
       fq(nx_ext,ny_ext,nz_ext),fu(nx_ext,ny_ext,nz_ext),               &
       fv(nx_ext,ny_ext,nz_ext),                                        &
       fs(nx_ext,ny_ext),fh(nx_ext,ny_ext),                             &
       fsoil(nx_ext,ny_ext),fgwet(nx_ext,ny_ext)

!
!  Original grid variables
!
  INTEGER :: iproj
  REAL :: scale,trlon,x0,y0
  REAL :: latnot(2)
!
!  Output external grid variables
!
  INTEGER :: iproj_ext
  REAL :: scale_ext,trlon_ext
  REAL :: latnot_ext(2)
  REAL :: x0_ext,y0_ext
  REAL :: x_ext(nx_ext),y_ext(ny_ext)
!
!  Output external variable arrays
!
  REAL :: lat_ext(nx_ext,ny_ext)
  REAL :: lon_ext(nx_ext,ny_ext)
  REAL :: p_ext(nx_ext,ny_ext,nz_ext)     ! (Pa)
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)   ! (m)
  REAL :: t_ext(nx_ext,ny_ext,nz_ext)     ! (K)
  REAL :: qv_ext(nx_ext,ny_ext,nz_ext)    ! (kg/kg)
  REAL :: u_ext(nx_ext,ny_ext,nz_ext)     ! (m/s)
  REAL :: v_ext(nx_ext,ny_ext,nz_ext)     ! (m/s)
  REAL :: qc_ext(nx_ext,ny_ext,nz_ext)    ! Cloud H2O mixing ratio (kg/kg)
  REAL :: qr_ext(nx_ext,ny_ext,nz_ext)    ! Rain  H2O mixing ratio (kg/kg)
  REAL :: qi_ext(nx_ext,ny_ext,nz_ext)    ! Ice   H2O mixing ratio (kg/kg)
  REAL :: qs_ext(nx_ext,ny_ext,nz_ext)    ! Snow  H2O mixing ratio (kg/kg)
  REAL :: qh_ext(nx_ext,ny_ext,nz_ext)    ! Hail  H2O mixing ratio (kg/kg)

  REAL :: tsfc_ext   (nx_ext,ny_ext)      ! Ground/sea-surface temp. (K)
  REAL :: tsoil_ext  (nx_ext,ny_ext)      ! Deep soil temperature (K)
  REAL :: wetsfc_ext (nx_ext,ny_ext)      ! Surface soil moisture (0-1)
  REAL :: wetdp_ext  (nx_ext,ny_ext)      ! Deep soil moisture (0-1)
  REAL :: wetcanp_ext(nx_ext,ny_ext)      ! Canopy water amount (0-1)


  INTEGER :: istatus
!
  INTEGER :: i,j,k,length_dir
  PARAMETER(len_dir=60)

  REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  ALLOCATE(utmp(nx_ext,ny_ext))
  ALLOCATE(vtmp(nx_ext,ny_ext))

!----------------------------------------------------------------------
!  Open and read the COAMPS file:
!  data--header file, fp--pressure (mb), ft--potential temp (K)
!  fq--water vapor mixing ratio (kg/kg), fu and fv -- wind (m/s)
!  fs--surface temp (K), fh--terrain height (m)
!----------------------------------------------------------------------
  length_dir = len_dir
  CALL strlnth(dir_extd, length_dir)
  cfile = dir_extd(1:length_dir)//extdinit//'+'//extdfcst
  PRINT *, 'Open coamps file= ', cfile
  OPEN(4,FILE=cfile, STATUS='old',FORM='unformatted')
  READ(4) DATA,fp,ft,fq,fu,fv,fs,fsoil,fgwet,fh
!
!-----------------------------------------------------------------------
!
!  Get the COAMPS grid configurations
!  In COAMPS 1=mercator, 2=lambert, 3=polar projection.
!
!-----------------------------------------------------------------------
!
  iproj_ext = INT(DATA(3))

  IF(iproj_ext == 1) THEN
    iproj_ext=3
  ELSE IF(iproj_ext == 3) THEN
    iproj_ext=1
  ELSE IF(iproj_ext > 3) THEN
    PRINT *, 'can not handle this projection', iproj_ext
    STOP
  END IF

  latnot_ext(1) = DATA(4)
  latnot_ext(2) = DATA(5)
  trlon_ext = DATA(6)-360.0
  scale_ext = 1.0
  swlat_ext = DATA(36)
  swlon_ext = DATA(37)-360.0
  dx_ext = DATA(9)
  dy_ext = DATA(10)
!
!-----------------------------------------------------------------------
!
!  Get the lat,lon of the COAMPS grid points
!
!-----------------------------------------------------------------------
!
  CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
  CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
  CALL lltoxy(1,1,swlat_ext,swlon_ext,x0_ext,y0_ext)
!
  DO i=1,nx_ext
    x_ext(i)=x0_ext+(i-1)*dx_ext
  END DO
  DO j=1,ny_ext
    y_ext(j)=y0_ext+(j-1)*dy_ext
  END DO

  CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
  PRINT *, ' maps point 34,17: ',lat_ext(34,17),lon_ext(34,17)
  PRINT *, ' maps point nx,ny: ',lat_ext(nx_ext,ny_ext),                &
                                 lon_ext(nx_ext,ny_ext)
!
!-----------------------------------------------------------------------
!
!  Get the COAMPS data
!
!-----------------------------------------------------------------------
!
  DO k=1,nz_ext
    nk=nz_ext-k+1
    DO j=1,ny_ext
      DO i=1,nx_ext
        hgt_ext(i,j,nk)= DATA(200+k)*(DATA(201)-fh(i,j))/DATA(201)      &
                        +fh(i,j)
      END DO
    END DO
  END DO

  DO k=1,nz_ext
    nk=nz_ext-k+1
    DO j=1,ny_ext
      DO i=1,nx_ext
        p_ext(i,j,nk)= fp(i,j,k)*100.0
        t_ext(i,j,nk)= ft(i,j,k)*(fp(i,j,k)/1000.0)**0.286
        qv_ext(i,j,nk) = fq(i,j,k)
        u_ext(i,j,nk)  = fu(i,j,k)
        v_ext(i,j,nk)  = fv(i,j,k)
      END DO
    END DO
  END DO

  DO j=1,ny_ext
    DO i=1,nx_ext
      tsfc_ext  (i,j) = fs(i,j)
      tsoil_ext (i,j) = fsoil(i,j)
      wetsfc_ext(i,j) = fgwet(i,j)
      wetdp_ext (i,j) = fgwet(i,j)  !assumed value!!!!!!!!!!!!!!!!
      wetcanp_ext(i,j)= 0.0         !assumed value!!!!!!!!!!!!!!!!
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Fill qc,qr,qi,qs,qh arrays with missing value.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz_ext
    DO j=1,ny_ext
      DO i=1,nx_ext
        qc_ext(i,j,k)=-999.
        qr_ext(i,j,k)=-999.
        qi_ext(i,j,k)=-999.
        qs_ext(i,j,k)=-999.
        qh_ext(i,j,k)=-999.
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Rotate winds to be relative to true north.
!  The COAMPS data are saved as grid-relative.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz_ext
!2001-05-16 GMB: Having umap & uear (or vmap & vear) point to
!the same array causes numerical errors when optimizing.
    CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k),               &
                 lon_ext,utmp,vtmp)
    u_ext(:,:,k) = utmp(:,:)
    v_ext(:,:,k) = vtmp(:,:)
  END DO
!
! test: print out
!
!  print *, 'iproj_ext,swlat_ext,swlon_ext,scale_ext,trlon_ext,
!    : latnot_ext(1),latnot_ext(2),lat_ext(1),lon_ext(1),dx_ext.dy_ext'
!  write(6,'(i3,10f12.4)') iproj_ext,swlat_ext,swlon_ext,
!    : scale_ext,trlon_ext,latnot_ext(1),latnot_ext(2),
!    : lat_ext(1,1),lon_ext(1,1),dx_ext,dy_ext
!  do 666 k=1,nz_ext
!  print *, 'k,p,t,q,u,v,s,h'
!  write(6,'(i3,6f12.4)')
!    : k,p_ext(20,20,k),t_ext(20,20,k),qv_ext(20,20,k),
!    : u_ext(20,20,k),v_ext(20,20,k),
!    : tsfc_ext(20,20)
!666   continue
!
!-----------------------------------------------------------------------
!
!  Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
  CALL setmapr(iproj,scale,latnot,trlon)
  CALL setorig(1,x0,y0)
!
!-----------------------------------------------------------------------
!
!  Set good status
!
!-----------------------------------------------------------------------
!
  istatus=1

  CLOSE(4)

  DEALLOCATE(utmp)
  DEALLOCATE(vtmp)

  RETURN
END SUBROUTINE getcoamps