!
! $Id: read_mask.f,v 1.1 2003/09/23 04:34:24 dimitri Exp $
!
!_ ---------------------------------------------------------------------
!_ $Log: read_mask.f,v $
!_ Revision 1.1  2003/09/23 04:34:24  dimitri
!_ o Mods and bug fixes to pkg/cal and pkg/exf needed for computation
!_   of tracer Green's fucntions for ocean inversion project.
!_
!_ Revision 1. 24/07/2003 Sara Fletcher
!_ Started Coding
!_
!***********************************************************************
!
! NAME: READ_MASK.F
!
! DESCRIPTION: This program demonstrates how to read a 3-D netCDF region mask 
!              file.
!
! COMPILING: 
!
!          f77 -o read_mask.x read_mask.f libnetcdf.a
!
! VARIABLES:
!       imt, jmt: indices of the region mask in the input file 
!       filename = Name of Net CDF file
!       ndyetrac  = Number of dye regions (30 or 16)
!       region_index = 2-D array containing a region index which assigns each 
!                      ocean region to a different index value
!       region_mask = 3-D (i,j,ndyetrac) array. Gridboxes included in each
!                     model region are set to 1 and those not included are set
!                     to 0.
!
!***********************************************************************

	PROGRAM read_mask
     
! No implicit declarations
        IMPLICIT NONE
        INCLUDE 'netcdf.inc'

! define variables
	CHARACTER*256 string256,filename
	INTEGER ndyetrac,imt,jmt, i,j,n 
        PARAMETER (ndyetrac=30)  !Equals 30 or 16
        PARAMETER (imt= 96, jmt=40)
        REAL*4 lon(imt)
        REAL*4 lat(jmt)
        CHARACTER*30, reg_name(ndyetrac)
        CHARACTER*1, string1
        CHARACTER*2, string2
        INTEGER region_index(imt,jmt)
        INTEGER mask_temp(imt,jmt), region_mask(imt,jmt,ndyetrac)
        

        INTEGER*4  NC_ID, STATUS
        INTEGER*4 idlon, idlat, idtrac, idindex, idmask(ndyetrac)



!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  Generate region names
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        
        do i=1,ndyetrac
         if (i.lt.10) then
          write(string1,'(I1.1)') i
          reg_name(i)='region'//string1
         else
          write(string2,'(I2.2)') i
          reg_name(i)='region'//string2
         endif 
       enddo


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NetCDF filename
! Either 30reg_regionmask.cdf or 16reg_regionmask.cdf 
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        
        if (ndyetrac.eq.30) then
          filename='30reg_regionmask.cdf'
        else if (ndyetrac.eq.16) then
          filename='16reg_regionmask.cdf'
        endif 

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Open the NetCDF file
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	STATUS=NF_OPEN(filename,NF_NOWRITE,NC_ID)
	IF (STATUS.NE.NF_NOERR) CALL HANDLE_ERRORS(STATUS)

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Get the netCDF variable dimensions 
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	
        STATUS=NF_INQ_DIMID(NC_ID,'tlon',idlon)
	IF (STATUS.NE.NF_NOERR) CALL HANDLE_ERRORS(STATUS)

	STATUS=NF_INQ_DIMID(NC_ID,'tlat',idlat)
         IF (STATUS.NE.NF_NOERR) CALL HANDLE_ERRORS(STATUS)
        
	STATUS=NF_INQ_VARID(NC_ID,'reg_indx',idindex)
	IF (STATUS.NE.NF_NOERR) CALL HANDLE_ERRORS(STATUS)

        do i=1,ndyetrac
	 STATUS=NF_INQ_VARID(NC_ID,reg_name(i),idmask(i))
	 IF (STATUS.NE.NF_NOERR) CALL HANDLE_ERRORS(STATUS)
        enddo

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! READ the netCDF DATA
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	STATUS=NF_GET_VAR_REAL(NC_ID,idlon, lon)
	IF (STATUS.NE.NF_NOERR) CALL HANDLE_ERRORS(STATUS)

	STATUS=NF_GET_VAR_REAL(NC_ID,idlat,lat)
	IF (STATUS.NE.NF_NOERR) CALL HANDLE_ERRORS(STATUS)

	STATUS=NF_GET_VAR_INT(NC_ID,idindex,region_index)
	IF (STATUS.NE.NF_NOERR) CALL HANDLE_ERRORS(STATUS)

        do n=1,ndyetrac
	 STATUS=NF_GET_VAR_INT(NC_ID,idmask(n),mask_temp)
	 IF (STATUS.NE.NF_NOERR) CALL HANDLE_ERRORS(STATUS)
         do i=1,imt
          do j=1,jmt
           region_mask(i,j,n)=0
           if(mask_temp(i,j).eq.1)region_mask(i,j,n)=1
           write(10,*),i,j,lon(i),lat(j),n,region_mask(i,j,n)
          enddo
         enddo
       enddo


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Close the NetCDF file
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	STATUS=NF_CLOSE(NC_ID)
	IF (STATUS.NE.NF_NOERR) CALL HANDLE_ERRORS(STATUS)

	END


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Error handling subroutine
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        SUBROUTINE HANDLE_ERRORS(STATUS)

        INCLUDE 'netcdf.inc'
        INTEGER*4 STATUS

        IF (STATUS .NE. nf_noerr) THEN
                PRINT *, nf_strerror(STATUS)
                STOP 'stopped'
        ENDIF

        END

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~