psmile_write_meta_byid.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2006-2010, SGI Germany, Munich, Germany.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: PSMILe_Write_meta_byid
!
! !INTERFACE:

      subroutine psmile_write_meta_byid (id_varid,id_taskid,ierror)
! !USES:
      use prism_constants
      use psmile, dummy_interface => psmile_write_meta_byid
#ifdef __PSMILE_WITH_IO
      use PRISM_calendar
      use psmile_io_utils
      use mpp_io_mod_oa, only: mpp_nullify_axistype_array
#endif
      implicit none
!
! !INPUT PARAMETERS:
!
      integer,intent(in)  :: id_varid
      integer,intent(in)  :: id_taskid
!
! !OUTPUT PARAMETERS:
!
      integer,intent(out) :: ierror
!
! !LOCAL VARIABLES
!
#ifdef __PSMILE_WITH_IO
      integer :: ierrp(2)
      integer :: il,ih,jl,jh,kl,kh,kk,ll,lh,ii,il_i,jj,ilbl,ic,jc,kc,iloc,jloc
      integer :: il_grid_id,il_method_id,il_grid_type,iicomp_id,blkid,kkj
      integer :: len,lenaz,nvcomp,vcomp_id,vmsk_id(3),il_gridparam(3)
      integer :: il_date(6)
      integer :: lenb,no_of_blocks,il_id
      integer :: myashape(2,5)
      integer :: myvshape(2,5)
      integer :: il_pack
      integer :: fullsize
      integer :: len_alloc(3),nof(ndim_3d),ncells(ndim_3d),nvertices(ndim_3d) 
                 ,noc,stride
      logical :: is_bundle,is_regular, is_regvrt, lbackward
      real ,pointer :: lonr(:),latr(:),vertr(:)
      real ,pointer :: lonrb(:),latrb(:)
      real          ,pointer :: vertrb(:,:)
      real ,pointer :: bound_lonr(:,:,:)
      real ,pointer :: bound_latr(:,:,:)
      real ,pointer :: bound_vertr(:,:,:,:)
      double precision :: dl_time
      double precision,pointer :: lond(:),latd(:)
      double precision,pointer :: londb(:),latdb(:)
      double precision,pointer :: vertd(:)
      double precision,pointer :: vertdb(:,:)
      double precision,pointer :: bound_lond(:,:,:)
      double precision,pointer :: bound_latd(:,:,:)
      double precision,pointer :: bound_vertd(:,:,:,:)
      character(len=max_name) :: cl_field(3),cl_latij(3) 
                                ,cl_lonij(3),cl_vertij(3) 
                                ,cl_mskij(3),cl_crnij(3) 
                                ,cl_subaij(3),cl_vert_longname(3) 
                                ,cl_grid_longname(3)
      character(len=max_name) :: cl_timeunit
      character(len=4)        :: label

      Type fds
         double precision,pointer::fdble(:)
      End Type fds
      Type(fds),pointer::fdbles(:,:,:)

      Type frs
         real,pointer::freal(:)
      End Type frs
      Type(frs),pointer::freals(:,:,:)
      Type flgs
          logical,pointer::flog(:)
      End Type flgs
      Type(flgs),pointer::flogs(:,:)

      Type(PRISM_Time_Struct) :: ljobstart
      Type (Corner_Block), pointer :: corner_pointer
#endif
!
! !DESCRIPTION:
! For a prognostic field identified by id_varid this routine accomplishes 
! the following:
!
! 1st) The routine declares the so-called axis of a NetCDF file.
! 2nd) Axis are assigned to field variables. Here, a field variable can
!      be a set of longitudes,latitudes, vertical axis information
!      masks and a prognostic  variable . The result is the 
!      NetCDF CF header information or the so-called metadata information.
! 3rd) All data associated with the metadata of latitudes, longitudes,
!      verticals and masks are written to the NetCDF CF file.
!
! The routine supports scalar and vector fields and bundles of fields.
! Component of vector fields are assigned to different NetCDF variables.
! Thus it is possible to plot vector fields - which are 4d in general -
! with 3d plotting facilities.
!
! For all of the types of prognostic variables mentioned above  GME-like
! blocked grids  are supported  where variables are assigned in a 
! round-robin fashion to blocks.
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 09.12.03    R. Vogelsang created
! 11.02.04    R. Vogelsang Modifications to use SMIOC related infos
! 04.03.01    R.Vogelsang  Added code to extract labels of bundles
! 23.3.04     R. Vogelsang Added the taskid to the interface and related code.
! 25.2.05     R. Vogelsang Corrrection of writing and extracting bounds
!                  Bounds stored all vertices(1), all vertices(2)....
! 02.10.05    R. Vogelsang Added writing of vertical bounds and handling
!              of Gaussian reduced grids.
!
! 28.02.06    S. Valckes   prism.inc removed, replaced by prism_constants
! 04.03.06    R.Vogelsang  Sequence of writing corner points changed
!                          to allow  reasonable plots. Nullify of
!                          pointers and pointer substructures.
!EOP
!----------------------------------------------------------------------
!
      Character(len=len_cvs_string), save :: mycvs = 
      '$Id: psmile_write_meta_byid.F90 2755 2010-11-19 15:19:52Z hanke $'

      ierror=0

#ifdef __PSMILE_WITH_IO
#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid: start'
      call psmile_flushstd

#endif
      
!
!     Not subject to I/O. Return!
!
      if (.not.associated(Fields(id_varid)%io_chan_infos)) return
      
      Nullify (fdbles)
      Nullify (freals)
      Nullify (flogs)
      Nullify (bound_lonr)
      Nullify (bound_latr)
      Nullify (bound_vertr)
      Nullify (bound_lond)
      Nullify (bound_latd)
      Nullify (bound_vertd)
!
!     IO infos are pointing to io_chan_infos(taskid)
!
      if(id_taskid .le. size(Fields(id_varid)%io_task_lookup)) then

        il_i=Fields(id_varid)%io_task_lookup(id_taskid)

      else

        ierror = PRISM_Error_IO_Meta
        ierrp (1) = id_taskid
        call psmile_error ( ierror &
                          , 'Task id out of range! ', &
                            ierrp, 1, __FILE__, __LINE__ )
        return

      endif

      if (il_i.gt.0) then

         Fields(id_varid)%io_infos =>  Fields(id_varid)%io_chan_infos(il_i)

      else
        ierror = PRISM_Error_IO_Meta
        ierrp (1) = id_taskid
        call psmile_error ( ierror &
                          , 'Negative task id! ', &
                            ierrp, 1, __FILE__, __LINE__ )
        return
      endif
!
!     The file is not considered for writing. Return!
!
      if ( Fields(id_varid)%io_infos%action .ne.MPP_OVERWR) return
!
!     For GME grid only the last varid of the set of related varids
!     does the writing of the meta data.
!
      if(associated(Fields(id_varid)%io_infos%related_ids)) then
        fullsize=size(Fields(id_varid)%io_infos%related_ids)
        if(fullsize.gt.1) then
          do il_i=1,fullsize
           if(id_varid.eq. Fields(id_varid)%io_infos%related_ids(il_i)) exit 
          enddo
         if(il_i.lt.fullsize) return
        endif
      endif
!     
!     Last check: Is the file open?
!
      if (.not. Fields(id_varid)%io_infos%opened) then
         ierror = PRISM_Error_IO_Meta
         ierrp (1) = ierror
         call psmile_error ( ierror &
                           , 'Attempt to write meta data to a non-open file', &
                             ierrp, 1, __FILE__, __LINE__ )
         return
      endif

!
!     Set the pelist for the current active component
!
      iicomp_id=Fields(id_varid)%comp_id
      call mpp_set_current_pelist(IO_Comps_infos(iicomp_id)%pelist)

!
!     Allocate axis descriptors
!
      
      if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%x_axis)) then
          allocate(Fields(id_varid)%io_infos%p_mpp_io%x_axis(2),stat=ierror)
          call mpp_nullify_axistype_array(Fields(id_varid)%io_infos%p_mpp_io%x_axis)
          if (ierror > 0) then
             ierrp (1) = ierror
             ierrp (2) = 1
             ierror = PRISM_Error_Alloc
             call psmile_error ( ierror, 'x_axis', &
                             ierrp, 2, __FILE__, __LINE__ )
             return
          endif
      endif


      if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%y_axis)) then
          allocate(Fields(id_varid)%io_infos%p_mpp_io%y_axis(1),stat=ierror)
          call mpp_nullify_axistype_array(Fields(id_varid)%io_infos%p_mpp_io%y_axis)
          if (ierror > 0) then
             ierrp (1) = ierror
             ierrp (2) = 1
             ierror = PRISM_Error_Alloc
             call psmile_error ( ierror, 'y_axis', &
                             ierrp, 2, __FILE__, __LINE__ )
             return
          endif
      endif
      
      if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%z_axis)) then
          allocate(Fields(id_varid)%io_infos%p_mpp_io%z_axis(1),stat=ierror)
          call mpp_nullify_axistype_array(Fields(id_varid)%io_infos%p_mpp_io%z_axis)
          if (ierror > 0) then
             ierrp (1) = ierror
             ierrp (2) = 1
             ierror = PRISM_Error_Alloc
             call psmile_error ( ierror, 'z_axis', &
                             ierrp, 2, __FILE__, __LINE__ )
             return
          endif
      endif

      if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%b_axis)) then
          allocate(Fields(id_varid)%io_infos%p_mpp_io%b_axis(1),stat=ierror)
          call mpp_nullify_axistype_array(Fields(id_varid)%io_infos%p_mpp_io%b_axis)
          if (ierror > 0) then
             ierrp (1) = ierror
             ierrp (2) = 1
             ierror = PRISM_Error_Alloc
             call psmile_error ( ierror, 'b_axis', &
                             ierrp, 2, __FILE__, __LINE__ )
             return
          endif
      endif

      if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%blk_axis)) then
          allocate(Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1),stat=ierror)
          call mpp_nullify_axistype_array(Fields(id_varid)%io_infos%p_mpp_io%blk_axis)
          if (ierror > 0) then
             ierrp (1) = ierror
             ierrp (2) = 1
             ierror = PRISM_Error_Alloc
             call psmile_error ( ierror, 'blk_axis', &
                             ierrp, 2, __FILE__, __LINE__ )
             return
          endif
      endif

      if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%t_axis)) then
          allocate(Fields(id_varid)%io_infos%p_mpp_io%t_axis(1),stat=ierror)
          call mpp_nullify_axistype_array(Fields(id_varid)%io_infos%p_mpp_io%t_axis)
          if (ierror > 0) then
             ierrp (1) = ierror
             ierrp (2) = 1
             ierror = PRISM_Error_Alloc
             call psmile_error ( ierror, 't_axis', &
                             ierrp, 2, __FILE__, __LINE__ )
             return
          endif
      endif

      if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%crnaxis)) then
          allocate(Fields(id_varid)%io_infos%p_mpp_io%crnaxis(3),stat=ierror)
          call mpp_nullify_axistype_array(Fields(id_varid)%io_infos%p_mpp_io%crnaxis)
          if (ierror > 0) then
             ierrp (1) = ierror
             ierrp (2) = 1
             ierror = PRISM_Error_Alloc
             call psmile_error ( ierror, 'crnaxis(3)', &
                             ierrp, 2, __FILE__, __LINE__ )
             return
          endif
      endif

      if(.not.associated(Fields(id_varid)%io_infos%p_mpp_io%domain)) then
         ierrp (1) = id_varid
         ierrp (2) = 1
         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'domain not initialized!', &
                         ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
!     Get the global partition
!
      call mpp_get_global_domain(Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
                                ,xbegin=il,xend=ih         &
                                ,ybegin=jl,yend=jh         )
!
!     Determine how many shape components of the grid where declared
!
      
      il_grid_id =Methods(Fields(id_varid)%method_id)%grid_id
      il_grid_type=Grids(il_grid_id)%grid_type
      if ( il_grid_type == PRISM_unstructlonlatvrt) then
         len = 1
      else if ( il_grid_type == PRISM_unstructlonlat_sigmavrt .or. &
                il_grid_type == PRISM_unstructlonlat_regvrt   ) then
         len = 2
      else if ( il_grid_type == PRISM_Gaussreduced_regvrt .or. &  
                il_grid_type == PRISM_Gaussreduced_sigmavrt ) then
         len = 3
      else
         len = Grids(il_grid_id)%n_dim
      endif
!------------------------------------------------------------------------------
!
!     Is the field a vector field or a scalar field ?
!
      if(associated(Methods(Fields(id_varid)%method_id)%vector_pointer)) then
!
!     Allocate field descriptors. Vector fields have 3 components
!     Anyway, use the statement below for the generality.
!
        nvcomp=size(Methods(Fields(id_varid)%method_id) &
                   %vector_pointer%array_of_point_ids)
       is_bundle=.false.
#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid: Field(' &
            ,id_varid,') is vector field! nvcomp = ',nvcomp
      call psmile_flushstd
#endif
      else ! Scalar field

        nvcomp=1
!     
!     Is the scalar field a bundle ?
!       
        lenb=Fields(id_varid)%var_shape(2,len+1) &
            -Fields(id_varid)%var_shape(1,len+1)+1
        if(lenb.gt.1) then
          is_bundle=.true.
#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid: Field(' &
            ,id_varid,') is bundle! lenb = ',lenb
      call psmile_flushstd
#endif
        else ! No bundle
          is_bundle=.false.
        endif
      endif

!
!     Allocate field descriptors for Fields(id_varid),lons,lats, verts 
!     bounds and masks.
!
        if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%field)) &
            allocate(Fields(id_varid)%io_infos%p_mpp_io%field(nvcomp),stat=ierror)
        if (ierror > 0) then
           ierrp (1) = ierror
           ierrp (2) = 1
           ierror = PRISM_Error_Alloc
           call psmile_error ( ierror, 'field', &
                           ierrp, 2, __FILE__, __LINE__ )
           return
        endif
        if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%latij)) &
            allocate(Fields(id_varid)%io_infos%p_mpp_io%latij(nvcomp),stat=ierror)
        if (ierror > 0) then
           ierrp (1) = ierror
           ierrp (2) = 1
           ierror = PRISM_Error_Alloc
           call psmile_error ( ierror, 'latij', &
                           ierrp, 2, __FILE__, __LINE__ )
           return
        endif
        if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%lonij)) &
            allocate(Fields(id_varid)%io_infos%p_mpp_io%lonij(nvcomp),stat=ierror)
        if (ierror > 0) then
           ierrp (1) = ierror
           ierrp (2) = 1
           ierror = PRISM_Error_Alloc
           call psmile_error ( ierror, 'lonij', &
                           ierrp, 2, __FILE__, __LINE__ )
           return
        endif
        if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%vertij)) &
            allocate(Fields(id_varid)%io_infos%p_mpp_io%vertij(nvcomp),stat=ierror)
        if (ierror > 0) then
           ierrp (1) = ierror
           ierrp (2) = 1
           ierror = PRISM_Error_Alloc
           call psmile_error ( ierror, 'vertij', &
                           ierrp, 2, __FILE__, __LINE__ )
           return
        endif

        if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%maskij)) &
            allocate(Fields(id_varid)%io_infos%p_mpp_io%maskij(nvcomp),stat=ierror)
        if (ierror > 0) then
           ierrp (1) = ierror
           ierrp (2) = 1
           ierror = PRISM_Error_Alloc
           call psmile_error ( ierror, 'maskij', &
                           ierrp, 2, __FILE__, __LINE__ )
           return
        endif

        if (.not.associated(Fields(id_varid)%io_infos%p_mpp_io%crnij)) &
            allocate(Fields(id_varid)%io_infos%p_mpp_io%crnij(3),stat=ierror)
        if (ierror > 0) then
           ierrp (1) = ierror
           ierrp (2) = 1
           ierror = PRISM_Error_Alloc
           call psmile_error ( ierror, 'crnij(3)', &
                           ierrp, 2, __FILE__, __LINE__ )
           return
        endif
!
!      Define axis. Assumption is here that each of the vector
!      components lives on the SAME grid indentified by grid_id.
!
        call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                           ,Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                           ,'X','1','global index space for x' &
                           ,domain=Fields(id_varid)%io_infos%p_mpp_io%xdom(1) &
                           ,data=(/(il_i-0.d0,il_i=il,ih)/))
!rr                        ,data=(/(il_i-1.d0,il_i=il,ih)/))
        call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                           ,Fields(id_varid)%io_infos%p_mpp_io%y_axis(1) &
                           ,'Y','1','global index space for y' &
                           ,domain=Fields(id_varid)%io_infos%p_mpp_io%ydom(1) &
                           ,data=(/(il_i-0.d0,il_i=jl,jh)/))
!rr                        ,data=(/(il_i-1.d0,il_i=jl,jh)/))
        if(len.eq.1) then
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                           ,Fields(id_varid)%io_infos%p_mpp_io%z_axis(1) &
                           ,'Z','1','global index space for z' &
                           ,data=(/1.d0/))
        else
          lenaz=Fields(id_varid)%var_shape(2,len) &
               -Fields(id_varid)%var_shape(1,len)+1
          if(lenaz.gt.1) then
!
!         The user has declared the third dimension for a field!
!
            kl=Grids(il_grid_id)%grid_shape(1,len)
            kh=Grids(il_grid_id)%grid_shape(2,len)
#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid: Field(' &
            ,id_varid,') has declared z-axis!',kl, kh
      call psmile_flushstd
#endif
          else
            kl=1
            kh=1
          endif
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                           ,Fields(id_varid)%io_infos%p_mpp_io%z_axis(1) &
                           ,'Z','1','global index space for z' &
                           ,data=(/(dble(real(il_i)),il_i=kl,kh)/))
          
        endif

        if(.not.is_bundle) then
!
!       Declare a 4th pseudo  axis for vector components. 
!       Maybe, I need that for a future desgin where a vector field
!       is written as 4d field.
!
!          ll=1
!          lh=size(Methods(Fields(id_varid)%method_id)%vector_pointer &
!                                                   %array_of_point_ids)
!          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
!                             ,Fields(id_varid)%io_infos%p_mpp_io%b_axis(1) &
!                             ,'B','1','index space for vector components' &
!                             ,data=(/(dble(il_i),il_i=ll,lh)/))
         else
!
!         Declare bundle axis
!
          ll=1
          lh=lenb

          if(size(Fields(id_varid)%io_infos%labels).ne.lenb) then
            
            ierror=PRISM_Error_IO_Meta
            ierrp(1)=lenb
            ierrp(2)=size(Fields(id_varid)%io_infos%labels)
            call psmile_error(ierror,'Number oflabels given by SMIOC '// &
                              'does not match number of bundles!', &
                               ierrp,2, __FILE__, __LINE__)
          endif

          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                             ,Fields(id_varid)%io_infos%p_mpp_io%b_axis(1) &
                             ,'B','1','labels of the bundle components' &
                             ,cdata=(/(Fields(id_varid)%io_infos%labels(il_i) &
                                      ,il_i=ll,lh)/))
         endif
!
!        Declare a block index axis for GME style grids
!
         no_of_blocks=1
         if(associated(Fields(id_varid)%io_infos%related_ids)) then
            no_of_blocks=size(Fields(id_varid)%io_infos%related_ids) 
            if(no_of_blocks .gt. 1) then
              ll=1
              lh=no_of_blocks
              call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                                 ,Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1) &
                                 ,'BLK','1','block indices' &
                                 ,data=(/(dble(real(il_i)),il_i=ll,lh)/))
            endif
         endif
!
!       Declare axis of bounds.
!
        Nullify ( corner_pointer )
        corner_pointer => Grids(il_grid_id)%corner_pointer

        if(associated(corner_pointer%corners_real(1)%vector)) then
          len_alloc(1)=size(corner_pointer%corners_real(1)%vector)
        else if(associated(corner_pointer%corners_dble(1)%vector)) then
          len_alloc(1)=size(corner_pointer%corners_dble(1)%vector)
        else
          ierror=PRISM_Error_IO_Meta
          call psmile_error(ierror,'No corners specified!', &
                           ierrp,0, __FILE__, __LINE__)
        endif

        if(associated(corner_pointer%corners_real(2)%vector)) then
          len_alloc(2)=size(corner_pointer%corners_real(2)%vector)
        else if(associated(corner_pointer%corners_dble(2)%vector)) then
          len_alloc(2)=size(corner_pointer%corners_dble(2)%vector)
        else
          ierror=PRISM_Error_IO_Meta
          call psmile_error(ierror,'No corners specified!', &
                           ierrp,0, __FILE__, __LINE__)
        endif

        if(associated(corner_pointer%corners_real(3)%vector)) then
          len_alloc(3)=size(corner_pointer%corners_real(3)%vector)
        else if(associated(corner_pointer%corners_dble(3)%vector)) then
          len_alloc(3)=size(corner_pointer%corners_dble(3)%vector)
        else
          ierror=PRISM_Error_IO_Meta
          call psmile_error(ierror,'No corners specified!', &
                           ierrp,0, __FILE__, __LINE__)
        endif

        noc=corner_pointer%nbr_corners

        print*,'psmile_write_meta_byid: len_alloc',len_alloc

        select case ( Grids(il_grid_id)%grid_type )
  
        case ( PRISM_Irrlonlat_sigmavrt)
  
           nof(1:2) = noc/2
           nof(ndim_3d) = 2

           do jj=1,3
             ncells(jj)=len_alloc(jj)/nof(jj)
             nvertices(jj)=noc
           enddo

           ncells(1)=ncells(1)/(corner_pointer%corner_shape(2,2) - &
                                corner_pointer%corner_shape(1,2)+1)
           ncells(2)=ncells(2)/(corner_pointer%corner_shape(2,1) - &
                                corner_pointer%corner_shape(1,1)+1)
           nvertices(3)=nof(3)

           is_regvrt=.false.
  
        case ( PRISM_Irrlonlat_regvrt )
  
           nof(1:2)     = noc/2
           nof(ndim_3d) = 2

           do jj=1,3
             ncells(jj)=len_alloc(jj)/nof(jj)
             nvertices(jj)=noc
           enddo

           ncells(1)=ncells(1)/(corner_pointer%corner_shape(2,2) - &
                                corner_pointer%corner_shape(1,2)+1)
           ncells(2)=ncells(2)/(corner_pointer%corner_shape(2,1) - &
                                corner_pointer%corner_shape(1,1)+1)
           nvertices(3)=nof(3)

           is_regvrt=.true.
  
        case ( PRISM_Reglonlat_sigmavrt )
  
           nof(1:ndim_3d) = 2

           do jj=1,3
             ncells(jj)=len_alloc(jj)/nof(jj)
             nvertices(jj)=noc
           enddo
           nvertices(3)=nof(3)


!RV           if(lenaz.eq.1) then
!RV             do jj=1,2
!RV               nvertices(jj)=max(nvertices(jj)/nof(jj),nof(jj))
!RV             enddo
!RV           endif

           is_regvrt=.false.
  
        case ( PRISM_Reglonlatvrt )
  
           nof(1:ndim_3d) = 2

           do jj=1,3
             ncells(jj)=len_alloc(jj)/nof(jj)
             nvertices(jj)=noc
           enddo
           nvertices(3)=nof(3)

!RV           if(lenaz.eq.1) then
!RV             do jj=1,2
!RV               nvertices(jj)=max(nvertices(jj)/nof(jj),nof(jj))
!RV             enddo
!RV           endif

           is_regvrt=.true.

        case ( PRISM_Irrlonlatvrt )
  
           nof(1:ndim_3d) = noc

           do jj=1,3
             ncells(jj)=len_alloc(jj)/nof(jj) &
                                     /(corner_pointer%corner_shape(2,3) - &
                                       corner_pointer%corner_shape(1,3)+1)
             nvertices(jj)=noc
           enddo
           ncells(1)=ncells(1)/(corner_pointer%corner_shape(2,2) - &
                                corner_pointer%corner_shape(1,2)+1)
           ncells(2)=ncells(2)/(corner_pointer%corner_shape(2,1) - &
                                corner_pointer%corner_shape(1,1)+1)
           ncells(3)=(corner_pointer%corner_shape(2,3) - &
                      corner_pointer%corner_shape(1,3)+1)
           is_regvrt=.false.
  
        case ( PRISM_Unstructlonlatvrt )
  
           nof(1:ndim_3d) = noc

           do jj=1,3
             ncells(jj)=len_alloc(jj)/nof(jj)
             nvertices(jj)=noc
           enddo
           is_regvrt=.false.
  
        case ( PRISM_Unstructlonlat_regvrt )
  
           nof(1:2) = noc/2
           nof(ndim_3d) = 2

           do jj=1,3
             ncells(jj)=len_alloc(jj)/nof(jj)
             nvertices(jj)=noc
           enddo
  
           nvertices(3)=nof(3)
           is_regvrt=.true.

           ierrp (1) = Grids(il_grid_id)%grid_type
  
           ierror = PRISM_Error_Internal
  
           call psmile_error(ierror,'grid generation type not supported yet', &
                             ierrp, 1, __FILE__, __LINE__ )
           return
  
        case ( PRISM_Unstructlonlat_sigmavrt )
  
           nof(1:2) = noc/2
           nof(ndim_3d) = 2

           do jj=1,3
             ncells(jj)=len_alloc(jj)/nof(jj)
             nvertices(jj)=noc
           enddo
           nvertices(3)=nof(3)
           is_regvrt=.false.
  
  
        case ( PRISM_Gaussreduced_regvrt )
  
       nof(1:2)     = 2
           nof(ndim_3d) = 2

           do jj=1,3
             ncells(jj)=len_alloc(jj)/nof(jj)
             nvertices(jj)=noc
           enddo
  
           nvertices(3)=nof(3)
           is_regvrt=.true.

  
        case ( PRISM_Gaussreduced_sigmavrt )
  
           nof(1:2)     = 2
           nof(ndim_3d) = 2

           do jj=1,3
             ncells(jj)=len_alloc(jj)/nof(jj)
             nvertices(jj)=noc
           enddo
           nvertices(3)=nof(3)
           is_regvrt=.false.
  
  
        case DEFAULT
  
           ierrp (1) = Grids(il_grid_id)%grid_type
  
           ierror = PRISM_Error_Internal
  
           call psmile_error(ierror, 'unknown grid generation type', &
                               ierrp, 1, __FILE__, __LINE__ )
           return
  
        end select

        call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                           ,Fields(id_varid)%io_infos%p_mpp_io%crnaxis(1) &
                          ,'NCORNERSLON','1','corner indices for longitudes' &
                          ,data=(/(dble(real(il_i)),il_i=1,nvertices(1))/))
        call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                           ,Fields(id_varid)%io_infos%p_mpp_io%crnaxis(2) &
                          ,'NCORNERSLAT','1','corner indices for latitudes ' &
                          ,data=(/(dble(real(il_i)),il_i=1,nvertices(2))/))
        call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                           ,Fields(id_varid)%io_infos%p_mpp_io%crnaxis(3) &
                          ,'NCORNERSZ','1','corner indices for z' &
                          ,data=(/(dble(real(il_i)),il_i=1,nvertices(3))/))

        
!
!       Time axis
!
        call psmile_io_get_jobstart_date(ljobstart,il_date,ierror) 
!
!       Remember the date as a Julian date in order to calculate the
!       the offset  of a Gregorian time  stamp with respect to the 
!       date of the job start.
!
        call psmile_date2ju (ljobstart &
                            ,Fields(id_varid)%io_infos%ju_start_day &
                            ,Fields(id_varid)%io_infos%ju_start_sec)
        call combine_with_date('second','since',il_date,cl_timeunit)

        call mpp_write_meta( Fields(id_varid)%io_infos%file_unit &
                           , Fields(id_varid)%io_infos%p_mpp_io%t_axis(1) &
                           , 'time',trim(cl_timeunit), 'Time' )
!--------------------------------------------------------------------------
!
!       NetCDF CF header informations
!
!
!       Assign field components, lats, lons, verts and mask to axis
!
        Select Case (il_grid_type)
        Case(PRISM_Reglonlatvrt)
          cl_vertij='regvert'
          cl_vert_longname='regular vertical'
          cl_grid_longname='PRISM_REGLONLATVRT'
          il_gridparam=PRISM_Reglonlatvrt
        Case(PRISM_Irrlonlat_regvrt)
          cl_vertij='regvert'
          cl_vert_longname='regular vertical'
          cl_grid_longname='PRISM_IRRLONLAT_REGVRT'
          il_gridparam=PRISM_Irrlonlat_regvrt
        Case(PRISM_Irrlonlatvrt)
          cl_vertij='irregvert'
          cl_vert_longname='irregular vertical'
          cl_grid_longname='PRISM_IRRLONLATVRT'
          il_gridparam=PRISM_Irrlonlat_regvrt
        Case(PRISM_Irrlonlat_sigmavrt)
          cl_vertij='sigmavert'
          cl_vert_longname='sigma vertical'
          cl_grid_longname='PRISM_IRRLONLAT_SIGMAVRT'
          il_gridparam=PRISM_Irrlonlat_sigmavrt
        Case(PRISM_Reglonlat_sigmavrt)
          cl_vertij='sigmavert'
          cl_vert_longname='sigma vertical'
          cl_grid_longname='PRISM_REGLONLAT_SIGMAVRT'
          il_gridparam=PRISM_Reglonlat_sigmavrt
        Case(PRISM_Unstructlonlat_regvrt)
          cl_vertij='regvert'
          cl_vert_longname='regular vertical'
          cl_grid_longname='PRISM_UNSTRUCTLONLAT_REGVRT'
          il_gridparam=PRISM_Unstructlonlat_regvrt
        Case(PRISM_Unstructlonlat_sigmavrt)
          cl_vertij='sigmavert'
          cl_vert_longname='sigma vertical'
          cl_grid_longname='PRISM_UNSTRUCTLONLAT_SIGMAVRT'
          il_gridparam=PRISM_Unstructlonlat_sigmavrt
        Case(PRISM_Unstructlonlatvrt)
          cl_vertij='unstructvert'
          cl_vert_longname='unstructured vertical'
          cl_grid_longname='PRISM_UNSTRUCTLONLATVRT'
          il_gridparam=PRISM_Unstructlonlatvrt
        Case(PRISM_Gaussreduced_regvrt)
          cl_vertij='regvert'
          cl_vert_longname='regular vertical'
          cl_grid_longname='PRISM_Gaussreduced_regvrt'
          il_gridparam=PRISM_Gaussreduced_regvrt
        Case(PRISM_Gaussreduced_sigmavrt)
          cl_vertij='sigmavrt'
          cl_vert_longname='sigma vertical'
          cl_grid_longname='PRISM_Gaussreduced_sigmavrt'
          il_gridparam=PRISM_Gaussreduced_sigmavrt
        Case DEFAULT

         ierror=PRISM_Error_IO_Meta
         ierrp(1)=il_grid_type
         call psmile_error(ierror,'Vector field on unkown grid type!', &
                           ierrp,1, __FILE__, __LINE__)
        End Select
!
!Header describing the grid type
!
        print*,'il_gridparam',il_gridparam
        call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                       ,Fields(id_varid)%io_infos%p_mpp_io%x_axis(2) &
                           ,'PRISM_GRIDTYPE','1' &
                           ,trim(cl_grid_longname(1)) &
                           ,data=dble(il_gridparam(1:1)))
!
!       Header of bounds of longitudes and latitudes
!
        if(associated(corner_pointer%corners_real(1)%vector)) then
          il_pack=2
        else if(associated(corner_pointer%corners_dble(1)%vector)) then
          il_pack=1
        endif

        if(il_grid_type.eq.PRISM_Unstructlonlatvrt .or. &
            il_grid_type .eq. PRISM_Gaussreduced_regvrt  .or. &  
            il_grid_type .eq. PRISM_Gaussreduced_sigmavrt  .or. &  
           il_grid_type.eq.PRISM_Unstructlonlat_sigmavrt ) then
          print*,'il_grid_type',il_grid_type

          is_regular=.FALSE.
!
!         For these grid types the assignment to the x axis
!         serves as the information about the domain decomposition
!         of the horizontal space.
!
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%crnij(1) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%crnaxis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%x_axis(1)    &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)/)   &
                        ,  'bounds_lon','degrees_east','longitude' &
                        ,  pack=il_pack)
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%crnij(2) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%crnaxis(2) &
                        ,Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)/) &
                        ,  'bounds_lat','degrees_north','latitude' &
                        ,  pack=il_pack)

          if(is_regvrt) then
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%crnaxis(3) &
                        ,Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)/) &
                        ,  'bounds_'//trim(cl_vertij(1)) &
                        ,  trim(Fields(id_varid)%io_infos%height_unit) &
                        ,  'bounds of '//trim(cl_vert_longname(1)) &
                        ,  pack=il_pack)
          else
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%crnaxis(3) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1) &
                        , Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)/) &
                        ,  'bounds_'//trim(cl_vertij(1)) &
                        ,  trim(Fields(id_varid)%io_infos%height_unit) &
                        ,  'bounds of '//trim(cl_vert_longname(1)) &
                        ,  pack=il_pack)
          endif 
        else if( il_grid_type.eq.PRISM_Irrlonlat_sigmavrt .or. &
                 il_grid_type.eq.PRISM_Irrlonlat_regvrt .or. &
                 il_grid_type.eq.PRISM_Irrlonlatvrt .or. &
                 il_grid_type.eq.PRISM_Reglonlat_sigmavrt .or. &
                 il_grid_type.eq.PRISM_Reglonlatvrt ) then

          is_regular=.TRUE.
!
!         For these grid types the rule is:
!         Number of cells for longitudes is equal  to 
!         the extent of the domain in x.
!         Number of cells for latitudes is equal  to 
!         the extent of the domain in y.
!

          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%crnij(1) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%crnaxis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)/) &
                        ,  'bounds_lon','degrees_east','longitude' &
                        ,  pack=il_pack)
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%crnij(2) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%crnaxis(2) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)/) &
                        ,  'bounds_lat','degrees_north','latitude' &
                        ,  pack=il_pack)

          if(is_regvrt) then
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                      ,(/Fields(id_varid)%io_infos%p_mpp_io%crnaxis(3) &
                      ,Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)/)   &
                        ,  'bounds_'//trim(cl_vertij(1)) &
                        ,  trim(Fields(id_varid)%io_infos%height_unit) &
                        ,  'bounds of '//trim(cl_vert_longname(1)) &
                        ,  pack=il_pack)
          else
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%crnaxis(3) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%x_axis(1)  &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)  &
                        ,Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)/)  &
                        ,  'bounds_'//trim(cl_vertij(1)) &
                        ,  trim(Fields(id_varid)%io_infos%height_unit) &
                        ,  'bounds of '//trim(cl_vert_longname(1)) &
                        ,  pack=il_pack)
           endif

         
        endif

        do ii=1,nvcomp

          vmsk_id(ii)=PRISM_UNDEFINED
          if(nvcomp.gt.1) then
            write(label,fmt='(i4.4)')ii
            cl_field(ii)=trim(Fields(id_varid)%io_infos%cfioname)//trim(label)
            cl_latij(ii)='lat'//trim(label)
            cl_lonij(ii)='lon'//trim(label)
            cl_vertij(ii)=trim(cl_vertij(ii))//trim(label)
            cl_vert_longname(ii)=trim(cl_vert_longname(ii))//' '//trim(label)
            cl_lonij(ii)='mask'//trim(label)
            vcomp_id=Methods(Fields(id_varid)%method_id)%vector_pointer &
                    %array_of_point_ids(ii)
!
          else
!
            cl_field(ii)=trim(Fields(id_varid)%io_infos%cfioname)
            cl_latij(ii)='lat'
            cl_lonij(ii)='lon'
            cl_mskij(ii)='mask'
            vcomp_id=Fields(id_varid)%method_id
            vmsk_id(ii)=Fields(id_varid)%mask_id
!
          endif

!
!         I assume that lons comes first, then lats and verts.
!         Pack modes are determined by the informations passed through
!         the PSMILe API.
!
!         Header of the Longitudes
!
          if(associated( &
             Methods(vcomp_id)%coords_pointer%coords_real(1)%vector)) then
             il_pack=2
          else
             il_pack=1
          endif

          if((il_grid_type .eq. PRISM_Unstructlonlatvrt) .or. &
             (il_grid_type .eq. PRISM_Irrlonlatvrt)) then
          if(no_of_blocks.eq.1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%lonij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)/)   &
                        ,  trim(cl_lonij(ii)),'degrees_east','longitude' &
                        ,  pack=il_pack)
          else if(no_of_blocks .gt. 1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%lonij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1)/)   &
                        ,  trim(cl_lonij(ii)),'degrees_east','longitude' &
                        ,  pack=il_pack)
          endif
          else
          if(no_of_blocks.eq.1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%lonij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)/)   &
                        ,  trim(cl_lonij(ii)),'degrees_east','longitude' &
                        ,  pack=il_pack)
          else if(no_of_blocks .gt. 1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%lonij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1)/)   &
                        ,  trim(cl_lonij(ii)),'degrees_east','longitude' &
                        ,  pack=il_pack)
          endif
          endif
!
!         Set the attribute bounds for the longitudes
!
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                    ,mpp_get_id(Fields(id_varid)%io_infos%p_mpp_io%lonij(ii)) &
                    ,'bounds' &
                    ,cval=trim('bounds_lon'))

!
!         Header of the Latitudes
!
          if(associated( &
             Methods(vcomp_id)%coords_pointer%coords_real(2)%vector)) then
             il_pack=2
          else
             il_pack=1
          endif
          if((il_grid_type .eq. PRISM_Unstructlonlatvrt) .or. &
             (il_grid_type .eq. PRISM_Irrlonlatvrt)) then
          if(no_of_blocks.eq.1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%latij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)/) &
                        ,  trim(cl_latij(ii)),'degrees_north','latitude' &
                        ,  pack=il_pack)
          else if(no_of_blocks .gt. 1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%latij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1)/) &
                        ,  trim(cl_latij(ii)),'degrees_north','latitude' &
                        ,  pack=il_pack)
          endif
          else
          if(no_of_blocks.eq.1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%latij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)/) &
                        ,  trim(cl_latij(ii)),'degrees_north','latitude' &
                        ,  pack=il_pack)
          else if(no_of_blocks .gt. 1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%latij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1)/) &
                        ,  trim(cl_latij(ii)),'degrees_north','latitude' &
                        ,  pack=il_pack)
          endif
          endif
!
!         Set the attribute bounds for the latitudes.
!
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                    ,mpp_get_id(Fields(id_varid)%io_infos%p_mpp_io%latij(ii)) &
                    ,'bounds' &
                    ,cval=trim('bounds_lat'))
!
!         Header of the Verticals
!
          if(associated( &
             Methods(vcomp_id)%coords_pointer%coords_real(3)%vector)) then
             il_pack=2
          else
             il_pack=1
          endif
          if((il_grid_type .eq. PRISM_Unstructlonlatvrt) .or. &
             (il_grid_type .eq. PRISM_Irrlonlatvrt)) then
          if(no_of_blocks.eq.1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)/)   &
                        ,  trim(cl_vertij(ii))                   &
                        ,  trim(Fields(id_varid)%io_infos%height_unit) &
                        ,  trim(cl_vert_longname(ii)) &
                        ,  pack=il_pack)
          else if(no_of_blocks .gt. 1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1)/)   &
                        ,  trim(cl_vertij(ii))                   &
                        ,  trim(Fields(id_varid)%io_infos%height_unit) &
                        ,  trim(cl_vert_longname(ii)) &
                        ,  pack=il_pack)
          endif
          else
          if(no_of_blocks.eq.1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)/) &
                        ,  trim(cl_vertij(ii))                   &
                        ,  trim(Fields(id_varid)%io_infos%height_unit) &
                        ,  trim(cl_vert_longname(ii)) &
                        ,  pack=il_pack)
          else if(no_of_blocks .gt. 1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%z_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1)/)   &
                        ,  trim(cl_vertij(ii))                   &
                        ,  trim(Fields(id_varid)%io_infos%height_unit) &
                        ,  trim(cl_vert_longname(ii)) &
                        ,  pack=il_pack)
          endif
          endif
!
!         Add the standard name to the vertical
!
          if(trim(adjustl(Fields(id_varid)%io_infos%height_stdname)) &
             .ne. '') then
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                    ,mpp_get_id(Fields(id_varid)%io_infos%p_mpp_io%vertij(ii)) &
                    ,'standard_name' &
                    ,cval=trim(Fields(id_varid)%io_infos%height_stdname))
          endif
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                    ,mpp_get_id(Fields(id_varid)%io_infos%p_mpp_io%vertij(ii)) &
                    ,'positive' &
                    ,cval=trim(Fields(id_varid)%io_infos%positive))

!
!         Add the formular terms to the attributes of the verticals
!
          if(trim(adjustl(Fields(id_varid)%io_infos%height_formular)) &
             .ne. '' ) then

            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                    ,mpp_get_id(Fields(id_varid)%io_infos%p_mpp_io%vertij(ii)) &
                    ,'formula_terms' &
                    ,cval=trim(Fields(id_varid)%io_infos%height_formular))

          endif
!
!         Set the attribute bounds for the vertical.
!
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                    ,mpp_get_id(Fields(id_varid)%io_infos%p_mpp_io%vertij(ii)) &
                    ,'bounds' &
                    ,cval='bounds_'//trim(cl_vertij(ii)))
!
!         Header of Masks
!
!         Masks are always 3d arrays. Assign all 3 spatial axis to masks.
!         I am writing mask as FLOAT  0. = .FALSE., 1.=.TRUE.
!
          if(vmsk_id(ii).ne.PRISM_UNDEFINED) then
          il_pack=2
          if(no_of_blocks.eq.1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%maskij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)/)   &
                        ,  trim(cl_mskij(ii)),'1','mask' &
                        ,  pack=il_pack)
          else if(no_of_blocks .gt. 1) then
            call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                        ,Fields(id_varid)%io_infos%p_mpp_io%maskij(ii) &
                        ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)   &
                        ,  Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1)/)   &
                        ,  trim(cl_mskij(ii)),'1','mask' &
                        ,  pack=il_pack)
          endif
          endif

!
!        Header of the Field
!
!rv
          il_pack=Fields(id_varid)%io_infos%pack_mode
!rv       The following coherence is performed in psmile_io_get_attr.F90
!rv          il_pack=1
!rv          select case(Fields(id_varid)%datatype)
!rv          case(PRISM_Integer)
!rv            il_pack=2
!rv          case(PRISM_Logical)
!rv            il_pack=2
!rv          case(PRISM_Real)
!rv            il_pack=2
!rv            if(Fields(id_varid)%io_infos%pack_mode .gt. il_pack) then
!rv              print*,trim(ch_id),' :  psmile_write_meta_byid: Warning:' &
!rv              ,'SMIOC pack mode not consistent with PRISM_Real'
!rv            endif
!rv          case(PRISM_Double_Precision)
!rv            il_pack=1
!rv            if(Fields(id_varid)%io_infos%pack_mode .gt. il_pack) then
!rv              print*,trim(ch_id),' :  psmile_write_meta_byid: Warning:' &
!rv              ,'SMIOC pack mode not consistent with PRISM_Double_Precision'
!rv            endif
!rv          case DEFAULT
!rv           ierrp (1) = Fields(id_varid)%datatype
!rv           ierrp (2) = 1
!rv           ierror =  PRISM_Error_IO_Meta
!rv           call PSMILe_Error ( ierror, 'Data type not yet supported!', &
!rv                           ierrp, 1, __FILE__, __LINE__ )
!rv           return
!rv
!rv          end select

          if(.not.is_bundle) then
            if(no_of_blocks.eq.1) then
              call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                          ,Fields(id_varid)%io_infos%p_mpp_io%field(ii) &
                          ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)   &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%t_axis(1)/) &
                          ,  trim(cl_field(ii))                    &
                          ,  trim(Fields(id_varid)%io_infos%units) &
                          ,  trim(Fields(id_varid)%io_infos%long_name) &
                          ,  min=Fields(id_varid)%io_infos%valid_min &
                          ,  max=Fields(id_varid)%io_infos%valid_max &
                          ,  missing=Fields(id_varid)%io_infos%missing_value &
                          ,  fill=Fields(id_varid)%io_infos%fill_value &
                          ,  scale=Fields(id_varid)%io_infos%scale &
                          ,  add=Fields(id_varid)%io_infos%add &
                          ,  pack=il_pack )
             call psmile_flushstd
            else if(no_of_blocks .gt. 1) then
              call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                          ,Fields(id_varid)%io_infos%p_mpp_io%field(ii) &
                          ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)   &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1)   &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%t_axis(1)/) &
                          ,  trim(cl_field(ii))                    &
                          ,  trim(Fields(id_varid)%io_infos%units) &
                          ,  trim(Fields(id_varid)%io_infos%long_name) &
                          ,  min=Fields(id_varid)%io_infos%valid_min &
                          ,  max=Fields(id_varid)%io_infos%valid_max &
                          ,  missing=Fields(id_varid)%io_infos%missing_value &
                          ,  fill=Fields(id_varid)%io_infos%fill_value &
                          ,  scale=Fields(id_varid)%io_infos%scale &
                          ,  add=Fields(id_varid)%io_infos%add &
                          ,  pack=il_pack)
            endif
          else
            if(no_of_blocks.eq.1) then
              call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                          ,Fields(id_varid)%io_infos%p_mpp_io%field(ii) &
                          ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)   &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%b_axis(1) &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%t_axis(1)/) &
                          ,  trim(cl_field(ii))                    &
                          ,  trim(Fields(id_varid)%io_infos%units) &
                          ,  trim(Fields(id_varid)%io_infos%long_name) &
                          ,  min=Fields(id_varid)%io_infos%valid_min &
                          ,  max=Fields(id_varid)%io_infos%valid_max &
                          ,  missing=Fields(id_varid)%io_infos%missing_value &
                          ,  fill=Fields(id_varid)%io_infos%fill_value &
                          ,  scale=Fields(id_varid)%io_infos%scale &
                          ,  add=Fields(id_varid)%io_infos%add &
                          ,  pack=Fields(id_varid)%io_infos%pack_mode )
            else if(no_of_blocks .gt. 1) then
              call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                          ,Fields(id_varid)%io_infos%p_mpp_io%field(ii) &
                          ,(/Fields(id_varid)%io_infos%p_mpp_io%x_axis(1) &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%y_axis(1)   &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%z_axis(1)   &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%b_axis(1) &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1) &
                          ,  Fields(id_varid)%io_infos%p_mpp_io%t_axis(1)/) &
                          ,  trim(cl_field(ii))                    &
                          ,  trim(Fields(id_varid)%io_infos%units) &
                          ,  trim(Fields(id_varid)%io_infos%long_name) &
                          ,  min=Fields(id_varid)%io_infos%valid_min &
                          ,  max=Fields(id_varid)%io_infos%valid_max &
                          ,  missing=Fields(id_varid)%io_infos%missing_value &
                          ,  fill=Fields(id_varid)%io_infos%fill_value &
                          ,  scale=Fields(id_varid)%io_infos%scale &
                          ,  add=Fields(id_varid)%io_infos%add &
                          ,  pack=Fields(id_varid)%io_infos%pack_mode )
            endif
          endif
          if(nvcomp.eq.1) then
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                    ,mpp_get_id(Fields(id_varid)%io_infos%p_mpp_io%field(ii)) &
                             ,'standard_name' &
                          ,cval=trim(Fields(id_varid)%io_infos%standard_name))
          else
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                    ,mpp_get_id(Fields(id_varid)%io_infos%p_mpp_io%field(ii)) &
                             ,'standard_name' &
                          ,cval=trim(Fields(id_varid)%io_infos%vcmp_names(ii)))
          endif
          call mpp_write_meta(Fields(id_varid)%io_infos%file_unit &
                    ,mpp_get_id(Fields(id_varid)%io_infos%p_mpp_io%field(ii)) &
                             ,'coordinates' &
                             ,cval=trim(cl_lonij(ii))//' '//trim(cl_latij(ii)) &
                                  //' '//trim(cl_vertij(ii)) )
          
          
        enddo
        

        
!
!-----------------------------------------------------------------------------
!    Here is the end of the NetCDF header section.
!
!    Now dump the x_axis,y_axis,z_axis,lons, lats, the verts &
!    and masks - if any - to the file.
!
!
!    The common axis
!
!       X
        call mpp_write(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%x_axis(1))
!       Y
        call mpp_write(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%y_axis(1))
!       Z
        call mpp_write(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%z_axis(1))
!       NVERTICES
        call mpp_write(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%crnaxis(1))
        call mpp_write(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%crnaxis(2))
        call mpp_write(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%crnaxis(3))
!       PRISM_GRIDTYPE
        call mpp_write(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%x_axis(2))
        
!       B
        if(is_bundle) then
        call mpp_write(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%b_axis(1))
        endif
!      BLOCK
       if(no_of_blocks.gt.1) then
         call mpp_write(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%blk_axis(1))
         myashape(1,4)=1
         myashape(2,4)=no_of_blocks
         myvshape(1,4)=1
         myvshape(2,4)=no_of_blocks
!
!        Stitching together the blocks of lons, lats and verts
!
         allocate(fdbles(no_of_blocks,nvcomp,3),stat=ierror)
         if (ierror > 0) then
            ierrp (1) = ierror
            ierrp (2) = 1
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'fdbles', &
                            ierrp, 2, __FILE__, __LINE__ )
            return
         endif
!
         allocate(freals(no_of_blocks,nvcomp,3),stat=ierror)
         if (ierror > 0) then
            ierrp (1) = ierror
            ierrp (2) = 1
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'freals', &
                            ierrp, 2, __FILE__, __LINE__ )
            return
         endif

         allocate(flogs(no_of_blocks,nvcomp),stat=ierror)
         if (ierror > 0) then
            ierrp (1) = ierror
            ierrp (2) = 1
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'flogs', &
                            ierrp, 2, __FILE__, __LINE__ )
            return
         endif
!
         do ii=1,nvcomp
           do jj=1,no_of_blocks
             il_id=Fields(id_varid)%io_infos%related_ids(jj)
             blkid=Fields(il_id)%io_infos%block_id

             vmsk_id(ii)=PRISM_UNDEFINED
             if(nvcomp.gt.1) then
               vcomp_id=Methods(Fields(il_id)%method_id)%vector_pointer &
                     %array_of_point_ids(ii)
             else
               vcomp_id=Fields(il_id)%method_id
               vmsk_id(ii) =Fields(il_id)%mask_id
             endif

#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid: Stitching blocks:' &
            ,'vector component, block, block_id :',ii,jj,blkid
      call psmile_flushstd
#endif
             if(associated( &
                Methods(vcomp_id)%coords_pointer%coords_dble(1)%vector))then 
                fdbles(blkid,ii,1)%fdble => Methods(vcomp_id)%coords_pointer &
                                                          %coords_dble(1)%vector
             else if(associated( &
                Methods(vcomp_id)%coords_pointer%coords_real(1)%vector))then
                freals(blkid,ii,1)%freal => Methods(vcomp_id)%coords_pointer &
                                                          %coords_real(1)%vector
             endif
             if(associated( &
                Methods(vcomp_id)%coords_pointer%coords_dble(2)%vector))then 
                fdbles(blkid,ii,2)%fdble => Methods(vcomp_id)%coords_pointer &
                                                          %coords_dble(2)%vector
             else if(associated( &
                Methods(vcomp_id)%coords_pointer%coords_real(2)%vector))then
                freals(blkid,ii,2)%freal => Methods(vcomp_id)%coords_pointer &
                                                          %coords_real(2)%vector
             endif
             if(associated( &
                Methods(vcomp_id)%coords_pointer%coords_dble(3)%vector))then 
                fdbles(blkid,ii,3)%fdble => Methods(vcomp_id)%coords_pointer &
                                                          %coords_dble(3)%vector
             else if(associated( &
                Methods(vcomp_id)%coords_pointer%coords_real(3)%vector))then
                freals(blkid,ii,3)%freal => Methods(vcomp_id)%coords_pointer &
                                                          %coords_real(3)%vector
             endif
             if(vmsk_id(ii).ne.PRISM_UNDEFINED) then
               flogs(blkid,ii)%flog => Masks(vmsk_id(ii))%mask_array
             endif
           enddo
         enddo
       endif
!
!      Writing the bounds for longitudes and latitudes
!
       if(associated(corner_pointer%corners_real(1)%vector)) then

         if(is_regular) then
           allocate(bound_lonr(nvertices(1),ncells(1),ncells(2)),stat=ierror)
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_lonr', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif
  
           allocate(bound_latr(nvertices(2),ncells(1),ncells(2)),stat=ierror)
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_latr', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif

           if (lenaz .gt. 1) then
             if(is_regvrt) then
               allocate(bound_vertr(nvertices(3),ncells(3),1,1) &
                   ,stat=ierror)
              else
                allocate(bound_vertr(nvertices(3),ncells(1),ncells(2) &
                        ,ncells(3)) ,stat=ierror)
              endif
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_latr', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif
           endif
!
         else
!
           allocate(bound_lonr(nvertices(1),ncells(1),1),stat=ierror)
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_lonr', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif
  
           allocate(bound_latr(nvertices(2),ncells(2),1),stat=ierror)
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_latr', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif

           if (lenaz .gt. 1) then
           allocate(bound_vertr(nvertices(3),ncells(3),1,1) &
                   ,stat=ierror)
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_latr', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif
           endif
!
         endif

#ifdef DEBUG
     print*,trim(ch_id),' :  psmile_write_meta_byid: ','is_regular ',is_regular
     print*,trim(ch_id),' :  psmile_write_meta_byid: '  &
              ,'ih-il+1,jh-jl+1,nvertices(1),ncells(1),nvertices(2),ncells(2)' &
              ,ih-il+1,jh-jl+1,nvertices(1),ncells(1),nvertices(2),ncells(2)
     print*,trim(ch_id),' :  psmile_write_meta_byid: ' &
                       ,'corner_pointer%corner_shape ' &
                       ,corner_pointer%corner_shape
      call psmile_flushstd
#endif
!
!        Branch for all logical rectangular grids
!
         if(is_regular) then

           if( il_grid_type == PRISM_Reglonlatvrt .or.  &
               il_grid_type == PRISM_Reglonlat_sigmavrt ) then

             lbackward=.false.
             do iloc=1,ncells(2)
               do ii=1,ncells(1)
                 ll=0
               do jj=1,nvertices(1)/nof(1)
                 lbackward=mod(jj-1,2).eq.0
                 if(.not.lbackward)then
                   do kk=1,nof(1)
                     ll=ll+1
                     bound_lonr(ll,ii,iloc)= &
                     corner_pointer%corners_real(1)%vector((kk-1)*ncells(1)+ii)
                   enddo
                 else
                   do kk=nof(1),1,-1
                     ll=ll+1
                     bound_lonr(ll,ii,iloc)= &
                     corner_pointer%corners_real(1)%vector((kk-1)*ncells(1)+ii)
                   enddo
                 endif
               enddo
              enddo
             enddo
  
             kkj=0
             do iloc=1,ncells(1)
               do ii=1,ncells(2)
                 ll=0
                 do jj=1,nvertices(2)/nof(2)
                   kkj=kkj+1
                   do kk=1,nof(2)
                     ll=ll+1
                     bound_latr(ll,iloc,ii)= &
                     corner_pointer%corners_real(2)%vector((kkj-1)*ncells(2)+ii)
                   enddo
                   if(kkj.eq.nof(2))kkj=0
                 enddo

               enddo
             enddo

           else

             stride = ncells(1)*ncells(2)

             ii = 0
             do jl=1,ncells(2)
                do il=1,ncells(1)
                   ii = ii + 1
                   ll=0
                   do jj=1,nvertices(1)/nof(1)
                      do kk=1,nof(1)
                         ll=ll+1
                         bound_lonr(ll,il,jl)= &
                         corner_pointer%corners_real(1)%vector((kk-1)*stride+ii)
                      enddo
                   enddo
                enddo
             enddo

             stride = ncells(1)*ncells(2)

             ii = 0
             do jl=1,ncells(2)
                do il=1,ncells(1)
                   ii = ii + 1
                   ll = 0
                   do jj=1,nvertices(2)/nof(2)
                      do kk=1,nof(2)
                         ll=ll+1
                         bound_latr(ll,il,jl)= &
                         corner_pointer%corners_real(2)%vector((kk-1)*stride+ii)
                      enddo
                   enddo
                enddo
             enddo


           endif ! il_grid_type == PRISM_Reglonlatvrt ||  &
                 ! il_grid_type == PRISM_Reglonlat_sigmavrt 
!
!  Extraction of height data for grids with logical rectangular
!  horizontal layers.
!
           if (lenaz .gt. 1 ) then
           if( il_grid_type == PRISM_Reglonlatvrt .or.  &
               il_grid_type == PRISM_Irrlonlat_regvrt ) then

             do ii=1,ncells(3)
                   ll=0
                   do jj=1,nvertices(3)/nof(3)
                     do kk=1,nof(3)
                       ll=ll+1
                       bound_vertr(ll,ii,1,1)= &
                       corner_pointer%corners_real(3) &
                       %vector((kk-1)*ncells(3)+ii)
                     enddo
                   enddo
             enddo

           else

             stride = ncells(1)*ncells(2)*ncells(3)

             ii = 0
             do kl=1,ncells(3)
               do jl=1,ncells(2)
                 do il=1,ncells(1)
                   ii = ii + 1
                   ll = 0
                   do jj=1,nvertices(3)/nof(3)
                      do kk=1,nof(3)
                         ll=ll+1
                         bound_vertr(ll,il,jl,kl)= &
                         corner_pointer%corners_real(3)%vector((kk-1)*stride+ii)
                      enddo
                   enddo
                 enddo
               enddo
             enddo

           endif
           endif

!
!       Branch for all grids which are not (logical) rectangular.
!
        else

           !Rene, do something here for reordering!!!

           if (il_grid_type == PRISM_Gaussreduced_regvrt ) then

              do ii=1,ncells(1)
                 ll=0
                 do jj=1,nof(1)
                    bound_lonr(ll+1,ii,1) = &
                         corner_pointer%corners_real(1)%vector((1-1)*ncells(1)+ii)
                    bound_lonr(ll+2,ii,1) = &
                         corner_pointer%corners_real(1)%vector((2-1)*ncells(1)+ii)
                    bound_lonr(ll+3,ii,1) = &
                         corner_pointer%corners_real(1)%vector((2-1)*ncells(1)+ii)
                    bound_lonr(ll+4,ii,1) = &
                         corner_pointer%corners_real(1)%vector((1-1)*ncells(1)+ii)
                    ll=ll+4
                 enddo
              enddo

              do ii=1,ncells(2)
                 ll=0
                 do jj=1,nof(2)
                    bound_latr(ll+1,ii,1)= &
                         corner_pointer%corners_real(2)%vector((1-1)*ncells(2)+ii)
                    bound_latr(ll+2,ii,1)= &
                         corner_pointer%corners_real(2)%vector((1-1)*ncells(2)+ii)
                    bound_latr(ll+3,ii,1)= &
                         corner_pointer%corners_real(2)%vector((2-1)*ncells(2)+ii)
                    bound_latr(ll+4,ii,1)= &
                         corner_pointer%corners_real(2)%vector((2-1)*ncells(2)+ii)
                    ll=ll+4
                 enddo
              enddo

           else

              do ii=1,ncells(1)
                 ll=0
                 do jj=1,nvertices(1)/nof(1)
                    do kk=1,nof(1)
                       ll=ll+1
                       bound_lonr(ll,ii,1)= &
                            corner_pointer%corners_real(1)%vector((kk-1)*ncells(1)+ii)
                    enddo
                 enddo
              enddo

              do ii=1,ncells(2)
                 ll=0
                 do jj=1,nvertices(2)/nof(2)
                    do kk=1,nof(2)
                       ll=ll+1
                       bound_latr(ll,ii,1)= &
                            corner_pointer%corners_real(2)%vector((kk-1)*ncells(2)+ii)
                    enddo
                 enddo
              enddo

           endif

           if (lenaz .gt. 1 ) then
           if (il_grid_type == PRISM_Unstructlonlat_regvrt .or. &
               il_grid_type == PRISM_Gaussreduced_regvrt   ) then

             do ii=1,ncells(3)
               ll=0
               do jj=1,nvertices(3)/nof(3)
                 do kk=1,nof(3)
                   ll=ll+1
                   bound_vertr(ll,ii,1,1)= &
                   corner_pointer%corners_real(3)%vector((kk-1)*ncells(3)+ii)
                 enddo
               enddo
             enddo

           elseif  (il_grid_type == PRISM_Unstructlonlat_sigmavrt .or. &
                    il_grid_type == PRISM_Gaussreduced_sigmavrt ) then

             stride = ncells(2)*ncells(3)

             ii = 0
             do kl=1,ncells(3)
               do jl=1,ncells(2)
                   ii = ii + 1
                   ll = 0
                   do jj=1,nvertices(3)/nof(3)
                      do kk=1,nof(3)
                         ll=ll+1
                         bound_vertr(ll,jl,1,kl)= &
                         corner_pointer%corners_real(3)%vector((kk-1)*stride+ii)
                      enddo
                   enddo
               enddo
             enddo

           else

             do ii=1,ncells(3)
               ll=0
               do jj=1,nvertices(3)/nof(3)
                 do kk=1,nof(3)
                   ll=ll+1
                   bound_vertr(ll,ii,1,1)= &
                   corner_pointer%corners_real(3)%vector((kk-1)*ncells(3)+ii)
                 enddo
               enddo
             enddo

           endif
           endif

         endif

!
!        corner_shape must match grid's valid shape
!
         if(is_regular) then
           myashape(1,1)=nvertices(1)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(1)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=corner_pointer%corner_shape(1:2,2)

           myvshape(1,1)=nvertices(1)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(1)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=Grids(il_grid_id)%grid_shape(1:2,2)

           call psmile_write_2d_real(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(1) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(2) &
               ,bound_lonr &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )

           myashape(1,1)=nvertices(2)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(2)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=corner_pointer%corner_shape(1:2,2)

           myvshape(1,1)=nvertices(2)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(2)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=Grids(il_grid_id)%grid_shape(1:2,2)

           call psmile_write_2d_real(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(2) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(2) &
               ,bound_latr &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
           if(lenaz.gt.1) then
           if(is_regvrt) then
           myashape(1,1)=nvertices(3)*(corner_pointer%corner_shape(1,3)-1)+1
           myashape(2,1)=nvertices(3)*corner_pointer%corner_shape(2,3)
           myvshape(1,1)=nvertices(3)*(Grids(il_grid_id)%grid_shape(1,3)-1)+1
           myvshape(2,1)=nvertices(3)*Grids(il_grid_id)%grid_shape(2,3)
           call psmile_write_1d_real(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(3) &
               ,bound_vertr &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
           else
           myashape(1,1)=nvertices(3)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(3)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=corner_pointer%corner_shape(1:2,2)
           myashape(1:2,3)=corner_pointer%corner_shape(1:2,3)

           myvshape(1,1)=nvertices(3)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(3)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=Grids(il_grid_id)%grid_shape(1:2,2)
           myvshape(1:2,3)=Grids(il_grid_id)%grid_shape(1:2,3)

           call psmile_write_3d_real(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(3) &
               ,bound_vertr &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
           endif
           endif
         else

           myashape(1,1)=nvertices(1)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(1)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=1

           myvshape(1,1)=nvertices(1)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(1)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=1

           call psmile_write_2d_real(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(1) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(2) &
               ,bound_lonr &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )

           myashape(1,1)=nvertices(2)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(2)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=1

           myvshape(1,1)=nvertices(2)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(2)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=1

           call psmile_write_2d_real(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(2) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(2) &
               ,bound_latr &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
           if (lenaz .gt. 1) then
           if(il_grid_type == PRISM_Unstructlonlat_regvrt .or. &
              il_grid_type == PRISM_Gaussreduced_regvrt ) then
           myashape(1,1)=nvertices(3)*(corner_pointer%corner_shape(1,3)-1)+1
           myashape(2,1)=nvertices(3)*corner_pointer%corner_shape(2,3)
           myvshape(1,1)=nvertices(3)*(Grids(il_grid_id)%grid_shape(1,3)-1)+1
           myvshape(2,1)=nvertices(3)*Grids(il_grid_id)%grid_shape(2,3)
           call psmile_write_1d_real(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(3) &
               ,bound_vertr &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
           
           else if( il_grid_type == PRISM_Unstructlonlat_sigmavrt .or. &
                    il_grid_type == PRISM_Gaussreduced_sigmavrt ) then

           myashape(1,1)=nvertices(3)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(3)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=1
           myashape(1:2,3)=corner_pointer%corner_shape(1:2,2)

           myvshape(1,1)=nvertices(3)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(3)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=1
           myvshape(1:2,3)=corner_pointer%corner_shape(1:2,2)
           call psmile_write_3d_real(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(3) &
               ,bound_vertr &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )

           else
           myashape(1,1)=nvertices(3)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(3)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=1
           myashape(1:2,3)=1

           myvshape(1,1)=nvertices(3)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(3)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=1
           myvshape(1:2,3)=1
           call psmile_write_3d_real(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(3) &
               ,bound_vertr &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
           endif
           endif
         endif

         if(associated(bound_lonr)) deallocate(bound_lonr,stat=ierror)
         if (ierror > 0) then
            ierrp (1) = ierror
            ierrp (2) = 1
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'bound_lonr', &
                            ierrp, 2, __FILE__, __LINE__ )
            return
         endif

         if(associated(bound_latr)) deallocate(bound_latr,stat=ierror)
         if (ierror > 0) then
            ierrp (1) = ierror
            ierrp (2) = 1
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'bound_latr', &
                            ierrp, 2, __FILE__, __LINE__ )
            return
         endif

         if(associated(bound_vertr)) deallocate(bound_vertr,stat=ierror)
         if (ierror > 0) then
            ierrp (1) = ierror
            ierrp (2) = 1
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'bound_vertr', &
                            ierrp, 2, __FILE__, __LINE__ )
            return
         endif
!
       else if(associated(corner_pointer%corners_dble(1)%vector)) then
!
         if(is_regular) then
           allocate(bound_lond(nvertices(1),ncells(1),ncells(2)),stat=ierror)
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_lond', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif
  
           allocate(bound_latd(nvertices(2),ncells(1),ncells(2)),stat=ierror)
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_latd', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif

           if (lenaz .gt. 1 ) then
             if(is_regvrt) then
               allocate(bound_vertd(nvertices(3),ncells(3),1,1) &
                   ,stat=ierror)
              else
                allocate(bound_vertd(nvertices(3),ncells(1),ncells(2) &
                        ,ncells(3)) ,stat=ierror)
              endif
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_latd', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif
           endif
!
         else
!
           allocate(bound_lond(nvertices(1),ncells(1),1),stat=ierror)
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_lond', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif
  
           allocate(bound_latd(nvertices(2),ncells(2),1),stat=ierror)
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_latd', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif

           if (lenaz .gt. 1 ) then
           allocate(bound_vertd(nvertices(3),ncells(3),1,1),stat=ierror)
           if (ierror > 0) then
              ierrp (1) = ierror
              ierrp (2) = 1
              ierror = PRISM_Error_Alloc
              call psmile_error ( ierror, 'bound_latd', &
                              ierrp, 2, __FILE__, __LINE__ )
              return
           endif
           endif
!
         endif

#ifdef DEBUG
     print*,trim(ch_id),' :  psmile_write_meta_byid: ','is_regular ',is_regular
     print*,trim(ch_id),' :  psmile_write_meta_byid: '  &
              ,'ih-il+1,jh-jl+1,nvertices(1),ncells(1),nvertices(2),ncells(2)' &
              ,ih-il+1,jh-jl+1,nvertices(1),ncells(1),nvertices(2),ncells(2)
     print*,trim(ch_id),' :  psmile_write_meta_byid: ' &
                       ,'corner_pointer%corner_shape ' &
                       ,corner_pointer%corner_shape
      call psmile_flushstd
#endif

!
!        Branch for all logical rectangular grids
!
         if(is_regular) then

           if( il_grid_type == PRISM_Reglonlatvrt .or.  &
               il_grid_type == PRISM_Reglonlat_sigmavrt ) then

             lbackward=.false.
             do iloc=1,ncells(2)
               do ii=1,ncells(1)
                 ll=0
               do jj=1,nvertices(1)/nof(1)
                 lbackward=mod(jj-1,2).eq.0
                 if(.not.lbackward)then
                 do kk=1,nof(1)
                   ll=ll+1
                   bound_lond(ll,ii,iloc)= &
                    corner_pointer%corners_dble(1)%vector((kk-1)*ncells(1)+ii)
                 enddo
                 else
                 do kk=nof(1),1,-1
                   ll=ll+1
                   bound_lond(ll,ii,iloc)= &
                    corner_pointer%corners_dble(1)%vector((kk-1)*ncells(1)+ii)
                 enddo
                 endif
               enddo
              enddo
             enddo
  
             kkj=0
             do iloc=1,ncells(1)
               do ii=1,ncells(2)
                 ll=0
                 do jj=1,nvertices(2)/nof(2)
                   kkj=kkj+1
                   do kk=1,nof(2)
                     ll=ll+1
                     bound_latd(ll,iloc,ii)= &
                     corner_pointer%corners_dble(2)%vector((kkj-1)*ncells(2)+ii)
                   enddo
                   if(kkj.eq.nof(2))kkj=0
                 enddo
               enddo
             enddo


           else

             stride = ncells(1)*ncells(2)

             ii = 0
             do jl=1,ncells(2)
                do il=1,ncells(1)
                   ii = ii + 1
                   ll=0
                   do jj=1,nvertices(1)/nof(1)
                      do kk=1,nof(1)
                         ll=ll+1
                         bound_lond(ll,il,jl)= &
                         corner_pointer%corners_dble(1)%vector((kk-1)*stride+ii)
                      enddo
                   enddo
                enddo
             enddo

             stride = ncells(1)*ncells(2)

             ii = 0
             do jl=1,ncells(2)
                do il=1,ncells(1)
                   ii = ii + 1
                   ll = 0
                   do jj=1,nvertices(2)/nof(2)
                      do kk=1,nof(2)
                         ll=ll+1
                         bound_latd(ll,il,jl)= &
                         corner_pointer%corners_dble(2)%vector((kk-1)*stride+ii)
                      enddo
                   enddo
                enddo
             enddo

           endif ! il_grid_type == PRISM_Reglonlatvrt ||  &
                 ! il_grid_type == PRISM_Reglonlat_sigmavrt 

!
!  Extraction of height data for grids with logical rectangular 
!  horizontal layers.
!
           if (lenaz .gt. 1 ) then
           if( il_grid_type == PRISM_Reglonlatvrt .or.  &
               il_grid_type == PRISM_Irrlonlat_regvrt ) then

             do ii=1,ncells(3)
                   ll=0
                   do jj=1,nvertices(3)/nof(3)
                     do kk=1,nof(3)
                       ll=ll+1
                       bound_vertd(ll,ii,1,1)= &
                       corner_pointer%corners_dble(3) &
                       %vector((kk-1)*ncells(3)+ii)
                     enddo
                   enddo
             enddo

           else

             stride = ncells(1)*ncells(2)*ncells(3)

             ii = 0
             do kl=1,ncells(3)
               do jl=1,ncells(2)
                 do il=1,ncells(1)
                   ii = ii + 1
                   ll = 0
                   do jj=1,nvertices(3)/nof(3)
                      do kk=1,nof(3)
                         ll=ll+1
                         bound_vertd(ll,il,jl,kl)= &
                         corner_pointer%corners_dble(3)%vector((kk-1)*stride+ii)
                      enddo
                   enddo
                 enddo
               enddo
             enddo

           endif
           endif

!
!       Branch for all grids which are not (logical) rectangular.
!
         else

           if (il_grid_type == PRISM_Gaussreduced_regvrt ) then

              do ii=1,ncells(1)
                 ll=0
                 do jj=1,nof(1)
                    bound_lond(ll+1,ii,1) = &
                         corner_pointer%corners_dble(1)%vector((1-1)*ncells(1)+ii)
                    bound_lond(ll+2,ii,1) = &
                         corner_pointer%corners_dble(1)%vector((2-1)*ncells(1)+ii)
                    bound_lond(ll+3,ii,1) = &
                         corner_pointer%corners_dble(1)%vector((2-1)*ncells(1)+ii)
                    bound_lond(ll+4,ii,1) = &
                         corner_pointer%corners_dble(1)%vector((1-1)*ncells(1)+ii)
                    ll=ll+4
                 enddo
              enddo

              do ii=1,ncells(2)
                 ll=0
                 do jj=1,nof(2)
                    bound_latd(ll+1,ii,1)= &
                         corner_pointer%corners_dble(2)%vector((1-1)*ncells(2)+ii)
                    bound_latd(ll+2,ii,1)= &
                         corner_pointer%corners_dble(2)%vector((1-1)*ncells(2)+ii)
                    bound_latd(ll+3,ii,1)= &
                         corner_pointer%corners_dble(2)%vector((2-1)*ncells(2)+ii)
                    bound_latd(ll+4,ii,1)= &
                         corner_pointer%corners_dble(2)%vector((2-1)*ncells(2)+ii)
                    ll=ll+4
                 enddo
              enddo

           else

              do ii=1,ncells(1)
                 ll=0
                 do jj=1,nvertices(1)/nof(1)
                    do kk=1,nof(1)
                       ll=ll+1
                       bound_lond(ll,ii,1)= &
                            corner_pointer%corners_dble(1)%vector((kk-1)*ncells(1)+ii)
                    enddo
                 enddo
              enddo

              do ii=1,ncells(2)
                 ll=0
                 do jj=1,nvertices(2)/nof(2)
                    do kk=1,nof(2)
                       ll=ll+1
                       bound_latd(ll,ii,1)= &
                            corner_pointer%corners_dble(2)%vector((kk-1)*ncells(2)+ii)
                    enddo
                 enddo
              enddo

           endif

           if (lenaz .gt. 1 ) then
           if (il_grid_type == PRISM_Unstructlonlat_regvrt .or. &
               il_grid_type == PRISM_Gaussreduced_regvrt ) then

             do ii=1,ncells(3)
               ll=0
               do jj=1,nvertices(3)/nof(3)
                 do kk=1,nof(3)
                   ll=ll+1
                   bound_vertd(ll,ii,1,1)= &
                   corner_pointer%corners_dble(3)%vector((kk-1)*ncells(3)+ii)
                 enddo
               enddo
             enddo

           elseif  (il_grid_type == PRISM_Unstructlonlat_sigmavrt .or. &
                    il_grid_type == PRISM_Gaussreduced_sigmavrt ) then

             stride = ncells(2)*ncells(3)

             ii = 0
             do kl=1,ncells(3)
               do jl=1,ncells(2)
                   ii = ii + 1
                   ll = 0
                   do jj=1,nvertices(3)/nof(3)
                      do kk=1,nof(3)
                         ll=ll+1
                         bound_vertd(ll,jl,1,kl)= &
                         corner_pointer%corners_dble(3)%vector((kk-1)*stride+ii)
                      enddo
                   enddo
               enddo
             enddo

           else

             do ii=1,ncells(3)
               ll=0
               do jj=1,nvertices(3)/nof(3)
                 do kk=1,nof(3)
                   ll=ll+1
                   bound_vertd(ll,ii,1,1)= &
                   corner_pointer%corners_dble(3)%vector((kk-1)*ncells(3)+ii)
                 enddo
               enddo
             enddo

           endif
           endif

         endif

!
!        corner_shape must match grid's valid shape
!
         if(is_regular) then
           myashape(1,1)=nvertices(1)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(1)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=corner_pointer%corner_shape(1:2,2)

           myvshape(1,1)=nvertices(1)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(1)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=Grids(il_grid_id)%grid_shape(1:2,2)

           call psmile_write_2d_dble(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(1) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(2) &
               ,bound_lond &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )

           myashape(1,1)=nvertices(2)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(2)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=corner_pointer%corner_shape(1:2,2)

           myvshape(1,1)=nvertices(2)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(2)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=Grids(il_grid_id)%grid_shape(1:2,2)

           call psmile_write_2d_dble(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(2) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(2) &
               ,bound_latd &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )

           if(lenaz.gt.1) then
           if(is_regvrt) then
           myashape(1,1)=nvertices(3)*(corner_pointer%corner_shape(1,3)-1)+1
           myashape(2,1)=nvertices(3)*corner_pointer%corner_shape(2,3)
           myvshape(1,1)=nvertices(3)*(Grids(il_grid_id)%grid_shape(1,3)-1)+1
           myvshape(2,1)=nvertices(3)*Grids(il_grid_id)%grid_shape(2,3)
           
           call psmile_write_1d_dble(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(3) &
               ,bound_vertd &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
           else
           myashape(1,1)=nvertices(3)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(3)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=corner_pointer%corner_shape(1:2,2)
           myashape(1:2,3)=corner_pointer%corner_shape(1:2,3)

           myvshape(1,1)=nvertices(3)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(3)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=Grids(il_grid_id)%grid_shape(1:2,2)
           myvshape(1:2,3)=Grids(il_grid_id)%grid_shape(1:2,3)

           call psmile_write_3d_dble(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(3) &
               ,bound_vertd &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
           endif
           endif

         else

           myashape(1,1)=nvertices(1)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(1)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=1

           myvshape(1,1)=nvertices(1)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(1)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=1

           call psmile_write_2d_dble(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(1) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(2) &
               ,bound_lond &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )

           myashape(1,1)=nvertices(2)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(2)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=1

           myvshape(1,1)=nvertices(2)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(2)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=1

           call psmile_write_2d_dble(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(2) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(2) &
               ,bound_latd &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
           if (lenaz .gt. 1) then
           if(il_grid_type == PRISM_Unstructlonlat_regvrt .or. &
              il_grid_type == PRISM_Gaussreduced_regvrt ) then
           myashape(1,1)=nvertices(3)*(corner_pointer%corner_shape(1,3)-1)+1
           myashape(2,1)=nvertices(3)*corner_pointer%corner_shape(2,3)
           myvshape(1,1)=nvertices(3)*(Grids(il_grid_id)%grid_shape(1,3)-1)+1
           myvshape(2,1)=nvertices(3)*Grids(il_grid_id)%grid_shape(2,3)
           call psmile_write_1d_dble(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(3) &
               ,bound_vertd &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
           
           else if( il_grid_type == PRISM_Unstructlonlat_sigmavrt .or. &
                    il_grid_type == PRISM_Gaussreduced_sigmavrt) then
           myashape(1,1)=nvertices(3)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(3)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=1
           myashape(1:2,3)=corner_pointer%corner_shape(1:2,2)

           myvshape(1,1)=nvertices(3)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(3)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=1
           myvshape(1:2,3)=corner_pointer%corner_shape(1:2,2)
           call psmile_write_3d_dble(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(3) &
               ,bound_vertd &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )

           else
           myashape(1,1)=nvertices(3)*(corner_pointer%corner_shape(1,1)-1)+1
           myashape(2,1)=nvertices(3)*corner_pointer%corner_shape(2,1)
           myashape(1:2,2)=1
           myashape(1:2,3)=1

           myvshape(1,1)=nvertices(3)*(Grids(il_grid_id)%grid_shape(1,1)-1)+1
           myvshape(2,1)=nvertices(3)*Grids(il_grid_id)%grid_shape(2,1)
           myvshape(1:2,2)=1
           myvshape(1:2,3)=1
           call psmile_write_3d_dble(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%crnij(3) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(3) &
               ,bound_vertd &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
           endif
           endif
         endif

         if(associated(bound_lond)) deallocate(bound_lond,stat=ierror)
         if (ierror > 0) then
            ierrp (1) = ierror
            ierrp (2) = 1
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'bound_lond', &
                            ierrp, 2, __FILE__, __LINE__ )
            return
         endif

         if(associated(bound_latd)) deallocate(bound_latd,stat=ierror)
         if (ierror > 0) then
            ierrp (1) = ierror
            ierrp (2) = 1
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'bound_latd', &
                            ierrp, 2, __FILE__, __LINE__ )
            return
         endif

         if(associated(bound_vertd))deallocate(bound_vertd,stat=ierror)
         if (ierror > 0) then
            ierrp (1) = ierror
            ierrp (2) = 1
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'bound_latd', &
                            ierrp, 2, __FILE__, __LINE__ )
            return
         endif

       endif


!
!       Loop over components of vector fields
!
      do ii=1,nvcomp
!
          if(nvcomp.gt.1) then
            vcomp_id=Methods(Fields(id_varid)%method_id)%vector_pointer &
                     %array_of_point_ids(ii)
          else
            vcomp_id=Fields(id_varid)%method_id
          endif
          if(il_grid_type == PRISM_Unstructlonlatvrt) then

            myashape(1:2,1)=Methods(vcomp_id)%coords_pointer%actual_shape(1:2,1)
            myashape(1:2,2)=1
            myashape(1:2,3)=1
!
            myvshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
            myvshape(1:2,2)=1
            myvshape(1:2,3)=1
#ifdef VERBOSE
    print*,trim(ch_id),' : psmile_write_meta_byid: PRISM_Unstructlonlatvrt: ' &
            ,'nvcomp, vcomp_id: myashape, myvshape:' &
            ,nvcomp,vcomp_id,':',myashape(1:2,1:3), myvshape(1:2,1:3)
    call psmile_flushstd
       
#endif
          else if((il_grid_type .eq. PRISM_Unstructlonlat_regvrt) .or. &
                  (il_grid_type .eq. PRISM_Gaussreduced_regvrt) .or. &  
                  ( il_grid_type .eq. PRISM_Gaussreduced_sigmavrt) .or. &  
                 (il_grid_type .eq. PRISM_Unstructlonlat_sigmavrt)) then
!           Lon and lat are a 1d array of length of number of pairs lat-lon
            myashape(1:2,1)=Methods(vcomp_id)%coords_pointer%actual_shape(1:2,1)
            myashape(1:2,2)=1
!           Vert is a 1d array.
            myashape(1:2,3)=Methods(vcomp_id)%coords_pointer%actual_shape(1:2,3)
!
            myvshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
            myvshape(1:2,2)=1
            myvshape(1:2,3)=Grids(il_grid_id)%grid_shape(1:2,3)
#ifdef VERBOSE
    print*,trim(ch_id),' : psmile_write_meta_byid:' &
                      ,' PRISM_Unstructlonlat_vrt|PRISM_Gaussreduced_vrt: ' &
            ,'nvcomp, vcomp_id: myashape, myvshape:' &
            ,nvcomp,vcomp_id,':',myashape(1:2,1:3), myvshape(1:2,1:3)
    call psmile_flushstd
       
#endif

          else if(il_grid_type .eq. PRISM_Irrlonlatvrt) then
!           Lon Lats Verts are arrays in I,j and K, means the grid is logical
!           rectangular in 3 dimensions with respect to the indices:
!           y=f(i,j,k) where y is lat,lon or vert
!
            myashape(1:2,1)=Methods(vcomp_id)%coords_pointer%actual_shape(1:2,1)
            myashape(1:2,2)=Methods(vcomp_id)%coords_pointer%actual_shape(1:2,2)
            myashape(1:2,3)=Methods(vcomp_id)%coords_pointer%actual_shape(1:2,3)
!
            myvshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
            myvshape(1:2,2)=Grids(il_grid_id)%grid_shape(1:2,2)
            myvshape(1:2,3)=Grids(il_grid_id)%grid_shape(1:2,3)
#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid:  PRISM_Irrlonlatvrt: ' &
            ,'nvcomp, vcomp_id: myashape, myvshape:' &
            ,nvcomp,vcomp_id,':',myashape(1:2,1:3), myvshape(1:2,1:3)
      call psmile_flushstd
       
#endif

          else if((il_grid_type .eq. PRISM_Irrlonlat_sigmavrt) .or. &
                  (il_grid_type .eq. PRISM_Irrlonlat_regvrt)) then

!           Lon and Lats are matrices in I and J
            myashape(1:2,1)=Methods(vcomp_id)%coords_pointer%coords_shape(1:2,1)
            myashape(1:2,2)=Methods(vcomp_id)%coords_pointer%coords_shape(1:2,2)
!           Vertical is a simple 1d array
            myashape(1:2,3)=Methods(vcomp_id)%coords_pointer%coords_shape(1:2,3)
!
            myvshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
            myvshape(1:2,2)=Grids(il_grid_id)%grid_shape(1:2,2)
            myvshape(1:2,3)=Grids(il_grid_id)%grid_shape(1:2,3)
#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid:  PRISM_Irrlonlat_vrt: ' &
            ,'nvcomp, vcomp_id: myashape, myvshape:' &
            ,nvcomp,vcomp_id,':',myashape(1:2,1:3), myvshape(1:2,1:3)
      call psmile_flushstd
       
#endif

          else if(il_grid_type .eq. PRISM_Reglonlatvrt) then
            myashape(1:2,1)=Methods(vcomp_id)%coords_pointer%actual_shape(1:2,1)
            myashape(1:2,2)=Methods(vcomp_id)%coords_pointer%actual_shape(1:2,2)
            myashape(1:2,3)=Methods(vcomp_id)%coords_pointer%actual_shape(1:2,3)
!
            myvshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
            myvshape(1:2,2)=Grids(il_grid_id)%grid_shape(1:2,2)
            myvshape(1:2,3)=Grids(il_grid_id)%grid_shape(1:2,3)
#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid:  PRISM_Reglonlatvrt: ' &
            ,'nvcomp, vcomp_id: myashape, myvshape:' &
            ,nvcomp,vcomp_id,':',myashape(1:2,1:3), myvshape(1:2,1:3)
      call psmile_flushstd
       
#endif
          endif
!
!    Longitudes
!
      if(no_of_blocks.eq.1) then
!
        if(associated( &
                  Methods(vcomp_id)%coords_pointer%coords_dble(1)%vector))then 
            if(il_grid_type .eq. PRISM_Reglonlatvrt) then
              allocate(lond((myashape(2,1)-myashape(1,1)+1)* &
                        (myashape(2,2)-myashape(1,2)+1)))
              iloc=1
              do jc=myashape(1,2),myashape(2,2)
!rv                do ic=myashape(1,1),myashape(2,1)
                do ic=1,myashape(2,1)-myashape(1,1)+1
                  lond(iloc)=Methods(vcomp_id) &
                             %coords_pointer%coords_dble(1)%vector(ic)
                  iloc=iloc+1
                enddo
              enddo
            else
              lond => Methods(vcomp_id) &
                           %coords_pointer%coords_dble(1)%vector
            endif

#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid:  Writing lond: ' &
            ,'nvcomp, vcomp_id: myashape, myvshape:' &
            ,nvcomp,vcomp_id,':',myashape(1:2,1:3), myvshape(1:2,1:3)
      call psmile_flushstd
       
#endif
            call psmile_write_2d_dble(Fields(id_varid)%io_infos%file_unit &
                                ,Fields(id_varid)%io_infos%p_mpp_io%lonij(ii) &
                                ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               ,lond &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
            if(il_grid_type .eq. PRISM_Reglonlatvrt) deallocate(lond)
        else
            if(il_grid_type .eq. PRISM_Reglonlatvrt) then
              allocate(lonr((myashape(2,1)-myashape(1,1)+1)* &
                        (myashape(2,2)-myashape(1,2)+1)))
              iloc=1
              do jc=myashape(1,2),myashape(2,2)
!rv                do ic=myashape(1,1),myashape(2,1)
                do ic=1,myashape(2,1)-myashape(1,1)+1
                  lonr(iloc)=Methods(vcomp_id) &
                             %coords_pointer%coords_real(1)%vector(ic)
                  iloc=iloc+1
                enddo
              enddo
            else
              lonr => Methods(vcomp_id) &
                           %coords_pointer%coords_real(1)%vector
            endif
#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid:  Writing lonr: ' &
            ,'nvcomp, vcomp_id: myashape, myvshape:' &
            ,nvcomp,vcomp_id,':',myashape(1:2,1:3), myvshape(1:2,1:3)
      call psmile_flushstd
       
#endif
            call psmile_write_2d_real(Fields(id_varid)%io_infos%file_unit &
                                 ,Fields(id_varid)%io_infos%p_mpp_io%lonij(ii) &
                                ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               ,lonr &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid:  Writing lonr done. ' 
      call psmile_flushstd
       
#endif
            if(il_grid_type .eq. PRISM_Reglonlatvrt) deallocate(lonr)
        endif
!
      else if(no_of_blocks.gt.1) then
!
        if(associated( &
             Methods(vcomp_id)%coords_pointer%coords_dble(1)%vector))then 

            allocate(londb((myashape(2,1)-myashape(1,1)+1) &
            *(myashape(2,2)-myashape(1,2)+1)*(myashape(2,4)-myashape(1,4)+1)))

            if(il_grid_type .eq. PRISM_Reglonlatvrt) then
              iloc=1
              do kc=myashape(1,4),myashape(2,4)
                do jc=myashape(1,2),myashape(2,2)
                  do ic=myashape(1,1),myashape(2,1)
                    londb(iloc)=fdbles(kc,ii,1)%fdble(ic)
                    iloc=iloc+1
                  enddo
                enddo
              enddo
            else
              londb = (/(fdbles(kc,ii,1)%fdble,kc=1,no_of_blocks)/)
            endif
            

            call psmile_write_3d_dble(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%lonij(ii) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               , londb &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )

            deallocate(londb)

        else
            allocate(lonrb((myashape(2,1)-myashape(1,1)+1) &
            *(myashape(2,2)-myashape(1,2)+1)*(myashape(2,4)-myashape(1,4)+1)))
            if(il_grid_type .eq. PRISM_Reglonlatvrt) then
              iloc=1
              do kc=myashape(1,4),myashape(2,4)
                do jc=myashape(1,2),myashape(2,2)
                  do ic=myashape(1,1),myashape(2,1)
                    lonrb(iloc)=freals(kc,ii,1)%freal(ic)
                    iloc=iloc+1
                  enddo
                enddo
              enddo
            else
              lonrb = (/(freals(kc,ii,1)%freal,kc=1,no_of_blocks)/)
            endif
            call psmile_write_3d_real(Fields(id_varid)%io_infos%file_unit &
                                 ,Fields(id_varid)%io_infos%p_mpp_io%lonij(ii) &
                                 ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               , lonrb &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
            deallocate(lonrb)
        endif
!
      endif
!
!    Latitudes
!
      if(no_of_blocks.eq.1) then
!
        if(associated( &
           Methods(vcomp_id)%coords_pointer%coords_dble(2)%vector))then 

            if(il_grid_type .eq. PRISM_Reglonlatvrt) then
              allocate(latd((myashape(2,1)-myashape(1,1)+1)* &
                        (myashape(2,2)-myashape(1,2)+1)))
              iloc=1
!rv              do jc=myashape(1,2),myashape(2,2)
              do jc=1,myashape(2,2)-myashape(1,2)+1
                do ic=myashape(1,1),myashape(2,1)
                  latd(iloc)=Methods(vcomp_id) &
                             %coords_pointer%coords_dble(2)%vector(jc)
                  iloc=iloc+1
                enddo
              enddo
            else
              latd => Methods(vcomp_id) &
                           %coords_pointer%coords_dble(2)%vector
            endif

                call psmile_write_2d_dble(Fields(id_varid)%io_infos%file_unit &
                                 ,Fields(id_varid)%io_infos%p_mpp_io%latij(ii) &
                                 ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               , latd &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )

            if(il_grid_type .eq. PRISM_Reglonlatvrt) deallocate(latd)

        else

            if(il_grid_type .eq. PRISM_Reglonlatvrt) then
              allocate(latr((myashape(2,1)-myashape(1,1)+1)* &
                        (myashape(2,2)-myashape(1,2)+1)))
              iloc=1
!rv              do jc=myashape(1,2),myashape(2,2)
              do jc=1,myashape(2,2)-myashape(1,2)+1
                do ic=myashape(1,1),myashape(2,1)
                  latr(iloc)=Methods(vcomp_id) &
                             %coords_pointer%coords_real(2)%vector(jc)
                  iloc=iloc+1
                enddo
              enddo
            else
              latr => Methods(vcomp_id) &
                           %coords_pointer%coords_real(2)%vector
            endif
            call psmile_write_2d_real(Fields(id_varid)%io_infos%file_unit &
                             ,Fields(id_varid)%io_infos%p_mpp_io%latij(ii) &
                             ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               ,latr &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )

            if(il_grid_type .eq. PRISM_Reglonlatvrt) deallocate(latr)

        endif
!
      else if(no_of_blocks.gt.1) then
!
        if(associated( &
           Methods(vcomp_id)%coords_pointer%coords_dble(2)%vector))then 

            allocate(latdb((myashape(2,1)-myashape(1,1)+1) &
            *(myashape(2,2)-myashape(1,2)+1)*(myashape(2,4)-myashape(1,4)+1)))
            if(il_grid_type .eq. PRISM_Reglonlatvrt) then
              iloc=1
              do kc=myashape(1,4),myashape(2,4)
                do jc=myashape(1,2),myashape(2,2)
                  do ic=myashape(1,1),myashape(2,1)
                    latdb(iloc)=fdbles(kc,ii,2)%fdble(ic)
                    iloc=iloc+1
                  enddo
                enddo
              enddo
            else
              latdb = (/(fdbles(kc,ii,2)%fdble,kc=1,no_of_blocks)/)
            endif

                call psmile_write_3d_dble(Fields(id_varid)%io_infos%file_unit &
                                 ,Fields(id_varid)%io_infos%p_mpp_io%latij(ii) &
                                 ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               ,(/(fdbles(il_i,ii,2)%fdble,il_i=1,no_of_blocks)/) &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
            deallocate(latdb)

        else
            allocate(latrb((myashape(2,1)-myashape(1,1)+1) &
            *(myashape(2,2)-myashape(1,2)+1)*(myashape(2,4)-myashape(1,4)+1)))
            if(il_grid_type .eq. PRISM_Reglonlatvrt) then
              iloc=1
              do kc=myashape(1,4),myashape(2,4)
                do jc=myashape(1,2),myashape(2,2)
                  do ic=myashape(1,1),myashape(2,1)
                    latrb(iloc)=freals(kc,ii,2)%freal(ic)
                    iloc=iloc+1
                  enddo
                enddo
              enddo
            else
              latrb = (/(freals(kc,ii,2)%freal,kc=1,no_of_blocks)/)
            endif

            call psmile_write_3d_real(Fields(id_varid)%io_infos%file_unit &
                                 ,Fields(id_varid)%io_infos%p_mpp_io%latij(ii) &
                                 ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               ,latrb &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
            deallocate(latrb)
        endif
      endif
!
!     Vertical
!
      if((il_grid_type .eq.PRISM_Unstructlonlatvrt) .or. &
         (il_grid_type .eq.PRISM_Irrlonlatvrt)) then
        if(no_of_blocks.eq.1) then
!
          if(associated( &
             Methods(vcomp_id)%coords_pointer%coords_dble(3)%vector))then 

                call psmile_write_3d_dble(Fields(id_varid)%io_infos%file_unit &
                                ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
                                ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               ,Methods(vcomp_id)%coords_pointer%coords_dble(3)%vector &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
          else
                call psmile_write_3d_real(Fields(id_varid)%io_infos%file_unit &
                           ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
                           ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               ,Methods(vcomp_id)%coords_pointer%coords_real(3)%vector &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
          endif
!
        else if(no_of_blocks.gt.1) then
!
          if(associated( &
             Methods(vcomp_id)%coords_pointer%coords_dble(3)%vector))then 

                call psmile_write_4d_dble(Fields(id_varid)%io_infos%file_unit &
                           ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
                           ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               ,(/(fdbles(il_i,ii,3)%fdble,il_i=1,no_of_blocks)/) &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
          else
                call psmile_write_4d_real(Fields(id_varid)%io_infos%file_unit &
                           ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
                           ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
               ,(/(freals(il_i,ii,3)%freal,il_i=1,no_of_blocks)/) &
               ,myashape &
               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
          endif
!
        endif
!
      else
!
!          We have a separate axis for z.
!
        if(no_of_blocks.eq.1) then
!
          if(associated( &
             Methods(vcomp_id)%coords_pointer%coords_dble(3)%vector))then 
            call psmile_write_1d_dble(Fields(id_varid)%io_infos%file_unit &
                       ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
                       ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
                       ,Methods(vcomp_id)%coords_pointer%coords_dble(3) &
                                  %vector &
               ,myashape(1,3) &
               ,myvshape(1,3),dl_time, .FALSE.,ilbl,.FALSE.,ierror )

          else
            call psmile_write_1d_real(Fields(id_varid)%io_infos%file_unit &
                       ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
                       ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
                       ,Methods(vcomp_id)%coords_pointer%coords_real(3) &
                                  %vector &
               ,myashape(1,3) &
               ,myvshape(1,3),dl_time, .FALSE.,ilbl,.FALSE.,ierror )

          endif
!
        else if(no_of_blocks.gt.1) then
!
          if(associated( &
             Methods(vcomp_id)%coords_pointer%coords_dble(3)%vector))then 
            allocate(vertdb(myashape(1,3):myashape(2,3) &
                           ,myashape(1,4):myashape(2,4)))

            do kc=myashape(1,4),myashape(2,4)
              vertdb(myashape(1,3):myashape(2,3),kc)=fdbles(kc,ii,3)%fdble
            enddo
            call mpp_write(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
            ,vertdb(myvshape(1,3):myvshape(2,3),myvshape(1,4):myvshape(2,4)))
            deallocate(vertdb)
          else
            allocate(vertdb(myashape(1,3):myashape(2,3) &
                           ,myashape(1,4):myashape(2,4)))
            do kc=myashape(1,4),myashape(2,4)
              vertdb(myashape(1,3):myashape(2,3),kc)=freals(kc,ii,3)%freal
            enddo
            call mpp_write(Fields(id_varid)%io_infos%file_unit &
                      ,Fields(id_varid)%io_infos%p_mpp_io%vertij(ii) &
            ,vertdb(myvshape(1,3):myvshape(2,3),myvshape(1,4):myvshape(2,4)))
            deallocate(vertdb)
          endif
!
        endif

      endif
!
!     Masks
!
      if(vmsk_id(ii).ne.PRISM_UNDEFINED) then
        myashape = 1
        myashape(1:2,1:ndim_3d) = Masks(vmsk_id(ii))%mask_shape
        if(no_of_blocks.eq.1) then
          call psmile_write_3d_log(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%maskij(ii) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(1)  &
                               ,Masks(vmsk_id(ii))%mask_array &
                               ,myashape &
                               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )

        else if(no_of_blocks.gt.1) then
          call psmile_write_4d_log(Fields(id_varid)%io_infos%file_unit &
                               ,Fields(id_varid)%io_infos%p_mpp_io%maskij(ii) &
                               ,Fields(id_varid)%io_infos%p_mpp_io%domain(1) &
                               ,(/(flogs(il_i,ii)%flog,il_i=1,no_of_blocks)/) &
                               ,myashape &
                               ,myvshape,dl_time, .FALSE.,ilbl,.FALSE.,ierror )
        endif
      endif
!
!       End of loop over components of vector fields
!
          
      enddo

#ifdef DEBUG
        ! Flush header information for debugging session
        call mpp_flush(Fields(id_varid)%io_infos%file_unit)
#endif           

!-----------------------------------------------------------------------------
!    End of dumping the grid informations for a vector field
!
!
     if(associated(fdbles)) deallocate(fdbles)
     if(associated(freals)) deallocate(freals)
     if(associated(flogs)) deallocate(flogs)
        


#ifdef VERBOSE
      print*,trim(ch_id),' :  psmile_write_meta_byid: end'
      call psmile_flushstd

#endif
#endif  /*   __PSMILE_WITH_IO  */
      end subroutine psmile_write_meta_byid