Changeset 2285 in lmdz_wrf for trunk/tools

Timestamp:

Jan 24, 2019, 3:01:33 PM (6 years ago)

Author:

lfita

Message:

Fixing `compute_cellbnds', right values at the x-axis borders

Location:

trunk/tools

Files:

• diagnostics.py (modified) (1 diff)
• module_ForDiagnostics.f90 (modified) (6 diffs)

trunk/tools/diagnostics.py

@@ -956 +956 @@
 
         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

trunk/tools/module_ForDiagnostics.f90

@@ -1069 +1069 @@
     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 +1116 @@
     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 +1125 @@
     indices(4,2,2,2) = -1
 
-    DO i=1,dx
+    DO i=2,dx-1
       DO j=1,dy
         DO iv=1,4
@@ -1130 +1131 @@
           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 +1144 @@
           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 +1152 @@
           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