global_cycle - Support NoahMP and fractional grid

sorc/global_cycle.fd/cycle.f90

@@ -109,12 +109,12 @@ PROGRAM SFC_DRV
  INTEGER :: ISOT, IVEGSRC, LENSFC, ZSEA1_MM, ZSEA2_MM, IERR
  INTEGER :: NPROCS, MYRANK, NUM_THREADS, NUM_PARTHDS, MAX_TASKS
  REAL    :: FH, DELTSFC, ZSEA1, ZSEA2
- LOGICAL :: USE_UFO, DO_NSST, DO_LNDINC, DO_SFCCYCLE
+ LOGICAL :: USE_UFO, DO_NSST, DO_LNDINC, DO_SFCCYCLE, FRAC_GRID
 !
  NAMELIST/NAMCYC/ IDIM,JDIM,LSOIL,LUGB,IY,IM,ID,IH,FH,&
                   DELTSFC,IALB,USE_UFO,DONST,             &
                   DO_SFCCYCLE,ISOT,IVEGSRC,ZSEA1_MM,      &
-                  ZSEA2_MM, MAX_TASKS, DO_LNDINC
+                  ZSEA2_MM, MAX_TASKS, DO_LNDINC, FRAC_GRID
 !
  DATA IDIM,JDIM,LSOIL/96,96,4/
  DATA IY,IM,ID,IH,FH/1997,8,2,0,0./
@@ -138,13 +138,14 @@ PROGRAM SFC_DRV
  DONST   = "NO"
  DO_LNDINC   = .FALSE.
  DO_SFCCYCLE = .TRUE.
+ FRAC_GRID = .FALSE.
 
  PRINT*
  PRINT*,"READ NAMCYC NAMELIST."
 
  CALL BAOPENR(36, "fort.36", IERR)
  READ(36, NML=NAMCYC)
- !IF (MYRANK==0) WRITE(6,NAMCYC)
+!IF (MYRANK==0) WRITE(6,NAMCYC)
 
  IF (MAX_TASKS < 99999 .AND. MYRANK > (MAX_TASKS - 1)) THEN
    PRINT*,"USER SPECIFIED MAX NUMBER OF TASKS: ", MAX_TASKS
@@ -170,7 +171,7 @@ PROGRAM SFC_DRV
  CALL SFCDRV(LUGB,IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
              IY,IM,ID,IH,FH,IALB,                  &
              USE_UFO,DO_NSST,DO_SFCCYCLE,DO_LNDINC, &
-             ZSEA1,ZSEA2,ISOT,IVEGSRC,MYRANK)
+             FRAC_GRID,ZSEA1,ZSEA2,ISOT,IVEGSRC,MYRANK)
 
  PRINT*
  PRINT*,'CYCLE PROGRAM COMPLETED NORMALLY ON RANK: ', MYRANK
@@ -286,6 +287,7 @@ END PROGRAM SFC_DRV
  !! @param[in] DO_SFCCYCLE Call sfccycle routine to update surface fields
  !! @param[in] DO_LNDINC Read in land increment files, and add increments to
  !!            requested states.
+ !! @param[in] FRAC_GRID When true, run with fractional grid.
  !! @param[in] ZSEA1 When running NSST model, this is the lower bound
  !!            of depth of sea temperature.  In whole mm.
  !! @param[in] ZSEA2 When running NSST model, this is the upper bound
@@ -297,7 +299,7 @@ END PROGRAM SFC_DRV
  SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
                    IY,IM,ID,IH,FH,IALB,                  &
                    USE_UFO,DO_NSST,DO_SFCCYCLE,DO_LNDINC,&
-                   ZSEA1,ZSEA2,ISOT,IVEGSRC,MYRANK)
+                   FRAC_GRID,ZSEA1,ZSEA2,ISOT,IVEGSRC,MYRANK)
 !
  USE READ_WRITE_DATA
  use machine
@@ -316,7 +318,7 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
  INTEGER, INTENT(IN) :: ISOT, IVEGSRC, MYRANK
 
  LOGICAL, INTENT(IN) :: USE_UFO, DO_NSST,DO_SFCCYCLE
- LOGICAL, INTENT(IN) :: DO_LNDINC
+ LOGICAL, INTENT(IN) :: DO_LNDINC, FRAC_GRID
 
  REAL, INTENT(IN)    :: FH, DELTSFC, ZSEA1, ZSEA2
 
@@ -332,6 +334,7 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
  INTEGER             :: I_INDEX(LENSFC), J_INDEX(LENSFC)
  INTEGER             :: IDUM(IDIM,JDIM)
  integer             :: num_parthds, num_threads
+! integer :: ichk, jchk, ijchk
 
  LOGICAL             :: IS_NOAHMP
  INTEGER             :: LSM
@@ -366,7 +369,7 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
  REAL, ALLOCATABLE   :: SLIFCS_FG(:)
  INTEGER, ALLOCATABLE :: LANDINC_MASK_FG(:), LANDINC_MASK(:)
  REAL, ALLOCATABLE   :: SND_BCK(:), SND_INC(:), SWE_BCK(:)
- REAL(KIND=KIND_IO8), ALLOCATABLE :: SLMASKL(:), SLMASKW(:)
+ REAL(KIND=KIND_IO8), ALLOCATABLE :: SLMASKL(:), SLMASKW(:), LANDFRAC(:)
 
  TYPE(NSST_DATA)     :: NSST
  real, dimension(idim,jdim) :: tf_clm,tf_trd,sal_clm
@@ -389,7 +392,6 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
  DO_SNO_INC = .FALSE.
  DO_SOI_INC = .FALSE.
 
-
  SIG1T = 0.0            ! Not a dead start!
 
  INPUT_NML_FILE = "NULL"
@@ -404,7 +406,13 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
 ! READ THE OROGRAPHY AND GRID POINT LAT/LONS FOR THE CUBED-SPHERE TILE.
 !--------------------------------------------------------------------------------
 
- CALL READ_LAT_LON_OROG(RLA,RLO,OROG,OROG_UF,TILE_NUM,IDIM,JDIM,LENSFC)
+ IF(FRAC_GRID) THEN
+   PRINT*,'- RUNNING WITH FRACTIONAL GRID.'
+   ALLOCATE(LANDFRAC(LENSFC))
+   CALL READ_LAT_LON_OROG(RLA,RLO,OROG,OROG_UF,TILE_NUM,IDIM,JDIM,LENSFC,LANDFRAC=LANDFRAC)
+ ELSE
+   CALL READ_LAT_LON_OROG(RLA,RLO,OROG,OROG_UF,TILE_NUM,IDIM,JDIM,LENSFC)
+ ENDIF
 
  DO I = 1, IDIM
    IDUM(I,:) = I
@@ -444,6 +452,10 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
  ENDIF
 
 IF (DO_LNDINC) THEN
+   IF (FRAC_GRID) THEN
+     print *, 'FATAL ERROR: land increment update does not work with fractional grids.'
+     call MPI_ABORT(MPI_COMM_WORLD, 17, IERR)
+   ENDIF
    ! identify variables to be updates, and allocate arrays.
    IF  (TRIM(LND_SOI_FILE) .NE. "NULL") THEN
        DO_SOI_INC = .TRUE.
@@ -486,6 +498,28 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
                 ABSFCS=ABSFCS,T2M=T2M      ,Q2M=Q2M      ,SLMASK=SLMASK,  &
                 ZSOIL=ZSOIL,   NSST=NSST)
 
+ print*,'is noahmp/fract grid ',is_noahmp, frac_grid
+
+ IF (FRAC_GRID .AND. .NOT. IS_NOAHMP) THEN
+   print *, 'FATAL ERROR: NOAH lsm update does not work with fractional grids.'
+   call MPI_ABORT(MPI_COMM_WORLD, 18, IERR)
+ ENDIF
+
+!ichk = 282
+!jchk = 362
+!ijchk = (jchk-1) * idim + ichk
+!print*,'chk point after read_data ',tile_num,landfrac(ijchk),slmask(ijchk),vetfcs(ijchk),vegfcs(ijchk)
+
+!print*,'vegfcs check ',maxval(vegfcs),minval(vegfcs)
+!do i=1,lensfc
+!  if(landfrac(i) > 0.0_kind_io8 .and. vegfcs(i) == 0.0) then
+!    print*,'bad point 1 ',tile_num,i,landfrac(i),vegfcs(i)
+!  endif
+!  if(landfrac(i) == 0.0_kind_io8 .and. vegfcs(i) > 0.0) then
+!    print*,'bad point 2 ',tile_num,i,landfrac(i),vegfcs(i)
+!  endif
+!enddo
+
  IF (IS_NOAHMP) THEN 
         LSM=LSM_NOAHMP
  ELSE
@@ -532,21 +566,56 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
 
  IF (DO_SFCCYCLE) THEN
    ALLOCATE(SLMASKL(LENSFC), SLMASKW(LENSFC))
+   IF (FRAC_GRID) THEN
+
+     DO I=1,LENSFC
+       IF(LANDFRAC(I) > 0.0_KIND_IO8) THEN
+         SLMASKL(I) = CEILING(LANDFRAC(I)-1.0E-6_KIND_IO8)
+         SLMASKW(I) =   FLOOR(LANDFRAC(I)+1.0E-6_KIND_IO8)
+       ELSE
+         IF(NINT(SLMASK(I)) == 1) THEN ! if landfrac is zero, this should not happen.
+                                       ! is this a band aid?
+            print*,'first if'
+           SLMASKL(I) = 1.0_KIND_io8
+           SLMASKW(I) = 1.0_KIND_io8
+         ELSE ! if landfrac is zero, this is expected.
+!           print*,'second if'
+           SLMASKL(I) = 0.0_KIND_io8
+           SLMASKW(I) = 0.0_KIND_io8
+         ENDIF
+       ENDIF
+
+       if(nint(slmaskl(i)) == 1 .and. vegfcs(i) == 0.0) then
+        print*,'bad point 3 ',tile_num,i,landfrac(i),slmaskl(i),vegfcs(i)
+       endif
+       if(nint(slmaskl(i)) == 0 .and. vegfcs(i) > 0.0) then
+        print*,'bad point 4 ',tile_num,i,landfrac(i),slmaskw(i),vegfcs(i)
+       endif
+
+     ENDDO
+
+   ELSE
 ! for running uncoupled (non-fractional grid)
+     DO I=1,LENSFC
+       IF(NINT(SLMASK(I)) == 1) THEN
+         SLMASKL(I) = 1.0_KIND_io8
+         SLMASKW(I) = 1.0_KIND_io8
+       ELSE
+         SLMASKL(I) = 0.0_KIND_io8
+         SLMASKW(I) = 0.0_KIND_io8
+       ENDIF
+     ENDDO  
+
+   ENDIF ! frac_grid
+
    DO I=1,LENSFC
-     IF(NINT(SLMASK(I)) == 1) THEN
-       SLMASKL(I) = 1.0_KIND_io8
-       SLMASKW(I) = 1.0_KIND_io8
-     ELSE
-       SLMASKL(I) = 0.0_KIND_io8
-       SLMASKW(I) = 0.0_KIND_io8
-     ENDIF
      if(nint(slmask(i)) == 0) then
        min_ice(i) = 0.15_KIND_io8
      else
        min_ice(i) = 0.0_KIND_io8
      endif
-   ENDDO  
+   ENDDO
+
    num_threads = num_parthds()
    PRINT*
    PRINT*,"CALL SFCCYCLE TO UPDATE SURFACE FIELDS."
@@ -562,9 +631,44 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC,  &
                SZ_NML, INPUT_NML_FILE,                   &
                min_ice, &
                IALB,ISOT,IVEGSRC,TILE_NUM,I_INDEX,J_INDEX)
+
+!  print*,'after call to sfccycle '
+
+!  print*,'chk point after sfccycle ',tile_num,landfrac(ijchk),slmask(ijchk),vetfcs(ijchk),vegfcs(ijchk)
+
+    goto 20
+
+     DO I=1,LENSFC
+       if(nint(slmaskl(i)) == 1 .and. vegfcs(i) == 0.0) then
+          if (nint(vetfcs(i)) /= 15) then
+        print*,'bad point 5 ',tile_num,i,landfrac(i),slmaskl(i),vegfcs(i)
+          endif
+       endif
+        print*,'in loop after first if'
+       if(nint(slmaskl(i)) == 0 .and. vegfcs(i) > 0.0) then
+        print*,'bad point 6 ',tile_num,i,landfrac(i),slmaskw(i),vegfcs(i)
+       endif
+       if(landfrac(i) > 0.0_kind_io8 .and. vegfcs(i) == 0.0) then
+          if (nint(vetfcs(i)) /= 15) then
+        print*,'bad point 7 ',tile_num,i,landfrac(i),vegfcs(i)
+          endif
+       endif
+       if(landfrac(i) == 0.0_kind_io8 .and. vegfcs(i) > 0.0) then
+        print*,'bad point 8 ',tile_num,i,landfrac(i),vegfcs(i)
+       endif
+       if(nint(slmaskw(i)) == 1 .and. sicfcs(i) > 0.0) then
+        print*,'bad ice point ',tile_num,i,landfrac(i),sicfcs(i)
+       endif
+
+     ENDDO
+
+ 20 continue
+
    DEALLOCATE(SLMASKL, SLMASKW)
  ENDIF
 
+ IF(FRAC_GRID) DEALLOCATE(LANDFRAC)
+
 !--------------------------------------------------------------------------------
 ! IF RUNNING WITH NSST, READ IN GSI FILE WITH THE UPDATED INCREMENTS (ON THE
 ! GAUSSIAN GRID), INTERPOLATE INCREMENTS TO THE CUBED-SPHERE TILE, AND PERFORM

sorc/global_cycle.fd/read_write_data.f90

@@ -677,9 +677,10 @@ end subroutine remove_checksum
  !! @param[out] TILE_NUM Cubed-sphere tile number
  !! @author George Gayno NOAA/EMC
  SUBROUTINE READ_LAT_LON_OROG(RLA,RLO,OROG,OROG_UF,&
-                              TILE_NUM,IDIM,JDIM,IJDIM)
+                              TILE_NUM,IDIM,JDIM,IJDIM,LANDFRAC)
 
  USE MPI
+ USE MACHINE
 
  IMPLICIT NONE
 
@@ -689,6 +690,7 @@ SUBROUTINE READ_LAT_LON_OROG(RLA,RLO,OROG,OROG_UF,&
 
  REAL, INTENT(OUT)      :: RLA(IJDIM),RLO(IJDIM)
  REAL, INTENT(OUT)      :: OROG(IJDIM),OROG_UF(IJDIM)
+ REAL(KIND=KIND_IO8), INTENT(OUT), OPTIONAL  :: LANDFRAC(IJDIM)
 
  CHARACTER(LEN=50)      :: FNOROG, FNGRID
  CHARACTER(LEN=3)       :: RANKCH
@@ -798,6 +800,14 @@ SUBROUTINE READ_LAT_LON_OROG(RLA,RLO,OROG,OROG_UF,&
  CALL NETCDF_ERR(ERROR, 'ERROR READING orog_filt RECORD' )
  OROG = RESHAPE(DUMMY4, (/IJDIM/))
 
+ IF(PRESENT(LANDFRAC))THEN
+   ERROR=NF90_INQ_VARID(NCID_OROG, 'land_frac', ID_VAR)
+   CALL NETCDF_ERR(ERROR, 'ERROR READING land_frac ID' )
+   ERROR=NF90_GET_VAR(NCID_OROG, ID_VAR, DUMMY4)
+   CALL NETCDF_ERR(ERROR, 'ERROR READING land_frac RECORD' )
+   LANDFRAC = RESHAPE(DUMMY4, (/IJDIM/))
+ ENDIF
+
  DEALLOCATE(DUMMY4)
 
  ERROR = NF90_CLOSE(NCID_OROG)

ush/global_cycle.sh

@@ -272,6 +272,7 @@ DO_SNO_INC=${DO_SNO_INC:-.false.}
 zsea1=${zsea1:-0}
 zsea2=${zsea2:-0}
 MAX_TASKS_CY=${MAX_TASKS_CY:-99999}
+FRAC_GRID=${FRAC_GRID:-.false.}
 
 FNGLAC=${FNGLAC:-${FIXam}/global_glacier.2x2.grb}
 FNMXIC=${FNMXIC:-${FIXam}/global_maxice.2x2.grb}
@@ -384,7 +385,8 @@ cat << EOF > fort.36
   iy=$iy, im=$im, id=$id, ih=$ih, fh=$FHOUR,
   deltsfc=$DELTSFC,ialb=$IALB,use_ufo=$use_ufo,donst=$DONST,
   do_sfccycle=$DO_SFCCYCLE,do_lndinc=$DO_LNDINC,isot=$ISOT,ivegsrc=$IVEGSRC,
-  zsea1_mm=$zsea1,zsea2_mm=$zsea2,MAX_TASKS=$MAX_TASKS_CY
+  zsea1_mm=$zsea1,zsea2_mm=$zsea2,MAX_TASKS=$MAX_TASKS_CY,
+  frac_grid=$FRAC_GRID
  /
 EOF
 

ush/global_cycle_driver.sh

@@ -53,6 +53,7 @@ export DO_SFCCYLE=${DO_SFCCYCLE:-".true."}
 export DO_LNDINC=${DO_LNDINC:-".false."}
 export LND_SOI_FILE=${LND_SOI_FILE:-"NULL"}
 export DO_SNO_INC=${DO_SNO_INC:-".false."}
+export FRAC_GRID=${FRAC_GRID:-".false."}
 
 CRES=$(echo $CASE | cut -c 2-)
 JCAP_CASE=$((2*CRES-2))
@@ -70,11 +71,16 @@ export DATA=${DATA:-$pwd/rundir$$}
 rm -fr $DATA
 mkdir -p $DATA
 
+# Make a copy of the input restart file as fnbgso. For fractional grids,
+# only a few fields will be updated by sfcsub.F. Only these records
+# will be overwritten in fnbgso.
+
 for n in $(seq 1 $ntiles); do
   ln -fs $COMIN/$PDY.${cyc}0000.sfc_data.tile${n}.nc      $DATA/fnbgsi.00$n
 
 # Make a copy of the input restart file in the working directory.
 # global_cycle will update the required records for noah-mp.
+
   cp $COMIN/$PDY.${cyc}0000.sfc_data.tile${n}.nc $COMOUT/$PDY.${cyc}0000.sfcanl_data.tile${n}.nc
   chmod 644  $COMOUT/$PDY.${cyc}0000.sfcanl_data.tile${n}.nc
   ln -fs $COMOUT/$PDY.${cyc}0000.sfcanl_data.tile${n}.nc  $DATA/fnbgso.00$n