Merge pull request NCAR#139 from SamuelTrahanNOAA/clm-unit-fix

clm lake fix: pressure is not density + NCAR#142 and NCAR#148

physics/GFS_surface_composites_pre.F90

@@ -241,8 +241,10 @@ subroutine GFS_surface_composites_pre_run (im, lkm, frac_grid, iopt_lake, iopt_l
         !mjz
           tsfcl(i) = huge
         endif
+        if (icy(i) .or. wet(i)) then ! init uustar_ice for all water/ice grids 
+           uustar_ice(i) = uustar(i)
+        endif
         if (icy(i)) then                   ! Ice
-          uustar_ice(i) = uustar(i)
           is_clm = lkm>0 .and. iopt_lake==iopt_lake_clm .and. use_lake_model(i)>0
           if(lsm /= lsm_ruc .and. .not.is_clm) then
             weasd_ice(i) = weasd(i)

physics/clm_lake.f90

@@ -238,7 +238,10 @@ subroutine calculate_z_dz_lake(i,input_lakedepth,clm_lakedepth,z_lake,dz_lake)
       real(kind_lake) :: dz_lake(nlevlake) ! layer thickness for lake (m)
       real(kind_lake) :: depthratio
 
-      if (input_lakedepth(i) == spval) then
+      if (input_lakedepth(i) == spval .or. input_lakedepth(i) < 0.1) then
+        ! This is a safeguard against:
+        ! 1. missing in the lakedepth database (== spval)
+        ! 2. errors in model cycling or unexpected changes in the orography database (< 0.1)
         clm_lakedepth(i) = zlak(nlevlake) + 0.5_kind_lake*dzlak(nlevlake)
         z_lake(1:nlevlake) = zlak(1:nlevlake)
         dz_lake(1:nlevlake) = dzlak(1:nlevlake)
@@ -267,8 +270,8 @@ SUBROUTINE clm_lake_run( &
 
          ! Atmospheric model state inputs:
          tg3, pgr, zlvl, gt0, prsi, phii, qvcurr, gu0, gv0, xlat_d, xlon_d,       &
-         ch, cm, dlwsfci, dswsfci, oro_lakedepth, wind, rho0, tsfc,               &
-         flag_iter, ISLTYP, rainncprv, raincprv,                                  &
+         ch, cm, dlwsfci, dswsfci, oro_lakedepth, wind, tsfc,                     &
+         flag_iter, flag_lakefreeze, ISLTYP, rainncprv, raincprv,                 &
 
          ! Feedback to atmosphere:
          evap_wat,     evap_ice,   hflx_wat,    hflx_ice,  gflx_wat, gflx_ice,    &
@@ -283,7 +286,7 @@ SUBROUTINE clm_lake_run( &
 
          salty, savedtke12d, snowdp2d, h2osno2d, snl2d, t_grnd2d, t_lake3d,       &
          lake_icefrac3d, t_soisno3d, h2osoi_ice3d, h2osoi_liq3d, h2osoi_vol3d,    &
-         z3d, dz3d, zi3d,                                                         &
+         z3d, dz3d, zi3d, t1, qv1, prsl1,                                         &
                    input_lakedepth, clm_lakedepth, cannot_freeze,                 &
 
          ! Error reporting:
@@ -321,10 +324,12 @@ SUBROUTINE clm_lake_run( &
     !
     REAL(KIND_PHYS), DIMENSION(:), INTENT(IN):: &
          tg3, pgr, zlvl, qvcurr, xlat_d, xlon_d, ch, cm, &
-         dlwsfci, dswsfci, oro_lakedepth, wind, rho0, &
-         rainncprv, raincprv
+         dlwsfci, dswsfci, oro_lakedepth, wind, &
+         rainncprv, raincprv, t1, qv1, prsl1
     REAL(KIND_PHYS), DIMENSION(:,:), INTENT(in) :: gu0, gv0, prsi, gt0, phii
     LOGICAL, DIMENSION(:), INTENT(IN) :: flag_iter
+    LOGICAL, DIMENSION(:), INTENT(INOUT) :: flag_lakefreeze
+
     INTEGER, DIMENSION(:), INTENT(IN) :: ISLTYP
 
     !
@@ -450,6 +455,7 @@ SUBROUTINE clm_lake_run( &
       logical, parameter :: feedback_to_atmosphere = .true. ! FIXME: REMOVE
 
       real(kind_lake) :: to_radians, lat_d, lon_d, qss, tkm, bd
+      real(kind_lake) :: rho0                    ! lowest model level air density
 
       integer :: month,num1,num2,day_of_month,isl
       real(kind_lake) :: wght1,wght2,Tclim,depthratio
@@ -693,12 +699,13 @@ SUBROUTINE clm_lake_run( &
 
                 !-- The CLM output is combined for fractional ice and water
                 if( t_grnd(c) >= tfrz ) then
-                  qfx         = eflx_lh_tot(c)*invhvap
+                  qfx           = eflx_lh_tot(c)*invhvap
                 else
-                  qfx         = eflx_lh_tot(c)*invhsub      ! heat flux (W/m^2)=>mass flux(kg/(sm^2))
+                  qfx           = eflx_lh_tot(c)*invhsub      ! heat flux (W/m^2)=>mass flux(kg/(sm^2))
                 endif
-                evap_wat(i) = qfx/rho0(i)                   ! kinematic_surface_upward_latent_heat_flux_over_water
-                hflx_wat(i)=eflx_sh_tot(c)/(rho0(i)*cpair)  ! kinematic_surface_upward_sensible_heat_flux_over_water
+                rho0            = prsl1(i) / (rair*t1(i)*(1.0 + con_fvirt*qv1(i)))
+                evap_wat(i)     = qfx/rho0                    ! kinematic_surface_upward_latent_heat_flux_over_water
+                hflx_wat(i)     = eflx_sh_tot(c)/(rho0*cpair) ! kinematic_surface_upward_sensible_heat_flux_over_water
                 gflx_wat(I)     = eflx_gnet(c)              ![W/m/m]   upward_heat_flux_in_soil_over_water
                 ep1d_water(i)   = eflx_lh_tot(c)            ![W/m/m]   surface_upward_potential_latent_heat_flux_over_water
                 tsurf_water(I)  = t_grnd(c)                 ![K]       surface skin temperature after iteration over water
@@ -754,6 +761,11 @@ SUBROUTINE clm_lake_run( &
                   weasd(i)      = weasdi(i)
                   snowd(i)      = snodi(c)                  ! surface_snow_thickness_water_equivalent_over_ice
 
+
+                  if (.not. icy(i)) then
+                     flag_lakefreeze(i)=.true.
+                  end if
+
                   ! Ice points are icy:
                   icy(i)=.true.                             ! flag_nonzero_sea_ice_surface_fraction
                   ice_points = ice_points+1

physics/clm_lake.meta

@@ -305,14 +305,6 @@
   type = real
   kind = kind_phys
   intent = in
-[rho0]
-  standard_name = air_pressure_at_surface_adjacent_layer
-  long_name = mean pressure at lowest model layer
-  units = Pa
-  dimensions = (horizontal_loop_extent)
-  type = real
-  kind = kind_phys
-  intent = in
 [tsfc]
   standard_name = surface_skin_temperature
   long_name = surface skin temperature
@@ -328,6 +320,13 @@
   dimensions = (horizontal_loop_extent)
   type = logical
   intent = in
+[flag_lakefreeze]
+  standard_name = flag_for_lake_water_freeze
+  long_name = flag for lake water freeze
+  units = flag
+  dimensions = (horizontal_loop_extent)
+  type = logical
+  intent = inout
 [isltyp]
   standard_name = soil_type_classification
   long_name = soil type at each grid cell
@@ -732,6 +731,30 @@
   type = real
   kind = kind_phys
   intent = in
+[t1]
+  standard_name = air_temperature_at_surface_adjacent_layer
+  long_name = mean temperature at lowest model layer
+  units = K
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[qv1]
+  standard_name = specific_humidity_at_surface_adjacent_layer
+  long_name = water vapor specific humidity at lowest model layer
+  units = kg kg-1
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[prsl1]
+  standard_name = air_pressure_at_surface_adjacent_layer
+  long_name = mean pressure at lowest model layer
+  units = Pa
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP

physics/module_sf_ruclsm.F90

@@ -1703,7 +1703,8 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j,         & !--- input varia
        IF (NEWSN > zero .and. snowfracnewsn > 0.99_kind_phys .and. rhosnfall < 450._kind_phys) THEN
        ! new snow
              KEEP_SNOW_ALBEDO = one
-             !snow_mosaic=0.  ! ???
+             ! turn off separate treatment of snow covered and snow-free portions of the grid cell
+             snow_mosaic=0.  ! ???
       ENDIF
 
     IF (debug_print ) THEN
@@ -2076,7 +2077,6 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j,         & !--- input varia
           hfx = hfxs*(one-snowfrac) + hfx*snowfrac
           s = ss*(one-snowfrac) + s*snowfrac
           evapl = evapls*(one-snowfrac)
-          sublim = sublim*snowfrac
           prcpl = prcpls*(one-snowfrac) + prcpl*snowfrac
           fltot = fltots*(one-snowfrac) + fltot*snowfrac
           ALB   = MAX(keep_snow_albedo*alb,              &
@@ -2088,10 +2088,6 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j,         & !--- input varia
 
           runoff1 = runoff1s*(one-snowfrac) + runoff1*snowfrac
           runoff2 = runoff2s*(one-snowfrac) + runoff2*snowfrac
-          smelt = smelt * snowfrac
-          snoh = snoh * snowfrac
-          snflx = snflx * snowfrac
-          snom = snom * snowfrac
           mavail = mavails*(one-snowfrac) + one*snowfrac
           infiltr = infiltrs*(one-snowfrac) + infiltr*snowfrac
 
@@ -2115,7 +2111,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j,         & !--- input varia
           qvg = qvgs*(one-snowfrac) + qvg*snowfrac
           qsg = qsgs*(one-snowfrac) + qsg*snowfrac
           qcg = qcgs*(one-snowfrac) + qcg*snowfrac
-          sublim = eeta*snowfrac
+          sublim = eeta
           eeta = eetas*(one-snowfrac) + eeta*snowfrac
           qfx = qfxs*(one-snowfrac) + qfx*snowfrac
           hfx = hfxs*(one-snowfrac) + hfx*snowfrac
@@ -2129,10 +2125,6 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j,         & !--- input varia
               (emissn - emiss_snowfree) * snowfrac), emissn))
           runoff1 = runoff1s*(one-snowfrac) + runoff1*snowfrac
           runoff2 = runoff2s*(one-snowfrac) + runoff2*snowfrac
-          smelt = smelt * snowfrac
-          snoh = snoh * snowfrac
-          snflx = snflx * snowfrac
-          snom = snom * snowfrac
     IF (debug_print ) THEN
       print *,'SOILT combined on ice', soilt
     ENDIF
@@ -2215,15 +2207,13 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j,         & !--- input varia
        IF (debug_print ) then
        !if (abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)then
          print *,'Snowfallac xlat, xlon',xlat,xlon
-         print *,'newsn,rhonewsn,newsnowratio=',newsn,rhonewsn,newsnowratio
+         print *,'newsn [m],rhonewsn,newsnowratio=',newsn,rhonewsn,newsnowratio
          print *,'Time-step newsn depth [m], swe [m]',newsn,newsn*rhonewsn
          print *,'Time-step smelt: swe [m]' ,smelt*delt
          print *,'Time-step sublim: swe,[kg m-2]',sublim*delt
        endif
 
-      snowfallac = snowfallac + max(zero,(newsn*rhonewsn -                               & ! source of snow (swe) [m]
-                                       (smelt+sublim*1.e-3_kind_phys)*delt*newsnowratio) & ! sink: melting and sublimation, (swe) [m]
-                                       /rhonewsn)*rhowater ! snow accumulation in snow depth [mm]
+      snowfallac = snowfallac + newsn * 1.e3_kind_phys    ! accumulated snow depth [mm], using variable snow density
 
        IF (debug_print ) THEN
        !if (abs(xlat-testptlat).lt.0.2 .and. abs(xlon-testptlon).lt.0.2)then
@@ -5596,7 +5586,7 @@ SUBROUTINE SNOWTEMP( debug_print,xlat,xlon,             &
         nmelt = 1
         soiltfrac=snowfrac*tfrz+(one-snowfrac)*SOILT
         QSG=min(QSG, QSN(soiltfrac,TBQ)/PP)
-        qvg=qsg
+        qvg=snowfrac*qsg+(one-snowfrac)*qvg
         T3      = STBOLT*TN*TN*TN
         UPFLUX  = T3 * 0.5_kind_phys*(TN + SOILTfrac)
         XINET   = EMISS*(GLW-UPFLUX)

physics/sfc_diff.f

@@ -60,6 +60,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav,                &  !intent(in)
      &                    sigmaf,vegtype,shdmax,ivegsrc,                &  !intent(in)
      &                    z0pert,ztpert,                                &  ! mg, sfc-perts !intent(in)
      &                    flag_iter,redrag,                             &  !intent(in)
+     &                    flag_lakefreeze,                              &  !intent(in)             
      &                    u10m,v10m,sfc_z0_type,                        &  !hafs,z0 type !intent(in)
      &                    wet,dry,icy,                                  &  !intent(in)
      &                    thsfc_loc,                                    &  !intent(in)
@@ -90,6 +91,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav,                &  !intent(in)
 
       logical, intent(in) :: redrag ! reduced drag coeff. flag for high wind over sea (j.han)
       logical, dimension(:), intent(in) :: flag_iter, dry, icy
+      logical, dimension(:), intent(in) :: flag_lakefreeze 
       logical, dimension(:), intent(inout) :: wet
 
       logical, intent(in) :: thsfc_loc ! Flag for reference pressure in theta calculation
@@ -168,7 +170,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav,                &  !intent(in)
 !       write(0,*)'in sfc_diff, sfc_z0_type=',sfc_z0_type
 
       do i=1,im
-        if(flag_iter(i)) then
+        if(flag_iter(i) .or. flag_lakefreeze(i)) then
 
           ! Need to initialize ztmax arrays
           ztmax_lnd(i) = 1. ! log(1) = 0

physics/sfc_diff.meta

@@ -194,6 +194,13 @@
   dimensions = ()
   type = logical
   intent = in
+[flag_lakefreeze]
+  standard_name = flag_for_lake_water_freeze
+  long_name = flag for lake water freeze
+  units = flag
+  dimensions = (horizontal_loop_extent)
+  type = logical
+  intent = in
 [u10m]
   standard_name = x_wind_at_10m
   long_name = 10 meter u wind speed