   MODULE utilities_netCDF
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  Purpose: A general-purpose module to write output in NetCDF format.
!           Each call to write_netCDF will:
!               Create a file associated with the field namelist
!               if it has not yet been created.
!               Write the field to the file associated.
!           A link list is used to keep track of everything associated
!               with the field.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   USE netcdf

   IMPLICIT NONE
   PRIVATE
   SAVE

   INTEGER, PARAMETER :: & ! the maximum rank of an array
      rank_max = 7         ! in fortran

   TYPE netCDF_node
      CHARACTER (LEN=128) :: &
         field_strng,path_strng,file_strng
      LOGICAL :: &
         l_file_open
      INTEGER :: &
         rank,index_max(rank_max+1),start_unlimited
      CHARACTER (LEN=40) :: &
         index_strng(rank_max+1),unit_strng
      INTEGER :: &
         file_ncid,index_dim(rank_max+1),field_varid
      TYPE (netCDF_node),POINTER :: &
         next
   END TYPE netCDF_node

   TYPE (netCDF_node),POINTER :: &
      netCDF_head

   LOGICAL :: &
      l_allocate_netCDF_head = .TRUE.

   PUBLIC :: write_netCDF
   
   CONTAINS
!======================================================================
!  BEGIN write_netCDF
!======================================================================
   SUBROUTINE write_netCDF (field,path,      &
                            index_max,index_strng,unit_strng, &
                            rk0,rk1,rk2,rk3,rk4,rk5,rk6,rk7)
!.......................................................................
!  INTENT IN
!.......................................................................
   CHARACTER (LEN=*),INTENT(IN) :: &
      field,path
   INTEGER, OPTIONAL :: &
      index_max(:)
   CHARACTER (LEN=40),INTENT(IN),OPTIONAL :: &
      index_strng(:),unit_strng
   REAL, INTENT(IN),OPTIONAL :: &
      rk0,rk1(:),rk2(:,:),rk3(:,:,:),rk4(:,:,:,:),rk5(:,:,:,:,:), &
      rk6(:,:,:,:,:,:),rk7(:,:,:,:,:,:,:)
!.......................................................................
!  LOCAL
!.......................................................................
   LOGICAL :: &
      l_add_new_node,l_array_passed
   INTEGER :: &
      n,rank_tmpry,index_tmpry(rank_max+1), &
      start(rank_max+1),count(rank_max+1),status
   CHARACTER (LEN=  1) :: &
      none_strng(8)
   CHARACTER (LEN=23) :: &
      time_strng
   TYPE (netCDF_node),POINTER :: &
      ptr
!.:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:.
!-----------------------------------------------------------------------
!  determine if an array is passed to subroutine.  if so, determine its
!  rank and shape.
!-----------------------------------------------------------------------
      l_array_passed=.FALSE.; rank_tmpry=-1; index_tmpry(:)=-1;
      IF (PRESENT (rk0)) THEN
         l_array_passed=.TRUE.; rank_tmpry=0;
      ENDIF
      IF (PRESENT (rk1)) THEN
         l_array_passed=.TRUE.; rank_tmpry=1; index_tmpry(1:1)=SHAPE (rk1);
      ENDIF
      IF (PRESENT (rk2)) THEN
         l_array_passed=.TRUE.; rank_tmpry=2; index_tmpry(1:2)=SHAPE (rk2);
      ENDIF
      IF (PRESENT (rk3)) THEN
         l_array_passed=.TRUE.; rank_tmpry=3; index_tmpry(1:3)=SHAPE (rk3);
      ENDIF
      IF (PRESENT (rk4)) THEN
         l_array_passed=.TRUE.; rank_tmpry=4; index_tmpry(1:4)=SHAPE (rk4);
      ENDIF
      IF (PRESENT (rk5)) THEN
         l_array_passed=.TRUE.; rank_tmpry=5; index_tmpry(1:5)=SHAPE (rk5);
      ENDIF
      IF (PRESENT (rk6)) THEN
         l_array_passed=.TRUE.; rank_tmpry=6; index_tmpry(1:6)=SHAPE (rk6);
      ENDIF
      IF (PRESENT (rk7)) THEN
         l_array_passed=.TRUE.; rank_tmpry=7; index_tmpry(1:7)=SHAPE (rk7);
      ENDIF
!-----------------------------------------------------------------------
!  the first time write_netCDF is called allocate the head node
!-----------------------------------------------------------------------
   l_add_new_node = .FALSE.
   IF (l_allocate_netCDF_head) THEN
      l_allocate_netCDF_head = .FALSE.
      l_add_new_node = .TRUE.
      ALLOCATE (netCDF_head); ptr => netCDF_head
      NULLIFY (ptr%next)
   ELSE
!-----------------------------------------------------------------------
!  look through the list to find a pre-existing field
!-----------------------------------------------------------------------
      ptr => netCDF_head
      DO 
         IF (TRIM (ptr%field_strng)==TRIM (field)) EXIT
         IF (.NOT.ASSOCIATED (ptr%next)) THEN
            l_add_new_node = .TRUE.
            EXIT
         ENDIF
         ptr => ptr%next
      ENDDO
   ENDIF
!-----------------------------------------------------------------------
!  if pre-existing field not found add a new node to the end of list
!-----------------------------------------------------------------------
   IF (l_add_new_node) THEN
      PRINT *," write_netCDF :: add new node :: field = ",TRIM (field)
      ALLOCATE (ptr%next)
      ptr => ptr%next
      NULLIFY (ptr%next)
   ENDIF
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  create file and set attributes of a new node
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   IF (l_add_new_node) THEN
!-----------------------------------------------------------------------
!  create a netCDF file and add timestamp attribute
!-----------------------------------------------------------------------
      ptr%field_strng = TRIM (field)
      ptr%path_strng  = TRIM (path)
      ptr%file_strng  = TRIM (path)//"/"//TRIM (field)//".nc"
      status = NF90_CREATE (TRIM (ptr%file_strng),NF90_CLOBBER,ptr%file_ncid)
      CALL NF_handle_error (status," write_netCDF :: NF90_CREATE")
      time_strng = time_stamp ()
      status = NF90_PUT_ATT (ptr%file_ncid,NF90_GLOBAL, &
                          "time_strng_::_time_of_file_creation",time_strng)
      CALL NF_handle_error (status," write_netCDF :: NF90_PUT_ATT time_strng")

!-----------------------------------------------------------------------
!  if index_max is passed to subroutine it determines the rank and shape
!  of the field.  verify that it is consistent with the passed array (if any).  
!  if index_max is not passed use the passed array to determine the rank and 
!  shape of the field.
!-----------------------------------------------------------------------
      IF (PRESENT (index_max)) THEN
         ptr%rank = SIZE (index_max(:),DIM=1)
         ptr%index_max(1:ptr%rank) = index_max(1:ptr%rank)

         IF (l_array_passed) THEN
            IF ((ptr%rank/=rank_tmpry).OR. &
                (ANY (ptr%index_max(1:ptr%rank)/=index_tmpry(1:ptr%rank)))) THEN
               PRINT "(A34,A128)"," netCDF ERROR :: ptr%file_strng = ", &
                                                         TRIM (ptr%file_strng)
               PRINT "(A45)"," index_max and passed array and inconsistent "
               STOP
            ENDIF
         ENDIF
      ELSE 
         IF (l_array_passed) THEN
            ptr%rank         = rank_tmpry
            ptr%index_max(:) = index_tmpry(:)
         ELSE
            PRINT "(A34,A128)"," netCDF ERROR :: ptr%file_strng = ", &
                                                         TRIM (ptr%file_strng)
            PRINT "(A37)"," cannot determine rank and index_max "
         ENDIF
      ENDIF
!-----------------------------------------------------------------------
! determine index_strng
!-----------------------------------------------------------------------
      none_strng = (/ "1","2","3","4","5","6","7","8" /)
      ptr%index_strng(:) = "undefined"
      IF (PRESENT (index_strng)) THEN
         DO n = 1,ptr%rank
            ptr%index_strng(n) = index_strng(n)
         ENDDO
      ELSE
         DO n = 1,ptr%rank
            ptr%index_strng(n) = "none_"//none_strng(n)
         ENDDO
      ENDIF
      ptr%index_strng(ptr%rank+1) = "unlimited"
!-----------------------------------------------------------------------
! define dimensions. the last dimension is defined UNLIMITED
!-----------------------------------------------------------------------
      DO n = 1,ptr%rank
         status = NF90_DEF_DIM (ptr%file_ncid,TRIM(ptr%index_strng(n)), &
                                    ptr%index_max(n),ptr%index_dim(n))
         CALL NF_handle_error (status," write_netCDF :: NF90_DEF_DIM loop ")
      ENDDO
      status = NF90_DEF_DIM (ptr%file_ncid,TRIM(ptr%index_strng(ptr%rank+1)), &
                                 NF90_UNLIMITED,ptr%index_dim(ptr%rank+1))
      CALL NF_handle_error (status," write_netCDF :: NF90_DEF_DIM UNLIMITED ")
!-----------------------------------------------------------------------
! define field
!-----------------------------------------------------------------------
      status = NF90_DEF_VAR (ptr%file_ncid,TRIM(ptr%field_strng),NF90_FLOAT, &
                       ptr%index_dim(1:ptr%rank+1),ptr%field_varid)
      CALL NF_handle_error (status," write_netCDF :: NF90_DEF_VAR ")
!-----------------------------------------------------------------------
! end definition
!-----------------------------------------------------------------------
      status = NF90_ENDDEF (ptr%file_ncid)
      CALL NF_handle_error (status," write_netCDF :: NF90_ENDDEF ")

      ptr%start_unlimited = 0

   ENDIF ! (l_add_new_node)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  write field
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   IF (l_array_passed) THEN
      status = NF90_OPEN (TRIM (ptr%file_strng),NF90_WRITE,ptr%file_ncid)
      CALL NF_handle_error (status," write_netCDF :: NF90_OPEN ")

      ptr%start_unlimited = ptr%start_unlimited + 1
      start(1:ptr%rank  ) = 1
      start(  ptr%rank+1) = ptr%start_unlimited

      count(1:ptr%rank  ) = ptr%index_max(1:ptr%rank)
      count(  ptr%rank+1) = 1

      SELECT CASE (ptr%rank)
      CASE (01)
         status = NF90_PUT_VAR (ptr%file_ncid,ptr%field_varid, &
                                rk1,start(1:2),count(1:2))
      CASE (02)
         status = NF90_PUT_VAR (ptr%file_ncid,ptr%field_varid, &
                                rk2,start(1:3),count(1:3))
      CASE (03)
         status = NF90_PUT_VAR (ptr%file_ncid,ptr%field_varid, &
                                rk3,start(1:4),count(1:4))
      CASE (04)
         status = NF90_PUT_VAR (ptr%file_ncid,ptr%field_varid, &
                                rk4,start(1:5),count(1:5))
      CASE (05)
         status = NF90_PUT_VAR (ptr%file_ncid,ptr%field_varid, &
                                rk5,start(1:6),count(1:6))
      CASE (06)
         status = NF90_PUT_VAR (ptr%file_ncid,ptr%field_varid, &
                                rk6,start(1:7),count(1:7))
      CASE (07)
         status = NF90_PUT_VAR (ptr%file_ncid,ptr%field_varid, &
                                rk7,start(1:8),count(1:8))
      END SELECT
      CALL NF_handle_error (status," write_netCDF :: NF90_PUT_VAR ")

      status = NF90_CLOSE (ptr%file_ncid)
      CALL NF_handle_error (status," write_netCDF :: NF90_CLOSE ")
   ENDIF

   END SUBROUTINE write_netCDF
!=======================================================================
!  END write_netCDF
!=======================================================================

!=======================================================================
!  BEGIN time_stamp
!=======================================================================
   FUNCTION time_stamp () RESULT (time_strng)
!.......................................................................
!  INTENT OUT
!.......................................................................
   CHARACTER (LEN=23) :: &
      time_strng
!.......................................................................
!  LOCAL
!.......................................................................
    INTEGER :: &
      date_time(8)
    CHARACTER (LEN=4) :: &
      s1
    CHARACTER (LEN=2) :: &
      s2,s3,s4,s5,s6
!.:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:.

   CALL DATE_AND_TIME (VALUES=date_time)

   WRITE (UNIT=s1,FMT="(I04)") date_time(1)
   WRITE (UNIT=s2,FMT="(I02)") date_time(2)
   WRITE (UNIT=s3,FMT="(I02)") date_time(3)
   WRITE (UNIT=s4,FMT="(I02)") date_time(5)
   WRITE (UNIT=s5,FMT="(I02)") date_time(6)
   WRITE (UNIT=s6,FMT="(I02)") date_time(7)

   time_strng = "("//s1//" "//s2//" "//s3//") ("//s4//":"//s5//":"//s6//")"

   END FUNCTION time_stamp
!=======================================================================
!  END time_stamp
!=======================================================================

!======================================================================
!  BEGIN NF_handle_error
!======================================================================
   SUBROUTINE NF_handle_error (status,message)
!.......................................................................
!  INTENT IN
!.......................................................................
   INTEGER :: status
   CHARACTER (LEN=*),INTENT(IN),OPTIONAL :: &
      message
!.:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:.

   IF (status==NF90_NOERR) RETURN

   IF (PRESENT (message)) THEN
      PRINT "(A27,A128)"," netCDF ERROR :: message = ",TRIM (message)
   ELSE
      PRINT "(A32,A128)"," netCDF ERROR :: NF_STRERROR = ",NF90_STRERROR (STATUS)
   ENDIF

   STOP

   END SUBROUTINE NF_handle_error
!=======================================================================
!  END NF_handle_error
!=======================================================================

   END MODULE utilities_netCDF
!cccccccc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5cccgtgccc6ccccccccc7cc

   PROGRAM main
   USE utilities_netCDF
   IMPLICIT NONE

   REAL, DIMENSION(10,10) :: temp, temp2
   INTEGER :: i, j, n
   REAL, PARAMETER :: pi = 4.*ATAN(1.)
   CHARACTER(LEN=40) :: xstring, ystring, tstring

   xstring = 'x'
   ystring = 'y'
   tstring = 'time'

   DO n = 1,3
      DO j = 1,10
         DO i= 1,10
            temp(i,j) = 273. + 5.*COS( 2.*pi*(i+j+n)/10. )
            temp2(i,j) = 273. + 5.*COS( 2.*pi*(i-j+n)/10. )
         ENDDO
      ENDDO

      CALL write_netCDF ('TEMP','./',      &
                         index_strng=(/xstring,ystring,tstring/) ,unit_strng='K', &
                         rk2=temp)
!      CALL write_netCDF ('TEMP','./',      &
!                         unit_strng='K', &
!                         rk2=temp)
      CALL write_netCDF ('TEMP2','./',      &
                         unit_strng='K', &
                         rk2=temp2)

   ENDDO

   END PROGRAM main
