diff --git a/physics/GFS_surface_composites_pre.F90 b/physics/GFS_surface_composites_pre.F90 index 98b9fecd2..fd16dea59 100644 --- a/physics/GFS_surface_composites_pre.F90 +++ b/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) diff --git a/physics/clm_lake.f90 b/physics/clm_lake.f90 index 620f79a96..91e8c71b7 100644 --- a/physics/clm_lake.f90 +++ b/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 diff --git a/physics/clm_lake.meta b/physics/clm_lake.meta index 11a44286a..373bfc308 100644 --- a/physics/clm_lake.meta +++ b/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 diff --git a/physics/module_sf_ruclsm.F90 b/physics/module_sf_ruclsm.F90 index 52fbc8123..b15592052 100644 --- a/physics/module_sf_ruclsm.F90 +++ b/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) diff --git a/physics/sfc_diff.f b/physics/sfc_diff.f index 6e834537a..c5ed8bfa6 100644 --- a/physics/sfc_diff.f +++ b/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 diff --git a/physics/sfc_diff.meta b/physics/sfc_diff.meta index eb30b8c50..1aaad7239 100644 --- a/physics/sfc_diff.meta +++ b/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