Index: /trunk/tools/diagnostics.py
===================================================================
--- /trunk/tools/diagnostics.py	(revision 2284)
+++ /trunk/tools/diagnostics.py	(revision 2285)
@@ -956,10 +956,12 @@
 
         ncvar.insert_variable(ncobj, 'lon_bnds', cellbndsx, diagoutd, diagoutvd, newnc)
+        newnc.sync()
         ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc)
+        newnc.sync()
 
         ovar = newnc.variables['XLONG']
-        ovar.setncattr('bounds', 'lon_bnds, lat_bnds')
+        ovar.setncattr('bounds', 'lon_bnds lat_bnds')
         ovar = newnc.variables['XLAT']
-        ovar.setncattr('bounds', 'lon_bnds, lat_bnds')
+        ovar.setncattr('bounds', 'lon_bnds lat_bnds')
 
 # mrso: total soil moisture SMOIS, DZS
Index: /trunk/tools/module_ForDiagnostics.f90
===================================================================
--- /trunk/tools/module_ForDiagnostics.f90	(revision 2284)
+++ /trunk/tools/module_ForDiagnostics.f90	(revision 2285)
@@ -1069,4 +1069,5 @@
     INTEGER                                              :: i,j,iv
     INTEGER                                              :: ix,ex,iy,ey
+    REAL(r_k)                                            :: tmpval1, tmpval2
     CHARACTER(len=2), DIMENSION(4)                       :: Svertex
     INTEGER, DIMENSION(4,2,2,2)                          :: indices
@@ -1115,5 +1116,5 @@
     indices(3,2,2,2) = 1
     ! SE
-    indices(4,1,1,1) = 1
+    indices(4,1,1,1) = 1 
     indices(4,1,1,2) = 1
     indices(4,1,2,1) = -1
@@ -1124,5 +1125,5 @@
     indices(4,2,2,2) = -1
 
-    DO i=1,dx
+    DO i=2,dx-1
       DO j=1,dy
         DO iv=1,4
@@ -1130,5 +1131,5 @@
           ix = MAX(i+indices(iv,1,1,1),1)
           !ex = MIN(i+indices(iv,1,1,2),dx)
-          ex = i+indices(iv,1,1,2)
+          ex = MAX(i+indices(iv,1,1,2),1)
           iy = MAX(j+indices(iv,1,2,1),1)
           ey = MIN(j+indices(iv,1,2,2),dy)
@@ -1143,5 +1144,5 @@
           iy = MAX(j+indices(iv,2,2,1),1)
           !ey = MIN(i+indices(iv,2,2,2),dy)
-          ey = j+indices(iv,2,2,2)
+          ey = MAX(j+indices(iv,2,2,2),1)
           paral(1,1) = vlon(ix,iy)
           paral(1,2) = vlat(ix,iy)
@@ -1151,10 +1152,106 @@
           CALL intersection_2Dlines(merid, paral, intsct, ptintsct)
           IF (.NOT.intsct) THEN
-            msg = 'not interection found for ' // Svertex(iv) // ' vertex'
+            msg = 'not intersection found for ' // Svertex(iv) // ' vertex'
             CALL ErrMsg(msg, fname, -1)
           END IF
           xbnds(i,j,iv) = ptintsct(1)
           ybnds(i,j,iv) = ptintsct(2)
-        END DO
+
+        END DO
+      END DO
+    END DO
+
+    ! Dealing with the boundary values
+    i = 1
+    DO j=1,dy
+      DO iv=1,4
+
+        ix = MAX(i+indices(iv,1,1,1),1)
+        !ex = MIN(i+indices(iv,1,1,2),dx)
+        ex = MAX(i+indices(iv,1,1,2),1)
+        iy = MAX(j+indices(iv,1,2,1),1)
+        ey = MIN(j+indices(iv,1,2,2),dy)
+        merid(1,1) = ulon(ix,iy)
+        merid(1,2) = ulat(ix,iy)
+        merid(2,1) = ulon(ex,ey)
+        merid(2,2) = ulat(ex,ey) 
+
+        ix = MAX(i+indices(iv,2,1,1),1)
+        ex = MIN(i+indices(iv,2,1,2),dx)
+        iy = MAX(j+indices(iv,2,2,1),1)
+        !ey = MIN(i+indices(iv,2,2,2),dy)
+        ey = MAX(j+indices(iv,2,2,2),1)
+        IF (iv == 1 .OR. iv == 2) THEN
+          ! Projecting values using dx from next grid point
+          tmpval1 = vlon(2,iy)
+          paral(2,1) = vlon(ex,ey)
+          tmpval2 = tmpval1 - paral(2,1)
+          paral(1,1) = paral(2,1) - tmpval2
+          tmpval1 = vlat(2,iy)
+          paral(2,2) = vlat(ex,ey)
+          tmpval2 = tmpval1 - paral(2,2)
+          paral(1,2) = paral(2,2) - tmpval2
+        ELSE
+          paral(1,1) = vlon(ix,iy)
+          paral(1,2) = vlat(ix,iy)
+          paral(2,1) = vlon(ex,ey)
+          paral(2,2) = vlat(ex,ey)
+        END IF
+
+        CALL intersection_2Dlines(merid, paral, intsct, ptintsct)
+        IF (.NOT.intsct) THEN
+          msg = 'not intersection found for ' // Svertex(iv) // ' vertex'
+          CALL ErrMsg(msg, fname, -1)
+        END IF
+        xbnds(i,j,iv) = ptintsct(1)
+        ybnds(i,j,iv) = ptintsct(2)
+
+      END DO
+    END DO
+
+    i = dx
+    DO j=1,dy
+      DO iv=1,4
+
+        ix = MAX(i+indices(iv,1,1,1),1)
+        !ex = MIN(i+indices(iv,1,1,2),dx)
+        ex = MAX(i+indices(iv,1,1,2),1)
+        iy = MAX(j+indices(iv,1,2,1),1)
+        ey = MIN(j+indices(iv,1,2,2),dy)
+        merid(1,1) = ulon(ix,iy)
+        merid(1,2) = ulat(ix,iy)
+        merid(2,1) = ulon(ex,ey)
+        merid(2,2) = ulat(ex,ey) 
+
+        ix = MAX(i+indices(iv,2,1,1),1)
+        ex = MIN(i+indices(iv,2,1,2),dx)
+        iy = MAX(j+indices(iv,2,2,1),1)
+        !ey = MIN(i+indices(iv,2,2,2),dy)
+        ey = MAX(j+indices(iv,2,2,2),1)
+        IF (iv == 3 .OR. iv == 4) THEN
+          ! Projecting values using dx from previous grid point
+          tmpval1 = vlon(dx-1,iy)
+          paral(2,1) = vlon(ex,ey)
+          tmpval2 = tmpval1 - paral(2,1)
+          paral(1,1) = paral(2,1) - tmpval2
+          tmpval1 = vlat(dx-1,iy)
+          paral(2,2) = vlat(ex,ey)
+          tmpval2 = tmpval1 - paral(2,2)
+          paral(1,2) = paral(2,2) - tmpval2
+        ELSE
+          paral(1,1) = vlon(ix,iy)
+          paral(1,2) = vlat(ix,iy)
+          paral(2,1) = vlon(ex,ey)
+          paral(2,2) = vlat(ex,ey)
+        END IF
+
+        CALL intersection_2Dlines(merid, paral, intsct, ptintsct)
+        IF (.NOT.intsct) THEN
+          msg = 'not intersection found for ' // Svertex(iv) // ' vertex'
+          CALL ErrMsg(msg, fname, -1)
+        END IF
+        xbnds(i,j,iv) = ptintsct(1)
+        ybnds(i,j,iv) = ptintsct(2)
+
       END DO
     END DO
