Fix mismatch isltyp and landmask (#746)

TYPE: bug fix

KEYWORDS: surface_input_source, isltyp, landmask

SOURCE: internal (reported by many users, including Priscilla Mooney from Norway)

DESCRIPTION OF CHANGES:
To accommodate users who would like to change LU_INDEX after running geogrid program, we changed namelist surface_input_source from default value of 1 (recomputing dominant categories in real) to 3 (using dominant categories from geogrid) in v3.8.

Unfortunately, we neglected to consider in the v3.8 fix that when dominant categories for land and soil are recomputed in real, the real program also does mismatch checks: ensuring that the lower resolution ISLTYP matches with LANDMASK data. This check in real has two two effects:

the actual ISLTYP, IVGTYP, and XLAND can be different from that from surface_input_source = 1 (a check later in real program effectively matched IVGTYP based on ISLTYP);
for certain cases, the real program would stop due to un-matched ISLTYP and P_.
This PR fixes the problem of this missing consistency check in the real program when surface_input_source = 3. The ISLTYP, IVGTYP, and XLAND are now identical to those coming from real after setting surface_input_source = 1.

LIST OF MODIFIED FILES:
M dyn_em/module_initialize_real.F

TESTS CONDUCTED:

Ran real before and after the change and results are correct.
Ran real for a case that failed before and now it succeeds.

RELEASE NOTE:
In V3.8, we changed default surface_input_source option value from 1 to 3, which uses the dominant categories computed in geogrid. But we neglected that when dominant categories for land and soil are recomputed in real, it also does mismatch checks, making sure the lower resolution isltyp matches with landmask data. This has two two effects: one is the actual isltyp, ivgtyp and xland can be different from that from surface_input_source = 1 (a check later in real program effectively matched ivgtyp based on isltyp); and the other effect is that for certain cases, the real program would stop due to un-matched isltyp and ivgtyp. This change fixes this problem, and isltyp, ivgtyp, and xland are now identical to those coming from surface_input_source = 1.

dyn_em/module_initialize_real.F

@@ -165,6 +165,7 @@ SUBROUTINE init_domain_rk ( grid &
       LOGICAL :: vnest !T if using vertical nesting with vet_refine_method=2, otherwise F
 
       INTEGER :: j_save
+      INTEGER :: change_soil, change_soilw, iforce
 
       LOGICAL :: wif_upside_down = .FALSE.
 
@@ -2995,6 +2996,33 @@ SUBROUTINE init_domain_rk ( grid &
             CALL wrf_error_fatal ( a_message )
          END IF
 
+         !  Need to match isltyp to landmask
+
+         iforce = 0
+         change_soil = 0
+         change_soilw = 0
+         DO j = jts, MIN(jde-1,jte)
+            DO i = its, MIN(ide-1,ite)
+               IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
+                  IF ( grid%landmask(i,j) .GT. 0.5 .AND. grid%isltyp(i,j) .EQ. grid%isoilwater ) THEN
+                     grid%isltyp(i,j) = 8
+                     change_soilw =  change_soilw + 1
+                     iforce =  iforce + 1
+                  ELSE IF ( grid%landmask(i,j) .LT. 0.5 .AND. grid%isltyp(i,j) .NE. grid%isoilwater ) THEN
+                     grid%isltyp(i,j) = grid%isoilwater
+                     change_soil =  change_soil + 1
+                     iforce =  iforce + 1
+                  END IF
+               END IF
+            END DO
+         END DO
+         IF ( change_soilw .GT. 0 .OR. change_soil .GT. 0 ) THEN
+            WRITE(a_message,FMT='(A,I4,A,I6)' ) &
+                  'forcing artificial silty clay loam at ',iforce,' points, out of ',&
+                  (MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
+                  CALL wrf_debug(0,a_message)
+         END IF
+
       END IF
 
       !  Split NUDAPT  Urban Parameters