Land surface upgrades for HR2 (#663)

* adding soil color data to NoahMP

* add iopt_diag for 2m t/q diagnostic option

* fixed a sfc_name2 error

* fix a counting error in Sfc_io_calculate_indices

ccpp/data/GFS_typedefs.F90

@@ -283,8 +283,10 @@ module GFS_typedefs
     real (kind=kind_phys), pointer :: vfrac  (:)   => null()  !< vegetation fraction
     integer,               pointer :: vtype  (:)   => null()  !< vegetation type
     integer,               pointer :: stype  (:)   => null()  !< soil type
+    integer,               pointer :: scolor  (:)   => null()  !< soil color
     integer,               pointer :: vtype_save (:) => null()!< vegetation type save
     integer,               pointer :: stype_save (:) => null()!< soil type save
+    integer,               pointer :: scolor_save (:) => null()!< soil color save
     real (kind=kind_phys), pointer :: uustar (:)   => null()  !< boundary layer parameter
     real (kind=kind_phys), pointer :: oro    (:)   => null()  !< orography
     real (kind=kind_phys), pointer :: oro_uf (:)   => null()  !< unfiltered orography
@@ -1045,6 +1047,8 @@ module GFS_typedefs
     integer              :: iopt_tbot !lower boundary of soil temperature (1->zero-flux; 2->noah)
     integer              :: iopt_stc  !snow/soil temperature time scheme (only layer 1)
     integer              :: iopt_trs  !thermal roughness scheme (1-z0h=z0m; 2-czil; 3-ec;4-kb inversed)
+    integer              :: iopt_diag !2m t/q diagnostic approach (1->external GFS sfc_diag 2->original NoahMP 2-title 3->NoahMP 
+                                      !2-title + internal GFS sfc_diag  )
 
     ! -- RUC LSM options
     integer              :: mosaic_lu=0     !< control for use of fractional landuse in RUC land surface model
@@ -2379,6 +2383,8 @@ subroutine sfcprop_create (Sfcprop, IM, Model)
     allocate (Sfcprop%vtype_save (IM))
     allocate (Sfcprop%stype      (IM))
     allocate (Sfcprop%stype_save (IM))
+    allocate (Sfcprop%scolor     (IM))
+    allocate (Sfcprop%scolor_save(IM))
     allocate (Sfcprop%uustar     (IM))
     allocate (Sfcprop%oro        (IM))
     allocate (Sfcprop%oro_uf     (IM))
@@ -2397,6 +2403,8 @@ subroutine sfcprop_create (Sfcprop, IM, Model)
     Sfcprop%vtype_save = zero
     Sfcprop%stype      = zero
     Sfcprop%stype_save = zero
+    Sfcprop%scolor      = zero
+    Sfcprop%scolor_save = zero
     Sfcprop%uustar     = clear_val
     Sfcprop%oro        = clear_val
     Sfcprop%oro_uf     = clear_val
@@ -3493,6 +3501,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     integer              :: iopt_tbot      =  2  !lower boundary of soil temperature (1->zero-flux; 2->noah)
     integer              :: iopt_stc       =  1  !snow/soil temperature time scheme (only layer 1)
     integer              :: iopt_trs       =  2  !thermal roughness scheme (1-z0h=z0m; 2-czil; 3-ec;4-kb reversed)
+    integer              :: iopt_diag      =  2  !2m t/q diagnostic approach (1->external GFS sfc_diag 2->original NoahMP 2-title
+                                                 !3->NoahMP 2-title + internal GFS sfc_diag  )
 
     integer              :: mosaic_lu      =  0  ! 1 - used of fractional landuse in RUC lsm
     integer              :: mosaic_soil    =  0  ! 1 - used of fractional soil in RUC lsm
@@ -3897,7 +3907,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
                           !    Noah MP options
                                iopt_dveg,iopt_crs,iopt_btr,iopt_run,iopt_sfc, iopt_frz,     &
                                iopt_inf, iopt_rad,iopt_alb,iopt_snf,iopt_tbot,iopt_stc,     &
-                               iopt_trs,                                                    &
+                               iopt_trs, iopt_diag,                                         &
                           !    RUC lsm options
                                mosaic_lu, mosaic_soil, isncond_opt, isncovr_opt,            &
                           !    GFDL surface layer options
@@ -4659,6 +4669,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     Model%iopt_tbot        = iopt_tbot
     Model%iopt_stc         = iopt_stc
     Model%iopt_trs         = iopt_trs
+    Model%iopt_diag        = iopt_diag
 
 ! RUC lsm options
     Model%mosaic_lu        = mosaic_lu
@@ -5597,6 +5608,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
         print *,'iopt_tbot   =  ',Model%iopt_tbot
         print *,'iopt_stc   =  ', Model%iopt_stc
         print *,'iopt_trs   =  ', Model%iopt_trs
+        print *,'iopt_diag  =  ', Model%iopt_diag
       elseif (Model%lsm == Model%lsm_ruc) then
         print *,' RUC Land Surface Model used'
         print *, 'The Physics options are'
@@ -6513,6 +6525,7 @@ subroutine control_print(Model)
         print *, ' iopt_tbot         : ', Model%iopt_tbot
         print *, ' iopt_stc          : ', Model%iopt_stc
         print *, ' iopt_trs          : ', Model%iopt_trs
+        print *, ' iopt_diag         : ', Model%iopt_diag
       elseif (Model%lsm == Model%lsm_ruc) then
         print *,' RUC Land Surface Model used'
         print *, 'The Physics options are'

ccpp/data/GFS_typedefs.meta

@@ -1053,6 +1053,18 @@
   units = index
   dimensions = (horizontal_loop_extent)
   type = integer
+[scolor]
+  standard_name = soil_color_classification
+  long_name = soil color for lsm
+  units = index
+  dimensions = (horizontal_loop_extent)
+  type = integer
+[scolor_save]
+  standard_name = soil_color_classification_save
+  long_name = soil color for lsm save
+  units = index
+  dimensions = (horizontal_loop_extent)
+  type = integer
 [uustar]
   standard_name = surface_friction_velocity
   long_name = boundary layer parameter
@@ -4905,6 +4917,12 @@
   units = index
   dimensions = ()
   type = integer
+[iopt_diag]
+  standard_name = control_for_land_surface_scheme_surface_diagnose_approach
+  long_name = choice for surface diagnose approach option (see noahmp module for definition)
+  units = index
+  dimensions = ()
+  type = integer
 [use_ufo]
   standard_name = flag_for_gcycle_surface_option
   long_name = flag for gcycle surface option

ccpp/driver/GFS_diagnostics.F90

@@ -3838,6 +3838,18 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
       ExtDiag(idx)%data(nb)%int2 => Sfcprop(nb)%stype(:)
     enddo
 
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'scolor'
+    ExtDiag(idx)%desc = 'soil color in integer 1-20' 
+    ExtDiag(idx)%unit = 'number'
+    ExtDiag(idx)%mod_name = 'gfs_sfc'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%int2 => Sfcprop(nb)%scolor(:)
+    enddo
+
+
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'lfrac'

io/fv3atm_restart_io.F90

@@ -162,7 +162,7 @@ subroutine fv3atm_checksum (Model, GFS_Data, Atm_block)
 
     ntr = size(GFS_Data(1)%Statein%qgrs,3)
 
-    nsfcprop2d = 93
+    nsfcprop2d = 94
     if (Model%lsm == Model%lsm_noahmp) then
       nsfcprop2d = nsfcprop2d + 49
       if (Model%use_cice_alb) then
@@ -228,6 +228,7 @@ subroutine fv3atm_checksum (Model, GFS_Data, Atm_block)
       call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,temp2d,GFS_Data(nb)%Sfcprop%vfrac)
       call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,temp2d,GFS_Data(nb)%Sfcprop%vtype)
       call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,temp2d,GFS_Data(nb)%Sfcprop%stype)
+      call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,temp2d,GFS_Data(nb)%Sfcprop%scolor)
       call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,temp2d,GFS_Data(nb)%Sfcprop%uustar)
       call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,temp2d,GFS_Data(nb)%Sfcprop%oro)
       call copy_from_GFS_Data(ii1,jj1,isc,jsc,nt,temp2d,GFS_Data(nb)%Sfcprop%oro_uf)

io/fv3atm_sfc_io.F90

@@ -107,7 +107,7 @@ function Sfc_io_calculate_indices(sfc, Model, reading, warm_start)
     integer :: nvar2m, nvar2o, nvar3, nvar2r, nvar2mp, nvar3mp, nvar2l
     integer :: nvar_before_lake
 
-    nvar2m = 48
+    nvar2m = 49
     if (Model%use_cice_alb .or. Model%lsm == Model%lsm_ruc) then
       nvar2m = nvar2m + 4
       !nvar2m = nvar2m + 5
@@ -438,6 +438,7 @@ subroutine Sfc_io_fill_2d_names(sfc,Model,warm_start)
     nt=nt+1 ; sfc%name2(nt) = 'slope'
     nt=nt+1 ; sfc%name2(nt) = 'snoalb'
     !--- variables below here are optional
+    nt=nt+1 ; sfc%name2(nt) = 'scolor'
     nt=nt+1 ; sfc%name2(nt) = 'sncovr'
     nt=nt+1 ; sfc%name2(nt) = 'snodl' !snowd on land portion of a cell
     nt=nt+1 ; sfc%name2(nt) = 'weasdl'!weasd on land portion of a cell
@@ -593,7 +594,7 @@ subroutine Sfc_io_register_2d_fields(sfc,Model,Sfc_restart,reading,warm_start)
            .or. trim(sfc%name2(num)) == 'albdirvis_ice' .or. trim(sfc%name2(num)) == 'albdirnir_ice' &
            .or. trim(sfc%name2(num)) == 'albdifvis_ice' .or. trim(sfc%name2(num)) == 'albdifnir_ice' &
            .or. trim(sfc%name2(num)) == 'emis_lnd'      .or. trim(sfc%name2(num)) == 'emis_ice'      &
-           .or. trim(sfc%name2(num)) == 'sncovr_ice') then
+           .or. trim(sfc%name2(num)) == 'sncovr_ice'    .or. trim(sfc%name2(num)) == 'scolor') then
         if(reading .and. sfc%is_lsoil) then
           call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=(/'lat','lon'/), is_optional=.true.)
         else
@@ -829,6 +830,7 @@ subroutine Sfc_io_transfer(sfc, reading, Model, Atm_block, Sfcprop, warm_start,
       call GFS_Data_transfer(reading,ii1,jj1,isc,jsc,nt,sfc%var2,Sfcprop(nb)%shdmax)  !--- shdmax
       call GFS_Data_transfer(reading,ii1,jj1,isc,jsc,nt,sfc%var2,Sfcprop(nb)%slope)   !--- slope
       call GFS_Data_transfer(reading,ii1,jj1,isc,jsc,nt,sfc%var2,Sfcprop(nb)%snoalb)  !--- snoalb
+      call GFS_Data_transfer(reading,ii1,jj1,isc,jsc,nt,sfc%var2,Sfcprop(nb)%scolor)  !--- scolor
       call GFS_Data_transfer(reading,ii1,jj1,isc,jsc,nt,sfc%var2,Sfcprop(nb)%sncovr)  !--- sncovr
       call GFS_Data_transfer(reading,ii1,jj1,isc,jsc,nt,sfc%var2,Sfcprop(nb)%snodl)   !--- snodl (snowd on land  portion of a cell)
       call GFS_Data_transfer(reading,ii1,jj1,isc,jsc,nt,sfc%var2,Sfcprop(nb)%weasdl)  !--- weasdl (weasd on land  portion of a cell)
@@ -1277,7 +1279,21 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
     i = Atm_block%index(1)%ii(1) - isc + 1
     j = Atm_block%index(1)%jj(1) - jsc + 1
 
-    if (sfc%var2(i,j,33) < -9990.0_kind_phys) then
+    if (sfc%var2(i,j,32) < -9990.0_kind_phys) then
+      if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - set init soil color')
+      !$omp parallel do default(shared) private(nb, ix)
+      do nb = 1, Atm_block%nblks
+        do ix = 1, Atm_block%blksz(nb)
+          if ( nint (Sfcprop(nb)%slmsk(ix)) == 1 ) then  !including glacier
+            Sfcprop(nb)%scolor(ix)  = 4
+          else
+            Sfcprop(nb)%scolor(ix)  = zero
+          endif
+        enddo
+      enddo
+    endif
+
+    if (sfc%var2(i,j,34) < -9990.0_kind_phys) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing snodl')
       !$omp parallel do default(shared) private(nb, ix, tem)
       do nb = 1, Atm_block%nblks
@@ -1292,7 +1308,7 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
       enddo
     endif
 
-    if (sfc%var2(i,j,34) < -9990.0_kind_phys) then
+    if (sfc%var2(i,j,35) < -9990.0_kind_phys) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing weasdl')
       !$omp parallel do default(shared) private(nb, ix, tem)
       do nb = 1, Atm_block%nblks
@@ -1307,7 +1323,7 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
       enddo
     endif
 
-    if (sfc%var2(i,j,36) < -9990.0_kind_phys) then
+    if (sfc%var2(i,j,37) < -9990.0_kind_phys) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing tsfcl')
       !$omp parallel do default(shared) private(nb, ix)
       do nb = 1, Atm_block%nblks
@@ -1317,7 +1333,7 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
       enddo
     endif
 
-    if (sfc%var2(i,j,37) < -9990.0_kind_phys) then
+    if (sfc%var2(i,j,38) < -9990.0_kind_phys) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorlw')
       !$omp parallel do default(shared) private(nb, ix)
       do nb = 1, Atm_block%nblks
@@ -1329,7 +1345,7 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
       enddo
     endif
 
-    if (sfc%var2(i,j,38) < -9990.0_kind_phys) then
+    if (sfc%var2(i,j,39) < -9990.0_kind_phys) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorll')
       !$omp parallel do default(shared) private(nb, ix)
       do nb = 1, Atm_block%nblks
@@ -1339,7 +1355,7 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
       enddo
     endif
 
-    if (sfc%var2(i,j,39) < -9990.0_kind_phys) then
+    if (sfc%var2(i,j,40) < -9990.0_kind_phys) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorli')
       !$omp parallel do default(shared) private(nb, ix)
       do nb = 1, Atm_block%nblks
@@ -1351,7 +1367,7 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
       enddo
     endif
 
-    if (sfc%var2(i,j,45) < -9990.0_kind_phys) then
+    if (sfc%var2(i,j,46) < -9990.0_kind_phys) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing emis_ice')
       !$omp parallel do default(shared) private(nb, ix)
       do nb = 1, Atm_block%nblks
@@ -1361,7 +1377,7 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
       enddo
     endif
 
-    if (sfc%var2(i,j,46) < -9990.0_kind_phys .and. Model%lsm /= Model%lsm_ruc) then
+    if (sfc%var2(i,j,47) < -9990.0_kind_phys .and. Model%lsm /= Model%lsm_ruc) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing sncovr_ice')
       !$omp parallel do default(shared) private(nb, ix)
       do nb = 1, Atm_block%nblks
@@ -1372,7 +1388,7 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
       enddo
     endif
 
-    if (sfc%var2(i,j,47) < -9990.0_kind_phys) then
+    if (sfc%var2(i,j,48) < -9990.0_kind_phys) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing snodi')
       !$omp parallel do default(shared) private(nb, ix, tem)
       do nb = 1, Atm_block%nblks
@@ -1387,7 +1403,7 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
       enddo
     endif
 
-    if (sfc%var2(i,j,48) < -9990.0_kind_phys) then
+    if (sfc%var2(i,j,49) < -9990.0_kind_phys) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing weasdi')
       !$omp parallel do default(shared) private(nb, ix, tem)
       do nb = 1, Atm_block%nblks
@@ -1403,7 +1419,7 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
     endif
 
     if (Model%use_cice_alb) then
-      if (sfc%var2(i,j,49) < -9990.0_kind_phys) then
+      if (sfc%var2(i,j,50) < -9990.0_kind_phys) then
         !$omp parallel do default(shared) private(nb, ix)
         do nb = 1, Atm_block%nblks
           do ix = 1, Atm_block%blksz(nb)