dtc/develop: final PR before merging to develop or master (fractional landmask changes, support for hera.gnu)

atmos_model.F90

@@ -1724,7 +1724,6 @@ subroutine assign_importdata(rc)
                   if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then
                     IPD_Data(nb)%Coupling%tseain_cpl(ix) = datar8(i,j)
                     IPD_Data(nb)%Sfcprop%tsfco(ix)       = datar8(i,j)
-!                   IPD_Data(nb)%Sfcprop%tsfc(ix)        = datar8(i,j)
                   endif
                 enddo
               enddo
@@ -1747,13 +1746,17 @@ subroutine assign_importdata(rc)
                   IPD_Data(nb)%Coupling%ficein_cpl(ix)   = zero
                   IPD_Data(nb)%Coupling%slimskin_cpl(ix) = IPD_Data(nb)%Sfcprop%slmsk(ix)
                   if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then
-                    if (datar8(i,j) >= IPD_control%min_seaice*IPD_Data(nb)%Sfcprop%oceanfrac(ix)) then
-                      IPD_Data(nb)%Coupling%ficein_cpl(ix) = max(zero, min(datar8(i,j),one))
-                      IPD_Data(nb)%Sfcprop%slmsk(ix)         = 2.                                        !slmsk=2 crashes in gcycle on partial land points
+                    IPD_Data(nb)%Coupling%ficein_cpl(ix) = max(zero, min(one, datar8(i,j)/IPD_Data(nb)%Sfcprop%oceanfrac(ix))) !LHS: ice frac wrt water area
+                    if (IPD_Data(nb)%Coupling%ficein_cpl(ix) > one-epsln) IPD_Data(nb)%Coupling%ficein_cpl(ix)=one
+                    if (IPD_Data(nb)%Coupling%ficein_cpl(ix) >= IPD_control%min_seaice) then
+                      if (abs(one-IPD_Data(nb)%Sfcprop%oceanfrac(ix)) < epsln) IPD_Data(nb)%Sfcprop%slmsk(ix) = 2. !slmsk=2 crashes in gcycle on partial land points
                       IPD_Data(nb)%Coupling%slimskin_cpl(ix) = 4.
-                    elseif (abs(one-IPD_Data(nb)%Sfcprop%oceanfrac(ix)) < epsln) then
-                      IPD_Data(nb)%Sfcprop%slmsk(ix)         = zero
-                      IPD_Data(nb)%Coupling%slimskin_cpl(ix) = zero
+                    else
+                      IPD_Data(nb)%Coupling%ficein_cpl(ix)   = zero
+                      if (abs(one-IPD_Data(nb)%Sfcprop%oceanfrac(ix)) < epsln) then
+                        IPD_Data(nb)%Sfcprop%slmsk(ix)         = zero
+                        IPD_Data(nb)%Coupling%slimskin_cpl(ix) = zero
+                      end if
                     endif
                   endif
                 enddo
@@ -1924,7 +1927,7 @@ subroutine assign_importdata(rc)
           ix = Atm_block%ixp(i,j)
           if (IPD_Data(nb)%Sfcprop%oceanfrac(ix) > zero) then
 !if it is ocean or ice get surface temperature from mediator
-            if(IPD_Data(nb)%Coupling%ficein_cpl(ix) >= IPD_control%min_seaice*IPD_Data(nb)%Sfcprop%oceanfrac(ix)) then
+            if(IPD_Data(nb)%Coupling%ficein_cpl(ix) >= IPD_control%min_seaice) then
               IPD_Data(nb)%Sfcprop%tisfc(ix) = IPD_Data(nb)%Coupling%tisfcin_cpl(ix)
               IPD_Data(nb)%Sfcprop%fice(ix)  = IPD_Data(nb)%Coupling%ficein_cpl(ix)
               IPD_Data(nb)%Sfcprop%hice(ix)  = IPD_Data(nb)%Coupling%hicein_cpl(ix)

ccpp/build_ccpp.sh

@@ -5,7 +5,7 @@ set -eu
 
 # List of valid/tested machines
 VALID_MACHINES=( wcoss_cray wcoss_dell_p3 gaea.intel jet.intel \
-                 hera.intel \
+                 hera.intel hera.gnu \
                  cheyenne.intel cheyenne.intel-impi cheyenne.gnu cheyenne.pgi endeavor.intel \
                  stampede.intel supermuc_phase2.intel macosx.gnu \
                  linux.intel linux.gnu linux.pgi )
@@ -128,6 +128,12 @@ fi
 if [[ "${MAKE_OPT}" == *"STATIC=Y"* ]]; then
   CCPP_CMAKE_FLAGS="${CCPP_CMAKE_FLAGS} -DSTATIC=ON"
 else
+  # hera.gnu uses the NCEPLIBS-external/NCEPLIBS umbrella build libraries,
+  # which cannot be linked dynamically at this point (missing -fPIC flag)
+  if [[ "${MACHINE_ID}" == "hera.gnu" ]]; then
+    echo "Dynamic CCPP build not supported on hera.gnu at this time."
+    exit 1
+  fi
   # Dynamic builds require linking the NCEPlibs, provide path to them
   CCPP_CMAKE_FLAGS="${CCPP_CMAKE_FLAGS} -DSTATIC=OFF -DBACIO_LIB4=${BACIO_LIB4} -DSP_LIBd=${SP_LIBd} -DW3NCO_LIBd=${W3NCO_LIBd}"
 fi

ccpp/set_compilers.sh

@@ -39,6 +39,13 @@ case "$MACHINE_ID" in
             export F77=mpiifort
             export F90=mpiifort
             ;;
+        hera.gnu)
+            export CC=mpicc
+            export CXX=mpicxx
+            export FC=mpif90
+            export F77=mpif77
+            export F90=mpif90
+            ;;
         cheyenne.intel)
             export CC=mpicc
             export CXX=mpicxx

gfsphysics/GFS_layer/GFS_physics_driver.F90

@@ -1069,7 +1069,7 @@ subroutine GFS_physics_driver                                 &
 !  --- ...  xw: transfer ice thickness & concentration from global to local variables
 !## CCPP ## global to local variable transfer not necessary for these two
         zice(i)    = Sfcprop%hice(i)
-        fice(i)    = Sfcprop%fice(i) ! ice fraction of lake/ocean wrt whole cell
+        fice(i)    = Sfcprop%fice(i)
 !*## CCPP ##* 
 !## CCPP ##* GFS_surface_composites.F90/GFS_surface_composites_pre_run
         tice(i)    = Sfcprop%tisfc(i)
@@ -1126,35 +1126,31 @@ subroutine GFS_physics_driver                                 &
 !*## CCPP ##
 
 !## CCPP ##* GFS_surface_composites.F90/GFS_surface_composites_pre
-      if (Model%frac_grid) then  ! here Sfcprop%fice is fraction of the whole grid that is ice
+      if (Model%frac_grid) then
         do i = 1, IM
           frland(i) = Sfcprop%landfrac(i)
           if (frland(i) > zero) dry(i) = .true.
-          tem = one - frland(i)
-          if (tem > epsln) then
+          if (frland(i) < one) then
             if (flag_cice(i)) then
-              if (fice(i) >= Model%min_seaice*tem) then
+              if (fice(i) >= Model%min_seaice) then
                 icy(i)  = .true.
               else
                 fice(i) = zero
               endif
             else
-              if (fice(i) >= Model%min_lakeice*tem) then
+              if (fice(i) >= Model%min_lakeice) then
                 icy(i) = .true.
-                fice(i) = fice(i)/tem  ! fice is fraction of ocean/lake
               else
                 fice(i) = zero
               endif
             endif
+            if (fice(i) < one) then
+              wet(i)=.true. !some open ocean/lake water exists
+              if (.not. Model%cplflx) Sfcprop%tsfco(i) = max(Sfcprop%tsfco(i), Sfcprop%tisfc(i), tgice)
+
+            end if
           else
             fice(i) = zero
-          endif
-                                        ! ocean/lake area that is not frozen
-          if (tem-fice(i) > epsln) then
-            wet(i) = .true.             ! there is some open water!
-            if (.not. Model%cplflx) Sfcprop%tsfco(i) = max(Sfcprop%tsfco(i), Sfcprop%tisfc(i), tgice)
-!           if (icy(i)) Sfcprop%tsfco(i) = max(Sfcprop%tsfco(i), tgice)
-!           if (icy(i)) Sfcprop%tsfco(i) = max(Sfcprop%tisfc(i), tgice)
           endif
         enddo
       else
@@ -1524,32 +1520,32 @@ subroutine GFS_physics_driver                                 &
 
         if (Model%frac_grid) then
           do i=1,im
-            tem = one - Sfcprop%fice(i) - frland(i)
+            tem = (one - frland(i)) * fice(i) ! tem = ice fraction wrt whole cell
             if (flag_cice(i)) then
               adjsfculw(i) = adjsfculw3(i,1) * frland(i)                 &
-                           + Coupling%ulwsfcin_cpl(i) * Sfcprop%fice(i)  &
-                           + adjsfculw3(i,3) * tem
+                           + Coupling%ulwsfcin_cpl(i) * tem	         &
+                           + adjsfculw3(i,3) * (one - frland(i) - tem)
             else
               adjsfculw(i) = adjsfculw3(i,1) * frland(i)                 &
-                           + adjsfculw3(i,2) * Sfcprop%fice(i)           &
-                           + adjsfculw3(i,3) * tem
+                           + adjsfculw3(i,2) * tem                       &
+                           + adjsfculw3(i,3) * (one - frland(i) - tem)
             endif
           enddo
         else
           do i=1,im
             if (dry(i)) then                     ! all land
               adjsfculw(i) = adjsfculw3(i,1)
             elseif (icy(i)) then                 ! ice (and water)
-              tem = one - Sfcprop%fice(i)
+              tem = one - fice(i)
               if (flag_cice(i)) then
                 if (wet(i) .and. adjsfculw3(i,3) /= huge) then
-                  adjsfculw(i) = Coupling%ulwsfcin_cpl(i)*Sfcprop%fice(i) + adjsfculw3(i,3)*tem
+                  adjsfculw(i) = Coupling%ulwsfcin_cpl(i)*fice(i) + adjsfculw3(i,3)*tem
                 else
                   adjsfculw(i) = Coupling%ulwsfcin_cpl(i)
                 endif
               else
                 if (wet(i) .and. adjsfculw3(i,3) /= huge) then
-                  adjsfculw(i) = adjsfculw3(i,2)*Sfcprop%fice(i) + adjsfculw3(i,3)*tem
+                  adjsfculw(i) = adjsfculw3(i,2)*fice(i) + adjsfculw3(i,3)*tem
                 else
                   adjsfculw(i) = adjsfculw3(i,2)
                 endif
@@ -2002,9 +1998,9 @@ subroutine GFS_physics_driver                                 &
         do i=1, im
 !
 ! Three-way composites (fields from sfc_diff)
-          txl = Sfcprop%landfrac(i)
-          txi = Sfcprop%fice(i)                 ! here Sfcprop%fice is grid fraction that is ice
-          txo = one - txl - txi
+          txl = frland(i)
+          txi = fice(i)*(one - frland(i)) ! txi = ice fraction wrt whole cell
+          txo = max(zero, one - txl - txi)
           Sfcprop%zorl(i)   = txl*zorl3(i,1)   + txi*zorl3(i,2)   + txo*zorl3(i,3)
           cd(i)             = txl*cd3(i,1)     + txi*cd3(i,2)     + txo*cd3(i,3)
           cdq(i)            = txl*cdq3(i,1)    + txi*cdq3(i,2)    + txo*cdq3(i,3)
@@ -2062,8 +2058,7 @@ subroutine GFS_physics_driver                                 &
           if (.not. flag_cice(i)) then
             if (islmsk(i) == 2) then                           ! return updated lake ice thickness & concentration to global array
               Sfcprop%hice(i)  = zice(i)
-!             Sfcprop%fice(i)  = fice(i) * Sfcprop%lakefrac(i) ! fice is fraction of lake area that is frozen
-              Sfcprop%fice(i)  = fice(i) * (one-Sfcprop%landfrac(i)) ! fice is fraction of wet area that is frozen
+              Sfcprop%fice(i)  = fice(i) 
               Sfcprop%tisfc(i) = tice(i)
             else                                               ! this would be over open ocean or land (no ice fraction)
               Sfcprop%hice(i)  = zero
@@ -2121,7 +2116,7 @@ subroutine GFS_physics_driver                                 &
           Sfcprop%zorlo(i)  = zorl3(i,3)
 
           if (flag_cice(i) .and. wet(i)) then    ! this was already done for lake ice in sfc_sice
-            txi = Sfcprop%fice(i)
+            txi = fice(i)
             txo = one - txi
             evap(i)         = txi * evap3(i,2) + txo * evap3(i,3)
             hflx(i)         = txi * hflx3(i,2) + txo * hflx3(i,3)
@@ -2843,34 +2838,31 @@ subroutine GFS_physics_driver                                 &
       if (Model%cplflx) then
         do i=1,im
           if (Sfcprop%oceanfrac(i) > zero) then               ! Ocean only, NO LAKES
-            if (fice(i) == Sfcprop%oceanfrac(i)) then ! use results from CICE
+            if (Sfcprop%fice(i) > one - epsln) then ! no open water, thus use results from CICE
               Coupling%dusfci_cpl(i) = dusfc_cice(i)
               Coupling%dvsfci_cpl(i) = dvsfc_cice(i)
               Coupling%dtsfci_cpl(i) = dtsfc_cice(i)
               Coupling%dqsfci_cpl(i) = dqsfc_cice(i)
-
-            elseif (wet(i)) then              ! use stress_ocean from sfc_diff for opw component at mixed point
-              if (icy(i) .or. dry(i)) then
-                tem1 = max(Diag%q1(i), 1.e-8)
-                rho = Statein%prsl(i,1) / (con_rd*Diag%t1(i)*(one+con_fvirt*tem1))
-                if (wind(i) > zero) then
-                  tem = - rho * stress3(i,3) / wind(i)
-                  Coupling%dusfci_cpl(i) = tem * Statein%ugrs(i,1)   ! U-momentum flux
-                  Coupling%dvsfci_cpl(i) = tem * Statein%vgrs(i,1)   ! V-momentum flux
-                else
-                  Coupling%dusfci_cpl(i) = zero
-                  Coupling%dvsfci_cpl(i) = zero
-                endif
-                Coupling%dtsfci_cpl(i) = con_cp   * rho * hflx3(i,3) ! sensible heat flux over open ocean
-                Coupling%dqsfci_cpl(i) = con_hvap * rho * evap3(i,3) ! latent heat flux over open ocean
+            elseif (icy(i) .or. dry(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point
+              tem1 = max(Diag%q1(i), 1.e-8)
+              rho = Statein%prsl(i,1) / (con_rd*Diag%t1(i)*(one+con_fvirt*tem1))
+              if (wind(i) > zero) then
+                tem = - rho * stress3(i,3) / wind(i)
+                Coupling%dusfci_cpl(i) = tem * Statein%ugrs(i,1)   ! U-momentum flux
+                Coupling%dvsfci_cpl(i) = tem * Statein%vgrs(i,1)   ! V-momentum flux
+              else
+                Coupling%dusfci_cpl(i) = zero
+                Coupling%dvsfci_cpl(i) = zero
+              endif
+              Coupling%dtsfci_cpl(i) = con_cp   * rho * hflx3(i,3) ! sensible heat flux over open ocean
+              Coupling%dqsfci_cpl(i) = con_hvap * rho * evap3(i,3) ! latent heat flux over open ocean
 !     if (lprnt .and. i == ipr) write(0,*)' hflx33=',hflx3(i,3),' evap33=',evap3(i,3), &
 !    ' con_cp=',con_cp,' rho=',rho,' con_hvap=',con_hvap
-              else                                                   ! use results from PBL scheme for 100% open ocean
-                Coupling%dusfci_cpl(i) = dusfc1(i)
-                Coupling%dvsfci_cpl(i) = dvsfc1(i)
-                Coupling%dtsfci_cpl(i) = dtsfc1(i)
-                Coupling%dqsfci_cpl(i) = dqsfc1(i)
-              endif
+            else                                                   ! use results from PBL scheme for 100% open ocean
+              Coupling%dusfci_cpl(i) = dusfc1(i)
+              Coupling%dvsfci_cpl(i) = dvsfc1(i)
+              Coupling%dtsfci_cpl(i) = dtsfc1(i)
+              Coupling%dqsfci_cpl(i) = dqsfc1(i)
             endif
 
             Coupling%dusfc_cpl (i) = Coupling%dusfc_cpl(i) + Coupling%dusfci_cpl(i) * dtf

gfsphysics/GFS_layer/GFS_typedefs.F90

@@ -6528,7 +6528,7 @@ subroutine interstitial_phys_reset (Interstitial, Model)
     ! Set flag for resetting maximum hourly output fields
     Interstitial%reset = mod(Model%kdt-1, nint(Model%avg_max_length/Model%dtp)) == 0
     ! Set flag for resetting radar reflectivity calculation
-    if (Interstitial%radar_reset<0) then
+    if (Model%nsradar_reset<0) then
       Interstitial%radar_reset = .true.
     else
       Interstitial%radar_reset = mod(Model%kdt-1, nint(Model%nsradar_reset/Model%dtp)) == 0