diff --git a/sorc/ncep_post.fd/AVIATION.f b/sorc/ncep_post.fd/AVIATION.f
index 7c8219ec0..4de228f92 100644
--- a/sorc/ncep_post.fd/AVIATION.f
+++ b/sorc/ncep_post.fd/AVIATION.f
@@ -5,6 +5,7 @@
 !! SUBPROGRAM:    CALLLWS       COMPUTES Low Level Wind Shear (0-2000feet) 
 !!   PRGRMMR: Binbin Zhou      /NCEP/EMC  DATE: 2005-08-16       
 !!   19-10-30  Bo CUI          - REMOVE "GOTO" STATEMENT
+!!   21-04-01  Jesse Meng      - computation on defined points only
 !!     
 !! ABSTRACT:  
 !!    This program computes the low level wind shear(LLWS) over 0-2000 feet (0-609.5m)
@@ -83,7 +84,7 @@ SUBROUTINE CALLLWS(U,V,H,LLWS)
 !
       USE vrbls2d, only: fis, u10, v10
       use params_mod, only: gi
-      use ctlblk_mod, only: jsta, jend, im, jm, lsm
+      use ctlblk_mod, only: jsta, jend, im, jm, lsm, spval
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       implicit none
 !
@@ -137,9 +138,10 @@ SUBROUTINE CALLLWS(U,V,H,LLWS)
             END IF
 
 !computer vector difference
-610       LLWS(I,J)=SQRT((U2-U10(I,J))**2+(V2-V10(I,J))**2)/     &
+         LLWS(I,J) = spval
+         if(U10(I,J)<spval.and.V10(I,J)<spval)                   &
+          LLWS(I,J)=SQRT((U2-U10(I,J))**2+(V2-V10(I,J))**2)/     &
                     609.6 * 1.943*609.6                         !unit: knot/2000ft
-
         ENDDO
  
 100   CONTINUE     
@@ -194,7 +196,7 @@ SUBROUTINE CALICING (T1,RH,OMGA, ICING)
 !     MACHINE : BLUE AT NCEP
 !$$$  
 !
-      use ctlblk_mod, only: jsta, jend, im
+      use ctlblk_mod, only: jsta, jend, im, spval
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       implicit none
 !     
@@ -208,7 +210,7 @@ SUBROUTINE CALICING (T1,RH,OMGA, ICING)
 !
       DO J=JSTA,JEND
         DO I=1,IM
-
+        IF(OMGA(I,J)<SPVAL.AND.T1(I,J)<SPVAL.AND.RH(I,J)<SPVAL) THEN
          IF(OMGA(I,J) < 0.0 .AND.                       &
             (T1(I,J) <= 273.0 .AND. T1(I,J) >= 251.0)   &
               .AND. RH(I,J) >= 70.0) THEN
@@ -217,6 +219,9 @@ SUBROUTINE CALICING (T1,RH,OMGA, ICING)
          ELSE
            ICING(I,J) = 0.0
          END IF
+        ELSE
+           ICING(I,J) = SPVAL
+        ENDIF
         ENDDO
       ENDDO
 
@@ -335,6 +340,8 @@ SUBROUTINE CALCAT(U,V,H,U_OLD,V_OLD,H_OLD,CAT)
         DO I=ISTART,ISTOP
 !
           IF(GRIDTYPE=='B')THEN
+           IF(U(I,J)<spval.and.U(I,J-1)<spval.and.U(I-1,J)<spval.and.U(I-1,J-1)<spval.and.&
+              V(I,J)<spval.and.V(I,J-1)<spval.and.V(I-1,J)<spval.and.V(I-1,J-1)<spval)THEN
 !dsh=dv/dx+du/dy 
            DSH=(0.5*(V(I,J)+V(I,J-1))-0.5*(V(I-1,J)+V(I-1,J-1)))*10000./DX(I,J) &
 	      +(0.5*(U(I,J)+U(I-1,J))-0.5*(U(I,J-1)+U(I-1,J-1)))*10000./DY(I,J)
@@ -346,7 +353,13 @@ SUBROUTINE CALCAT(U,V,H,U_OLD,V_OLD,H_OLD,CAT)
 !cvg=-(du/dx+dv/dy)
            CVG = -((0.5*(U(I,J)+U(I,J-1))-0.5*(U(I-1,J)+U(I-1,J-1)))*10000./DX(I,J) &
                 +(0.5*(V(I,J)+V(I-1,J))-0.5*(V(I,J-1)+V(I-1,J-1)))*10000./DY(I,J))	   
+           ELSE
+            DEF = SPVAL
+            CVG = SPVAL
+           ENDIF
           ELSE
+           IF(U(I,J+1)<spval.and.U(I,J-1)<spval.and.U(I+IHE(J),J)<spval.and.U(I+IHW(J),J)<spval.and.&
+              V(I,J+1)<spval.and.V(I,J-1)<spval.and.V(I+IHE(J),J)<spval.and.V(I+IHW(J),J)<spval)THEN
 !dsh=dv/dx+du/dy           
            DSH = (V(I+IHE(J),J) - V(I+IHW(J),J))*10000./(2*DX(I,J))   &  
               + (U(I,J+1) - U(I,J-1))*10000./(2*DY(I,J))
@@ -360,42 +373,106 @@ SUBROUTINE CALCAT(U,V,H,U_OLD,V_OLD,H_OLD,CAT)
 !cvg=-(du/dx+dv/dy)
            CVG = -( (U(I+IHE(J),J) - U(I+IHW(J),J))*10000./(2*DX(I,J)) &
                   +(V(I,J+1) - V(I,J-1))*10000./(2*DY(I,J)) )
+           ELSE
+            DEF = SPVAL
+            CVG = SPVAL
+           ENDIF
           END IF
 	  
           IF(GRIDTYPE == 'A')THEN
 !vws=d|U|/dz
+           IF(U_OLD(I,J)<spval.and.U(I,J)<spval.and.&
+              V_OLD(I,J)<spval.and.V(I,J)<spval.and.&
+              H_OLD(I,J)<spval.and.H(I,J)<spval)THEN
            VWS = ( SQRT(U_OLD(I,J)**2+V_OLD(I,J)**2 ) -               &
                   SQRT(U(I,J)**2+V(I,J)**2 )   ) *                    &
                   1000.0/(H_OLD(I,J) - H(I,J))
+           ELSE
+            VWS = SPVAL
+           ENDIF
           else IF(GRIDTYPE == 'E')THEN
 !vws=d|U|/dz
+           IF(U_OLD(I+IHE(J),J)<spval.and.U(I+IHE(J),J)<spval.and.&
+              V_OLD(I+IHE(J),J)<spval.and.V(I+IHE(J),J)<spval)THEN
+            
 	   VWS1 = ( SQRT(U_OLD(I+IHE(J),J)**2+V_OLD(I+IHE(J),J)**2 ) -&
                   SQRT(U(I+IHE(J),J)**2+V(I+IHE(J),J)**2 )   ) 
+           ELSE
+            VWS1 = SPVAL
+           ENDIF
 !vws=d|U|/dz
+           IF(U_OLD(I+IHW(J),J)<spval.and.U(I+IHW(J),J)<spval.and.&
+              V_OLD(I+IHW(J),J)<spval.and.V(I+IHW(J),J)<spval)THEN
            VWS2 = ( SQRT(U_OLD(I+IHW(J),J)**2+V_OLD(I+IHW(J),J)**2 ) -&   
                   SQRT(U(I+IHW(J),J)**2+V(I+IHW(J),J)**2 )   ) 
+           ELSE
+            VWS2 = SPVAL
+           ENDIF
 !vws=d|U|/dz
+           IF(U_OLD(I,J-1)<spval.and.U(I,J-1)<spval.and.&
+              V_OLD(I,J-1)<spval.and.V(I,J-1)<spval)THEN
            VWS3 = ( SQRT(U_OLD(I,J-1)**2+V_OLD(I,J-1)**2 ) -          & 
                   SQRT(U(I,J-1)**2+V(I,J-1)**2 )   ) 
+           ELSE
+            VWS3 = SPVAL
+           ENDIF
 !vws=d|U|/dz
+           IF(U_OLD(I,J+1)<spval.and.U(I,J+1)<spval.and.&
+              V_OLD(I,J+1)<spval.and.V(I,J+1)<spval)THEN
            VWS4 = ( SQRT(U_OLD(I,J+1)**2+V_OLD(I,J+1)**2 ) -          & 
                   SQRT(U(I,J+1)**2+V(I,J+1)**2 )   ) 
+           ELSE
+            VWS4 = SPVAL
+           ENDIF
+
+           IF(VWS1<spval.and.VWS2<spval.and.VWS3<spval.and.VWS4<spval.and.&
+              H_OLD(I,J)<spval.and.H(I,J)<spval)THEN
            VWS=1000.0*(VWS1+VWS2+VWS3+VWS4)/4.0/(H_OLD(I,J) - H(I,J))
+           ELSE
+            VWS = SPVAL
+           ENDIF
 	  ELSE IF(GRIDTYPE == 'B')THEN
+           IF(U_OLD(I+IHE(J),J)<spval.and.U(I+IHE(J),J)<spval.and.&
+              V_OLD(I+IHE(J),J)<spval.and.V(I+IHE(J),J)<spval)THEN
 	   VWS1 = ( SQRT(U_OLD(I+IHE(J),J)**2+V_OLD(I+IHE(J),J)**2 ) -&
                   SQRT(U(I+IHE(J),J)**2+V(I+IHE(J),J)**2 )   ) 
+           ELSE
+            VWS1 = SPVAL
+           ENDIF
 !vws=d|U|/dz
+           IF(U_OLD(I+IHW(J),J)<spval.and.U(I+IHW(J),J)<spval.and.&
+              V_OLD(I+IHW(J),J)<spval.and.V(I+IHW(J),J)<spval)THEN
            VWS2 = ( SQRT(U_OLD(I+IHW(J),J)**2+V_OLD(I+IHW(J),J)**2 ) -&   
                   SQRT(U(I+IHW(J),J)**2+V(I+IHW(J),J)**2 )   ) 
+           ELSE
+            VWS2 = SPVAL
+           ENDIF
 !vws=d|U|/dz
+           IF(U_OLD(I,J-1)<spval.and.U(I,J-1)<spval.and.&
+              V_OLD(I,J-1)<spval.and.V(I,J-1)<spval)THEN
            VWS3 = ( SQRT(U_OLD(I,J-1)**2+V_OLD(I,J-1)**2 ) -          & 
                   SQRT(U(I,J-1)**2+V(I,J-1)**2 )   ) 
+           ELSE
+            VWS3 = SPVAL
+           ENDIF
 !vws=d|U|/dz
+           IF(U_OLD(I-1,J-1)<spval.and.U(I-1,J-1)<spval.and.&
+              V_OLD(I-1,J-1)<spval.and.V(I-1,J-1)<spval)THEN
            VWS4 = ( SQRT(U_OLD(I-1,J-1)**2+V_OLD(I-1,J-1)**2 ) -          & 
                   SQRT(U(I-1,J-1)**2+V(I-1,J-1)**2 )   ) 
+           ELSE
+            VWS4 = SPVAL
+           ENDIF
+
+           IF(VWS1<spval.and.VWS2<spval.and.VWS3<spval.and.VWS4<spval.and.&
+              H_OLD(I,J)<spval.and.H(I,J)<spval)THEN
            VWS=1000.0*(VWS1+VWS2+VWS3+VWS4)/4.0/(H_OLD(I,J) - H(I,J)) 
+           ELSE
+            VWS=SPVAL
+           ENDIF
 	  END IF  
-           
+          
+         IF(VWS<spval.and.DEF<spval.and.CVG<spval)THEN 
           TRBINDX = ABS(VWS)*(DEF + ABS(CVG))
 	  
           IF(TRBINDX<=4.) THEN
@@ -407,7 +484,9 @@ SUBROUTINE CALCAT(U,V,H,U_OLD,V_OLD,H_OLD,CAT)
           ELSE
             CAT(I,J)=3.0
           END IF        
- 
+         ELSE
+          CAT(I,J)=SPVAL
+         ENDIF
         ENDDO
  
 100   CONTINUE     
@@ -526,7 +605,7 @@ SUBROUTINE CALFLTCND (CEILING,FLTCND)
 !$$$  
 !
       use vrbls2d, only: vis
-      use ctlblk_mod, only: jsta, jend, im
+      use ctlblk_mod, only: jsta, jend, im, spval
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       implicit none
 !     
@@ -543,6 +622,7 @@ SUBROUTINE CALFLTCND (CEILING,FLTCND)
       DO J=JSTA,JEND
         DO I=1,IM
  
+        IF(CEILING(I,J)<spval.and.VIS(I,J)<spval)THEN
           CEIL = CEILING(I,J) * 3.2808               !from m -> feet
           VISI = VIS(I,J) / 1609.0                   !from m -> miles       
 
@@ -561,7 +641,9 @@ SUBROUTINE CALFLTCND (CEILING,FLTCND)
              FLTCND(I,J) = 4.0
 
           END IF
-
+        ELSE
+          FLTCND(I,J) = SPVAL
+        ENDIF
         ENDDO
       ENDDO
 
diff --git a/sorc/ncep_post.fd/BNDLYR.f b/sorc/ncep_post.fd/BNDLYR.f
index e06c53c7c..f78c2f08b 100644
--- a/sorc/ncep_post.fd/BNDLYR.f
+++ b/sorc/ncep_post.fd/BNDLYR.f
@@ -31,6 +31,7 @@
 !!   00-01-04  JIM TUCCILLO - MPI VERSION 
 !!   02-01-15  MIKE BALDWIN - WRF VERSION
 !!   20-11-10  JESSE MENG   - USE UPP_PHYSICS MODULE
+!!   21-04-01  JESSE MENG   - COMPUTATION ON DEFINED POINTS ONLY
 !!     
 !!     USAGE:    CALL BNDLYR(PBND,TBND,QBND,RHBND,UBND,VBND,
 !!                            WBND,OMGBND,PWTBND,QCNVBND)
@@ -273,6 +274,8 @@ SUBROUTINE BNDLYR(PBND,TBND,QBND,RHBND,UBND,VBND,       &
         DO J=JSTA,JEND
           DO I=1,IM
             IF(PSUM(I,J,LBND)/=0.)THEN
+            IF(T(I,J,LBND)<spval.and.Q(I,J,LBND)<spval.and.&
+               UH(I,J,LBND)<spval.and.VH(I,J,LBND)<spval) THEN
               RPSUM           = 1./PSUM(I,J,LBND)
               LVLBND(I,J,LBND)= LVLBND(I,J,LBND)/NSUM(I,J,LBND)
               PBND(I,J,LBND)  = (PBINT(I,J,LBND)+PBINT(I,J,LBND+1))*0.5
@@ -286,6 +289,17 @@ SUBROUTINE BNDLYR(PBND,TBND,QBND,RHBND,UBND,VBND,       &
               END IF 
               WBND(I,J,LBND)     = WBND(I,J,LBND)*RPSUM
               QCNVBND(I,J,LBND)  = QCNVBND(I,J,LBND)*RPSUM
+            ELSE
+              LVLBND(I,J,LBND)= spval
+              PBND(I,J,LBND)  = spval
+              TBND(I,J,LBND)  = spval
+              QBND(I,J,LBND)  = spval
+              OMGBND(I,J,LBND)= spval
+              UBND(I,J,LBND)  = spval
+              VBND(I,J,LBND)  = spval
+              WBND(I,J,LBND)  = spval
+              QCNVBND(I,J,LBND)= spval
+            ENDIF
             ENDIF
           ENDDO
         ENDDO
@@ -294,9 +308,14 @@ SUBROUTINE BNDLYR(PBND,TBND,QBND,RHBND,UBND,VBND,       &
           DO J=JSTA_M,JEND_M
             DO I=2,IM-1
               IF(PVSUM(I,J,LBND)/=0.)THEN
+              IF(UBND(I,J,LBND)<spval.and.VBND(I,J,LBND)<spval.and.PVSUM(I,J,LBND)<spval)THEN
                 RPVSUM         = 1./PVSUM(I,J,LBND)
                 UBND(I,J,LBND) = UBND(I,J,LBND)*RPVSUM
                 VBND(I,J,LBND) = VBND(I,J,LBND)*RPVSUM
+              ELSE
+                UBND(I,J,LBND) = spval
+                VBND(I,J,LBND) = spval
+              ENDIF
               ENDIF
             ENDDO
           ENDDO
@@ -336,6 +355,7 @@ SUBROUTINE BNDLYR(PBND,TBND,QBND,RHBND,UBND,VBND,       &
                 VBND(I,J,LBND) = VH(I,J,L)
               END IF 
               WBND(I,J,LBND)    = WH(I,J,L)
+            IF(T(I,J,LBND)<spval.and.Q(I,J,LBND)<spval)THEN
               QCNVBND(I,J,LBND) = QCNVG(I,J,L)
               IF(MODELNAME == 'GFS' .OR. MODELNAME == 'FV3R')THEN
                 ES   = FPVSNEW(T(I,J,L))
@@ -347,10 +367,17 @@ SUBROUTINE BNDLYR(PBND,TBND,QBND,RHBND,UBND,VBND,       &
               QSBND(I,J,LBND)  = QSAT
               OMGBND(I,J,LBND) = OMGA(I,J,L)
               PWTBND(I,J,LBND) = (Q(I,J,L)+CWM(I,J,L))*DP*GI
+            ELSE
+              QCNVBND(I,J,LBND)= spval
+              QSBND(I,J,LBND)  = spval
+              OMGBND(I,J,LBND) = spval
+              PWTBND(I,J,LBND) = spval
+            ENDIF
             ENDIF
 !
 !   RH, BOUNDS CHECK
 !
+          IF(T(I,J,LBND)<spval.and.Q(I,J,LBND)<spval)THEN            
             RHBND(I,J,LBND) = QBND(I,J,LBND)/QSBND(I,J,LBND)
             IF (RHBND(I,J,LBND)>1.0) THEN
               RHBND(I,J,LBND) = 1.0
@@ -360,6 +387,10 @@ SUBROUTINE BNDLYR(PBND,TBND,QBND,RHBND,UBND,VBND,       &
               RHBND(I,J,LBND) = 0.01
               QBND(I,J,LBND)  = RHBND(I,J,LBND)*QSBND(I,J,LBND)
             ENDIF
+          ELSE
+              RHBND(I,J,LBND) = spval
+              QBND(I,J,LBND)  = spval
+          ENDIF
           ENDDO
         ENDDO
 !
diff --git a/sorc/ncep_post.fd/CALHEL.f b/sorc/ncep_post.fd/CALHEL.f
index 84d5fc841..8520bf5cd 100644
--- a/sorc/ncep_post.fd/CALHEL.f
+++ b/sorc/ncep_post.fd/CALHEL.f
@@ -84,7 +84,7 @@ SUBROUTINE CALHEL(DEPTH,UST,VST,HELI,USHR1,VSHR1,USHR6,VSHR6)
       use params_mod, only: g
       use lookup_mod, only: ITB,JTB,ITBQ,JTBQ
       use ctlblk_mod, only: jsta, jend, jsta_m, jend_m, jsta_2l, jend_2u, &
-                            lm, im, jm, me
+                            lm, im, jm, me, spval
       use gridspec_mod, only: gridtype
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       implicit none
@@ -444,6 +444,10 @@ SUBROUTINE CALHEL(DEPTH,UST,VST,HELI,USHR1,VSHR1,USHR6,VSHR6)
                 DU2 = UH(I,J,L)-UH(I,J,L-1)
                 DV1 = VH(I,J,L+1)-VH(I,J,L)
                 DV2 = VH(I,J,L)-VH(I,J,L-1)
+               IF( VH(I,J,L)  <spval.and.UH(I,J,L)  <spval.and.         &
+                   VH(I,J,L+1)<spval.and.UH(I,J,L+1)<spval.and.         &
+                   VH(I,J,L-1)<spval.and.UH(I,J,L-1)<spval.and.         &
+                   VST(I,J)   <spval.and.UST(I,J)   <spval)             &
                 HELI(I,J,N) = ((VH(I,J,L)-VST(I,J))*                      &
                                (DZ2*(DU1/DZ1)+DZ1*(DU2/DZ2))              &
                             -  (UH(I,J,L)-UST(I,J))*                      &
diff --git a/sorc/ncep_post.fd/CALHEL2.f b/sorc/ncep_post.fd/CALHEL2.f
index e1560bd49..c7678b1c9 100644
--- a/sorc/ncep_post.fd/CALHEL2.f
+++ b/sorc/ncep_post.fd/CALHEL2.f
@@ -88,7 +88,7 @@ SUBROUTINE CALHEL2(LLOW,LUPP,DEPTH,UST,VST,HELI,CANGLE)
       use params_mod, only: g
       use lookup_mod, only: ITB,JTB,ITBQ,JTBQ
       use ctlblk_mod, only: jsta, jend, jsta_m, jend_m, jsta_2l, jend_2u, &
-                            lm, im, jm, me
+                            lm, im, jm, me, spval
       use gridspec_mod, only: gridtype
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       implicit none
@@ -461,6 +461,10 @@ SUBROUTINE CALHEL2(LLOW,LUPP,DEPTH,UST,VST,HELI,CANGLE)
                 DV1 = VH(I,J,L+1)-VH(I,J,L)
                 DV2 = VH(I,J,L)-VH(I,J,L-1)
               IF( L >= LUPP(I,J) .AND. L <= LLOW(I,J) ) THEN
+               IF( VH(I,J,L)  <spval.and.UH(I,J,L)  <spval.and.         &
+                   VH(I,J,L+1)<spval.and.UH(I,J,L+1)<spval.and.         &
+                   VH(I,J,L-1)<spval.and.UH(I,J,L-1)<spval.and.         &
+                   VST(I,J)   <spval.and.UST(I,J)   <spval)             &
                 HELI(I,J,N) = ((VH(I,J,L)-VST(I,J))*                      &
                                (DZ2*(DU1/DZ1)+DZ1*(DU2/DZ2))              &
                             -  (UH(I,J,L)-UST(I,J))*                      &
@@ -484,11 +488,16 @@ SUBROUTINE CALHEL2(LLOW,LUPP,DEPTH,UST,VST,HELI,CANGLE)
 
           DO J=JSTART,JSTOP
             DO I=ISTART,ISTOP
+             IF(VSHR05(I,J)<spval.and.USHR05(I,J)<spval.and. &
+                VST(I,J)<spval.and.UST(I,J)<spval) THEN
                CANGLE(I,J)=ATAN2(VSHR05(I,J),USHR05(I,J))-ATAN2(VST(I,J),UST(I,J))
                CANGLE(I,J)=(CANGLE(I,J)/PI)*180.
                IF(CANGLE(I,J) >  180.) CANGLE(I,J)=360.-CANGLE(I,J)
                IF(CANGLE(I,J) < 0. .AND. CANGLE(I,J) >= -180.) CANGLE(I,J)=-CANGLE(I,J)
                IF(CANGLE(I,J) < -180.) CANGLE(I,J)=360.+CANGLE(I,J)
+             ELSE
+               CANGLE(I,J)=spval
+             ENDIF
             ENDDO
           ENDDO
 !
diff --git a/sorc/ncep_post.fd/CALHEL3.f b/sorc/ncep_post.fd/CALHEL3.f
index 1bd11a0bc..976778bfb 100644
--- a/sorc/ncep_post.fd/CALHEL3.f
+++ b/sorc/ncep_post.fd/CALHEL3.f
@@ -86,7 +86,7 @@ SUBROUTINE CALHEL3(LLOW,LUPP,UST,VST,HELI)
       use params_mod, only: g
       use lookup_mod, only: ITB,JTB,ITBQ,JTBQ
       use ctlblk_mod, only: jsta, jend, jsta_m, jend_m, jsta_2l, jend_2u, &
-                            lm, im, jm, me
+                            lm, im, jm, me, spval
       use gridspec_mod, only: gridtype
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       implicit none
@@ -458,6 +458,10 @@ SUBROUTINE CALHEL3(LLOW,LUPP,UST,VST,HELI)
                 DV1 = VH(I,J,L+1)-VH(I,J,L)
                 DV2 = VH(I,J,L)-VH(I,J,L-1)
               IF( L >= LUPP(I,J) .AND. L <= LLOW(I,J) ) THEN
+               IF( VH(I,J,L)  <spval.and.UH(I,J,L)  <spval.and.         &
+                   VH(I,J,L+1)<spval.and.UH(I,J,L+1)<spval.and.         &
+                   VH(I,J,L-1)<spval.and.UH(I,J,L-1)<spval.and.         &
+                   VST(I,J)   <spval.and.UST(I,J)   <spval)             &
                 HELI(I,J) = ((VH(I,J,L)-VST(I,J))*                      &
                                (DZ2*(DU1/DZ1)+DZ1*(DU2/DZ2))              &
                             -  (UH(I,J,L)-UST(I,J))*                      &
diff --git a/sorc/ncep_post.fd/CALMCVG.f b/sorc/ncep_post.fd/CALMCVG.f
index 2e96e10d3..9edfcdc13 100644
--- a/sorc/ncep_post.fd/CALMCVG.f
+++ b/sorc/ncep_post.fd/CALMCVG.f
@@ -21,6 +21,7 @@
 !!   02-04-23  MIKE BALDWIN - WRF C-GRID VERSION     
 !!   05-07-07  BINBIN ZHOU - ADD RSM A GRID
 !!   06-04-25  H CHUANG - BUG FIXES TO CORECTLY COMPUTE MC AT BOUNDARIES 
+!!   21-04-01  J MENG   - COMPUTATION ON DEFINED POINTS ONLY
 !!   
 !! USAGE:    CALL CALMCVG(Q1D,U1D,V1D,QCNVG)
 !!   INPUT ARGUMENT LIST:
@@ -102,6 +103,8 @@ SUBROUTINE CALMCVG(Q1D,U1D,V1D,QCNVG)
              QUDX   = (Q1D(I+1,J)*UWND(I+1,J)-Q1D(I-1,J)*UWND(I-1,J))*R2DX
              QVDY   = (Q1D(I,J+1)*VWND(I,J+1)-Q1D(I,J-1)*VWND(I,J-1))*R2DY
              QCNVG(I,J) = -(QUDX + QVDY)
+           ELSE
+             QCNVG(I,J) = SPVAL
            ENDIF
          ENDDO
        ENDDO
@@ -119,8 +122,13 @@ SUBROUTINE CALMCVG(Q1D,U1D,V1D,QCNVG)
          ISTA = 1+MOD(J+1,2)
          IEND = IM-MOD(J,2)
          DO I=ISTA,IEND
+          IF(Q1D(I,J-1)<SPVAL.AND.Q1D(I+IVW(J),J)<SPVAL.AND.&
+             Q1D(I+IVE(J),J)<SPVAL.AND.Q1D(I,J+1)<SPVAL) THEN
            QV(I,J) = D25*(Q1D(I,J-1)+Q1D(I+IVW(J),J)                   &
                          +Q1D(I+IVE(J),J)+Q1D(I,J+1))
+          ELSE
+           QV(I,J) = SPVAL
+          ENDIF
          END DO
        END DO
 
@@ -132,6 +140,8 @@ SUBROUTINE CALMCVG(Q1D,U1D,V1D,QCNVG)
        DO J=JSTA_M2,JEND_M2
          IEND = IM-1-MOD(J,2)
          DO I=2,IEND
+          IF(QV(I+IHE(J),J)<SPVAL.AND.UWND(I+IHE(J),J)<SPVAL.AND.&
+             QV(I+IHW(J),J)<SPVAL.AND.UWND(I+IHW(J),J)<SPVAL) THEN
            R2DX   = 1./(2.*DX(I,J))
            R2DY   = 1./(2.*DY(I,J))
            QUDX   = (QV(I+IHE(J),J)*UWND(I+IHE(J),J)                   &
@@ -139,6 +149,9 @@ SUBROUTINE CALMCVG(Q1D,U1D,V1D,QCNVG)
            QVDY   = (QV(I,J+1)*VWND(I,J+1)-QV(I,J-1)*VWND(I,J-1))*R2DY
 
            QCNVG(I,J) = -(QUDX + QVDY) * HBM2(I,J)
+          ELSE
+           QCNVG(I,J) = SPVAL
+          ENDIF
          ENDDO
        ENDDO
       ELSE IF(gridtype=='B')THEN
@@ -148,6 +161,12 @@ SUBROUTINE CALMCVG(Q1D,U1D,V1D,QCNVG)
 !$omp  parallel do private(i,j,qudx,qvdy,r2dx,r2dy)
        DO J=JSTA_M,JEND_M
         DO I=2,IM-1
+         IF(UWND(I,J)<SPVAL.AND.UWND(I,J-1)<SPVAL.AND.&
+            UWND(I-1,J)<SPVAL.AND.UWND(I-1,J-1)<SPVAL.AND.&
+            Q1D(I,J)<SPVAL.AND.Q1D(I+1,J)<SPVAL.AND.Q1D(I-1,J)<SPVAL.AND.&
+            VWND(I,J)<SPVAL.AND.VWND(I-1,J)<SPVAL.AND.&
+            VWND(I,J-1)<SPVAL.AND.VWND(I-1,J-1)<SPVAL.AND.&
+            Q1D(I,J+1)<SPVAL.AND.Q1D(I,J-1)<SPVAL) THEN
           R2DX   = 1./DX(I,J)
           R2DY   = 1./DY(I,J)
           QUDX=(0.5*(UWND(I,J)+UWND(I,J-1))*0.5*(Q1D(I,J)+Q1D(I+1,J))        &
@@ -156,6 +175,9 @@ SUBROUTINE CALMCVG(Q1D,U1D,V1D,QCNVG)
                -0.5*(VWND(I,J-1)+VWND(I-1,J-1))*0.5*(Q1D(I,J)+Q1D(I,J-1)))*R2DY
   
           QCNVG(I,J) = -(QUDX + QVDY)
+         ELSE
+          QCNVG(I,J) = SPVAL
+         ENDIF
 !	  print*,'mcvg=',i,j,r2dx,r2dy,QCNVG(I,J)
         ENDDO
        ENDDO
diff --git a/sorc/ncep_post.fd/CALVIS_GSD.f b/sorc/ncep_post.fd/CALVIS_GSD.f
index 95ff6efb6..1ed19a1d5 100644
--- a/sorc/ncep_post.fd/CALVIS_GSD.f
+++ b/sorc/ncep_post.fd/CALVIS_GSD.f
@@ -310,7 +310,7 @@ SUBROUTINE CALVIS_GSD(CZEN,VIS)
        end if
 
 !  Calculation of visibility based on hydrometeor and aerosols.  (RH effect not yet included.)
-        VIS(I,J)=MIN(90.,CONST1/BETAV+extcof55(i,j,lm))      ! max of 90km
+        VIS(I,J)=MIN(90.,CONST1/(BETAV+extcof55(i,j,lm)))      ! max of 90km
 
         if (vis(i,j)<vis_min) vis_min = vis(i,j)
         if (visrh<visrh_min) visrh_min = visrh
diff --git a/sorc/ncep_post.fd/CLDRAD.f b/sorc/ncep_post.fd/CLDRAD.f
index efc577db0..749e4f3bf 100644
--- a/sorc/ncep_post.fd/CLDRAD.f
+++ b/sorc/ncep_post.fd/CLDRAD.f
@@ -72,6 +72,7 @@
 !!   20-11-10  Jesse Meng - USE UPP_PHYSICS MODULE
 !!   21-02-08  Anning Cheng, read aod550, aod550_du/su/ss/oc/bc
 !!             directly from fv3gfs and output to grib2 by setting rdaod
+!!   21-04-01  Jesse Meng - COMPUTATION ON DEFINED POINTS ONLY
 !!     
 !! USAGE:    CALL CLDRAD
 !!   INPUT ARGUMENT LIST:
@@ -475,10 +476,12 @@ SUBROUTINE CLDRAD
 !     
 !     TOTAL COLUMN CLOUD WATER
       IF (IGET(200) > 0 .or. IGET(575) > 0) THEN 
+       GRID1 = spval
+       GRID2 = spval
        IF (MODELNAME == 'RAPR') THEN
           DO J=JSTA,JEND
             DO I=1,IM
-              GRID1(I,J) = LWP(I,J)/1000.0 ! use WRF-diagnosed value
+              IF(LWP(I,J) < SPVAL) GRID1(I,J) = LWP(I,J)/1000.0 ! use WRF-diagnosed value
             ENDDO
           ENDDO
        ELSE
@@ -489,7 +492,11 @@ SUBROUTINE CLDRAD
 !$omp parallel do private(i,j)
           DO J=JSTA,JEND
             DO I=1,IM
+             IF(GRID1(I,J)<SPVAL.and.GRID2(I,J)<SPVAL)THEN 
               GRID1(I,J) = GRID1(I,J) + GRID2(I,J)
+             ELSE
+              GRID1(I,J) = SPVAL
+             ENDIF
             ENDDO
           ENDDO
         END IF ! GFS
@@ -527,10 +534,11 @@ SUBROUTINE CLDRAD
 !
 !     TOTAL COLUMN CLOUD ICE
       IF (IGET(201) > 0) THEN
+       GRID1 = spval
        IF (MODELNAME == 'RAPR') THEN
           DO J=JSTA,JEND
             DO I=1,IM
-              GRID1(I,J) = IWP(I,J)/1000.0 ! use WRF-diagnosed value
+              IF(IWP(I,J) < SPVAL) GRID1(I,J) = IWP(I,J)/1000.0 ! use WRF-diagnosed value
             ENDDO
           ENDDO
        ELSE
@@ -697,7 +705,7 @@ SUBROUTINE CLDRAD
 !$omp  parallel do
          DO J=JSTA,JEND
            DO I=1,IM
-             GRID1(I,J) = GRID1(I,J)*RRNUM
+             IF(GRID1(I,J) < SPVAL) GRID1(I,J) = GRID1(I,J)*RRNUM
            ENDDO
          ENDDO
          ID(1:25)=0
@@ -747,7 +755,7 @@ SUBROUTINE CLDRAD
 !$omp  parallel do
          DO J=JSTA,JEND
          DO I=1,IM
-           GRID1(I,J) = GRID1(I,J)*RRNUM
+           IF(GRID1(I,J) < SPVAL) GRID1(I,J) = GRID1(I,J)*RRNUM
          ENDDO
          ENDDO
          ID(1:25)=0
@@ -1334,10 +1342,14 @@ SUBROUTINE CLDRAD
 !             ENDIF
 !ADDED BRAD'S MODIFICATION
               RSUM = D00
+             IF (NCFRST(I,J)<SPVAL.and.ACFRST(I,J)<SPVAL)THEN
               IF (NCFRST(I,J) > 0) RSUM=ACFRST(I,J)/NCFRST(I,J)
               IF (NCFRCV(I,J) > 0)                               &
                 RSUM=MAX(RSUM, ACFRCV(I,J)/NCFRCV(I,J))
               GRID1(I,J) = RSUM*100.
+             ELSE
+               GRID1(I,J) = spval
+             ENDIF
             ENDDO
           ENDDO
         END IF 
@@ -1389,11 +1401,15 @@ SUBROUTINE CLDRAD
            ELSE 
             DO J=JSTA,JEND
             DO I=1,IM
+              IF (NCFRST(I,J)<SPVAL.and.ACFRST(I,J)<SPVAL)THEN
                IF (NCFRST(I,J)>0.0) THEN
                   GRID1(I,J) = ACFRST(I,J)/NCFRST(I,J)*100.
                ELSE
                   GRID1(I,J) = D00
                ENDIF
+              ELSE
+                  GRID1(I,J) = spval
+              ENDIF     
             ENDDO
             ENDDO
            END IF 
@@ -1437,11 +1453,15 @@ SUBROUTINE CLDRAD
 	   ELSE  
             DO J=JSTA,JEND
             DO I=1,IM
+              IF (NCFRCV(I,J)<SPVAL.and.ACFRCV(I,J)<SPVAL)THEN
                IF (NCFRCV(I,J)>0.0) THEN
                   GRID1(I,J) = ACFRCV(I,J)/NCFRCV(I,J)*100.
                ELSE
                   GRID1(I,J) = D00
                ENDIF
+              ELSE
+                  GRID1(I,J) = spval
+              ENDIF
             ENDDO
             ENDDO
 	   END IF
@@ -3668,14 +3688,17 @@ SUBROUTINE CLDRAD
 !     
 !     CURRENT INCOMING SW RADIATION AT THE SURFACE.
       IF (IGET(156)>0) THEN
+         GRID1=SPVAL
          DO J=JSTA,JEND
          DO I=1,IM
+          IF(RSWIN(I,J)<SPVAL) THEN
            IF(CZMEAN(I,J)>1.E-6) THEN
              FACTRS=CZEN(I,J)/CZMEAN(I,J)
            ELSE
              FACTRS=0.0
            ENDIF
-           GRID1(I,J)=RSWIN(I,J)*FACTRS
+           IF(RSWIN(I,J)<SPVAL) GRID1(I,J)=RSWIN(I,J)*FACTRS
+          ENDIF
          ENDDO
          ENDDO
 !
@@ -3695,6 +3718,7 @@ SUBROUTINE CLDRAD
           IF(MODELNAME=='RSM' .OR. MODELNAME == 'RAPR') THEN      !add by Binbin: RSM has direct RLWIN output
            GRID1(I,J)=RLWIN(I,J)
           ELSE
+           IF(SIGT4(I,J)<SPVAL.and.T(I,J,NINT(LMH(I,J)))<SPVAL) THEN
            IF(SIGT4(I,J)>0.0) THEN
              LLMH=NINT(LMH(I,J))
              TLMH=T(I,J,LLMH)
@@ -3703,6 +3727,7 @@ SUBROUTINE CLDRAD
              FACTRL=0.0
            ENDIF
            IF(RLWIN(I,J) < spval) GRID1(I,J)=RLWIN(I,J)*FACTRL
+           ENDIF
           ENDIF
          ENDDO
          ENDDO
@@ -3716,15 +3741,18 @@ SUBROUTINE CLDRAD
 !     
 !     CURRENT OUTGOING SW RADIATION AT THE SURFACE.
       IF (IGET(141)>0) THEN
+        GRID1 = spval
 !$omp parallel do private(i,j)
         DO J=JSTA,JEND
           DO I=1,IM
+            IF(RSWOUT(I,J)<SPVAL) THEN
              IF(CZMEAN(I,J)>1.E-6) THEN
                FACTRS=CZEN(I,J)/CZMEAN(I,J)
              ELSE
                FACTRS=0.0
              ENDIF
-             GRID1(I,J)=RSWOUT(I,J)*FACTRS
+             IF(RSWOUT(I,J)<SPVAL) GRID1(I,J)=RSWOUT(I,J)*FACTRS
+            ENDIF
            ENDDO
          ENDDO
 !
@@ -3810,15 +3838,18 @@ SUBROUTINE CLDRAD
 
 !     CURRENT (instantaneous) INCOMING CLEARSKY SW RADIATION AT THE SURFACE.
       IF (IGET(262)>0) THEN
+         GRID1 = spval
 !$omp parallel do private(i,j)
          DO J=JSTA,JEND
            DO I=1,IM
+            IF(RSWINC(I,J)<SPVAL) THEN
              IF(CZMEAN(I,J)>1.E-6) THEN
                FACTRS=CZEN(I,J)/CZMEAN(I,J)
              ELSE
                FACTRS=0.0
              ENDIF
-             GRID1(I,J) = RSWINC(I,J)*FACTRS
+             IF(RSWINC(I,J)<SPVAL) GRID1(I,J) = RSWINC(I,J)*FACTRS
+            ENDIF
            ENDDO
          ENDDO
          if(grib=="grib2" )then
@@ -4971,15 +5002,18 @@ SUBROUTINE CLDRAD
 
 ! AER ASYM FACTOR AT 340 NM
           IF ( IGET(649) > 0 )  THEN
+          GRID1 = spval
 !$omp parallel do private(i,j)
           DO J=JSTA,JEND
           DO I=1,IM
+           IF(SCA2D(I,J)<SPVAL.and.ASY2D(I,J)<SPVAL) THEN
             IF ( SCA2D(I,J) > 0.0 ) THEN
              ASY2D(I,J) = ASY2D(I,J) / SCA2D(I,J)
             ELSE
              ASY2D(I,J) = 0.
             ENDIF
-          GRID1(I,J)=ASY2D(I,J)
+            IF(ASY2D(I,J)<SPVAL) GRID1(I,J)=ASY2D(I,J)
+           ENDIF
           ENDDO
           ENDDO
           CALL BOUND(GRID1,D00,H99999)
@@ -4992,15 +5026,18 @@ SUBROUTINE CLDRAD
 
 ! AER SINGLE SCATTER ALB AT 340 NM
           IF ( IGET(648) > 0 )  THEN
+          GRID1 = SPVAL
 !$omp parallel do private(i,j)
           DO J=JSTA,JEND
           DO I=1,IM
+           IF(AOD(I,J)<SPVAL.and.SCA2D(I,J)<SPVAL) THEN
             IF ( AOD(I,J) > 0.0 ) THEN
              SCA2D(I,J) = SCA2D(I,J) / AOD(I,J)
             ELSE
              SCA2D(I,J) = 1.0
             ENDIF
-             GRID1(I,J)=SCA2D(I,J)
+             IF(SCA2D(I,J)<SPVAL) GRID1(I,J)=SCA2D(I,J)
+           ENDIF
           ENDDO
           ENDDO
           CALL BOUND(GRID1,D00,H99999)
@@ -5118,8 +5155,9 @@ SUBROUTINE CLDRAD
 !$omp parallel do private(i,j)
          DO J = JSTA,JEND
             DO I = 1,IM
-               GRID1(I,J) = DUEM(I,J,1)*1.E-6
+               IF(DUEM(I,J,1)<SPVAL) GRID1(I,J) = DUEM(I,J,1)*1.E-6
                DO K=2,NBIN_DU
+               IF(DUEM(I,J,K)<SPVAL)&
                 GRID1(I,J) = GRID1(I,J) + DUEM(I,J,K)*1.E-6
                END DO
             END DO
@@ -5136,8 +5174,9 @@ SUBROUTINE CLDRAD
 !$omp parallel do private(i,j)
          DO J = JSTA,JEND
             DO I = 1,IM
-               GRID1(I,J) = DUSD(I,J,1)*1.E-6
+               IF(DUSD(I,J,1)<SPVAL) GRID1(I,J) = DUSD(I,J,1)*1.E-6
                DO K=2,NBIN_DU
+               IF(DUSD(I,J,K)<SPVAL)&
                 GRID1(I,J) = GRID1(I,J)+ DUSD(I,J,K)*1.E-6
                END DO
             END DO
@@ -5255,11 +5294,12 @@ SUBROUTINE CLDRAD
       ENDIF
 !! ADD TOTAL AEROSOL PM10 COLUMN DENSITY (kg/m2) !
       IF (IGET(621)>0 ) THEN
+         GRID1=SPVAL
 !$omp parallel do private(i,j)
          DO J = JSTA,JEND
             DO I = 1,IM
                !GRID1(I,J) = DUCMASS(I,J) * 1.E-6
-               GRID1(I,J) = DUCMASS(I,J) * 1.E-9
+               IF(DUCMASS(I,J)<SPVAL) GRID1(I,J) = DUCMASS(I,J) * 1.E-9
             END DO
          END DO
          if(grib=='grib2') then
@@ -5271,11 +5311,12 @@ SUBROUTINE CLDRAD
 
 !! ADD TOTAL AEROSOL PM2.5 COLUMN DENSITY (kg/m2)  
       IF (IGET(622)>0 ) THEN
+         GRID1=SPVAL
 !$omp parallel do private(i,j)
          DO J = JSTA,JEND
             DO I = 1,IM
                !GRID1(I,J) = DUCMASS25(I,J) * 1.E-6
-               GRID1(I,J) = DUCMASS25(I,J) * 1.E-9
+               IF(DUCMASS25(I,J)<SPVAL) GRID1(I,J) = DUCMASS25(I,J) * 1.E-9
             END DO
          END DO
          if(grib=='grib2') then
@@ -5287,10 +5328,11 @@ SUBROUTINE CLDRAD
 
 !! ADD DUST PM2.5 COLUMN DENSITY (kg/m2)  
       IF (IGET(646)>0 ) THEN
+         GRID1=SPVAL
 !$omp parallel do private(i,j)
          DO J = JSTA,JEND
             DO I = 1,IM
-               GRID1(I,J) = DUSTCB(I,J) * 1.E-9
+               IF(DUSTCB(I,J)<SPVAL) GRID1(I,J) = DUSTCB(I,J) * 1.E-9
             END DO
          END DO
          if(grib=='grib2') then
@@ -5302,10 +5344,11 @@ SUBROUTINE CLDRAD
 
 !! ADD SEA SALT PM2.5 COLUMN DENSITY (kg/m2)  
       IF (IGET(647)>0 ) THEN
+         GRID1=SPVAL
 !$omp parallel do private(i,j)
          DO J = JSTA,JEND
             DO I = 1,IM
-               GRID1(I,J) = SSCB(I,J) * 1.E-9
+               IF(SSCB(I,J)<SPVAL) GRID1(I,J) = SSCB(I,J) * 1.E-9
             END DO
          END DO
          if(grib=='grib2') then
@@ -5316,10 +5359,11 @@ SUBROUTINE CLDRAD
       ENDIF
 !! ADD BC COLUMN DENSITY (kg/m2)  
       IF (IGET(616)>0 ) THEN
+         GRID1=SPVAL
 !$omp parallel do private(i,j)
          DO J = JSTA,JEND
             DO I = 1,IM
-               GRID1(I,J) = BCCB(I,J) * 1.E-9
+               IF(BCCB(I,J)<SPVAL) GRID1(I,J) = BCCB(I,J) * 1.E-9
             END DO
          END DO
          if(grib=='grib2') then
@@ -5331,10 +5375,11 @@ SUBROUTINE CLDRAD
 
 !! ADD OC COLUMN DENSITY (kg/m2)  !
       IF (IGET(617)>0 ) THEN
+         GRID1=SPVAL
 !$omp parallel do private(i,j)
          DO J = JSTA,JEND
             DO I = 1,IM
-               GRID1(I,J) = OCCB(I,J) * 1.E-9
+               IF(OCCB(I,J)<SPVAL) GRID1(I,J) = OCCB(I,J) * 1.E-9
             END DO
          END DO
          if(grib=='grib2') then
@@ -5346,10 +5391,11 @@ SUBROUTINE CLDRAD
 
 !! ADD SULF COLUMN DENSITY (kg/m2)  !
       IF (IGET(618)>0 ) THEN
+         GRID1=SPVAL
 !$omp parallel do private(i,j)
          DO J = JSTA,JEND
             DO I = 1,IM
-               GRID1(I,J) = SULFCB(I,J) * 1.E-9
+               IF(SULFCB(I,J)<SPVAL) GRID1(I,J) = SULFCB(I,J) * 1.E-9
             END DO
          END DO
          if(grib=='grib2') then
@@ -5616,8 +5662,9 @@ subroutine wrt_aero_diag(igetfld,nbin,data)
 !$omp parallel do private(i,j)
       DO J = JSTA,JEND
         DO I = 1,IM
-          grid1(I,J) = data(I,J,1)
+          if(data(I,J,1)<spval) grid1(I,J) = data(I,J,1)
           DO K=2,NBIN
+           if(data(I,J,K)<spval)&
             GRID1(I,J) = GRID1(I,J)+ data(I,J,K)
           END DO
         END DO
diff --git a/sorc/ncep_post.fd/FIXED.f b/sorc/ncep_post.fd/FIXED.f
index 38b977b69..64c549d7c 100644
--- a/sorc/ncep_post.fd/FIXED.f
+++ b/sorc/ncep_post.fd/FIXED.f
@@ -15,6 +15,7 @@
 !!   02-06-19  MIKE BALDWIN - WRF VERSION
 !!   11-02-06  JUN WANG     - grib2 option
 !!   20-03-25  JESSE MENG   - remove grib1
+!!   21-04-01  JESSE MENG   - computation on defined points only
 !!     
 !! USAGE:    CALL FIXED
 !!   INPUT ARGUMENT LIST:
@@ -176,8 +177,11 @@ SUBROUTINE FIXED
 !           SNOFAC = AMIN1(SNOK*50.0,1.0)
 !           EGRID1(I,J)=ALB(I,J)+(1.-VEGFRC(I,J))*SNOFAC
 !     1                *(SNOALB-ALB(I,J))
-          IF(ABS(ALBEDO(I,J)-SPVAL)>SMALL)                   &
+          IF(ABS(ALBEDO(I,J)-SPVAL)>SMALL) THEN
            GRID1(I,J)=ALBEDO(I,J)
+          ELSE
+           GRID1(I,J) = SPVAL
+          ENDIF
          ENDDO
        ENDDO
 !       CALL E2OUT(150,000,GRID1,GRID2,GRID1,GRID2,IM,JM)
@@ -212,8 +216,11 @@ SUBROUTINE FIXED
 !$omp parallel do private(i,j)
             DO J=JSTA,JEND
               DO I=1,IM
-                IF(ABS(AVGALBEDO(I,J)-SPVAL)>SMALL)           &
+                IF(ABS(AVGALBEDO(I,J)-SPVAL)>SMALL) THEN
                   GRID1(I,J) = AVGALBEDO(I,J)*100.
+                ELSE
+                  GRID1(I,J) = SPVAL
+                ENDIF
               ENDDO
             ENDDO
        
@@ -261,6 +268,8 @@ SUBROUTINE FIXED
      &               (abs(SICE(I,J)-1.) < 1.0E-5) ) THEN
                MXSNAL(I,J)=0.60
              ENDIF
+             ELSE
+               MXSNAL(I,J)=SPVAL
              ENDIF
            ENDDO
          ENDDO
@@ -287,11 +296,13 @@ SUBROUTINE FIXED
 !$omp parallel do private(i,j)
          DO J=JSTA,JEND
            DO I=1,IM
+             GRID1(I,J) = SPVAL
              IF (MODELNAME == 'NMM') THEN
                IF( (abs(SM(I,J)-1.) < 1.0E-5) ) THEN
                  GRID1(I,J) = SST(I,J)
                ELSE
-                 GRID1(I,J) = THS(I,J)*(PINT(I,J,LM+1)/P1000)**CAPA
+                 IF(THS(I,J)<SPVAL.and.PINT(I,J,LM+1)<SPVAL)&
+                  GRID1(I,J) = THS(I,J)*(PINT(I,J,LM+1)/P1000)**CAPA
                END IF  
              ELSE
                GRID1(I,J) = SST(I,J)
diff --git a/sorc/ncep_post.fd/MDL2AGL.f b/sorc/ncep_post.fd/MDL2AGL.f
index f113acad6..8847fc148 100644
--- a/sorc/ncep_post.fd/MDL2AGL.f
+++ b/sorc/ncep_post.fd/MDL2AGL.f
@@ -14,6 +14,7 @@
 !!   19-10-30  B CUI - REMOVE "GOTO" STATEMENT
 !!   20-03-25  J MENG - remove grib1 
 !!   21-03-11  B Cui - change local arrays to dimension (im,jsta:jend)
+!!   21-04-01  J MENG - computation on defined points only
 !!     
 !! USAGE:    CALL MDL2P
 !!   INPUT ARGUMENT LIST:
@@ -1291,10 +1292,15 @@ SUBROUTINE MDL2AGL
           IF((IGET(411)>0) ) THEN
             DO J=JSTA,JEND
             DO I=1,IM
+             IF(QAGL(I,J)<SPVAL.and.PAGL(I,J)<SPVAL.and.TAGL(I,J)<SPVAL.and.&
+                 UAGL(I,J)<SPVAL.and.VAGL(I,J)<SPVAL)THEN
               QAGL(I,J)=QAGL(I,J)/1000.0
               PV=QAGL(I,J)*PAGL(I,J)/(EPS*(1-QAGL(I,J)) + QAGL(I,J))
               RHO=(1/TAGL(I,J))*(((PAGL(I,J)-PV)/RD) + PV/461.495)
               GRID1(I,J)=0.5*RHO*(SQRT(UAGL(I,J)**2+VAGL(I,J)**2))**3
+             ELSE
+              GRID1(I,J)=SPVAL
+             ENDIF
             ENDDO
             ENDDO
             if(grib=="grib2" )then
diff --git a/sorc/ncep_post.fd/MDL2P.f b/sorc/ncep_post.fd/MDL2P.f
index 15c07f87f..8f1d27d79 100644
--- a/sorc/ncep_post.fd/MDL2P.f
+++ b/sorc/ncep_post.fd/MDL2P.f
@@ -29,6 +29,7 @@
 !!   20-05-20  J MENG   - CALRH unification with NAM scheme
 !!   20-11-10  J MENG   - USE UPP_PHYSICS MODULE
 !!   21-03-11  B Cui - change local arrays to dimension (im,jsta:jend)
+!!   21-04-01  J MENG   - COMPUTATION ON DEFINED POINTS ONLY
 !!
 !! USAGE:    CALL MDL2P
 !!   INPUT ARGUMENT LIST:
@@ -332,8 +333,8 @@ SUBROUTINE MDL2P(iostatusD3D)
                  IF(Q(I,J,1) < SPVAL)   QSL(I,J) = Q(I,J,1)
 
                  IF(gridtype == 'A')THEN
-                   USL(I,J) = UH(I,J,1)
-                   VSL(I,J) = VH(I,J,1)
+                   IF(UH(I,J,1) < SPVAL) USL(I,J) = UH(I,J,1)
+                   IF(VH(I,J,1) < SPVAL) VSL(I,J) = VH(I,J,1)
                  END IF
 
 !                if ( J == JSTA.and. I == 1.and.me == 0)    &
@@ -1183,7 +1184,11 @@ SUBROUTINE MDL2P(iostatusD3D)
 !$omp parallel do private(i,j)
             DO J=JSTA,JEND
               DO I=1,IM
+              IF(TSL(I,J) < SPVAL .AND. QSL(I,J) < SPVAL) THEN
                 GRID1(I,J) = TSL(I,J)*(1.+0.608*QSL(I,J))
+              ELSE
+                GRID1(I,J) = SPVAL
+              ENDIF
               ENDDO
             ENDDO
 
@@ -1328,6 +1333,7 @@ SUBROUTINE MDL2P(iostatusD3D)
 !$omp  parallel do private(i,j)
             DO J=JSTA,JEND
               DO I=1,IM
+                GRID1(I,J) = SPVAL
                 CFRSL(I,J) = MIN(MAX(0.0,CFRSL(I,J)),1.0)
                 IF(abs(CFRSL(I,J)-SPVAL) > SMALL)                   &    
                       GRID1(I,J) = CFRSL(I,J)*H100
@@ -1621,7 +1627,9 @@ SUBROUTINE MDL2P(iostatusD3D)
 !$omp  parallel do private(i,j)
             DO J=JSTA,JEND
               DO I=1,IM
+              IF(FSL(I,J)<SPVAL)THEN
                 EGRID2(I,J) = FSL(I,J)*GI
+              ENDIF
               ENDDO
             ENDDO
             CALL CALSTRM(EGRID2(1,jsta),EGRID1(1,jsta))
@@ -1684,8 +1692,12 @@ SUBROUTINE MDL2P(iostatusD3D)
 !$omp  parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(QW1(I,J) < SPVAL .AND. QI1(I,J) < SPVAL) THEN
                    GRID1(I,J) = QW1(I,J) + QI1(I,J)
                    QI1(I,J) = spval
+                 ELSE
+                   GRID1(I,J) = SPVAL
+                 ENDIF
                  ENDDO
                ENDDO
              ELSE
@@ -1940,7 +1952,11 @@ SUBROUTINE MDL2P(iostatusD3D)
 !$omp  parallel do private(i,j)
               DO J=JSTA,JEND
                 DO I=1,IM
+                  IF(FSL(I,J)<spval)THEN
                   GRID1(I,J)  = FSL(I,J)*GI
+                  ELSE
+                  GRID1(I,J)  = SPVAL
+                  ENDIF
                   EGRID1(I,J) = SPVAL
                 ENDDO
               ENDDO
@@ -1977,7 +1993,11 @@ SUBROUTINE MDL2P(iostatusD3D)
           DO I=1,IM
             USL_OLD(I,J) = USL(I,J)
             VSL_OLD(I,J) = VSL(I,J)
+            IF(FSL(I,J)<spval)THEN
             FSL_OLD(I,J) = FSL(I,J)*GI
+            ELSE
+            FSL_OLD(I,J) = SPVAL
+            ENDIF
           ENDDO
         ENDDO
 !
@@ -2013,7 +2033,11 @@ SUBROUTINE MDL2P(iostatusD3D)
 !$omp  parallel do private(i,j)
              DO J=JSTA,JEND
                DO I=1,IM
+               IF(SMOKESL(I,J,1)<SPVAL.and.SPL(LP)<SPVAL.and.TSL(I,J)<SPVAL)THEN
                  GRID1(I,J) = (1./RD)*SMOKESL(I,J,1)*(SPL(LP)/TSL(I,J))
+               ELSE
+                 GRID1(I,J) = SPVAL
+               ENDIF
                ENDDO
              ENDDO
              if(grib == 'grib2')then
@@ -3526,7 +3550,11 @@ SUBROUTINE MDL2P(iostatusD3D)
                    ELSE
                      IMOIS = 3
                    END IF
+                 IF(TSL(I,J)<spval.and.TPRS(I,J,LUHI)<spval.and.TDSL(I,J)<spval)THEN
                    HAINES(I,J) = ISTA + IMOIS
+                 ELSE
+                   HAINES(I,J) = SPVAL
+                 ENDIF
 ! 	       if(i==570 .and. j==574)print*,'high hainesindex:',i,j,luhi,tsl(i,j) &
 ! 	       ,tprs(i,j,luhi),tdsl(i,j),ista,imois,spl(luhi),spl(lp),haines(i,j)
                  END IF 
@@ -3568,7 +3596,11 @@ SUBROUTINE MDL2P(iostatusD3D)
                  END IF
 ! 	       if(i==570 .and. j==574)print*,'mid haines index:',i,j,luhi,tsl(i,j) &
 ! 	       ,tprs(i,j,luhi),tdsl(i,j),ista,imois,spl(luhi),spl(lp),haines(i,j)
+               IF(TSL(I,J)<spval.and.TPRS(I,J,LUHI)<spval.and.TDSL(I,J)<spval)THEN
                  HAINES(I,J) = ISTA + IMOIS
+               ELSE
+                 HAINES(I,J) = SPVAL
+               ENDIF
                END IF 
              END DO
            END DO  
@@ -3608,7 +3640,11 @@ SUBROUTINE MDL2P(iostatusD3D)
                  END IF
 ! 	       if(i==570 .and. j==574)print*,'low haines index:',i,j,luhi,tsl(i,j) &
 ! 	       ,tprs(i,j,luhi),tdsl(i,j),ista,imois,spl(luhi),spl(lp),haines(i,j)
+               IF(TSL(I,J)<spval.and.TPRS(I,J,LUHI)<spval.and.TDSL(I,J)<spval)THEN
                  HAINES(I,J) = ISTA + IMOIS
+               ELSE
+                 HAINES(I,J) = SPVAL
+               ENDIF
                END IF 
              END DO
            END DO  
@@ -3808,7 +3844,11 @@ SUBROUTINE MDL2P(iostatusD3D)
 !$omp  parallel do private(i,j)
                  DO J=JSTA,JEND
                    DO I=1,IM
+                   IF(FSL(I,J)<SPVAL)THEN
                      GRID1(I,J) = FSL(I,J)*GI
+                   ELSE
+                     GRID1(I,J) = SPVAL
+                   ENDIF
                    ENDDO
                  ENDDO    
                ELSE
diff --git a/sorc/ncep_post.fd/MDLFLD.f b/sorc/ncep_post.fd/MDLFLD.f
index ba916b0ba..37e91ff4d 100644
--- a/sorc/ncep_post.fd/MDLFLD.f
+++ b/sorc/ncep_post.fd/MDLFLD.f
@@ -42,6 +42,7 @@
 !!   20-05-20  J MENG - CALRH unification with NAM scheme
 !!   20-11-10  J MENG - USE UPP_MATH MODULE
 !!   20-11-10  J MENG - USE UPP_PHYSICS MODULE
+!!   21-04-01  J MENG - COMPUTATION ON DEFINED POINTS ONLY
 !!
 !! USAGE:    CALL MDLFLD
 !!   INPUT ARGUMENT LIST:
@@ -365,6 +366,7 @@ SUBROUTINE MDLFLD
                 ! - Model output DBZ is present - proceed with calc
                 DO J=JSTA,JEND
                 DO I=1,IM
+                IF(P1D(I,J)<spval.and.T1D(I,J)<spval.and.Q1D(I,J)<spval)THEN
                   ze_nc=10.**(0.1*REF_10CM(I,J,L))
                   DBZ1(I,J)=10.*LOG10(max(Zmin,(ze_nc+CUREFL(I,J))))
                   DBZR1(I,J)=MIN(DBZR1(I,J), REF_10CM(I,J,L))
@@ -398,6 +400,10 @@ SUBROUTINE MDLFLD
                       DBZI1(I,J)=10.*LOG10(ze_s)
                     ENDIF  refl_adj
                   ENDIF    refl_comp
+                ELSE
+                      DBZR1(I,J)=DBZmin
+                      DBZI1(I,J)=DBZmin
+                ENDIF
                 ENDDO
                 ENDDO
               ELSE
@@ -427,8 +433,13 @@ SUBROUTINE MDLFLD
 !
            DO J=JSTA,JEND
            DO I=1,IM
+           IF(C1D(I,J)<spval.and.FI1D(I,J)<spval)THEN
               QI1(I,J)=C1D(I,J)*FI1D(I,J)
               QW1(I,J)=C1D(I,J)-QI1(I,J)
+           ELSE
+              QI1(I,J)=D00
+              QW1(I,J)=D00
+           ENDIF
               QR1(I,J)=D00
               QS1(I,J)=D00
               DBZ1(I,J)=DBZmin
@@ -539,7 +550,6 @@ SUBROUTINE MDLFLD
             DBZ(I,J,L)=MAX(DBZmin, DBZ(I,J,L))
             DBZR(I,J,L)=MAX(DBZmin, DBZR(I,J,L))
             DBZI(I,J,L)=MAX(DBZmin, DBZI(I,J,L))
-    
           ENDIF       !-- End IF (L > LMH(I,J)) ...
          ENDDO         !-- End DO I loop
         ENDDO
@@ -630,6 +640,7 @@ SUBROUTINE MDLFLD
                 DBZC(I,J,L)=CUREFL(I,J)
                ENDIF                   !-- End IF (CUREFL_S(I,J) > 0.)
 
+             IF(T(I,J,L)<spval) THEN
 !              IF(T(I,J,L)  <  1.0E-3) print*,'ZERO T'    
                IF(T(I,J,L)  >  1.0E-3)                            &
      &         DENS = PMID(I,J,L)/(RD*T(I,J,L)*(Q(I,J,L)*D608+1.0))      ! DENSITY
@@ -695,6 +706,12 @@ SUBROUTINE MDLFLD
                  DBZI(I,J,L) = MAX(DBZmin, DBZI(I,J,L))
                  DBZC(I,J,L) = MAX(DBZmin, DBZC(I,J,L))
                END IF 
+             ELSE
+                 DBZ(I,J,L)  = DBZmin
+                 DBZR(I,J,L) = DBZmin
+                 DBZI(I,J,L) = DBZmin
+                 DBZC(I,J,L) = DBZmin
+             ENDIF !(T(I,J,L)<spval)
              ENDDO
            ENDDO
          ENDDO
@@ -726,6 +743,7 @@ SUBROUTINE MDLFLD
 
         DO L=1,LM
           LL=LM-L+1
+          IF(T(I,J,LL)<spval)THEN
             IF(T(I,J,LL) < 1.0E-3)print*,'ZERO T'
             IF(T(I,J,LL) > 1.0E-3)                            &
              RHOD=PMID(I,J,LL)/                                  &
@@ -821,7 +839,11 @@ SUBROUTINE MDLFLD
              DBZ(i,j,ll) = ze_sum
              DBZR(i,j,ll) = ze_r*1.E18
              DBZI(i,j,ll) = (ze_s+ze_g)*1.E18
-
+           ELSE
+             DBZ(i,j,ll)  = DBZmin
+             DBZR(i,j,ll) = DBZmin
+             DBZI(i,j,ll) = DBZmin
+           ENDIF !T(I,J,LL)<spval
            ENDDO
 !         parameterized convection
 !         -------------------------
@@ -1439,7 +1461,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(T(I,J,LL)<spval.and.Q(I,J,LL)<spval)THEN
                    GRID1(I,J)=T(I,J,LL)*(1.+D608*Q(I,J,LL))
+                 ELSE
+                   GRID1(I,J)=spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2" )then
@@ -1503,7 +1529,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(P1D(I,J)<spval.and.T1D(I,J)<spval.and.Q(I,J,LL)<spval)THEN
                    GRID1(I,J) = EGRID3(I,J) * (1.+D608*Q(I,J,LL))
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -1541,9 +1571,15 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
               DO J=JSTA,JEND
                 DO I=1,IM
+                IF(P1D(I,J)<spval.and.T1D(I,J)<spval.and.Q1D(I,J)<spval)THEN
                   GRID1(I,J)   = EGRID4(I,J)*100.
                   RH3D(I,J,LL) = GRID1(I,J)
                   EGRID2(I,J)  = Q(I,J,LL)/max(1.e-8,EGRID4(I,J)) ! Revert QS to compute cloud cover later
+                ELSE
+                  GRID1(I,J)   = spval
+                  RH3D(I,J,LL) = spval
+                  EGRID2(I,J)  = spval
+                ENDIF
                 ENDDO
               ENDDO
               IF (item > 0) then
@@ -1579,7 +1615,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(P1D(I,J)<spval.and.T1D(I,J)<spval.and.Q1D(I,J)<spval)THEN
                    GRID1(I,J) = EGRID3(I,J)
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -1630,7 +1670,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(Q(I,J,LL)<spval)THEN
                    GRID1(I,J) = Q(I,J,LL) / (1.-Q(I,J,LL))
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                CALL BOUND(GRID1,H1M12,H99999)
@@ -1668,8 +1712,13 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(Q1D(I,J)<spval.and.EGRID1(I,J)<spval.and.EGRID2(I,J)<spval)THEN
                    GRID1(I,J)   = EGRID3(I,J)
                    MCVG(I,J,LL) = EGRID3(I,J)
+                 ELSE
+                   GRID1(I,J)   = spval
+                   MCVG(I,J,LL) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                IF(IGET(083)>0 .AND. LLL>0)THEN
@@ -1791,7 +1840,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(EGRID3(I,J)<spval)THEN
                    GRID1(I,J) = EGRID3(I,J)
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2003,7 +2056,11 @@ SUBROUTINE MDLFLD
 !$omp  parallel do
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(TRAIN(I,J,LL)<spval)THEN
                    GRID1(I,J) = TRAIN(I,J,LL)*RRNUM
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                ID(1:25) = 0
@@ -2055,7 +2112,11 @@ SUBROUTINE MDLFLD
 !$omp  parallel do
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(TCUCN(I,J,LL)<spval)THEN
                    GRID1(I,J) = TCUCN(I,J,LL)*RRNUM
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                ID(1:25) = 0
@@ -2128,7 +2189,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                DO I=1,IM
+               IF(PMID(I,J,LL)<spval.and.T(I,J,LL)<spval.and.SMOKE(I,J,LL,1)<spval)THEN
                  GRID1(I,J) = (1./RD)*(PMID(I,J,LL)/T(I,J,LL))*SMOKE(I,J,LL,1)
+               ELSE
+                 GRID1(I,J) = spval
+               ENDIF
                ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2153,8 +2218,12 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                DO I=1,IM
+               IF(DUST(I,J,LL,1)<spval.and.RHOMID(I,J,LL)<spval)THEN
                  !GRID1(I,J) = DUST(I,J,LL,1)
                  GRID1(I,J) = DUST(I,J,LL,1)*RHOMID(I,J,LL) !lzhang ug/kg-->ug/m3
+               ELSE
+                 GRID1(I,J) = spval
+               ENDIF               
                ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2179,8 +2248,12 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                DO I=1,IM
+               IF(DUST(I,J,LL,2)<spval.and.RHOMID(I,J,LL)<spval)THEN
                  !GRID1(I,J) = DUST(I,J,LL,2)
                  GRID1(I,J) = DUST(I,J,LL,2)*RHOMID(I,J,LL) !lzhang ug/kg-->ug/m3
+               ELSE
+                 GRID1(I,J) = spval
+               ENDIF
                ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2205,9 +2278,13 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(DUST(I,J,LL,3)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    !GRID1(I,J) = DUST(I,J,LL,3)
                    GRID1(I,J) = DUST(I,J,LL,3)*RHOMID(I,J,LL) !lzhang ug/kg-->ug/m3
-                 ENDDO
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
+               ENDDO
                ENDDO
                if(grib=="grib2") then
                  cfld=cfld+1
@@ -2231,8 +2308,12 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(DUST(I,J,LL,4)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    !GRID1(I,J) = DUST(I,J,LL,4)
                    GRID1(I,J) = DUST(I,J,LL,4)*RHOMID(I,J,LL) !lzhang ug/kg-->ug/m3
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2257,8 +2338,12 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(DUST(I,J,LL,5)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    !GRID1(I,J) = DUST(I,J,LL,5)
                    GRID1(I,J) = DUST(I,J,LL,5)*RHOMID(I,J,LL) !lzhang ug/kg-->ug/m3
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2283,7 +2368,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(SALT(I,J,LL,1)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    GRID1(I,J) = SALT(I,J,LL,1)*RHOMID(I,J,LL) !lzhang ug/kg-->ug/m3
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2308,7 +2397,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(SALT(I,J,LL,2)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    GRID1(I,J) = SALT(I,J,LL,2)*RHOMID(I,J,LL) !lzhang ug/kg-->ug/m3
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2333,7 +2426,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(SALT(I,J,LL,3)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    GRID1(I,J) = SALT(I,J,LL,3)*RHOMID(I,J,LL) !lzhang ug/kg-->ug/m3
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2358,7 +2455,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(SALT(I,J,LL,4)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    GRID1(I,J) = SALT(I,J,LL,4)*RHOMID(I,J,LL) !lzhang ug/kg-->ug/m3
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2383,7 +2484,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(SALT(I,J,LL,5)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    GRID1(I,J) = SALT(I,J,LL,5)*RHOMID(I,J,LL) !lzhang ug/kg-->ug/m3
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2408,8 +2513,12 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(SUSO(I,J,LL,1)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    !GRID1(I,J) = SUSO(I,J,LL,1)
                    GRID1(I,J) = SUSO(I,J,LL,1)*RHOMID(I,J,LL) !lzhang ug/kg-->ug/m3
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2434,8 +2543,12 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(WASO(I,J,LL,1)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    !GRID1(I,J) = WASO(I,J,LL,1)
                    GRID1(I,J) = WASO(I,J,LL,1)*RHOMID(I,J,LL) !lzhang
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2460,8 +2573,12 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(WASO(I,J,LL,2)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    !GRID1(I,J) = WASO(I,J,LL,2)
                    GRID1(I,J) = WASO(I,J,LL,2)*RHOMID(I,J,LL) !lzhang
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2486,8 +2603,12 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(SOOT(I,J,LL,1)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    !GRID1(I,J) = SOOT(I,J,LL,1)
                    GRID1(I,J) = SOOT(I,J,LL,1)*RHOMID(I,J,LL) !lzhang
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2512,8 +2633,12 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                 IF(SOOT(I,J,LL,2)<spval.and.RHOMID(I,J,LL)<spval)THEN
                    !GRID1(I,J) = SOOT(I,J,LL,2)
                    GRID1(I,J) = SOOT(I,J,LL,2)*RHOMID(I,J,LL) !lzhang
+                 ELSE
+                   GRID1(I,J) = spval
+                 ENDIF
                  ENDDO
                ENDDO
                if(grib=="grib2") then
@@ -2726,7 +2851,7 @@ SUBROUTINE MDLFLD
           DO I=1,IM
             GRID1(I,J)=0.0
             DO L=1,NINT(LMH(I,J))
-              if(zint(i,j,l) < spval) then
+              if(zint(i,j,l) < spval .and.zint(i,j,l+1)<spval.and.DBZ(I,J,L)<spval) then
                 GRID1(I,J)=GRID1(I,J)+0.00344* &
                 (10.**(DBZ(I,J,L)/10.))**0.57143*(ZINT(I,J,L)-ZINT(I,J,L+1))/1000.
               endif
@@ -2923,10 +3048,16 @@ SUBROUTINE MDLFLD
             DO I=1,IM
                GRID1(I,J)=0.0
                DO L=1,NINT(LMH(I,J))
+               IF(QQR(I,J,L)<spval.and.QQS(I,J,L)<spval.and.QQG(I,J,L)<spval.and.&
+                  ZINT(I,J,L)<spval.and.ZINT(I,J,L+1)<spval.and.&
+                  PMID(I,J,L)<spval.and.T(I,J,L)<spval.and.Q(I,J,L)<spval)THEN
                   GRID1(I,J)=GRID1(I,J) + (QQR(I,J,L) +      &
                                QQS(I,J,L) + QQG(I,J,L))*     &
                              (ZINT(I,J,L)-ZINT(I,J,L+1))*PMID(I,J,L)/  &
                              (RD*T(I,J,L)*(Q(I,J,L)*D608+1.0))
+               ELSE
+                  GRID1(I,J)=spval
+               ENDIF
                ENDDO
             ENDDO
          ENDDO
@@ -3079,7 +3210,11 @@ SUBROUTINE MDLFLD
 	  IF(abs(vis(i,j))>24135.1)print*,'bad visbility'     &
        , i,j,Q1D(i,j),QW1(i,j),QR1(i,j),QI1(i,j)                 &
        , QS1(i,j),T1D(i,j),P1D(i,j),vis(i,j)
-	  GRID1(I,J)=VIS(I,J)
+          IF(Q1D(I,J)<spval.and.T1D(I,J)<spval.and.P1D(I,J)<spval)THEN
+       	    GRID1(I,J)=VIS(I,J)
+          ELSE
+            GRID1(I,J)=spval
+          ENDIF
 	END DO
 	END DO  
         if(grib=="grib2") then
@@ -3097,7 +3232,11 @@ SUBROUTINE MDLFLD
         CALL CALVIS_GSD(CZEN,VIS)
         DO J=JSTA,JEND
         DO I=1,IM
-          GRID1(I,J)=VIS(I,J)
+          IF(Q1D(I,J)<spval.and.T1D(I,J)<spval.and.P1D(I,J)<spval)THEN
+            GRID1(I,J)=VIS(I,J)
+          ELSE
+            GRID1(I,J)=spval
+          ENDIF
         END DO
         END DO
         if(grib=="grib2") then
@@ -3399,7 +3538,11 @@ SUBROUTINE MDLFLD
 !$omp parallel do private(i,j)
               DO J=JSTA,JEND
                 DO I=1,IM
+                IF(PBLRI(I,J)<spval.and.ZINT(I,J,LM+1)<spval)THEN
                   EGRID3(I,J) = PBLRI(I,J) + ZINT(I,J,LM+1)
+                ELSE
+                  EGRID3(I,J) = spval
+                ENDIF
                 END DO
               END DO  
 ! compute U and V separately because they are on different locations for B grid
@@ -3418,7 +3561,9 @@ SUBROUTINE MDLFLD
                  HCOUNT=0
                  DO J=JSTA,JEND
                   DO I=1,IM
-
+                   if (EGRID4(I,J)<spval.and.EGRID5(I,J)<spval.and.&
+                       EGRID6(I,J)<spval.and.EGRID7(I,J)<spval.and.&
+                       UH(I,J,1)<spval)THEN
                    if (EGRID5(I,J)  <=  EGRID4(I,J)) then
 !       if (I == 50 .and. J == 50) then
 !        write(0,*) 'working with L : ', L
@@ -3430,6 +3575,7 @@ SUBROUTINE MDLFLD
 !                  else
 !                    exit vert_loopu
                    endif
+                   endif
                   end do
                 end do 
                  if(HCOUNT < 1 )exit vert_loopu
@@ -3463,6 +3609,9 @@ SUBROUTINE MDLFLD
 		 HCOUNT=0
                  DO J=JSTA,JEND
                   DO I=1,IM
+                   if (EGRID4(I,J)<spval.and.EGRID5(I,J)<spval.and.&
+                       EGRID6(I,J)<spval.and.EGRID7(I,J)<spval.and.&
+                       VH(I,J,1)<spval)THEN
                    if (EGRID5(I,J) <= EGRID4(I,J)) then
                      HCOUNT=HCOUNT+1
                      DP = EGRID6(I,J) - EGRID7(I,J)
@@ -3471,6 +3620,7 @@ SUBROUTINE MDLFLD
 !                  else
 !                    exit vert_loopu
                    endif
+                   endif 
                   end do
                  end do 
                  if(HCOUNT<1)exit vert_loopv
diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f
index c353f265e..3538cd00e 100644
--- a/sorc/ncep_post.fd/MISCLN.f
+++ b/sorc/ncep_post.fd/MISCLN.f
@@ -44,6 +44,7 @@
 !!   20-03-24  J Meng - remove grib1
 !!   20-11-10  J Meng - USE UPP_PHYSICS MODULE
 !!   21-03-25  E Colon - 3D-RTMA-specific SPC fields added as output
+!!   21-04-01  J Meng - computation on defined points only
 !!     
 !! USAGE:    CALL MISCLN
 !!   INPUT ARGUMENT LIST:
@@ -1489,11 +1490,12 @@ SUBROUTINE MISCLN
           END IF
 
 !        HIGHEST FREEZING LEVEL RELATIVE HUMIDITY
-          IF (IGET(350)>0)THEN  
+          IF (IGET(350)>0)THEN
+               GRID1=spval 
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                DO I=1,IM
-                 GRID1(I,J)=RH1D(I,J)*100.
+                 IF(RH1D(I,J) < spval) GRID1(I,J)=RH1D(I,J)*100.
                ENDDO
                ENDDO
             CALL BOUND (GRID1,H1,H100)
@@ -1562,10 +1564,11 @@ SUBROUTINE MISCLN
 
 !        HIGHEST -10C ISOTHERM RELATIVE HUMIDITY
           IF (IGET(777)>0)THEN  
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                DO I=1,IM
-                 GRID1(I,J)=RH1D(I,J)*100.
+                 IF(RH1D(I,J) < spval) GRID1(I,J)=RH1D(I,J)*100.
                ENDDO
                ENDDO
             CALL BOUND (GRID1,H1,H100)
@@ -1634,10 +1637,11 @@ SUBROUTINE MISCLN
 
 !        HIGHEST -20C ISOTHERM RELATIVE HUMIDITY
           IF (IGET(780)>0)THEN
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                DO I=1,IM
-                 GRID1(I,J)=RH1D(I,J)*100.
+                 IF(RH1D(I,J) < spval) GRID1(I,J)=RH1D(I,J)*100.
                ENDDO
                ENDDO
             CALL BOUND (GRID1,H1,H100)
@@ -2111,6 +2115,7 @@ SUBROUTINE MISCLN
 !$omp parallel do private(i,j)
            DO J=JSTA,JEND
              DO I=1,IM
+
                IF (EGRID1(I,J) > EGRID2(I,J)) THEN
                  EGRID2(I,J) = EGRID1(I,J)
                  LB2(I,J)    = LVLBND(I,J,LBND)
@@ -2212,10 +2217,11 @@ SUBROUTINE MISCLN
             CALL CALLCL(PBND(1,jsta,1),TBND(1,jsta,1),          &
                         QBND(1,jsta,1),EGRID1,EGRID2)
             IF (IGET(109)>0) THEN
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = EGRID2(I,J)
+                   IF(TBND(I,J,1) < spval) GRID1(I,J) = EGRID2(I,J)
                  ENDDO
                ENDDO
                if(grib=='grib2') then
@@ -2231,10 +2237,11 @@ SUBROUTINE MISCLN
                endif
             ENDIF
             IF (IGET(110)>0) THEN
+               GRID1=spval  
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = EGRID1(I,J)
+                   IF(TBND(I,J,1) < spval) GRID1(I,J) = EGRID1(I,J)
                  ENDDO
                ENDDO
                if(grib=='grib2') then
@@ -2368,10 +2375,11 @@ SUBROUTINE MISCLN
 !     
 !           SIGMA 0.89671 TEMPERATURE
              IF (IGET(097) > 0) THEN
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = T89671(I,J)
+                   IF(T(I,J,LM) < spval) GRID1(I,J) = T89671(I,J)
 !                  IF(T89671(I,J)>350.)PRINT*,'LARGE T89671 ',   &
 !                    I,J,T89671(I,J)
                  ENDDO
@@ -2392,10 +2400,11 @@ SUBROUTINE MISCLN
 !     
 !           SIGMA 0.78483 TEMPERATURE
              IF (IGET(098)>0) THEN
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = T78483(I,J)
+                   IF(T(I,J,LM) < spval) GRID1(I,J) = T78483(I,J)
                  ENDDO
                ENDDO
                if(grib=='grib2') then
@@ -2785,10 +2794,11 @@ SUBROUTINE MISCLN
 !           SIGMA 0.85000-1.00000 MOISTURE CONVERGENCE.
             IF (IGET(103)>0) THEN
 !           CONVERT TO DIVERGENCE FOR GRIB
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = -1.0*QM8510(I,J)
+                   IF(QM8510(I,J) < spval) GRID1(I,J) = -1.0*QM8510(I,J)
                  ENDDO
                ENDDO
                if(grib=='grib2') then
@@ -2817,10 +2827,11 @@ SUBROUTINE MISCLN
 !     
 !           SIGMA 0.44-1.00 MEAN RELATIVE HUMIIDITY.
             IF (IGET(318)>0) THEN
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = RH4410(I,J)*100.
+                   IF(RH4410(I,J) < spval) GRID1(I,J) = RH4410(I,J)*100.
                  ENDDO
                ENDDO
                CALL BOUND(GRID1,D00,H100)
@@ -2840,10 +2851,11 @@ SUBROUTINE MISCLN
 !     
 !           SIGMA 0.72-0.94 MEAN RELATIVE HUMIIDITY.
             IF (IGET(319)>0) THEN
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = RH7294(I,J)*100.
+                   IF(RH7294(I,J) < spval) GRID1(I,J) = RH7294(I,J)*100.
                  ENDDO
                ENDDO
                CALL BOUND(GRID1,D00,H100)
@@ -2863,10 +2875,11 @@ SUBROUTINE MISCLN
 !     
 !           SIGMA 0.44-0.72 MEAN RELATIVE HUMIIDITY.
             IF (IGET(320)>0) THEN
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J)=RH4472(I,J)*100.
+                   IF(RH4472(I,J) < spval) GRID1(I,J)=RH4472(I,J)*100.
                  ENDDO
                ENDDO
                CALL BOUND(GRID1,D00,H100)
@@ -2910,9 +2923,11 @@ SUBROUTINE MISCLN
 	   END DO
 ! Temperature	   
 	   IF (IGET(321)>0) THEN
+             GRID1=spval
 !$omp parallel do private(i,j)
              DO J=JSTA,JEND
                DO I=1,IM
+                IF(T(I,J,LM)<spval.and.T(I,J,LM-1)<spval.and.EGRID1(I,J)<spval)&
                  GRID1(I,J) = T(I,J,LM)+(T(I,J,LM-1)-T(I,J,LM)) &
                             * EGRID1(I,J)
                ENDDO
@@ -2934,9 +2949,11 @@ SUBROUTINE MISCLN
             ENDIF
 ! Potential Temperature	    
             IF (IGET(322)>0) THEN
+             GRID2=spval
 !$omp parallel do private(i,j)
              DO J=JSTA,JEND
                DO I=1,IM
+                IF(T(I,J,LM)<spval.and.T(I,J,LM-1)<spval.and.EGRID1(I,J)<spval)&
                  GRID2(I,J) = T(I,J,LM)+(T(I,J,LM-1)-T(I,J,LM))     &
                             * EGRID1(I,J)
                ENDDO
@@ -2957,9 +2974,12 @@ SUBROUTINE MISCLN
             ENDIF
 ! RH	    
             IF (IGET(323)>0) THEN
+                 GRID1=spval
 !$omp parallel do private(i,j,es1,qs1,rh1,es2,qs2,rh2)
              DO J=JSTA,JEND
                DO I=1,IM
+                IF(PMID(I,J,LM)<spval.and.PMID(I,J,LM-1)<spval.and.&
+                   Q(I,J,LM)<spval.and.Q(I,J,LM-1)<spval) THEN
                  ES1 = min(PMID(I,J,LM),FPVSNEW(T(I,J,LM)))
                  QS1 = CON_EPS*ES1/(PMID(I,J,LM)+CON_EPSM1*ES1)
                  RH1 = Q(I,J,LM)/QS1
@@ -2967,6 +2987,7 @@ SUBROUTINE MISCLN
                  QS2 = CON_EPS*ES2/(PMID(I,J,LM-1)+CON_EPSM1*ES2)
                  RH2 = Q(I,J,LM-1)/QS2
                  GRID1(I,J) = (RH1+(RH2-RH1)*EGRID1(I,J))*100.
+                ENDIF
                ENDDO
              ENDDO
              CALL BOUND(GRID1,D00,H100)
@@ -2985,9 +3006,11 @@ SUBROUTINE MISCLN
             ENDIF    
 ! U	   
             IF (IGET(324)>0) THEN
+             GRID1=spval
 !$omp parallel do private(i,j)
              DO J=JSTA,JEND
                DO I=1,IM
+                IF(UH(I,J,LM)<spval.and.UH(I,J,LM-1)<spval.and.EGRID1(I,J)<spval)&
                  GRID1(I,J) = UH(I,J,LM)+(UH(I,J,LM-1)-UH(I,J,LM))    &
                             * EGRID1(I,J)
                ENDDO
@@ -3007,9 +3030,11 @@ SUBROUTINE MISCLN
             ENDIF
 ! V	   
             IF (IGET(325)>0) THEN
+             GRID1=spval
 !$omp parallel do private(i,j)
              DO J=JSTA,JEND
                DO I=1,IM
+                IF(VH(I,J,LM)<spval.and.VH(I,J,LM-1)<spval.and.EGRID1(I,J)<spval)&
                  GRID1(I,J) = VH(I,J,LM)+(VH(I,J,LM-1)-VH(I,J,LM))    &
                             * EGRID1(I,J)
                ENDDO
@@ -3029,9 +3054,11 @@ SUBROUTINE MISCLN
            ENDIF
 ! OMGA	   
            IF (IGET(326)>0) THEN
+             GRID1=spval
 !$omp parallel do private(i,j)
              DO J=JSTA,JEND
                DO I=1,IM
+                IF(OMGA(I,J,LM)<spval.and.OMGA(I,J,LM-1)<spval.and.EGRID1(I,J)<spval)&
                  GRID1(I,J) = OMGA(I,J,LM)+(OMGA(I,J,LM-1)-OMGA(I,J,LM))  &
                             * EGRID1(I,J)
                ENDDO
@@ -3163,9 +3190,10 @@ SUBROUTINE MISCLN
          IF ( (IGET(109)>0).OR.(IGET(110)>0) ) THEN
             CALL CALLCL(P1D,T1D,Q1D,EGRID1,EGRID2)
             IF (IGET(109)>0) THEN
+               GRID1=spval
 	       DO J=JSTA,JEND
                DO I=1,IM
-                 GRID1(I,J)=EGRID2(I,J)
+                 IF(T1D(I,J) < spval) GRID1(I,J)=EGRID2(I,J)
                  IF (SUBMODELNAME == 'RTMA') MLLCL(I,J) = GRID1(I,J)
                ENDDO
                ENDDO
@@ -3294,10 +3322,11 @@ SUBROUTINE MISCLN
               
 !    EQUILLIBRIUM HEIGHT
            IF (IGET(443)>0) THEN
+               GRID1 = spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = EGRID4(I,J)
+                   IF(T1D(I,J) < spval) GRID1(I,J) = EGRID4(I,J)
                  ENDDO
                ENDDO
              if(grib=='grib2') then
@@ -3339,10 +3368,11 @@ SUBROUTINE MISCLN
 !      PRESSURE OF LEVEL FROM WHICH 300 MB MOST UNSTABLE CAPE
 !      PARCEL WAS LIFTED (eq. PRESSURE OF LEVEL OF HIGHEST THETA-E)
            IF (IGET(246)>0) THEN
+              GRID1 = spval
 !$omp parallel do private(i,j)
               DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = EGRID3(I,J)
+                   IF(T1D(I,J) < spval) GRID1(I,J) = EGRID3(I,J)
                  ENDDO
                ENDDO
                CALL BOUND(GRID1,D00,H99999)
@@ -3364,14 +3394,17 @@ SUBROUTINE MISCLN
 
 !    GENERAL THUNDER PARAMETER  ??? 458 ???
         IF (IGET(444)>0) THEN
+               GRID1 = spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
+                  IF(CPRATE(I,J) < spval) THEN
                    IF (CPRATE(I,J) > PTHRESH) THEN
                     GRID1(I,J) = EGRID5(I,J)
                    ELSE
                     GRID1(I,J) = 0
                    ENDIF
+                  ENDIF
                  ENDDO
                ENDDO
                CALL BOUND(GRID1,D00,H99999)
@@ -3505,10 +3538,11 @@ SUBROUTINE MISCLN
 !    LFC HEIGHT
 
             IF (IGET(952)>0) THEN
+              GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = EGRID3(I,J)
+                   IF(T1D(I,J) < spval) GRID1(I,J) = EGRID3(I,J)
                  ENDDO
                ENDDO
                CALL BOUND(GRID1,D00,H99999)
@@ -3641,8 +3675,10 @@ SUBROUTINE MISCLN
 
 !Height of effbot
             IF (IGET(979)>0) THEN
+             GRID1=spval
              DO J=JSTA,JEND
                DO I=1,IM
+                IF(ZINT(I,J,LLOW(I,J))<spval.and.HTSFC(I,J)<spval)&
                  GRID1(I,J) = ZINT(I,J,LLOW(I,J)) - HTSFC(I,J)
                ENDDO
              ENDDO
@@ -3661,8 +3697,10 @@ SUBROUTINE MISCLN
             ENDIF
 !Height of effbot
             IF (IGET(980)>0) THEN
+             GRID1=spval
              DO J=JSTA,JEND
                DO I=1,IM
+                IF(ZINT(I,J,LUPP(I,J))<spval.and.HTSFC(I,J)<spval)&
                  GRID1(I,J) = ZINT(I,J,LUPP(I,J)) - HTSFC(I,J)
                ENDDO
              ENDDO
@@ -3683,14 +3721,18 @@ SUBROUTINE MISCLN
 !U inflow based to 50% EL shear vector
 
             IF (IGET(983)>0) THEN
+             GRID1=spval
              DO J=JSTA,JEND
                DO I=1,IM
+                 IF(LLOW(I,J)<spval.and.LUPP(I,J)<spval.and.&
+                    UH(I,J,MIDCAL)<spval.and.UH(I,J,LLOW(I,J))<spval) THEN
                        MIDCAL=INT(LLOW(I,J)+D50*(LUPP(I,J)-LLOW(I,J)))       
                                                             !mid-layer 
                                                             !vertical
                                                             !index
                        UVECT(I,J)=UH(I,J,MIDCAL)-UH(I,J,LLOW(I,J))
                        GRID1(I,J)=UVECT(I,J)
+                 ENDIF
                ENDDO
              ENDDO
              if(grib=='grib2') then
@@ -3709,14 +3751,18 @@ SUBROUTINE MISCLN
 
 !V inflow based to 50% EL shear vector
             IF (IGET(984)>0) THEN
+             GRID1=spval
              DO J=JSTA,JEND
                DO I=1,IM
+                 IF(LLOW(I,J)<spval.and.LUPP(I,J)<spval.and.&
+                    VH(I,J,MIDCAL)<spval.and.VH(I,J,LLOW(I,J))<spval) THEN
                        MIDCAL=INT(LLOW(I,J)+D50*(LUPP(I,J)-LLOW(I,J)))
                                                             !mid-layer 
                                                             !vertical
                                                             !index
                        VVECT(I,J)=VH(I,J,MIDCAL)-VH(I,J,LLOW(I,J))
                        GRID1(I,J)=VVECT(I,J)
+                 ENDIF
                ENDDO
              ENDDO
              if(grib=='grib2') then
@@ -3735,12 +3781,15 @@ SUBROUTINE MISCLN
 
 !Inflow based (ESFC) to (50%) EL shear magnitude
             IF (IGET(985)>0) THEN
+             GRID1=spval
              DO J=JSTA,JEND
                DO I=1,IM
+                 IF(UVECT(I,J)<spval.and.VVECT(I,J)<spval) THEN
                        ESHR(I,J)=SQRT((UVECT(I,J)**2)+(VVECT(I,J))**2)
                                                                !effshear
                                                                !calc
                        GRID1(I,J)=ESHR(I,J) !Effective
+                 ENDIF
              ENDDO
              ENDDO
              if(grib=='grib2') then
@@ -3765,8 +3814,10 @@ SUBROUTINE MISCLN
 !
 
             IF (IGET(986)>0) THEN
+             GRID1=spval
              DO J=JSTA,JEND
                DO I=1,IM
+                  IF(LLOW(I,J)<spval.and.LUPP(I,J)<spval)&
                        GRID1(I,J)=EFFUST(I,J)
                ENDDO
              ENDDO
@@ -3786,8 +3837,10 @@ SUBROUTINE MISCLN
 
 !V Bunkers Effective right motion
             IF (IGET(987)>0) THEN
+             GRID1=spval
              DO J=JSTA,JEND
                DO I=1,IM
+                  IF(LLOW(I,J)<spval.and.LUPP(I,J)<spval)&
                        GRID1(I,J)=EFFVST(I,J)
                ENDDO
              ENDDO
@@ -3807,8 +3860,10 @@ SUBROUTINE MISCLN
 
 !Effective layer helicity
             IF (IGET(988)>0) THEN
+             GRID1=spval
              DO J=JSTA,JEND
                DO I=1,IM
+                  IF(LLOW(I,J)<spval.and.LUPP(I,J)<spval)&
                        GRID1(I,J)=ESRH(I,J)
                ENDDO
              ENDDO
@@ -3853,11 +3908,14 @@ SUBROUTINE MISCLN
                 ENDIF
                 STP=(MLCAPE(I,J)/D1500)*MLLCLtmp*(ESRH(I,J)/150.)*&
                         ESHRtmp*MLCINtmp
+                GRID1(I,J) = SPVAL
+                IF(LLOW(I,J)<spval.and.LUPP(I,J)<spval) THEN
                 IF (STP>0) THEN
                    GRID1(I,J)=STP
                 ELSE
                    GRID1(I,J)=D00
                 ENDIF
+                ENDIF
                ENDDO
             ENDDO
             if(grib=='grib2') then
@@ -3923,11 +3981,14 @@ SUBROUTINE MISCLN
                 ENDIF
                 STP=(EGRID1(I,J)/D1500)*SLCLtmp*(HELI(I,J,2)/150.)*&
                         FSHRtmp*SCINtmp
+                GRID1(I,J) = spval
+                IF(T1D(I,J) < spval) THEN
                 IF (STP>0) THEN
                    GRID1(I,J)=STP
                 ELSE
                    GRID1(I,J)=D00
                 ENDIF
+                ENDIF
                ENDDO
             ENDDO
             if(grib=='grib2') then
@@ -3962,11 +4023,14 @@ SUBROUTINE MISCLN
                 ENDIF
                 STP=(MUCAPE(I,J)/D1000)*(ESRH(I,J)/50.)*&
                         ESHRtmp*MUCINtmp
+                GRID1(I,J) = spval
+                IF(T1D(I,J) < spval) THEN
                 IF (STP>0) THEN
                    GRID1(I,J)=STP
                 ELSE
                    GRID1(I,J)=D00
                 ENDIF
+                ENDIF
                ENDDO
             ENDDO
             if(grib=='grib2') then
@@ -4013,10 +4077,10 @@ SUBROUTINE MISCLN
              CALL CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,LB2,            &
                            EGRID1,EGRID2,EGRID3,EGRID4,EGRID5,     &
                            EGRID6,EGRID7,EGRID8)
-
+             GRID1=spval
              DO J=JSTA,JEND
                DO I=1,IM
-                 GRID1(I,J) = EGRID3(I,J)
+                 IF(T1D(I,J) < spval) GRID1(I,J) = EGRID3(I,J)
                ENDDO
              ENDDO
              CALL BOUND(GRID1,D00,H99999)
@@ -4059,8 +4123,12 @@ SUBROUTINE MISCLN
 
 !Hail parameter
             IF (IGET(993)>0) THEN
+            GRID1=spval
             DO J=JSTA,JEND
                DO I=1,IM
+               IF(T700(I,J) < spval .and. T500(I,J) < spval .and.&
+                  Z700(I,J) < spval .and. Z500(I,J) < spval .and.&
+                  MUCAPE(I,J) < spval .and. MUQ1D(I,J) < spval .and. FSHR(I,J) < spval) THEN
                LAPSE=-((T700(I,J)-T500(I,J))/((Z700(I,J)-Z500(I,J))))
                 SHIP=(MUCAPE(I,J)*D1000*MUQ1D(I,J)*LAPSE*(T500(I,J)-K2C)*FSHR(I,J))/HCONST
                 IF (MUCAPE(I,J)<1300.)THEN
@@ -4073,6 +4141,7 @@ SUBROUTINE MISCLN
                    SHIP=SHIP*(FREEZELVL(I,J)/2400.)
                 ENDIF
                 GRID1(I,J)=SHIP
+               ENDIF
                ENDDO
             ENDDO
             if(grib=='grib2') then
@@ -4095,10 +4164,11 @@ SUBROUTINE MISCLN
 !    Critical Angle
 
             IF (IGET(957)>0) THEN
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                    GRID1(I,J) = CANGLE(I,J)   
+                   IF(T1D(I,J) < spval ) GRID1(I,J) = CANGLE(I,J)   
            !         IF(EGRID1(I,J)<100. .OR. EGRID2(I,J)>-250.) THEN
            !           GRID1(I,J) = 0.
            !         ENDIF
@@ -4121,10 +4191,11 @@ SUBROUTINE MISCLN
 !    Dendritic Layer Depth, -17C < T < -12C
 
             IF (IGET(955)>0) THEN
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = EGRID7(I,J)
+                   IF(T1D(I,J) < spval ) GRID1(I,J) = EGRID7(I,J)
                  ENDDO
                ENDDO
                CALL BOUND(GRID1,D00,H99999)
@@ -4145,10 +4216,11 @@ SUBROUTINE MISCLN
 !    Enhanced Stretching Potential
 
             IF (IGET(956)>0) THEN
+               GRID1=spval
 !$omp parallel do private(i,j)
                DO J=JSTA,JEND
                  DO I=1,IM
-                   GRID1(I,J) = EGRID8(I,J)
+                   IF(T1D(I,J) < spval ) GRID1(I,J) = EGRID8(I,J)
                  ENDDO
                ENDDO
                CALL BOUND(GRID1,D00,H99999)
diff --git a/sorc/ncep_post.fd/SURFCE.f b/sorc/ncep_post.fd/SURFCE.f
index 73d125375..3af081583 100644
--- a/sorc/ncep_post.fd/SURFCE.f
+++ b/sorc/ncep_post.fd/SURFCE.f
@@ -37,6 +37,7 @@
 !! -  20-05-20  J MENG    - CALRH unification with NAM scheme
 !! -  20-11-10  J MENG    - USE UPP_PHYSICS MODULE
 !! -  21-03-11  B Cui - change local arrays to dimension (im,jsta:jend)
+!! -  21-04-01  J MENG    - COMPUTATION ON DEFINED POINTS ONLY
 !!     
 !! USAGE:    CALL SURFCE
 !!   INPUT ARGUMENT LIST:
@@ -1252,7 +1253,7 @@ SUBROUTINE SURFCE
 !
 !        SHELTER LEVEL TEMPERATURE
         IF (IGET(106)>0) THEN
-!          GRID1=spval
+           GRID1=SPVAL
            DO J=JSTA,JEND
              DO I=1,IM
 !              GRID1(I,J)=TSHLTR(I,J)
@@ -1395,14 +1396,18 @@ SUBROUTINE SURFCE
 
 !
            IF ((IGET(547)>0).OR.(IGET(548)>0)) THEN
+            GRID1=SPVAL
+            GRID2=SPVAL
             DO J=JSTA,JEND
             DO I=1,IM
+            if(TSHLTR(I,J)/=spval.and.PSHLTR(I,J)/=spval.and.QSHLTR(I,J)/=spval) then
 ! DEWPOINT DEPRESSION in GRID1
              GRID1(i,j)=max(0.,TSHLTR(I,J)*(PSHLTR(I,J)*1.E-5)**CAPA-EGRID1(i,j))
 
 ! SURFACE EQIV POT TEMP in GRID2
              APE=(H10E5/PSHLTR(I,J))**CAPA
              GRID2(I,J)=TSHLTR(I,J)*EXP(ELOCP*QSHLTR(I,J)*APE/TSHLTR(I,J))
+            endif 
             ENDDO
             ENDDO
 !       print *,' MAX/MIN --> DEWPOINT DEPRESSION',maxval(grid1(1:im,jsta:jend)),&
@@ -1478,9 +1483,11 @@ SUBROUTINE SURFCE
            ENDIF
 
            IF(IGET(808)>0)THEN
+             GRID2=SPVAL
 !$omp parallel do private(i,j,dum1,dum2,dum3,dum216,dum1s,dum3s)
              DO J=JSTA,JEND
                DO I=1,IM
+               if(T1D(I,J)/=spval.and.U10H(I,J)/=spval.and.V10H(I,J)<spval) then
                  DUM1 = (T1D(I,J)-TFRZ)*1.8+32.
                  DUM2 = SQRT(U10H(I,J)**2.0+V10H(I,J)**2.0)/0.44704
                  DUM3 = EGRID1(I,J) * 100.0
@@ -1508,6 +1515,7 @@ SUBROUTINE SURFCE
                  END IF
 !                if(abs(gdlon(i,j)-120.)<1. .and. abs(gdlat(i,j))<1.) &
 !                 print*,'Debug AT: OUTPUT',Grid2(i,j)
+               endif
                ENDDO
              ENDDO
 
@@ -1647,9 +1655,10 @@ SUBROUTINE SURFCE
 !
 !        SHELTER LEVEL MAX RH.
          IF (IGET(347)>0) THEN       
+         GRID1=SPVAL
             DO J=JSTA,JEND
             DO I=1,IM
-             GRID1(I,J)=MAXRHSHLTR(I,J)*100.
+             if(MAXRHSHLTR(I,J)/=spval) GRID1(I,J)=MAXRHSHLTR(I,J)*100.
             ENDDO
             ENDDO
 	    ID(1:25) = 0
@@ -1698,9 +1707,10 @@ SUBROUTINE SURFCE
 !
 !        SHELTER LEVEL MIN RH.
          IF (IGET(348)>0) THEN       
+            GRID1=SPVAL
             DO J=JSTA,JEND
             DO I=1,IM
-             GRID1(I,J)=MINRHSHLTR(I,J)*100.
+             if(MINRHSHLTR(I,J)/=spval) GRID1(I,J)=MINRHSHLTR(I,J)*100.
             ENDDO
             ENDDO
 	    ID(1:25) = 0
@@ -1827,8 +1837,10 @@ SUBROUTINE SURFCE
 ! E. James - 12 Sep 2018: SMOKE from WRF-CHEM on lowest model level
 !
          IF (IGET(739)>0) THEN
+           GRID1=SPVAL
            DO J=JSTA,JEND
              DO I=1,IM
+             if(T(I,J,LM)/=spval.and.PMID(I,J,LM)/=spval.and.SMOKE(I,J,LM,1)/=spval)&
                GRID1(I,J) = (1./RD)*(PMID(I,J,LM)/T(I,J,LM))*SMOKE(I,J,LM,1)
              ENDDO
            ENDDO
@@ -2265,10 +2277,11 @@ SUBROUTINE SURFCE
       IF (IGET(249)>0) THEN
          RDTPHS=1000./DTQ2     !--- 1000 kg/m**3, density of liquid water
 !        RDTPHS=1000./(TRDLW*3600.)
+         GRID1=SPVAL
 !$omp parallel do private(i,j)
          DO J=JSTA,JEND
            DO I=1,IM
-             GRID1(I,J) = CPRATE(I,J)*RDTPHS
+            if(CPRATE(I,J)/=spval) GRID1(I,J) = CPRATE(I,J)*RDTPHS
 !             GRID1(I,J) = CUPPT(I,J)*RDTPHS
            ENDDO
          ENDDO
@@ -2290,14 +2303,17 @@ SUBROUTINE SURFCE
 !MEB need to get physics DT
          RDTPHS=1./(DTQ2) 
 !MEB need to get physics DT
+         GRID1=SPVAL
 !$omp parallel do private(i,j)
          DO J=JSTA,JEND
            DO I=1,IM
+           if(PREC(I,J)/=spval) then
              IF(MODELNAME /= 'RSM') THEN
               GRID1(I,J) = PREC(I,J)*RDTPHS*1000.
              ELSE        !Add by Binbin
               GRID1(I,J) = PREC(I,J)
              END IF
+           endif
            ENDDO
          ENDDO
          if(grib=='grib2') then
@@ -2316,9 +2332,10 @@ SUBROUTINE SURFCE
 ! MAXIMUM INSTANTANEOUS PRECIPITATION RATE.
       IF (IGET(508)>0) THEN
 !-- PRATE_MAX in units of mm/h from NMMB history files
+         GRID1=SPVAL
          DO J=JSTA,JEND
            DO I=1,IM
-             GRID1(I,J)=PRATE_MAX(I,J)*SEC2HR
+            if(PRATE_MAX(I,J)/=spval) GRID1(I,J)=PRATE_MAX(I,J)*SEC2HR
            ENDDO
          ENDDO
          if(grib=='grib2') then
@@ -2344,9 +2361,10 @@ SUBROUTINE SURFCE
 ! MAXIMUM INSTANTANEOUS *FROZEN* PRECIPITATION RATE.
       IF (IGET(509)>0) THEN
 !-- FPRATE_MAX in units of mm/h from NMMB history files
+         GRID1=SPVAL
          DO J=JSTA,JEND
            DO I=1,IM
-             GRID1(I,J)=FPRATE_MAX(I,J)*SEC2HR
+            if(FPRATE_MAX(I,J)/=spval) GRID1(I,J)=FPRATE_MAX(I,J)*SEC2HR
            ENDDO
          ENDDO
          if(grib=='grib2') then
@@ -2526,7 +2544,11 @@ SUBROUTINE SURFCE
 !$omp parallel do private(i,j)
            DO J=JSTA,JEND
              DO I=1,IM
+              IF(ACPREC(I,J) < SPVAL)THEN
                GRID1(I,J) = ACPREC(I,J)*1000.
+              ELSE
+               GRID1(I,J) = SPVAL
+              ENDIF
              ENDDO
            ENDDO
          END IF 
@@ -2672,7 +2694,11 @@ SUBROUTINE SURFCE
 !$omp parallel do private(i,j)
            DO J=JSTA,JEND
             DO I=1,IM
+            IF(CUPREC(I,J) < SPVAL)THEN
              GRID1(I,J) = CUPREC(I,J)*1000.
+            ELSE
+             GRID1(I,J) = SPVAL
+            ENDIF
             ENDDO
            ENDDO
          END IF 
@@ -2902,13 +2928,14 @@ SUBROUTINE SURFCE
 !     
 !     ACCUMULATED LAND SURFACE PRECIPITATION.
       IF (IGET(256)>0) THEN
+      GRID1=SPVAL
 !$omp parallel do private(i,j)
          DO J=JSTA,JEND
            DO I=1,IM
             IF(LSPA(I,J)<=-1.0E-6)THEN
-             GRID1(I,J) = ACPREC(I,J)*1000
+             if(ACPREC(I,J)/=spval) GRID1(I,J) = ACPREC(I,J)*1000
             ELSE
-             GRID1(I,J) = LSPA(I,J)*1000.
+             if(LSPA(I,J)/=spval) GRID1(I,J) = LSPA(I,J)*1000.
             END IF
            ENDDO
          ENDDO
@@ -4735,9 +4762,10 @@ SUBROUTINE SURFCE
             ELSE
               RRNUM=0.
             ENDIF
+            GRID1=SPVAL
             DO J=JSTA,JEND
             DO I=1,IM
-             GRID1(I,J) = SUBSHX(I,J)*RRNUM
+             if(SUBSHX(I,J)/=spval) GRID1(I,J) = SUBSHX(I,J)*RRNUM
             ENDDO
             ENDDO
             ID(1:25) = 0
@@ -4784,9 +4812,10 @@ SUBROUTINE SURFCE
             ELSE
               RRNUM=0.
             ENDIF
+            GRID1=SPVAL
             DO J=JSTA,JEND
             DO I=1,IM
-             GRID1(I,J) = SNOPCX(I,J)*RRNUM
+             if(SNOPCX(I,J)/=spval) GRID1(I,J) = SNOPCX(I,J)*RRNUM
             ENDDO
             ENDDO
             ID(1:25) = 0
@@ -4886,9 +4915,10 @@ SUBROUTINE SURFCE
             ELSE
               RRNUM=0.
             ENDIF
+            GRID1=SPVAL
             DO J=JSTA,JEND
             DO I=1,IM
-             GRID1(I,J) = SFCUX(I,J)*RRNUM
+             if(SFCUX(I,J)/=spval) GRID1(I,J) = SFCUX(I,J)*RRNUM
             ENDDO
             ENDDO
             ID(1:25) = 0
@@ -4935,9 +4965,10 @@ SUBROUTINE SURFCE
             ELSE
               RRNUM=0.
             ENDIF
+            GRID1=SPVAL
             DO J=JSTA,JEND
             DO I=1,IM
-             GRID1(I,J) = SFCVX(I,J)*RRNUM
+             if(SFCVX(I,J)/=spval) GRID1(I,J) = SFCVX(I,J)*RRNUM
             ENDDO
             ENDDO
             ID(1:25) = 0
@@ -4974,9 +5005,10 @@ SUBROUTINE SURFCE
 !     
 !     ACCUMULATED SURFACE EVAPORATION
          IF (IGET(047)>0) THEN
+            GRID1=SPVAL
             DO J=JSTA,JEND
               DO I=1,IM
-                GRID1(I,J) = SFCEVP(I,J)*1000.
+               if(SFCEVP(I,J)/=spval) GRID1(I,J) = SFCEVP(I,J)*1000.
               ENDDO
             ENDDO
             ID(1:25) = 0
@@ -5016,9 +5048,10 @@ SUBROUTINE SURFCE
 !     
 !     ACCUMULATED POTENTIAL EVAPORATION
          IF (IGET(137)>0) THEN
+            GRID1=SPVAL
             DO J=JSTA,JEND
               DO I=1,IM
-               GRID1(I,J) = POTEVP(I,J)*1000.
+               if(POTEVP(I,J)/=spval) GRID1(I,J) = POTEVP(I,J)*1000.
               ENDDO
             ENDDO
             ID(1:25) = 0
@@ -5351,9 +5384,10 @@ SUBROUTINE SURFCE
 !     
 !     GREEN VEG FRACTION
       IF (IGET(170)>0) THEN
+            GRID1=SPVAL
             DO J=JSTA,JEND
             DO I=1,IM
-             GRID1(I,J)=VEGFRC(I,J)*100.
+             if(VEGFRC(I,J)/=spval) GRID1(I,J)=VEGFRC(I,J)*100.
             ENDDO
             ENDDO
            if(grib=='grib2') then
@@ -5366,9 +5400,10 @@ SUBROUTINE SURFCE
 !
 !     MIN GREEN VEG FRACTION
       IF (IGET(726)>0) THEN
+            GRID1=SPVAL
             DO J=JSTA,JEND
             DO I=1,IM
-             GRID1(I,J)=shdmin(I,J)*100.
+             if(shdmin(I,J)/=spval) GRID1(I,J)=shdmin(I,J)*100.
             ENDDO
             ENDDO
           if(grib=='grib2') then
@@ -5380,9 +5415,10 @@ SUBROUTINE SURFCE
 !
 !     MAX GREEN VEG FRACTION
       IF (IGET(729)>0) THEN
+            GRID1=SPVAL
             DO J=JSTA,JEND
             DO I=1,IM
-             GRID1(I,J)=shdmax(I,J)*100.
+             if(shdmax(I,J)/=spval) GRID1(I,J)=shdmax(I,J)*100.
             ENDDO
             ENDDO
           if(grib=='grib2') then
diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f
index fc0a9c22d..b9ca9995d 100644
--- a/sorc/ncep_post.fd/UPP_PHYSICS.f
+++ b/sorc/ncep_post.fd/UPP_PHYSICS.f
@@ -284,7 +284,7 @@ SUBROUTINE CALRH_GSD(P1,T1,Q1,RHB)
 !------------------------------------------------------------------
 !
 
-      use ctlblk_mod, only: jsta, jend, im
+      use ctlblk_mod, only: jsta, jend, im, spval
 
       implicit none
 
@@ -295,7 +295,7 @@ SUBROUTINE CALRH_GSD(P1,T1,Q1,RHB)
 
       DO J=JSTA,JEND
         DO I=1,IM
-
+         IF (T1(I,J) < SPVAL .AND. P1(I,J) < SPVAL .AND. Q1(I,J) < SPVAL) THEN
 ! - compute relative humidity
           Tx=T1(I,J)-273.15
           POL = 0.99999683       + TX*(-0.90826951E-02 +    &
@@ -308,7 +308,9 @@ SUBROUTINE CALRH_GSD(P1,T1,Q1,RHB)
           ES = esx
           E = P1(I,J)/100.*Q1(I,J)/(0.62197+Q1(I,J)*0.37803)
           RHB(I,J) = MIN(1.,E/ES)
-
+         ELSE
+          RHB(I,J) = SPVAL
+         ENDIF
         ENDDO
       ENDDO
 
@@ -324,7 +326,7 @@ SUBROUTINE CALRH_PW(RHPW)
 
       use vrbls3d, only: q, pmid, t
       use params_mod, only: g
-      use ctlblk_mod, only: lm, jsta, jend, lm, im
+      use ctlblk_mod, only: lm, jsta, jend, lm, im, spval
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        implicit none
 
@@ -342,6 +344,7 @@ SUBROUTINE CALRH_PW(RHPW)
         k=lm-l+1
        DO J=JSTA,JEND
         DO I=1,IM
+          if(t(i,j,k)<spval.and.q(i,j,k)<spval) then
 ! -- use specific humidity for PW calculation
            sh = q(i,j,k)
            qv = sh/(1.-sh)
@@ -371,7 +374,9 @@ SUBROUTINE CALRH_PW(RHPW)
 
 !sgb - This IS RH w.r.t. PW-sat.
            RHPW (i,j) = min(1.,PW(i,j) / pw_sat(i,j)) * 100.
-
+          else
+           RHPW (i,j) = spval
+          endif
         ENDDO
        ENDDO
       ENDDO