WriteMosaicField Subroutine

public subroutine WriteMosaicField(field, inputfile, mosaic, rank, dims, rc)

Arguments

Type IntentOptional Attributes Name
type(ESMF_Field), intent(in) :: field
character(len=*), intent(in) :: inputfile
type(ESMF_Mosaic), intent(in) :: mosaic
integer, intent(in) :: rank
integer, intent(in) :: dims(:)
integer, intent(out), optional :: rc

Source Code

subroutine WriteMosaicField(field, inputfile, mosaic, rank, dims, rc)

   type(ESMF_Field), intent(in)   :: field
   character(*), intent(in)       :: inputfile
   type(ESMF_Mosaic), intent(in)  :: mosaic
   integer, intent(in)            :: rank
   integer, intent(in)            :: dims(:)
   integer, optional, intent(out) :: rc

   ! -- local variables
    integer :: localrc
    integer :: localDe, localDeCount
    integer :: de, deCount, dimCount, tile, tileCount
    integer, dimension(:), allocatable :: deToTileMap, localDeToDeMap
    integer, dimension(:,:), allocatable :: minIndexPDe, maxIndexPDe
    integer, dimension(:,:), allocatable :: minIndexPTile, maxIndexPTile
    type(ESMF_Grid)     :: grid
    type(ESMF_DistGrid) :: distgrid
    type(ESMF_Array) :: array
    type(ESMF_VM) :: vm
    integer       :: PetNo, PetCnt, localroot
    type(ESMF_StaggerLoc)         :: staggerloc
    character(len=ESMF_MAXPATHLEN):: fileName
    character(len=MAXNAMELEN):: fieldName
    real(ESMF_KIND_R8), pointer   :: fptr2d(:,:), fptr3d(:,:,:), fptr4d(:,:,:,:)
    real(ESMF_KIND_R8), allocatable   :: buff2d(:,:), buff3d(:,:,:), buff4d(:,:,:,:)
    real(ESMF_KIND_R8), allocatable   :: sndrcvbuffer(:)
    integer :: sndcount, rcvcount, xsize, ysize
    integer, allocatable :: tilenos(:), denos(:), tilearray(:)
    integer :: lncid, varId, ncStatus
    integer :: i, j, i1, i2, i3, i4, tiles(2)

#ifdef ESMF_NETCDF
    rc = ESMF_FAILURE

     call ESMF_VMGetCurrent(vm, rc=localrc)
     if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
     call ESMF_VMGet(vm, petCount = PetCnt, localPet = PetNo, rc=localrc)
     if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
     call ESMF_FieldGet(field, grid=grid, &
        staggerloc=staggerloc, localDeCount=localDeCount, rc=localrc)

       ! -- get domain decomposition
      call ESMF_GridGet(grid, staggerloc, distgrid=distgrid, rc=localrc)
      if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out

      call ESMF_DistGridGet(distgrid, deCount=deCount, dimCount=dimCount, &
        tileCount=tileCount, rc=localrc)
      if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out

      allocate(minIndexPDe(dimCount, deCount), maxIndexPDe(dimCount, deCount),  &
        minIndexPTile(dimCount, tileCount), maxIndexPTile(dimCount, tileCount), &
        deToTileMap(deCount), localDeToDeMap(localDeCount), stat=localrc)
      if (ESMF_LogFoundAllocError(statusToCheck=localrc, ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out

      call ESMF_DistGridGet(distgrid, &
        minIndexPDe=minIndexPDe, maxIndexPDe=maxIndexPDe, &
        minIndexPTile=minIndexPTile, maxIndexPTile=maxIndexPTile, &
        rc=localrc)
      if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out

      call ESMF_FieldGet(field, name = fieldName, array=array, rc=localrc)
      if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
           ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out

      call ESMF_ArrayGet(array, deToTileMap=deToTileMap, &
        localDeToDeMap=localDeToDeMap, rc=localrc)
      if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out

      ! only write out to file if iamroot is true
      if (deCount >= 12) then 
      ! More than one PETs per tile, need to gather data to one PET to write
      ! out to file, in this case, each PET's localDeCount is always <= 1
         if (localDeCount > 1) then
            call ESMF_LogSetError(rcToCheck=ESMF_FAILURE, &
                 msg="- localDeCount should not be more than 1 ", &
                 ESMF_CONTEXT, rcToReturn=rc)
         endif
         if (localDeCount > 0) then
            de   = localDeToDeMap(1) + 1
            tile = deToTileMap(de)
         else
            de = -1
            tile = -1
         endif
         allocate(tilearray(PetCnt*2))
         tiles(1) = de
         tiles(2) = tile
         call ESMF_VMAllGather(vm, tiles, tilearray, 2, rc=localrc)
         if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
              ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
         allocate(tilenos(PetCnt), denos(PetCnt))
         do i=1,PetCnt
	    tilenos(i)=tilearray(i*2)
            denos(i) = tilearray(i*2-1)
         enddo
         deallocate(tilearray)
         ! Find the root of my local tile
	 localDe = 0
         if (tile > 0) then
            ! get the fortran pointer from my own field first 
            if (rank==2) then
               call ESMF_FieldGet(field, localDe=localDe, farrayPtr=fptr2d, &
                    rc=localrc)
               if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
                    ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
               sndcount = size(fptr2d,1)*size(fptr2d,2)
               allocate(sndrcvbuffer(sndcount))
               sndrcvbuffer=reshape(fptr2d, (/sndcount/))
            elseif (rank==3) then
               call ESMF_FieldGet(field, localDe=localDe, farrayPtr=fptr3d, &
                    rc=localrc)
               if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
                    ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
               sndcount = size(fptr3d,1)*size(fptr3d,2)*dims(3)
               allocate(sndrcvbuffer(sndcount))
               sndrcvbuffer=reshape(fptr3d, (/sndcount/))
            elseif (rank==4) then
               call ESMF_FieldGet(field, localDe=localDe, farrayPtr=fptr4d, &
                    rc=localrc)
               if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
                    ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
               sndcount = size(fptr4d,1)*size(fptr4d,2)*dims(3)*dims(4)
               allocate(sndrcvbuffer(sndcount))
               sndrcvbuffer=reshape(fptr4d, (/sndcount/)) 
            endif
            ! Find the first PET that owns the same tile
	    do i=1, PetCnt
               if (tilenos(i)==tile) then
                  localroot = i-1
                  exit
               endif
            enddo
	    if (PetNo == localroot) then
              ! gather all the data from the PETs that own the tile before writing it out
              if (rank==2) then
                 allocate(buff2d(dims(1),dims(2)))
              elseif (rank==3) then
                 allocate(buff3d(dims(1),dims(2),dims(3)))
              elseif (rank==4) then
                 allocate(buff4d(dims(1),dims(2),dims(3),dims(4)))
              endif
              do i=PetNo+1, PetCnt
                 if (tilenos(i)==tile) then
                    !allocate array to receive the data
                    !Receive the data from PetNo
                    xsize = maxIndexPDe(1,denos(i))-minIndexPDe(1,denos(i))+1
                    ysize = maxIndexPDe(2,denos(i))-minIndexPDe(2,denos(i))+1
                    if (i /= PetNo+1) then 
                       if (rank==2) then
                          rcvcount = xsize*ysize
                       elseif (rank==3) then
                          rcvcount = xsize*ysize*dims(3)
                       elseif (rank==4) then
                          rcvcount = xsize*ysize*dims(3)*dims(4)
                       endif
                       ! Get the data from other PETs that own the same tile
                       allocate(sndrcvbuffer(rcvcount))
                       call ESMF_VMRecv(vm, sndrcvbuffer, rcvcount, i-1, rc=localrc)
                       if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
                            ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
                    endif
                    ! Put the block of data into the big array
                    if (rank==2) then
                       j=1
                       do i2=minIndexPDe(2,denos(i)),maxIndexPDe(2,denos(i))
                          do i1=minIndexPDe(1,denos(i)),maxIndexPDe(1,denos(i))
                             buff2d(i1,i2)=sndrcvbuffer(j)
                             j=j+1
                          enddo
                       enddo
                    elseif (rank==3) then
                       j=1
                       do i3 = 1, dims(3)
                          do i2=minIndexPDe(2,denos(i)),maxIndexPDe(2,denos(i))
                             do i1=minIndexPDe(1,denos(i)),maxIndexPDe(1,denos(i))
                                buff3d(i1,i2,i3)=sndrcvbuffer(j)
                                j=j+1
                             enddo
                          enddo
                       enddo
                    elseif (rank==4) then
                       j=1
                       do i4 = 1, dims(4)
                          do i3 = 1, dims(3)
                             do i2=minIndexPDe(2,denos(i)),maxIndexPDe(2,denos(i))
                                do i1=minIndexPDe(1,denos(i)),maxIndexPDe(1,denos(i))
                                   buff4d(i1,i2,i3,i4)=sndrcvbuffer(j)
                                   j=j+1
                                enddo
                             enddo
                          enddo
                       enddo
                    endif
		    deallocate(sndrcvbuffer)
                 endif
              enddo
              ! write it out
              filename = trim(mosaic%tileDirectory)//trim(inputfile)//"."//trim(mosaic%tilenames(tile))//".nc"
              ncStatus = nf90_open(path=trim(fileName), mode=NF90_WRITE, ncid=lncid)
              if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                   msg="Error opening file "//trim(fileName), &
                   ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out

              ncStatus = nf90_inq_varid(lncid, trim(fieldName), varId)
              if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                   msg="Error inquiring variable "//trim(fieldName)//" in "//trim(fileName), &
                   ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out

              if (rank==2) then 
                 ncStatus = nf90_put_var(lncid, varId, buff2d)
                 if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                      msg="Error reading "//trim(fieldName)//" in "//trim(fileName), &
                      ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
                 deallocate(buff2d)
              elseif (rank==3) then
                 ncStatus = nf90_put_var(lncid, varId, buff3d)
                 if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                      msg="Error reading "//trim(fieldName)//" in "//trim(fileName), &
                      ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
                 deallocate(buff3d)
              elseif (rank==4) then
                 ncStatus = nf90_put_var(lncid, varId, buff4d)
                 if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                      msg="Error reading "//trim(fieldName)//" in "//trim(fileName), &
                      ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
                 deallocate(buff4d)
              endif
              ncStatus = nf90_close(lncid)
              if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                   msg="Error closing NetCDF data set", &
                   ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
            else 
               ! If I am not a rootpet of a tile, just send my data to the rootpet
               ! get the fortran pointer from the field 
               call ESMF_VMSend(vm, sndrcvbuffer, sndcount, localroot, rc=localrc)
               if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
                    ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
               deallocate(sndrcvbuffer)
            endif
         endif
         deallocate(tilenos, denos)
      else !deCount < 12
         do localDe = 0, localDeCount-1
            de   = localDeToDeMap(localDe+1) + 1
            tile = deToTileMap(de)
            
            filename = trim(mosaic%tileDirectory)//trim(inputfile)//"."//trim(mosaic%tilenames(tile))//".nc"
            ncStatus = nf90_open(path=trim(fileName), mode=NF90_WRITE, ncid=lncid)
            if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                 msg="Error opening file "//trim(fileName), &
                 ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out

            ncStatus = nf90_inq_varid(lncid, trim(fieldName), varId)
            if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                 msg="Error inquiring variable "//trim(fieldName)//" in "//trim(fileName), &
                 ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out

            if (rank==2) then
               call ESMF_FieldGet(field, localDe=localDe, farrayPtr=fptr2d, &
                    rc=localrc)
               if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
                    ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
               ncStatus = nf90_put_var(lncid, varId, fptr2d)
               if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                    msg="Error reading "//trim(fieldName)//" in "//trim(fileName), &
                    ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
            elseif (rank==3) then
               call ESMF_FieldGet(field, localDe=localDe, farrayPtr=fptr3d, &
                    rc=localrc)
               if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
                    ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
               ncStatus = nf90_put_var(lncid, varId, fptr3d)
               if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                    msg="Error reading "//trim(fieldName)//" in "//trim(fileName), &
                    ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
            elseif (rank==4) then
               call ESMF_FieldGet(field, localDe=localDe, farrayPtr=fptr4d, &
                    rc=localrc)
               if (ESMF_LogFoundError(rcToCheck=localrc, ESMF_ERR_PASSTHRU, &
                    ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
               ncStatus = nf90_put_var(lncid, varId, fptr4d)
               if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                    msg="Error reading "//trim(fieldName)//" in "//trim(fileName), &
                    ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
            endif
            ncStatus = nf90_close(lncid)
            if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, &
                 msg="Error closing NetCDF data set", &
                 ESMF_CONTEXT, rcToReturn=rc)) return  ! bail out
         enddo
      endif
      
      rc=ESMF_SUCCESS
      return
#else
     call ESMF_LogSetError(ESMF_RC_LIB_NOT_PRESENT, &
                 msg="- ESMF_NETCDF not defined when lib was compiled", &
                 ESMF_CONTEXT, rcToReturn=rc)
     return
#endif
   end subroutine WriteMosaicField