From d2900143f7de8c3fe79e8135a0bab5f3267267b9 Mon Sep 17 00:00:00 2001 From: uturuncoglu Date: Tue, 10 Jan 2023 15:31:25 -0700 Subject: [PATCH 01/14] add sfc_land and allow sbs along with fully coupled --- CODEOWNERS | 2 + physics/noahmpdrv.F90 | 12 ++- physics/noahmpdrv.meta | 14 ++++ physics/sfc_land.f | 146 ++++++++++++++++++++++++++++++++ physics/sfc_land.meta | 186 +++++++++++++++++++++++++++++++++++++++++ 5 files changed, 358 insertions(+), 2 deletions(-) create mode 100644 physics/sfc_land.f create mode 100644 physics/sfc_land.meta diff --git a/CODEOWNERS b/CODEOWNERS index cf7a886aa..19e0eb2a5 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -126,6 +126,8 @@ physics/h2ophys.* @AlexBelochitski-NOAA physics/sfc_sice.* @wd20xw @climbfuji @SamuelTrahanNOAA @grantfirl @mzhangw @panll @mkavulich @ChunxiZhang-NOAA physics/sfc_cice.* @wd20xw @climbfuji @SamuelTrahanNOAA @grantfirl @mzhangw @panll @mkavulich @ChunxiZhang-NOAA +physics/sfc_land.* @uturuncoglu @barlage + ######################################################################## # Lines starting with '#' are comments. diff --git a/physics/noahmpdrv.F90 b/physics/noahmpdrv.F90 index fed823ead..d15a9e82a 100644 --- a/physics/noahmpdrv.F90 +++ b/physics/noahmpdrv.F90 @@ -140,7 +140,7 @@ subroutine noahmpdrv_run & iopt_stc, iopt_trs,xlatin, xcoszin, iyrlen, julian, garea, & rainn_mp, rainc_mp, snow_mp, graupel_mp, ice_mp, & con_hvap, con_cp, con_jcal, rhoh2o, con_eps, con_epsm1, & - con_fvirt, con_rd, con_hfus, thsfc_loc, & + con_fvirt, con_rd, con_hfus, thsfc_loc, cpllnd, cpllnd2atm,& ! --- in/outs: weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, & @@ -274,6 +274,9 @@ subroutine noahmpdrv_run & logical , intent(in) :: thsfc_loc ! Flag for reference pressure in theta calculation + logical , intent(in) :: cpllnd ! Flag for land coupling (atm->lnd) + logical , intent(in) :: cpllnd2atm ! Flag for land coupling (lnd->atm) + real(kind=kind_phys), dimension(:) , intent(inout) :: weasd ! water equivalent accumulated snow depth [mm] real(kind=kind_phys), dimension(:) , intent(inout) :: snwdph ! snow depth [mm] real(kind=kind_phys), dimension(:) , intent(inout) :: tskin ! ground surface skin temperature [K] @@ -601,7 +604,12 @@ subroutine noahmpdrv_run & errmsg = '' errflg = 0 -do i = 1, im +! +! --- Just return if external land component is activated for two-way interaction +! + if (cpllnd .and. cpllnd2atm) return + + do i = 1, im if (flag_iter(i) .and. dry(i)) then diff --git a/physics/noahmpdrv.meta b/physics/noahmpdrv.meta index 3235b7c90..643987d98 100644 --- a/physics/noahmpdrv.meta +++ b/physics/noahmpdrv.meta @@ -611,6 +611,20 @@ dimensions = () type = logical intent = in +[cpllnd] + standard_name = flag_for_land_coupling + long_name = flag controlling cpllnd collection (default off) + units = flag + dimensions = () + type = logical + intent = in +[cpllnd2atm] + standard_name = flag_for_one_way_land_coupling_to_atmosphere + long_name = flag controlling land coupling to the atmosphere (default off) + units = flag + dimensions = () + type = logical + intent = in [weasd] standard_name = water_equivalent_accumulated_snow_depth_over_land long_name = water equiv of acc snow depth over land diff --git a/physics/sfc_land.f b/physics/sfc_land.f new file mode 100644 index 000000000..0c3130bbe --- /dev/null +++ b/physics/sfc_land.f @@ -0,0 +1,146 @@ +!> \file sfc_land.f +!! This file contains the code for coupling to land component + +!> This module contains the CCPP-compliant GFS land post +!! interstitial codes, which returns updated surface +!! properties such as latent heat and sensible heat +!! provided by the component version of land model + +!> This module contains the CCPP-compliant GFS land scheme. + module sfc_land + + contains + +!> \defgroup sfc_land for coupling to land +!! @{ +!! \section diagram Calling Hierarchy Diagram +!! \section intraphysics Intraphysics Communication +!! +!> \brief Brief description of the subroutine +!! +!! \section arg_table_sfc_land_run Arguments +!! \htmlinclude sfc_land_run.html +!! + +!! +!! \section general General Algorithm +!! \section detailed Detailed Algorithm +!! @{ + + +!! use physcons, only : hvap => con_hvap, cp => con_cp, & +!! & rvrdm1 => con_fvirt, rd => con_rd +! +!----------------------------------- + subroutine sfc_land_run & +! --- inputs: + & ( im, cpllnd, cpllnd2atm, flag_iter, dry, & + & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, & + & ep_lnd, t2mmp_lnd, q2mp_lnd, & +! --- outputs: + & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, & + & errmsg, errflg, naux2d, aux2d + & ) + +! ===================================================================== ! +! description: ! +! Dec 2022 -- Ufuk Turuncoglu created for coupling to land ! +! ! +! usage: ! +! ! +! call sfc_land ! +! inputs: ! +! ( im, cpllnd, cpllnd2atm, flag_iter, dry, ! +! sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ! +! ep_lnd, t2mmp_lnd, q2mp_lnd, ! +! outputs: ! +! sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, ! +! errmsg, errflg) ! +! ! +! ==================== defination of variables ==================== ! +! ! +! inputs: +! im - integer, horiz dimension +! cpllnd - logical, flag for land coupling +! cpllnd2atm - logical, flag for land coupling (lnd->atm) +! flag_iter - logical, flag for iteration +! dry - logical, eq T if a point with any land +! sncovr1_lnd - real , surface snow area fraction +! qsurf_lnd - real , specific humidity at sfc +! evap_lnd - real , evaporation from latent heat +! hflx_lnd - real , sensible heat +! ep_lnd - real , surface upward potential latent heat flux +! t2mmp_lnd - real , 2m temperature +! q2mp_lnd - real , 2m specific humidity +! outputs: +! sncovr1 - real , snow cover over land +! qsurf - real , specific humidity at sfc +! evap - real , evaporation from latent heat +! hflx - real , sensible heat +! ep - real , potential evaporation +! t2mmp - real , temperature at 2m +! q2mp - real , specific humidity at 2m +! ==================== end of description ===================== ! +! +! + use machine , only : kind_phys + implicit none + +! --- inputs: + integer, intent(in) :: im + logical, intent(in) :: cpllnd, cpllnd2atm + logical, dimension(:), intent(in) :: flag_iter + logical, dimension(:), intent(in) :: dry + + real (kind=kind_phys), dimension(:), intent(in) :: & + & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ep_lnd, & + & t2mmp_lnd, q2mp_lnd + +! --- outputs: + real (kind=kind_phys), dimension(:), intent(out) :: & + & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp +! + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + integer, intent(in) :: naux2d + real(kind_phys), intent(out) :: aux2d(:,:) + +! --- locals: + + integer :: i + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 +! + if (.not. cpllnd2atm) return +! + do i = 1, im + !if (flag_iter(i) .and. dry(i)) then + !if (dry(i)) then + sncovr1(i) = sncovr1_lnd(i) + qsurf(i) = qsurf_lnd(i) + hflx(i) = hflx_lnd(i) + evap(i) = evap_lnd(i) + ep(i) = ep_lnd(i) + t2mmp(i) = t2mmp_lnd(i) + q2mp(i) = q2mp_lnd(i) + !end if + enddo + + aux2d(:,1) = dry(:) !sncovr1(:) + aux2d(:,2) = qsurf(:) + aux2d(:,3) = hflx(:) + aux2d(:,4) = evap(:) + aux2d(:,5) = ep(:) + aux2d(:,6) = qsurf_lnd(:) !t2mmp(:) + aux2d(:,7) = q2mp(:) + + return +!----------------------------------- + end subroutine sfc_land_run +!----------------------------------- + +!> @} + end module sfc_land diff --git a/physics/sfc_land.meta b/physics/sfc_land.meta new file mode 100644 index 000000000..f31d779ae --- /dev/null +++ b/physics/sfc_land.meta @@ -0,0 +1,186 @@ +[ccpp-table-properties] + name = sfc_land + type = scheme + dependencies = machine.F + +######################################################################## +[ccpp-arg-table] + name = sfc_land_run + type = scheme +[im] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in +[cpllnd] + standard_name = flag_for_land_coupling + long_name = flag controlling cpllnd collection (default off) + units = flag + dimensions = () + type = logical + intent = in +[cpllnd2atm] + standard_name = flag_for_one_way_land_coupling_to_atmosphere + long_name = flag controlling land coupling to the atmosphere (default off) + units = flag + dimensions = () + type = logical + intent = in +[flag_iter] + standard_name = flag_for_iteration + long_name = flag for iteration + units = flag + dimensions = (horizontal_loop_extent) + type = logical + intent = in +[dry] + standard_name = flag_nonzero_land_surface_fraction + long_name = flag indicating presence of some land surface area fraction + units = flag + dimensions = (horizontal_loop_extent) + type = logical + intent = in +[sncovr1_lnd] + standard_name = surface_snow_area_fraction_over_land_from_land + long_name = surface snow area fraction over land for coupling + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[qsurf_lnd] + standard_name = surface_specific_humidity_over_land_from_land + long_name = surface air saturation specific humidity over land + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[evap_lnd] + standard_name = surface_upward_latent_heat_flux_over_land_from_land + long_name = sfc latent heat flux input over land for coupling + units = kg kg-1 m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[hflx_lnd] + standard_name = surface_upward_sensible_heat_flux_over_land_from_land + long_name = sfc sensible heat flux input over land for coupling + units = K m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[ep_lnd] + standard_name = surface_upward_potential_latent_heat_flux_over_land_from_land + long_name = surface upward potential latent heat flux over land for coupling + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[t2mmp_lnd] + standard_name = temperature_at_2m_over_land_from_land + long_name = 2 meter temperature over land for coupling + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[q2mp_lnd] + standard_name = specific_humidity_at_2m_over_land_from_land + long_name = 2 meter specific humidity over land for coupling + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[sncovr1] + standard_name = surface_snow_area_fraction_over_land + long_name = surface snow area fraction + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[qsurf] + standard_name = surface_specific_humidity_over_land + long_name = surface air saturation specific humidity over land + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[evap] + standard_name = kinematic_surface_upward_latent_heat_flux_over_land + long_name = kinematic surface upward latent heat flux over land + units = kg kg-1 m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[hflx] + standard_name = kinematic_surface_upward_sensible_heat_flux_over_land + long_name = kinematic surface upward sensible heat flux over land + units = K m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[ep] + standard_name = surface_upward_potential_latent_heat_flux_over_land + long_name = surface upward potential latent heat flux over land + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[t2mmp] + standard_name = temperature_at_2m_from_noahmp + long_name = 2 meter temperature from noahmp + units = K + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[q2mp] + standard_name = specific_humidity_at_2m_from_noahmp + long_name = 2 meter specific humidity from noahmp + units = kg kg-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out +[errflg] + standard_name = ccpp_error_code + long_name = error code for error handling in CCPP + units = 1 + dimensions = () + type = integer + intent = out +[naux2d] + standard_name = number_of_xy_dimensioned_auxiliary_arrays + long_name = number of 2d auxiliary arrays to output (for debugging) + units = count + dimensions = () + type = integer + intent = in +[aux2d] + standard_name = auxiliary_2d_arrays + long_name = auxiliary 2d arrays to output (for debugging) + units = none + dimensions = (horizontal_loop_extent,number_of_xy_dimensioned_auxiliary_arrays) + type = real + kind = kind_phys + intent = out From 7e2954615d03ee085c2cbc91d6326ee2ecac88cb Mon Sep 17 00:00:00 2001 From: uturuncoglu Date: Tue, 7 Feb 2023 15:12:09 -0700 Subject: [PATCH 02/14] update sfc_land --- physics/sfc_land.f | 39 +++++++++++++++++++++------------------ physics/sfc_land.meta | 16 ++++++++++++++++ 2 files changed, 37 insertions(+), 18 deletions(-) diff --git a/physics/sfc_land.f b/physics/sfc_land.f index 0c3130bbe..ab0691251 100644 --- a/physics/sfc_land.f +++ b/physics/sfc_land.f @@ -36,9 +36,10 @@ subroutine sfc_land_run & ! --- inputs: & ( im, cpllnd, cpllnd2atm, flag_iter, dry, & & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, & - & ep_lnd, t2mmp_lnd, q2mp_lnd, & + & ep_lnd, t2mmp_lnd, q2mp_lnd, gflux_lnd, & ! --- outputs: & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, & + & gflux, & & errmsg, errflg, naux2d, aux2d & ) @@ -52,9 +53,10 @@ subroutine sfc_land_run & ! inputs: ! ! ( im, cpllnd, cpllnd2atm, flag_iter, dry, ! ! sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ! -! ep_lnd, t2mmp_lnd, q2mp_lnd, ! +! ep_lnd, t2mmp_lnd, q2mp_lnd, gflux_lnd, ! ! outputs: ! ! sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, ! +! gflux, ! ! errmsg, errflg) ! ! ! ! ==================== defination of variables ==================== ! @@ -71,7 +73,8 @@ subroutine sfc_land_run & ! hflx_lnd - real , sensible heat ! ep_lnd - real , surface upward potential latent heat flux ! t2mmp_lnd - real , 2m temperature -! q2mp_lnd - real , 2m specific humidity +! q2mp_lnd - real , 2m specific humidity +! gflux_lnd - real , soil heat flux over land ! outputs: ! sncovr1 - real , snow cover over land ! qsurf - real , specific humidity at sfc @@ -79,7 +82,8 @@ subroutine sfc_land_run & ! hflx - real , sensible heat ! ep - real , potential evaporation ! t2mmp - real , temperature at 2m -! q2mp - real , specific humidity at 2m +! q2mp - real , specific humidity at 2m +! gflux - real , soil heat flux over land ! ==================== end of description ===================== ! ! ! @@ -94,11 +98,11 @@ subroutine sfc_land_run & real (kind=kind_phys), dimension(:), intent(in) :: & & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ep_lnd, & - & t2mmp_lnd, q2mp_lnd + & t2mmp_lnd, q2mp_lnd, gflux_lnd ! --- outputs: real (kind=kind_phys), dimension(:), intent(out) :: & - & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp + & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, gflux ! character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -117,25 +121,24 @@ subroutine sfc_land_run & if (.not. cpllnd2atm) return ! do i = 1, im - !if (flag_iter(i) .and. dry(i)) then - !if (dry(i)) then - sncovr1(i) = sncovr1_lnd(i) - qsurf(i) = qsurf_lnd(i) - hflx(i) = hflx_lnd(i) - evap(i) = evap_lnd(i) - ep(i) = ep_lnd(i) - t2mmp(i) = t2mmp_lnd(i) - q2mp(i) = q2mp_lnd(i) - !end if + sncovr1(i) = sncovr1_lnd(i) + qsurf(i) = qsurf_lnd(i) + hflx(i) = hflx_lnd(i) + evap(i) = evap_lnd(i) + ep(i) = ep_lnd(i) + t2mmp(i) = t2mmp_lnd(i) + q2mp(i) = q2mp_lnd(i) + gflux(i) = gflux_lnd(i) enddo - aux2d(:,1) = dry(:) !sncovr1(:) + aux2d(:,1) = sncovr1(:) aux2d(:,2) = qsurf(:) aux2d(:,3) = hflx(:) aux2d(:,4) = evap(:) aux2d(:,5) = ep(:) - aux2d(:,6) = qsurf_lnd(:) !t2mmp(:) + aux2d(:,6) = t2mmp(:) aux2d(:,7) = q2mp(:) + aux2d(:,8) = gflux(:) return !----------------------------------- diff --git a/physics/sfc_land.meta b/physics/sfc_land.meta index f31d779ae..50ddecd46 100644 --- a/physics/sfc_land.meta +++ b/physics/sfc_land.meta @@ -98,6 +98,14 @@ type = real kind = kind_phys intent = in +[gflux_lnd] + standard_name = upward_heat_flux_in_soil_over_land_from_land + long_name = soil heat flux over land for coupling + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in [sncovr1] standard_name = surface_snow_area_fraction_over_land long_name = surface snow area fraction @@ -154,6 +162,14 @@ type = real kind = kind_phys intent = out +[gflux] + standard_name = upward_heat_flux_in_soil_over_land + long_name = soil heat flux over land + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 40e092d0d92d0c78848acc47ae83629908375a6c Mon Sep 17 00:00:00 2001 From: uturuncoglu Date: Thu, 23 Feb 2023 14:28:26 -0700 Subject: [PATCH 03/14] add runoff and drain to land coupling --- physics/sfc_land.f | 19 +++++++++++++++---- physics/sfc_land.meta | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 4 deletions(-) diff --git a/physics/sfc_land.f b/physics/sfc_land.f index ab0691251..f7aebe171 100644 --- a/physics/sfc_land.f +++ b/physics/sfc_land.f @@ -37,9 +37,10 @@ subroutine sfc_land_run & & ( im, cpllnd, cpllnd2atm, flag_iter, dry, & & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, & & ep_lnd, t2mmp_lnd, q2mp_lnd, gflux_lnd, & + & runoff_lnd, drain_lnd, ! --- outputs: & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, & - & gflux, & + & gflux, runoff, drain, & & errmsg, errflg, naux2d, aux2d & ) @@ -54,9 +55,10 @@ subroutine sfc_land_run & ! ( im, cpllnd, cpllnd2atm, flag_iter, dry, ! ! sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ! ! ep_lnd, t2mmp_lnd, q2mp_lnd, gflux_lnd, ! +! runoff_lnd, drain_lnd, ! ! outputs: ! ! sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, ! -! gflux, ! +! gflux, runoff, drain, ! ! errmsg, errflg) ! ! ! ! ==================== defination of variables ==================== ! @@ -75,6 +77,8 @@ subroutine sfc_land_run & ! t2mmp_lnd - real , 2m temperature ! q2mp_lnd - real , 2m specific humidity ! gflux_lnd - real , soil heat flux over land +! runoff_lnd - real , surface runoff +! drain_lnd - real , subsurface runoff ! outputs: ! sncovr1 - real , snow cover over land ! qsurf - real , specific humidity at sfc @@ -84,6 +88,8 @@ subroutine sfc_land_run & ! t2mmp - real , temperature at 2m ! q2mp - real , specific humidity at 2m ! gflux - real , soil heat flux over land +! runoff - real , surface runoff +! drain - real , subsurface runoff ! ==================== end of description ===================== ! ! ! @@ -98,11 +104,12 @@ subroutine sfc_land_run & real (kind=kind_phys), dimension(:), intent(in) :: & & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ep_lnd, & - & t2mmp_lnd, q2mp_lnd, gflux_lnd + & t2mmp_lnd, q2mp_lnd, gflux_lnd, runoff_lnd, drain_lnd ! --- outputs: real (kind=kind_phys), dimension(:), intent(out) :: & - & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, gflux + & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, gflux, & + & runoff, drain ! character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -129,6 +136,8 @@ subroutine sfc_land_run & t2mmp(i) = t2mmp_lnd(i) q2mp(i) = q2mp_lnd(i) gflux(i) = gflux_lnd(i) + drain(i) = drain_lnd(i) + runoff(i) = runoff_lnd(i) enddo aux2d(:,1) = sncovr1(:) @@ -139,6 +148,8 @@ subroutine sfc_land_run & aux2d(:,6) = t2mmp(:) aux2d(:,7) = q2mp(:) aux2d(:,8) = gflux(:) + aux2d(:,9) = drain(:) + aux2d(:,10) = runoff(:) return !----------------------------------- diff --git a/physics/sfc_land.meta b/physics/sfc_land.meta index 50ddecd46..60a853b89 100644 --- a/physics/sfc_land.meta +++ b/physics/sfc_land.meta @@ -106,6 +106,22 @@ type = real kind = kind_phys intent = in +[runoff_lnd] + standard_name = surface_runoff_flux_from_land + long_name = surface runoff flux over land for coupling + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in +[drain_lnd] + standard_name = subsurface_runoff_flux_from_land + long_name = subsurface runoff flux over land for coupling + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in [sncovr1] standard_name = surface_snow_area_fraction_over_land long_name = surface snow area fraction @@ -170,6 +186,22 @@ type = real kind = kind_phys intent = out +[runoff] + standard_name = surface_runoff_flux + long_name = surface runoff flux + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[drain] + standard_name = subsurface_runoff_flux + long_name = subsurface runoff flux + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 64fdd5aa53672132bf1295701a9d7844232f336b Mon Sep 17 00:00:00 2001 From: uturuncoglu Date: Tue, 7 Mar 2023 11:00:00 -0700 Subject: [PATCH 04/14] add exchange coefficents --- physics/sfc_land.f | 23 ++++++++++++++++------- physics/sfc_land.meta | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 48 insertions(+), 7 deletions(-) diff --git a/physics/sfc_land.f b/physics/sfc_land.f index f7aebe171..0436519a5 100644 --- a/physics/sfc_land.f +++ b/physics/sfc_land.f @@ -37,10 +37,10 @@ subroutine sfc_land_run & & ( im, cpllnd, cpllnd2atm, flag_iter, dry, & & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, & & ep_lnd, t2mmp_lnd, q2mp_lnd, gflux_lnd, & - & runoff_lnd, drain_lnd, + & runoff_lnd, drain_lnd, cmm_lnd, chh_lnd, ! --- outputs: & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, & - & gflux, runoff, drain, & + & gflux, runoff, drain, cmm, chh, & & errmsg, errflg, naux2d, aux2d & ) @@ -55,10 +55,10 @@ subroutine sfc_land_run & ! ( im, cpllnd, cpllnd2atm, flag_iter, dry, ! ! sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ! ! ep_lnd, t2mmp_lnd, q2mp_lnd, gflux_lnd, ! -! runoff_lnd, drain_lnd, ! +! runoff_lnd, drain_lnd, cmm_lnd, chh_lnd, ! ! outputs: ! ! sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, ! -! gflux, runoff, drain, ! +! gflux, runoff, drain, cmm, chh, ! ! errmsg, errflg) ! ! ! ! ==================== defination of variables ==================== ! @@ -79,6 +79,8 @@ subroutine sfc_land_run & ! gflux_lnd - real , soil heat flux over land ! runoff_lnd - real , surface runoff ! drain_lnd - real , subsurface runoff +! cmm_lnd - real , surface drag wind speed for momentum +! chh_lnd - real , surface drag mass flux for heat and moisture ! outputs: ! sncovr1 - real , snow cover over land ! qsurf - real , specific humidity at sfc @@ -89,7 +91,9 @@ subroutine sfc_land_run & ! q2mp - real , specific humidity at 2m ! gflux - real , soil heat flux over land ! runoff - real , surface runoff -! drain - real , subsurface runoff +! drain - real , subsurface runoff +! cmm - real , surface drag wind speed for momentum +! chh - real , surface drag mass flux for heat and moisture ! ==================== end of description ===================== ! ! ! @@ -104,12 +108,13 @@ subroutine sfc_land_run & real (kind=kind_phys), dimension(:), intent(in) :: & & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ep_lnd, & - & t2mmp_lnd, q2mp_lnd, gflux_lnd, runoff_lnd, drain_lnd + & t2mmp_lnd, q2mp_lnd, gflux_lnd, runoff_lnd, drain_lnd, & + & cmm_lnd, chh_lnd ! --- outputs: real (kind=kind_phys), dimension(:), intent(out) :: & & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, gflux, & - & runoff, drain + & runoff, drain, cmm, chh ! character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -138,6 +143,8 @@ subroutine sfc_land_run & gflux(i) = gflux_lnd(i) drain(i) = drain_lnd(i) runoff(i) = runoff_lnd(i) + cmm(i) = cmm_lnd(i) + chh(i) = chh_lnd(i) enddo aux2d(:,1) = sncovr1(:) @@ -150,6 +157,8 @@ subroutine sfc_land_run & aux2d(:,8) = gflux(:) aux2d(:,9) = drain(:) aux2d(:,10) = runoff(:) + aux2d(:,11) = cmm(:) + aux2d(:,12) = chh(:) return !----------------------------------- diff --git a/physics/sfc_land.meta b/physics/sfc_land.meta index 60a853b89..99a795c65 100644 --- a/physics/sfc_land.meta +++ b/physics/sfc_land.meta @@ -122,6 +122,22 @@ type = real kind = kind_phys intent = in +[cmm_lnd] + standard_name = surface_drag_wind_speed_for_momentum_in_air_over_land_from_land + long_name = momentum exchange coefficient over land for coupling + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[chh_lnd] + standard_name = surface_drag_mass_flux_for_heat_and_moisture_in_air_over_land_from_land + long_name = thermal exchange coefficient over land for coupling + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in [sncovr1] standard_name = surface_snow_area_fraction_over_land long_name = surface snow area fraction @@ -202,6 +218,22 @@ type = real kind = kind_phys intent = out +[cmm] + standard_name = surface_drag_wind_speed_for_momentum_in_air_over_land + long_name = momentum exchange coefficient over land + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[chh] + standard_name = surface_drag_mass_flux_for_heat_and_moisture_in_air_over_land + long_name = thermal exchange coefficient over land + units = kg m-2 s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 91a9b44eac42dfaf047cfc79718fa0332db0979d Mon Sep 17 00:00:00 2001 From: uturuncoglu Date: Tue, 7 Mar 2023 22:40:41 -0700 Subject: [PATCH 05/14] add zvfun to land coupling --- physics/sfc_land.f | 15 ++++++++++----- physics/sfc_land.meta | 16 ++++++++++++++++ 2 files changed, 26 insertions(+), 5 deletions(-) diff --git a/physics/sfc_land.f b/physics/sfc_land.f index 0436519a5..44a9b5a06 100644 --- a/physics/sfc_land.f +++ b/physics/sfc_land.f @@ -37,10 +37,10 @@ subroutine sfc_land_run & & ( im, cpllnd, cpllnd2atm, flag_iter, dry, & & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, & & ep_lnd, t2mmp_lnd, q2mp_lnd, gflux_lnd, & - & runoff_lnd, drain_lnd, cmm_lnd, chh_lnd, + & runoff_lnd, drain_lnd, cmm_lnd, chh_lnd, zvfun_lnd, & ! --- outputs: & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, & - & gflux, runoff, drain, cmm, chh, & + & gflux, runoff, drain, cmm, chh, zvfun, & & errmsg, errflg, naux2d, aux2d & ) @@ -56,9 +56,10 @@ subroutine sfc_land_run & ! sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ! ! ep_lnd, t2mmp_lnd, q2mp_lnd, gflux_lnd, ! ! runoff_lnd, drain_lnd, cmm_lnd, chh_lnd, ! +! zvfun_lnd, ! ! outputs: ! ! sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, ! -! gflux, runoff, drain, cmm, chh, ! +! gflux, runoff, drain, cmm, chh, zvfun, ! ! errmsg, errflg) ! ! ! ! ==================== defination of variables ==================== ! @@ -81,6 +82,7 @@ subroutine sfc_land_run & ! drain_lnd - real , subsurface runoff ! cmm_lnd - real , surface drag wind speed for momentum ! chh_lnd - real , surface drag mass flux for heat and moisture +! zvfun_lnd - real , function of surface roughness length and green vegetation fraction ! outputs: ! sncovr1 - real , snow cover over land ! qsurf - real , specific humidity at sfc @@ -94,6 +96,7 @@ subroutine sfc_land_run & ! drain - real , subsurface runoff ! cmm - real , surface drag wind speed for momentum ! chh - real , surface drag mass flux for heat and moisture +! zvfun - real , function of surface roughness length and green vegetation fraction ! ==================== end of description ===================== ! ! ! @@ -109,12 +112,12 @@ subroutine sfc_land_run & real (kind=kind_phys), dimension(:), intent(in) :: & & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ep_lnd, & & t2mmp_lnd, q2mp_lnd, gflux_lnd, runoff_lnd, drain_lnd, & - & cmm_lnd, chh_lnd + & cmm_lnd, chh_lnd, zvfun_lnd ! --- outputs: real (kind=kind_phys), dimension(:), intent(out) :: & & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, gflux, & - & runoff, drain, cmm, chh + & runoff, drain, cmm, chh, zvfun ! character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -145,6 +148,7 @@ subroutine sfc_land_run & runoff(i) = runoff_lnd(i) cmm(i) = cmm_lnd(i) chh(i) = chh_lnd(i) + zvfun(i) = zvfun_lnd(i) enddo aux2d(:,1) = sncovr1(:) @@ -159,6 +163,7 @@ subroutine sfc_land_run & aux2d(:,10) = runoff(:) aux2d(:,11) = cmm(:) aux2d(:,12) = chh(:) + aux2d(:,13) = zvfun(:) return !----------------------------------- diff --git a/physics/sfc_land.meta b/physics/sfc_land.meta index 99a795c65..a146dec3a 100644 --- a/physics/sfc_land.meta +++ b/physics/sfc_land.meta @@ -138,6 +138,14 @@ type = real kind = kind_phys intent = in +[zvfun_lnd] + standard_name = function_of_surface_roughness_length_and_green_vegetation_fraction_from_land + long_name = function of surface roughness length and green vegetation fraction + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in [sncovr1] standard_name = surface_snow_area_fraction_over_land long_name = surface snow area fraction @@ -234,6 +242,14 @@ type = real kind = kind_phys intent = out +[zvfun] + standard_name = function_of_surface_roughness_length_and_green_vegetation_fraction + long_name = function of surface roughness length and green vegetation fraction + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 8442a202769b4082d548c4d0d0b5d0c1cb7a80c1 Mon Sep 17 00:00:00 2001 From: uturuncoglu Date: Wed, 5 Apr 2023 15:49:56 -0600 Subject: [PATCH 06/14] clean sfc_land --- physics/sfc_land.f | 22 +--------------------- physics/sfc_land.meta | 15 --------------- 2 files changed, 1 insertion(+), 36 deletions(-) diff --git a/physics/sfc_land.f b/physics/sfc_land.f index 44a9b5a06..aec47ff77 100644 --- a/physics/sfc_land.f +++ b/physics/sfc_land.f @@ -27,9 +27,6 @@ module sfc_land !! \section detailed Detailed Algorithm !! @{ - -!! use physcons, only : hvap => con_hvap, cp => con_cp, & -!! & rvrdm1 => con_fvirt, rd => con_rd ! !----------------------------------- subroutine sfc_land_run & @@ -41,7 +38,7 @@ subroutine sfc_land_run & ! --- outputs: & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, & & gflux, runoff, drain, cmm, chh, zvfun, & - & errmsg, errflg, naux2d, aux2d + & errmsg, errflg & ) ! ===================================================================== ! @@ -122,9 +119,6 @@ subroutine sfc_land_run & character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg - integer, intent(in) :: naux2d - real(kind_phys), intent(out) :: aux2d(:,:) - ! --- locals: integer :: i @@ -150,20 +144,6 @@ subroutine sfc_land_run & chh(i) = chh_lnd(i) zvfun(i) = zvfun_lnd(i) enddo - - aux2d(:,1) = sncovr1(:) - aux2d(:,2) = qsurf(:) - aux2d(:,3) = hflx(:) - aux2d(:,4) = evap(:) - aux2d(:,5) = ep(:) - aux2d(:,6) = t2mmp(:) - aux2d(:,7) = q2mp(:) - aux2d(:,8) = gflux(:) - aux2d(:,9) = drain(:) - aux2d(:,10) = runoff(:) - aux2d(:,11) = cmm(:) - aux2d(:,12) = chh(:) - aux2d(:,13) = zvfun(:) return !----------------------------------- diff --git a/physics/sfc_land.meta b/physics/sfc_land.meta index a146dec3a..979cca377 100644 --- a/physics/sfc_land.meta +++ b/physics/sfc_land.meta @@ -265,18 +265,3 @@ dimensions = () type = integer intent = out -[naux2d] - standard_name = number_of_xy_dimensioned_auxiliary_arrays - long_name = number of 2d auxiliary arrays to output (for debugging) - units = count - dimensions = () - type = integer - intent = in -[aux2d] - standard_name = auxiliary_2d_arrays - long_name = auxiliary 2d arrays to output (for debugging) - units = none - dimensions = (horizontal_loop_extent,number_of_xy_dimensioned_auxiliary_arrays) - type = real - kind = kind_phys - intent = out From 9891fff7bd9eb8660a60ff203112298bf1349406 Mon Sep 17 00:00:00 2001 From: uturuncoglu Date: Mon, 31 Jul 2023 10:50:36 -0600 Subject: [PATCH 07/14] switch from in to inout for output variables --- physics/sfc_land.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/sfc_land.f b/physics/sfc_land.f index aec47ff77..d4e88c25a 100644 --- a/physics/sfc_land.f +++ b/physics/sfc_land.f @@ -112,7 +112,7 @@ subroutine sfc_land_run & & cmm_lnd, chh_lnd, zvfun_lnd ! --- outputs: - real (kind=kind_phys), dimension(:), intent(out) :: & + real (kind=kind_phys), dimension(:), intent(inout) :: & & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, gflux, & & runoff, drain, cmm, chh, zvfun ! From 70038f6f5f70572f09489732e3563f1d11066a1f Mon Sep 17 00:00:00 2001 From: uturuncoglu Date: Mon, 31 Jul 2023 11:50:02 -0600 Subject: [PATCH 08/14] update meta file for sfc_land too --- physics/sfc_land.meta | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/physics/sfc_land.meta b/physics/sfc_land.meta index 979cca377..493d2a70b 100644 --- a/physics/sfc_land.meta +++ b/physics/sfc_land.meta @@ -153,7 +153,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [qsurf] standard_name = surface_specific_humidity_over_land long_name = surface air saturation specific humidity over land @@ -161,7 +161,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [evap] standard_name = kinematic_surface_upward_latent_heat_flux_over_land long_name = kinematic surface upward latent heat flux over land @@ -169,7 +169,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [hflx] standard_name = kinematic_surface_upward_sensible_heat_flux_over_land long_name = kinematic surface upward sensible heat flux over land @@ -177,7 +177,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [ep] standard_name = surface_upward_potential_latent_heat_flux_over_land long_name = surface upward potential latent heat flux over land @@ -185,7 +185,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [t2mmp] standard_name = temperature_at_2m_from_noahmp long_name = 2 meter temperature from noahmp @@ -193,7 +193,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [q2mp] standard_name = specific_humidity_at_2m_from_noahmp long_name = 2 meter specific humidity from noahmp @@ -201,7 +201,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [gflux] standard_name = upward_heat_flux_in_soil_over_land long_name = soil heat flux over land @@ -209,7 +209,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [runoff] standard_name = surface_runoff_flux long_name = surface runoff flux @@ -217,7 +217,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [drain] standard_name = subsurface_runoff_flux long_name = subsurface runoff flux @@ -225,7 +225,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [cmm] standard_name = surface_drag_wind_speed_for_momentum_in_air_over_land long_name = momentum exchange coefficient over land @@ -233,7 +233,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [chh] standard_name = surface_drag_mass_flux_for_heat_and_moisture_in_air_over_land long_name = thermal exchange coefficient over land @@ -241,7 +241,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [zvfun] standard_name = function_of_surface_roughness_length_and_green_vegetation_fraction long_name = function of surface roughness length and green vegetation fraction @@ -249,7 +249,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 6ec72e918d2f2d04f78a9afaf5b58968370a44a4 Mon Sep 17 00:00:00 2001 From: uturuncoglu Date: Mon, 31 Jul 2023 16:21:27 -0600 Subject: [PATCH 09/14] update sfc_land --- physics/sfc_land.F90 | 108 ++++++++++++++++++++++++++++++ physics/sfc_land.f | 154 ------------------------------------------- 2 files changed, 108 insertions(+), 154 deletions(-) create mode 100644 physics/sfc_land.F90 delete mode 100644 physics/sfc_land.f diff --git a/physics/sfc_land.F90 b/physics/sfc_land.F90 new file mode 100644 index 000000000..2b0696ed8 --- /dev/null +++ b/physics/sfc_land.F90 @@ -0,0 +1,108 @@ +!> \file sfc_land.F90 +!! This file contains the code for coupling to land component + +!> This module contains the CCPP-compliant GFS land post +!! interstitial codes, which returns updated surface +!! properties such as latent heat and sensible heat +!! provided by the component version of land model + +!> This module contains the CCPP-compliant GFS land scheme. + module sfc_land + + use machine, only : kind_phys + + contains + +!> \defgroup sfc_land for coupling to land +!! @{ +!! \section diagram Calling Hierarchy Diagram +!! \section intraphysics Intraphysics Communication +!! +!> \brief Brief description of the subroutine +!! +!! \section arg_table_sfc_land_run Arguments +!! \htmlinclude sfc_land_run.html +!! + +!! +!! \section general General Algorithm +!! \section detailed Detailed Algorithm +!! @{ + subroutine sfc_land_run(im, cpllnd, cpllnd2atm, flag_iter, dry, & + sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, & + ep_lnd, t2mmp_lnd, q2mp_lnd, gflux_lnd, & + runoff_lnd, drain_lnd, cmm_lnd, chh_lnd, zvfun_lnd, & + sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, & + gflux, runoff, drain, cmm, chh, zvfun, & + errmsg, errflg) + + implicit none + + ! Inputs + integer , intent(in) :: im + logical , intent(in) :: cpllnd + logical , intent(in) :: cpllnd2atm + logical , intent(in) :: flag_iter(:) + logical , intent(in) :: dry(:) + real(kind=kind_phys), intent(in) :: sncovr1_lnd(:) + real(kind=kind_phys), intent(in) :: qsurf_lnd(:) + real(kind=kind_phys), intent(in) :: evap_lnd(:) + real(kind=kind_phys), intent(in) :: hflx_lnd(:) + real(kind=kind_phys), intent(in) :: ep_lnd(:) + real(kind=kind_phys), intent(in) :: t2mmp_lnd(:) + real(kind=kind_phys), intent(in) :: q2mp_lnd(:) + real(kind=kind_phys), intent(in) :: gflux_lnd(:) + real(kind=kind_phys), intent(in) :: runoff_lnd(:) + real(kind=kind_phys), intent(in) :: drain_lnd(:) + real(kind=kind_phys), intent(in) :: cmm_lnd(:) + real(kind=kind_phys), intent(in) :: chh_lnd(:) + real(kind=kind_phys), intent(in) :: zvfun_lnd(:) + ! Inputs/Outputs + real(kind=kind_phys), intent(inout) :: sncovr1(:) + real(kind=kind_phys), intent(inout) :: qsurf(:) + real(kind=kind_phys), intent(inout) :: evap(:) + real(kind=kind_phys), intent(inout) :: hflx(:) + real(kind=kind_phys), intent(inout) :: ep(:) + real(kind=kind_phys), intent(inout) :: t2mmp(:) + real(kind=kind_phys), intent(inout) :: q2mp(:) + real(kind=kind_phys), intent(inout) :: gflux(:) + real(kind=kind_phys), intent(inout) :: runoff(:) + real(kind=kind_phys), intent(inout) :: drain(:) + real(kind=kind_phys), intent(inout) :: cmm(:) + real(kind=kind_phys), intent(inout) :: chh(:) + real(kind=kind_phys), intent(inout) :: zvfun(:) + ! Outputs + character(len=*) , intent(out) :: errmsg + integer , intent(out) :: errflg + + ! Locals + integer :: i + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Check coupling from component land to atmosphere + if (.not. cpllnd2atm) return + + ! Fill variables + do i = 1, im + sncovr1(i) = sncovr1_lnd(i) + qsurf(i) = qsurf_lnd(i) + hflx(i) = hflx_lnd(i) + evap(i) = evap_lnd(i) + ep(i) = ep_lnd(i) + t2mmp(i) = t2mmp_lnd(i) + q2mp(i) = q2mp_lnd(i) + gflux(i) = gflux_lnd(i) + drain(i) = drain_lnd(i) + runoff(i) = runoff_lnd(i) + cmm(i) = cmm_lnd(i) + chh(i) = chh_lnd(i) + zvfun(i) = zvfun_lnd(i) + enddo + + end subroutine sfc_land_run + +!> @} + end module sfc_land diff --git a/physics/sfc_land.f b/physics/sfc_land.f deleted file mode 100644 index d4e88c25a..000000000 --- a/physics/sfc_land.f +++ /dev/null @@ -1,154 +0,0 @@ -!> \file sfc_land.f -!! This file contains the code for coupling to land component - -!> This module contains the CCPP-compliant GFS land post -!! interstitial codes, which returns updated surface -!! properties such as latent heat and sensible heat -!! provided by the component version of land model - -!> This module contains the CCPP-compliant GFS land scheme. - module sfc_land - - contains - -!> \defgroup sfc_land for coupling to land -!! @{ -!! \section diagram Calling Hierarchy Diagram -!! \section intraphysics Intraphysics Communication -!! -!> \brief Brief description of the subroutine -!! -!! \section arg_table_sfc_land_run Arguments -!! \htmlinclude sfc_land_run.html -!! - -!! -!! \section general General Algorithm -!! \section detailed Detailed Algorithm -!! @{ - -! -!----------------------------------- - subroutine sfc_land_run & -! --- inputs: - & ( im, cpllnd, cpllnd2atm, flag_iter, dry, & - & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, & - & ep_lnd, t2mmp_lnd, q2mp_lnd, gflux_lnd, & - & runoff_lnd, drain_lnd, cmm_lnd, chh_lnd, zvfun_lnd, & -! --- outputs: - & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, & - & gflux, runoff, drain, cmm, chh, zvfun, & - & errmsg, errflg - & ) - -! ===================================================================== ! -! description: ! -! Dec 2022 -- Ufuk Turuncoglu created for coupling to land ! -! ! -! usage: ! -! ! -! call sfc_land ! -! inputs: ! -! ( im, cpllnd, cpllnd2atm, flag_iter, dry, ! -! sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ! -! ep_lnd, t2mmp_lnd, q2mp_lnd, gflux_lnd, ! -! runoff_lnd, drain_lnd, cmm_lnd, chh_lnd, ! -! zvfun_lnd, ! -! outputs: ! -! sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, ! -! gflux, runoff, drain, cmm, chh, zvfun, ! -! errmsg, errflg) ! -! ! -! ==================== defination of variables ==================== ! -! ! -! inputs: -! im - integer, horiz dimension -! cpllnd - logical, flag for land coupling -! cpllnd2atm - logical, flag for land coupling (lnd->atm) -! flag_iter - logical, flag for iteration -! dry - logical, eq T if a point with any land -! sncovr1_lnd - real , surface snow area fraction -! qsurf_lnd - real , specific humidity at sfc -! evap_lnd - real , evaporation from latent heat -! hflx_lnd - real , sensible heat -! ep_lnd - real , surface upward potential latent heat flux -! t2mmp_lnd - real , 2m temperature -! q2mp_lnd - real , 2m specific humidity -! gflux_lnd - real , soil heat flux over land -! runoff_lnd - real , surface runoff -! drain_lnd - real , subsurface runoff -! cmm_lnd - real , surface drag wind speed for momentum -! chh_lnd - real , surface drag mass flux for heat and moisture -! zvfun_lnd - real , function of surface roughness length and green vegetation fraction -! outputs: -! sncovr1 - real , snow cover over land -! qsurf - real , specific humidity at sfc -! evap - real , evaporation from latent heat -! hflx - real , sensible heat -! ep - real , potential evaporation -! t2mmp - real , temperature at 2m -! q2mp - real , specific humidity at 2m -! gflux - real , soil heat flux over land -! runoff - real , surface runoff -! drain - real , subsurface runoff -! cmm - real , surface drag wind speed for momentum -! chh - real , surface drag mass flux for heat and moisture -! zvfun - real , function of surface roughness length and green vegetation fraction -! ==================== end of description ===================== ! -! -! - use machine , only : kind_phys - implicit none - -! --- inputs: - integer, intent(in) :: im - logical, intent(in) :: cpllnd, cpllnd2atm - logical, dimension(:), intent(in) :: flag_iter - logical, dimension(:), intent(in) :: dry - - real (kind=kind_phys), dimension(:), intent(in) :: & - & sncovr1_lnd, qsurf_lnd, evap_lnd, hflx_lnd, ep_lnd, & - & t2mmp_lnd, q2mp_lnd, gflux_lnd, runoff_lnd, drain_lnd, & - & cmm_lnd, chh_lnd, zvfun_lnd - -! --- outputs: - real (kind=kind_phys), dimension(:), intent(inout) :: & - & sncovr1, qsurf, evap, hflx, ep, t2mmp, q2mp, gflux, & - & runoff, drain, cmm, chh, zvfun -! - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - -! --- locals: - - integer :: i - - ! Initialize CCPP error handling variables - errmsg = '' - errflg = 0 -! - if (.not. cpllnd2atm) return -! - do i = 1, im - sncovr1(i) = sncovr1_lnd(i) - qsurf(i) = qsurf_lnd(i) - hflx(i) = hflx_lnd(i) - evap(i) = evap_lnd(i) - ep(i) = ep_lnd(i) - t2mmp(i) = t2mmp_lnd(i) - q2mp(i) = q2mp_lnd(i) - gflux(i) = gflux_lnd(i) - drain(i) = drain_lnd(i) - runoff(i) = runoff_lnd(i) - cmm(i) = cmm_lnd(i) - chh(i) = chh_lnd(i) - zvfun(i) = zvfun_lnd(i) - enddo - - return -!----------------------------------- - end subroutine sfc_land_run -!----------------------------------- - -!> @} - end module sfc_land From fe77e06ab907accb3b860857d050f990ba88662d Mon Sep 17 00:00:00 2001 From: Ufuk Turuncoglu Date: Mon, 8 Jan 2024 11:54:41 -0600 Subject: [PATCH 10/14] move sfc_land to new location --- physics/{ => SFC_Models/Land}/sfc_land.F90 | 0 physics/{ => SFC_Models/Land}/sfc_land.meta | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename physics/{ => SFC_Models/Land}/sfc_land.F90 (100%) rename physics/{ => SFC_Models/Land}/sfc_land.meta (100%) diff --git a/physics/sfc_land.F90 b/physics/SFC_Models/Land/sfc_land.F90 similarity index 100% rename from physics/sfc_land.F90 rename to physics/SFC_Models/Land/sfc_land.F90 diff --git a/physics/sfc_land.meta b/physics/SFC_Models/Land/sfc_land.meta similarity index 100% rename from physics/sfc_land.meta rename to physics/SFC_Models/Land/sfc_land.meta From 09b02350a7526f1eddc0885f8bcb26e9cac5c72d Mon Sep 17 00:00:00 2001 From: Ufuk Turuncoglu Date: Mon, 8 Jan 2024 11:55:57 -0600 Subject: [PATCH 11/14] fix meta file --- physics/SFC_Models/Land/sfc_land.meta | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/SFC_Models/Land/sfc_land.meta b/physics/SFC_Models/Land/sfc_land.meta index 493d2a70b..6a4bd8fbe 100644 --- a/physics/SFC_Models/Land/sfc_land.meta +++ b/physics/SFC_Models/Land/sfc_land.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = sfc_land type = scheme - dependencies = machine.F + dependencies = ../../hooks/machine.F ######################################################################## [ccpp-arg-table] From 02b3440378e5bb04ddc4ba50b44f2532eb7cab08 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Fri, 19 Jan 2024 18:24:37 +0000 Subject: [PATCH 12/14] "Supplementary physics updates for RRFS code freeze" --- physics/CONV/Grell_Freitas/cu_gf_deep.F90 | 5 +- physics/CONV/Grell_Freitas/cu_gf_driver.F90 | 8 +- physics/CONV/Grell_Freitas/cu_gf_driver.meta | 7 + physics/PBL/MYNN_EDMF/module_bl_mynn.F90 | 146 ++++++++---------- .../SFC_Models/Land/RUC/module_sf_ruclsm.F90 | 6 +- physics/smoke_dust/dep_dry_mod_emerson.F90 | 58 ++++--- physics/smoke_dust/dust_fengsha_mod.F90 | 11 +- physics/smoke_dust/module_add_emiss_burn.F90 | 75 +++++---- physics/smoke_dust/module_plumerise.F90 | 61 +++++--- physics/smoke_dust/module_smoke_plumerise.F90 | 40 ++--- physics/smoke_dust/module_wetdep_ls.F90 | 13 +- physics/smoke_dust/plume_data_mod.F90 | 51 ------ physics/smoke_dust/rrfs_smoke_config.F90 | 3 +- physics/smoke_dust/rrfs_smoke_postpbl.meta | 3 +- physics/smoke_dust/rrfs_smoke_wrapper.F90 | 145 ++++++++++------- physics/smoke_dust/rrfs_smoke_wrapper.meta | 85 +++++++++- 16 files changed, 424 insertions(+), 293 deletions(-) delete mode 100755 physics/smoke_dust/plume_data_mod.F90 diff --git a/physics/CONV/Grell_Freitas/cu_gf_deep.F90 b/physics/CONV/Grell_Freitas/cu_gf_deep.F90 index a1bca36c9..8a2c73600 100644 --- a/physics/CONV/Grell_Freitas/cu_gf_deep.F90 +++ b/physics/CONV/Grell_Freitas/cu_gf_deep.F90 @@ -142,13 +142,13 @@ subroutine cu_gf_deep_run( & !! betwee -1 and +1 ,do_capsuppress,cap_suppress_j & ! ,k22 & ! - ,jmin,tropics) ! + ,jmin,kdt,tropics) ! implicit none integer & ,intent (in ) :: & - nranflag,itf,ktf,its,ite, kts,kte,ipr,imid + nranflag,itf,ktf,its,ite, kts,kte,ipr,imid,kdt integer, intent (in ) :: & ichoice,nchem real(kind=kind_phys), dimension (its:ite,4) & @@ -591,6 +591,7 @@ subroutine cu_gf_deep_run( & sig(i)=(1.-frh)**2 !frh_out(i) = frh if(forcing(i,7).eq.0.)sig(i)=1. + if(kdt.le.(3600./dtime))sig(i)=1. frh_out(i) = frh*sig(i) enddo !$acc end kernels diff --git a/physics/CONV/Grell_Freitas/cu_gf_driver.F90 b/physics/CONV/Grell_Freitas/cu_gf_driver.F90 index 92f8760b0..54a23ca74 100644 --- a/physics/CONV/Grell_Freitas/cu_gf_driver.F90 +++ b/physics/CONV/Grell_Freitas/cu_gf_driver.F90 @@ -68,7 +68,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& dfi_radar_max_intervals,ldiag3d,qci_conv,do_cap_suppress, & maxupmf,maxMF,do_mynnedmf,ichoice_in,ichoicem_in,ichoice_s_in, & spp_cu_deep,spp_wts_cu_deep,nchem,chem3d,fscav,wetdpc_deep, & - do_smoke_transport,errmsg,errflg) + do_smoke_transport,kdt,errmsg,errflg) !------------------------------------------------------------- implicit none integer, parameter :: maxiens=1 @@ -95,7 +95,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& integer :: ishallow_g3 ! depend on imfshalcnv !------------------------------------------------------------- integer :: its,ite, jts,jte, kts,kte - integer, intent(in ) :: im,km,ntracer, nchem + integer, intent(in ) :: im,km,ntracer,nchem,kdt integer, intent(in ) :: ichoice_in,ichoicem_in,ichoice_s_in logical, intent(in ) :: flag_init, flag_restart, do_mynnedmf logical, intent(in ) :: flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend @@ -766,7 +766,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& ! betwee -1 and +1 ,do_cap_suppress_here,cap_suppress_j & ,k22m & - ,jminm,tropics) + ,jminm,kdt,tropics) !$acc kernels do i=its,itf do k=kts,ktf @@ -853,7 +853,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& ! betwee -1 and +1 ,do_cap_suppress_here,cap_suppress_j & ,k22 & - ,jmin,tropics) + ,jmin,kdt,tropics) jpr=0 ipr=0 !$acc kernels diff --git a/physics/CONV/Grell_Freitas/cu_gf_driver.meta b/physics/CONV/Grell_Freitas/cu_gf_driver.meta index fe9b4c375..d0b661fd8 100644 --- a/physics/CONV/Grell_Freitas/cu_gf_driver.meta +++ b/physics/CONV/Grell_Freitas/cu_gf_driver.meta @@ -651,6 +651,13 @@ type = real kind = kind_phys intent = inout +[kdt] + standard_name = index_of_timestep + long_name = current forecast iteration + units = index + dimensions = () + type = integer + intent = in [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 b/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 index 6840f80bf..cc7a47ce6 100644 --- a/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 +++ b/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 @@ -2001,9 +2001,9 @@ SUBROUTINE mym_length ( & uonset= 15. wt_u = (1.0 - min(max(ugrid - uonset, 0.0)/30.0, 0.5)) cns = 2.7 !was 3.5 - alp1 = 0.22 + alp1 = 0.23 alp2 = 0.3 - alp3 = 2.0 * wt_u !taper off bouyancy enhancement in shear-driven pbls + alp3 = 2.5 * wt_u !taper off bouyancy enhancement in shear-driven pbls alp4 = 5.0 alp5 = 0.3 alp6 = 50. @@ -2059,12 +2059,12 @@ SUBROUTINE mym_length ( & ! ** Length scale limited by the buoyancy effect ** IF ( dtv(k) .GT. 0.0 ) THEN - bv = max( sqrt( gtr*dtv(k) ), 0.001) + bv = max( sqrt( gtr*dtv(k) ), 0.0001) elb = MAX(alp2*qkw(k), & & alp6*edmf_a1(k-1)*edmf_w1(k-1)) / bv & & *( 1.0 + alp3*SQRT( vsc/(bv*elt) ) ) elb = MIN(elb, zwk) - elf = 0.80 * qkw(k)/bv + elf = 1.0 * qkw(k)/bv elBLavg(k) = MAX(elBLavg(k), alp6*edmf_a1(k-1)*edmf_w1(k-1)/bv) ELSE elb = 1.0e10 @@ -2084,8 +2084,10 @@ SUBROUTINE mym_length ( & !add blending to use BouLac mixing length in free atmos; !defined relative to the PBLH (zi) + transition layer (h1) !el(k) = MIN(elb/( elb/elt+elb/els+1.0 ),elf) - !try squared-blending - el(k) = SQRT( els**2/(1. + (els**2/elt**2) +(els**2/elb**2))) + !try squared-blending - but take out elb (makes it underdiffusive) + !el(k) = SQRT( els**2/(1. + (els**2/elt**2) +(els**2/elb**2))) + el(k) = sqrt( els**2/(1. + (els**2/elt**2))) + el(k) = min(el(k), elb) el(k) = MIN (el(k), elf) el(k) = el(k)*(1.-wt) + alp5*elBLavg(k)*wt @@ -3633,13 +3635,13 @@ SUBROUTINE mym_condensation (kts,kte, & real(kind_phys):: qsl,esat,qsat,dqsl,cld0,q1k,qlk,eq1,qll, & &q2p,pt,rac,qt,t,xl,rsl,cpm,Fng,qww,alpha,beta,bb, & - &ls,wt,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, & + &ls,wt,wt2,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, & &qmq,qsat_tk,q1_rh,rh_hack,dzm1,zsl,maxqc real(kind_phys), parameter :: qpct_sfc=0.025 real(kind_phys), parameter :: qpct_pbl=0.030 real(kind_phys), parameter :: qpct_trp=0.040 real(kind_phys), parameter :: rhcrit =0.83 !for cloudpdf = 2 - real(kind_phys), parameter :: rhmax =1.01 !for cloudpdf = 2 + real(kind_phys), parameter :: rhmax =1.02 !for cloudpdf = 2 integer :: i,j,k real(kind_phys):: erf @@ -3864,25 +3866,18 @@ SUBROUTINE mym_condensation (kts,kte, & !Add condition for falling/settling into low-RH layers, so at least !some cloud fraction is applied for all qc, qs, and qi. rh_hack= rh(k) + wt2 = min(max( zagl - pblh2, 0.0 )/300., 1.0) !ensure adequate RH & q1 when qi is at least 1e-9 (above the PBLH) - if (qi(k)>1.e-9 .and. zagl .gt. pblh2) then - rh_hack =min(rhmax, rhcrit + 0.07*(9.0 + log10(qi(k)))) + if ((qi(k)+qs(k))>1.e-9 .and. (zagl .gt. pblh2)) then + rh_hack =min(rhmax, rhcrit + wt2*0.045*(9.0 + log10(qi(k)+qs(k)))) rh(k) =max(rh(k), rh_hack) !add rh-based q1 q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit) q1(k) =max(q1_rh, q1(k) ) endif - !ensure adequate RH & q1 when qc is at least 1e-6 - if (qc(k)>1.e-6) then - rh_hack =min(rhmax, rhcrit + 0.09*(6.0 + log10(qc(k)))) - rh(k) =max(rh(k), rh_hack) - !add rh-based q1 - q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit) - q1(k) =max(q1_rh, q1(k) ) - endif - !ensure adequate RH & q1 when qs is at least 1e-8 (above the PBLH) - if (qs(k)>1.e-8 .and. zagl .gt. pblh2) then - rh_hack =min(rhmax, rhcrit + 0.07*(8.0 + log10(qs(k)))) + !ensure adequate rh & q1 when qc is at least 1e-6 (above the PBLH) + if (qc(k)>1.e-6 .and. (zagl .gt. pblh2)) then + rh_hack =min(rhmax, rhcrit + wt2*0.08*(6.0 + log10(qc(k)))) rh(k) =max(rh(k), rh_hack) !add rh-based q1 q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit) @@ -3994,7 +3989,7 @@ SUBROUTINE mym_condensation (kts,kte, & fac_damp = min(zagl * 0.0025, 1.0) !cld_factor = 1.0 + fac_damp*MAX(0.0, ( RH(k) - 0.75 ) / 0.26 )**1.9 !HRRRv4 !cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.25 )**2, 0.3) - cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.145)**2, 0.35) + cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.145)**2, 0.37) cldfra_bl1D(K) = min( 1., cld_factor*cldfra_bl1D(K) ) enddo @@ -4181,38 +4176,33 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & k=kts -!original approach (drag in b-vector): -! a(1)=0. -! b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff -! c(1)=-dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff -! d(1)=u(k) + dtz(k)*uoce*ust**2/wspd - dtz(k)*s_awu(k+1)*onoff + & -! sub_u(k)*delt + det_u(k)*delt - !rho-weighted (drag in b-vector): a(k)= -dtz(k)*kmdz(k)*rhoinv(k) - b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k) & - & - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k) & + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) & - & - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff - d(k)=u(k) + dtz(k)*uoce*ust**2/wspd - dtz(k)*s_awu(k+1)*onoff - & - & dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff + sub_u(k)*delt + det_u(k)*delt - -!rho-weighted with drag term moved out of b-array -! a(k)= -dtz(k)*kmdz(k)*rhoinv(k) -! b(k)=1.+dtz(k)*(kmdz(k+1))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff -! c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff -! d(k)=u(k)*(1.-ust**2/wspd*dtz(k)*rhosfc/rho(k)) + dtz(k)*uoce*ust**2/wspd - & -! !!!d(k)=u(k)*(1.-ust**2/wspd*dtz(k)) + dtz(k)*uoce*ust**2/wspd - & -! & dtz(k)*rhoinv(k)*s_awu(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff + sub_u(k)*delt + det_u(k)*delt - - DO k=kts+1,kte-1 - a(k)= -dtz(k)*kmdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff - b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + & - & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff - c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff - d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff + dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff + & - & sub_u(k)*delt + det_u(k)*delt - ENDDO + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + d(k)=u(k) + dtz(k)*uoce*ust**2/wspd & + & - dtz(k)*rhoinv(k)*s_awu(k+1)*onoff & + & + dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff & + & + sub_u(k)*delt + det_u(k)*delt + + do k=kts+1,kte-1 + a(k)= -dtz(k)*kmdz(k)*rhoinv(k) & + & + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff & + & + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff + b(k)=1.+ dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) & + & + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff & + & + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff + c(k)= - dtz(k)*kmdz(k+1)*rhoinv(k) & + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff & + & - dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff & + & + sub_u(k)*delt + det_u(k)*delt + enddo !! no flux at the top ! a(kte)=-1. @@ -4247,37 +4237,33 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & k=kts -!original approach (drag in b-vector): -! a(1)=0. -! b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff -! c(1)= - dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff -! d(1)=v(k) + dtz(k)*voce*ust**2/wspd - dtz(k)*s_awv(k+1)*onoff + & -! sub_v(k)*delt + det_v(k)*delt - !rho-weighted (drag in b-vector): a(k)= -dtz(k)*kmdz(k)*rhoinv(k) - b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k) & - & - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff - c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff - d(k)=v(k) + dtz(k)*voce*ust**2/wspd - dtz(k)*s_awv(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff + & - & sub_v(k)*delt + det_v(k)*delt - -!rho-weighted with drag term moved out of b-array -! a(k)= -dtz(k)*kmdz(k)*rhoinv(k) -! b(k)=1.+dtz(k)*(kmdz(k+1))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff -! c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff -! d(k)=v(k)*(1.-ust**2/wspd*dtz(k)*rhosfc/rho(k)) + dtz(k)*voce*ust**2/wspd - & -! !!!d(k)=v(k)*(1.-ust**2/wspd*dtz(k)) + dtz(k)*voce*ust**2/wspd - & -! & dtz(k)*rhoinv(k)*s_awv(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff + sub_v(k)*delt + det_v(k)*delt - - DO k=kts+1,kte-1 - a(k)= -dtz(k)*kmdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff - b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + & - & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff - c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff - d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff + dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff + & - & sub_v(k)*delt + det_v(k)*delt - ENDDO + b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k) & + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) & + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + d(k)=v(k) + dtz(k)*voce*ust**2/wspd & + & - dtz(k)*rhoinv(k)*s_awv(k+1)*onoff & + & + dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff & + & + sub_v(k)*delt + det_v(k)*delt + + do k=kts+1,kte-1 + a(k)= -dtz(k)*kmdz(k)*rhoinv(k) & + & + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff & + & + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff + b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) & + & + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff & + & + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff + c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) & + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff & + & - dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff & + & + sub_v(k)*delt + det_v(k)*delt + enddo !! no flux at the top ! a(kte)=-1. diff --git a/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 b/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 index b15592052..2d01f96c9 100644 --- a/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 +++ b/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 @@ -1687,7 +1687,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia endif if(newsn > zero ) then - SNOWFRACnewsn=MIN(one,SNHEI/SNHEI_CRIT_newsn) + SNOWFRACnewsn=MIN(one,snowfallac*1.e-3_kind_phys/SNHEI_CRIT_newsn) endif !-- due to steep slopes and blown snow, limit snow fraction in the @@ -1700,7 +1700,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia if(snowfrac < 0.75_kind_phys) snow_mosaic = one KEEP_SNOW_ALBEDO = zero - IF (NEWSN > zero .and. snowfracnewsn > 0.99_kind_phys .and. rhosnfall < 450._kind_phys) THEN + IF (snowfracnewsn > 0.99_kind_phys .and. rhosnfall < 450._kind_phys) THEN ! new snow KEEP_SNOW_ALBEDO = one ! turn off separate treatment of snow covered and snow-free portions of the grid cell @@ -1735,7 +1735,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia ! hwlps with these biases.. if( snow_mosaic == one) then ALBsn=alb_snow - if(newsn > zero .and. KEEP_SNOW_ALBEDO > 0.9_kind_phys .and. albsn < 0.4_kind_phys) then + if(KEEP_SNOW_ALBEDO > 0.9_kind_phys .and. albsn < 0.4_kind_phys) then !-- Albedo correction with fresh snow and deep snow pack !-- will reduce warm bias in western Canada !-- and US West coast, where max snow albedo is low (0.3-0.5). diff --git a/physics/smoke_dust/dep_dry_mod_emerson.F90 b/physics/smoke_dust/dep_dry_mod_emerson.F90 index 76fdc2411..771801c44 100755 --- a/physics/smoke_dust/dep_dry_mod_emerson.F90 +++ b/physics/smoke_dust/dep_dry_mod_emerson.F90 @@ -9,7 +9,7 @@ module dep_dry_emerson_mod use machine , only : kind_phys use dep_data_mod ! JLS - use rrfs_smoke_config, only : num_chem, p_smoke, p_dust_1, p_coarse_pm + use rrfs_smoke_config, only : num_chem, p_smoke, p_dust_1, p_coarse_pm, n_dbg_lines implicit none @@ -23,7 +23,7 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & settling_flag,drydep_flux,settling_flux,dbg_opt, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte ) + its,ite, jts,jte, kts,kte, curr_secs, mpiid, xlat, xlong ) ! ! compute dry deposition velocity for aerosol particles ! Based on Emerson et al. (2020), PNAS, @@ -37,6 +37,9 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte + + REAL(kind_phys) :: curr_secs + REAL(kind_phys), DIMENSION( ims:ime , jms:jme ) , & INTENT(IN) :: ustar, rmol, znt, snowh REAL(kind_phys), DIMENSION( ims:ime , kms:kme , jms:jme ), & @@ -80,6 +83,9 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & real(kind_phys), dimension( kts:kte, ndvel ) :: cblk_col, vg_col integer, dimension(ndvel) :: ndt_settl integer :: i, j, k, ntdt, nv + integer :: icall=0 + integer, INTENT(IN) :: mpiid + real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: xlat,xlong ! chem pointers (p_*) are not sequentially numbered, need to define so nv loops work integer, dimension(ndvel) :: chem_pointers !> -- Gas constant @@ -87,11 +93,15 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & chem_pointers(1) = p_smoke chem_pointers(2) = p_dust_1 chem_pointers(3) = p_coarse_pm - + growth_fac = 1.0 conver=1.e-9 converi=1.e9 + if (mod(int(curr_secs),1800) .eq. 0) then + icall = 0 + endif + do j = jts, jte do i = its, ite aer_res(i,j) = 0.0 @@ -116,7 +126,7 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & aer_res(i,j) = max(aer_res(i,j)/100._kind_phys,0._kind_phys) ! Air kinematic viscosity (cm^2/s) airkinvisc = ( 1.8325e-4 * ( 416.16 / ( t_phy(i,k,j) + 120.0 ) ) * & - ( ( t_phy(i,k,j) / 296.16 )**1.5 ) ) / ( rho_phy(i,k,j) / 28.966e3 ) ! Convert density to mol/cm^3 + ( ( t_phy(i,k,j) / 296.16 )**1.5 ) ) / ( rho_phy(i,k,j) / 1.e3 ) ! Convert density to mol/cm^3 ! Air molecular freepath (cm) ! Check against XLM from above freepath = 7.39758e-4 * airkinvisc / sqrt( t_phy(i,k,j) ) do nv = 1, ndvel @@ -141,11 +151,11 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & amu_corrected = amu / Cc ! Gravitational Settling vg = aerodens * dp * dp * gravity * 100. * Cc / & ! Convert gravity to cm/s^2 - ( 18. * airkinvisc * ( rho_phy(i,k,j) / 28.966e3 ) ) ! Convert density to mol/cm^3 + ( 18. * airkinvisc * ( rho_phy(i,k,j) / 1.e3 ) ) ! Convert density to mol/cm^3 ! -- Rest of loop for the surface when deposition velocity needs to be cacluated if ( k == kts ) then ! Brownian Diffusion - DDp = ( boltzmann * t_phy(i,k,j) ) * Cc / (3. * pi * airkinvisc * ( rho_phy(i,k,j) / 28.966e3 ) * dp) ! Convert density to mol/cm^3 + DDp = ( boltzmann * t_phy(i,k,j) ) * Cc / (3. * pi * airkinvisc * ( rho_phy(i,k,j) / 1.e3 ) * dp) ! Convert density to mol/cm^3 ! Schmit number Sc = airkinvisc / DDp ! Brownian Diffusion @@ -179,13 +189,17 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & Rs = 1. / ( ( ustar(i,j) * 100.) * ( Eb + Eim + Ein) * eps0 ) ! Convert ustar to cm/s ! Compute final ddvel = aer_res + RS, set max at max_dep_vel in dep_data_mod.F[ m/s] ! The /100. term converts from cm/s to m/s, required for MYNN. - ddvel(i,j,nv) = min( (1. / (aer_res(i,j) + Rs ))/100., max_dep_vel) - if ( dbg_opt ) then - WRITE(6,*) 'dry_dep_mod_emerson: i,j,nv',i,j,nv - WRITE(6,*) 'dry_dep_mod_emerson: deposition velocity (m/s) ',ddvel(i,j,nv) + if ( settling_flag == 1 ) then + ddvel(i,j,nv) = max(min( ( vg + 1./(aer_res(i,j)+Rs) )/100., max_dep_vel),0._kind_phys) + else + ddvel(i,j,nv) = max(min( ( 1./(aer_res(i,j)+Rs) )/100., max_dep_vel),0._kind_phys) endif - drydep_flux(i,j,nv) = chem(i,kts,j,chem_pointers(nv))*p_phy(i,kts,j) / & - (RSI*t_phy(i,kts,j))*ddvel(i,j,nv)*dt*1.E-6 + if ( dbg_opt .and. (icall .le. n_dbg_lines) ) then + WRITE(1000+mpiid,*) 'dry_dep_mod_emer:xlat,xlong,curr_secs,nv',xlat(i,j),xlong(i,j),int(curr_secs),nv + WRITE(1000+mpiid,*) 'dry_dep_mod_emer:xlat,xlong,curr_secs,deposition velocity (m/s)',xlat(i,j),xlong(i,j),int(curr_secs),ddvel(i,j,nv) + icall = icall + 1 + endif + drydep_flux(i,j,nv) = chem(i,kts,j,chem_pointers(nv))*rho_phy(i,k,j)*ddvel(i,j,nv)/100.0*dt endif ! k == kts vgrav(i,k,j,nv) = vg ! Fill column variables @@ -220,25 +234,25 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & IF (ndt_settl(nv) > 12) ndt_settl(nv) = 12 dt_settl(nv) = REAL(ntdt,kind=kind_phys) /REAL(ndt_settl(nv),kind=kind_phys) enddo - do nv = 1, ndvel - chem_before(nv) = 0._kind_phys - do k = kts, kte - chem_before(nv) = chem_before(nv) + (cblk_col(k,nv) * rho_phy(i,k,j) * delz(i,k,j) ) ! ug/m2 - enddo - enddo + !do nv = 1, ndvel + ! chem_before(nv) = 0._kind_phys + ! do k = kts, kte + ! chem_before(nv) = chem_before(nv) + (cblk_col(k,nv) * rho_phy(i,k,j) * delz(i,k,j) ) ! ug/m2 + ! enddo + !enddo ! Perform gravitational settling if desired if ( settling_flag == 1 ) then call particle_settling(cblk_col,rho_col,delz_col,vg_col,dt_settl,ndt_settl,ndvel,kts,kte) endif ! Put cblk back into chem array do nv= 1, ndvel - chem_after(nv) = 0._kind_phys - settling_flux(i,j,nv) = 0._kind_phys + !chem_after(nv) = 0._kind_phys + !settling_flux(i,j,nv) = 0._kind_phys do k = kts, kte chem(i,k,j,chem_pointers(nv)) = cblk_col(k,nv) - chem_after(nv) = chem_after(nv) + (cblk_col(k,nv) * rho_phy(i,k,j) * delz(i,k,j) ) ! ug/m2 + !chem_after(nv) = chem_after(nv) + (cblk_col(k,nv) * rho_phy(i,k,j) * delz(i,k,j) ) ! ug/m2 enddo ! k - settling_flux(i,j,nv) = settling_flux(i,j,nv) + (chem_before(nv) - chem_after(nv)) ! ug/m2 + !settling_flux(i,j,nv) = settling_flux(i,j,nv) + (chem_before(nv) - chem_after(nv)) ! ug/m2 enddo ! nv end do ! j end do ! i diff --git a/physics/smoke_dust/dust_fengsha_mod.F90 b/physics/smoke_dust/dust_fengsha_mod.F90 index 54e66712d..6ec8f8d4a 100755 --- a/physics/smoke_dust/dust_fengsha_mod.F90 +++ b/physics/smoke_dust/dust_fengsha_mod.F90 @@ -21,8 +21,8 @@ module dust_fengsha_mod contains subroutine gocart_dust_fengsha_driver(dt, & - chem,rho_phy,smois,p8w,ssm, & - isltyp,snowh,xland,area,g,emis_dust, & + chem,rho_phy,smois,stemp,p8w,ssm, & + isltyp,snowh,xland,area,g,emis_dust, & ust,znt,clay,sand,rdrag,uthr, & num_emis_dust,num_chem,num_soil_layers, & ids,ide, jds,jde, kds,kde, & @@ -54,7 +54,7 @@ subroutine gocart_dust_fengsha_driver(dt, & REAL(kind_phys), DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN) :: rho_phy REAL(kind_phys), DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), INTENT(INOUT) :: chem REAL(kind_phys), DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL, INTENT(INOUT) :: emis_dust - REAL(kind_phys), DIMENSION( ims:ime, num_soil_layers, jms:jme ), INTENT(IN) :: smois + REAL(kind_phys), DIMENSION( ims:ime, num_soil_layers, jms:jme ), INTENT(IN) :: smois, stemp !0d input variables REAL(kind_phys), INTENT(IN) :: dt ! time step @@ -146,6 +146,11 @@ subroutine gocart_dust_fengsha_driver(dt, & ilwi = 0 endif + ! Don't emit over frozen soil + if (stemp(i,1,j) < 268.0) then ! -5C + ilwi = 0 + endif + ! Do not allow areas with bedrock, lava, or land-ice to loft IF (isltyp(i,j) .eq. 15 .or. isltyp(i,j) .eq. 16. .or. & diff --git a/physics/smoke_dust/module_add_emiss_burn.F90 b/physics/smoke_dust/module_add_emiss_burn.F90 index 95005b973..0a22fcfd7 100755 --- a/physics/smoke_dust/module_add_emiss_burn.F90 +++ b/physics/smoke_dust/module_add_emiss_burn.F90 @@ -9,11 +9,11 @@ module module_add_emiss_burn subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & chem,julday,gmt,xlat,xlong, & fire_end_hr, peak_hr,time_int, & - coef_bb_dc, fhist, hwp, hwp_prevd, & + coef_bb_dc, fire_hist, hwp, hwp_prevd, & swdown,ebb_dcycle, ebu_in, ebu,fire_type,& ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte ) + its,ite, jts,jte, kts,kte,mpiid ) IMPLICIT NONE @@ -22,6 +22,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte + INTEGER, INTENT(IN) :: mpiid real(kind_phys), DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem ! shall we set num_chem=1 here? @@ -29,7 +30,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & INTENT(INOUT ) :: ebu real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: xlat,xlong, swdown - real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: hwp, peak_hr, fire_end_hr, ebu_in !, vfrac + real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: hwp, peak_hr, fire_end_hr, ebu_in !RAR: Shall we make fire_end integer? real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: coef_bb_dc ! RAR: real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: hwp_prevd @@ -38,17 +39,17 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & real(kind_phys), INTENT(IN) :: time_int,pi ! RAR: time in seconds since start of simulation INTEGER, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: fire_type integer, INTENT(IN) :: ebb_dcycle ! RAR: this is going to be namelist dependent, ebb_dcycle=means - real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fhist + real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fire_hist !>--local integer :: i,j,k,n,m + integer :: icall=0 real(kind_phys) :: conv_rho, conv, dm_smoke, dc_hwp, dc_gp, dc_fn !daero_num_wfa, daero_num_ifa !, lu_sum1_5, lu_sum12_14 INTEGER, PARAMETER :: kfire_max=51 ! max vertical level for BB plume rise - ! real, parameter :: cx = 2.184936 * 3600., rinti = 2.1813936e-8 , ax = 2000.6038 - ! bx_bburn = 20.041288 * 3600., RAR: this depends on the vegetation class, location (local time) etc. real(kind_phys) :: timeq, fire_age, age_hr, dt1,dt2,dtm ! For BB emis. diurnal cycle calculation - ! For Gaussian diurnal cycle + +! For Gaussian diurnal cycle real(kind_phys), PARAMETER :: sc_factor=1. ! to scale up the wildfire emissions, TBD later real(kind_phys), PARAMETER :: rinti=2.1813936e-8, ax2=3400., const2=130., & coef2=10.6712963e-4, cx2=7200., timeq_max=3600.*24. @@ -90,32 +91,39 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & ! Constants for the fire diurnal cycle calculation do j=jts,jte do i=its,ite - fire_age= time_int + (fire_end_hr(i,j))*3600. + fire_age= time_int + (fire_end_hr(i,j)-1._kind_phys)*3600._kind_phys !One hour delay is due to the latency of the RAVE files + fire_age= MAX(0._kind_phys,fire_age) SELECT CASE ( fire_type(i,j) ) !Ag, urban fires, bare land etc. CASE (1) ! these fires will have exponentially decreasing diurnal cycle, - ! We assume 1hr latency in ingesting the sat. data - coef_bb_dc(i,j) = 1._kind_phys/((2*pi)**0.5_kind_phys * sigma_fire_dur(1) *fire_age) * & - exp(- ( log(fire_age) - avg_fire_dur(1))**2 /(2*sigma_fire_dur(1)**2 )) + coef_bb_dc(i,j) = 1._kind_phys/((2._kind_phys*pi)**0.5_kind_phys * sigma_fire_dur(1) *fire_age) * & + exp(- ( log(fire_age) - avg_fire_dur(1))**2 /(2._kind_phys*sigma_fire_dur(1)**2 )) + + ! IF ( dbg_opt .AND. time_int<5000.) then + ! WRITE(6,*) 'i,j,peak_hr(i,j) ',i,j,peak_hr(i,j) + ! WRITE(6,*) 'coef_bb_dc(i,j) ',coef_bb_dc(i,j) + ! END IF + CASE (3) age_hr= fire_age/3600._kind_phys - IF (swdown(i,j)<.1 .AND. age_hr> 12. .AND. fhist(i,j)>0.75) THEN - fhist(i,j)= 0.75_kind_phys + IF (swdown(i,j)<.1 .AND. age_hr> 12. .AND. fire_hist(i,j)>0.75) THEN + fire_hist(i,j)= 0.75_kind_phys ENDIF - IF (swdown(i,j)<.1 .AND. age_hr> 24. .AND. fhist(i,j)>0.5) THEN - fhist(i,j)= 0.5_kind_phys + IF (swdown(i,j)<.1 .AND. age_hr> 24. .AND. fire_hist(i,j)>0.5) THEN + fire_hist(i,j)= 0.5_kind_phys ENDIF - IF (swdown(i,j)<.1 .AND. age_hr> 48. .AND. fhist(i,j)>0.25) THEN - fhist(i,j)= 0.25_kind_phys + IF (swdown(i,j)<.1 .AND. age_hr> 48. .AND. fire_hist(i,j)>0.25) THEN + fire_hist(i,j)= 0.25_kind_phys ENDIF ! this is based on hwp, hourly or instantenous TBD - dc_hwp= ebu_in(i,j)* hwp(i,j)/ MAX(1._kind_phys,hwp_prevd(i,j)) + dc_hwp= hwp(i,j)/ MAX(5._kind_phys,hwp_prevd(i,j)) dc_hwp= MAX(0._kind_phys,dc_hwp) + dc_hwp= MIN(25._kind_phys,dc_hwp) - !coef_bb_dc(i,j)= sc_factor* fhist(i,j)* rate_ebb2(i,j)* (1. + log( + !coef_bb_dc(i,j)= sc_factor* fire_hist(i,j)* rate_ebb2(i,j)* (1. + log( !hwp_(i,j)/ hwp_day_avg(i,j))) ! RAR: Gaussian profile for wildfires @@ -125,17 +133,30 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & dc_gp = rinti*( ax2 * exp(- dtm**2/(2._kind_phys*cx2**2) ) + const2 - coef2*timeq ) dc_gp = MAX(0._kind_phys,dc_gp) - dc_fn = MAX(dc_hwp/dc_gp,3._kind_phys) - coef_bb_dc(i,j) = fhist(i,j)* dc_fn + dc_fn = MIN(dc_hwp/dc_gp,3._kind_phys) + !coef_bb_dc(i,j) = fire_hist(i,j)* dc_fn + coef_bb_dc(i,j) = fire_hist(i,j)* dc_hwp + + ! IF ( dbg_opt .AND. time_int<5000.) then + ! WRITE(6,*) 'i,j,fire_hist(i,j),peak_hr(i,j) ', i,j,fire_hist(i,j),peak_hr(i,j) + ! WRITE(6,*) 'dc_gp,dc_hwp,dc_fn,coef_bb_dc(i,j) ',dc_gp,dc_hwp,dc_fn,coef_bb_dc(i,j) + ! END IF + CASE DEFAULT END SELECT enddo enddo endif + if (mod(int(time_int),1800) .eq. 0) then + icall = 0 + endif + do j=jts,jte do i=its,ite do k=kts,kfire_max + if (ebu(i,k,j)<0.001_kind_phys) cycle + if (ebb_dcycle==1) then conv= dtstep/(rho_phy(i,k,j)* dz8w(i,k,j)) elseif (ebb_dcycle==2) then @@ -143,14 +164,14 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & endif dm_smoke= conv*ebu(i,k,j) chem(i,k,j,p_smoke) = chem(i,k,j,p_smoke) + dm_smoke - chem(i,k,j,p_smoke) = MIN(chem(i,k,j,p_smoke),5.e+3) + chem(i,k,j,p_smoke) = MIN(MAX(chem(i,k,j,p_smoke),0._kind_phys),5.e+3_kind_phys) - if ( dbg_opt .and. (k==kts .OR. k==kfire_max) ) then - WRITE(6,*) 'add_emiss_burn: i,j,k ',i,j,k - WRITE(6,*) 'add_emiss_burn: rho_phy(i,k,j),dz8w(i,k,j),conv',rho_phy(i,k,j),dz8w(i,k,j),conv - WRITE(6,*) 'add_emiss_burn: ebu(i,k,j),dm_smoke ', ebu(i,k,j),dm_smoke - endif + if ( dbg_opt .and. (k==kts .OR. k==kfire_max) .and. (icall .le. n_dbg_lines) ) then + WRITE(1000+mpiid,*) 'add_emiss_burn:xlat,xlong,curr_secs,fire_type,fire_hist,peak_hr', xlat(i,j),xlong(i,j),int(time_int),fire_type(i,j),fire_hist(i,j),peak_hr(i,j) + WRITE(1000+mpiid,*) 'add_emiss_burn:xlat,xlong,curr_secs,coef_bb_dc,ebu',xlat(i,j),xlong(i,j),int(time_int),coef_bb_dc(i,j),ebu(i,k,j) + endif enddo + icall = icall + 1 enddo enddo diff --git a/physics/smoke_dust/module_plumerise.F90 b/physics/smoke_dust/module_plumerise.F90 index 8a1d6ab25..5f7ef2a0e 100755 --- a/physics/smoke_dust/module_plumerise.F90 +++ b/physics/smoke_dust/module_plumerise.F90 @@ -24,10 +24,11 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & kpbl_thetav, & ! SRB: added kpbl_thetav ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte, errmsg, errflg) + its,ite, jts,jte, kts,kte, errmsg, errflg,curr_secs, & + xlat, xlong , uspdavg2, hpbl_thetav2, mpiid) use rrfs_smoke_config - use plume_data_mod + !use plume_data_mod USE module_zero_plumegen_coms USE module_smoke_plumerise IMPLICIT NONE @@ -40,6 +41,8 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & real(kind=kind_phys), DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: frp_inst ! RAR: FRP array + real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: xlat,xlong ! SRB + real(kind_phys), DIMENSION(ims:ime, jms:jme), INTENT(IN) :: kpbl_thetav ! SRB character(*), intent(inout) :: errmsg @@ -47,6 +50,7 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte + real(kind_phys) :: curr_secs INTEGER, INTENT(IN ) :: wind_eff_opt real(kind=kind_phys), DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: ebu real(kind=kind_phys), INTENT(IN ) :: g, con_cp, con_rd @@ -57,23 +61,32 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & ! Local variables... INTEGER :: nv, i, j, k, kp1, kp2 + INTEGER :: icall=0 INTEGER, DIMENSION(ims:ime, jms:jme), INTENT (OUT) :: k_min, k_max ! Min and max ver. levels for BB injection spread + REAL, DIMENSION(ims:ime, jms:jme), INTENT (OUT) :: uspdavg2, hpbl_thetav2 ! SRB real(kind_phys), dimension (kte) :: u_in ,v_in ,w_in ,theta_in ,pi_in, rho_phyin ,qv_in ,zmid, z_lev, uspd ! SRB real(kind=kind_phys) :: dz_plume, cpor, con_rocp, uspdavg ! SRB +! MPI variables + INTEGER, INTENT(IN) :: mpiid + cpor =con_cp/con_rd con_rocp=con_rd/con_cp - IF ( dbg_opt ) then - WRITE(*,*) 'module_plumerise: its,ite,jts,jte ', its,ite,jts,jte - WRITE(*,*) 'module_plumerise: ims,ime,jms,jme ', ims,ime,jms,jme - WRITE(*,*) 'module_plumerise: maxval(ebu(:,kts,:)) ', maxval(ebu(:,kts,:)) + if (mod(int(curr_secs),1800) .eq. 0) then + icall = 0 + endif + + IF ( dbg_opt .and. icall .le. n_dbg_lines) then + WRITE(1000+mpiid,*) 'module_plumerise: its,ite,jts,jte ', its,ite,jts,jte + WRITE(1000+mpiid,*) 'module_plumerise: ims,ime,jms,jme ', ims,ime,jms,jme + WRITE(1000+mpiid,*) 'module_plumerise: maxval(ebu(:,kts,:)) ', maxval(ebu(:,kts,:)) END IF ! RAR: setting to zero the ebu emissions at the levels k>1, this is necessary when the plumerise is called, so the emissions at k>1 are updated !do nv=1,num_ebu do j=jts,jte - do k=kts+1,kte + do k=kts,kte do i=its,ite ebu(i,k,j)=0. enddo @@ -112,12 +125,10 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & uspd(k)= wind_phy(i,k,j) ! SRB enddo - IF (dbg_opt) then - WRITE(*,*) 'module_plumerise: i,j ',i,j - WRITE(*,*) 'module_plumerise: frp_inst(i,j) ',frp_inst(i,j) - WRITE(*,*) 'module_plumerise: ebu(i,kts,j) ',ebu(i,kts,j) - WRITE(*,*) 'module_plumerise: u_in(10),v_in(10),w_in(kte),qv_in(10),pi_in(10) ',u_in(10),v_in(10),w_in(kte),qv_in(10),pi_in(10) - WRITE(*,*) 'module_plumerise: zmid(kte),z_lev(kte),rho_phyin(kte),theta_in(kte) ',zmid(kte),z_lev(kte),rho_phyin(kte),theta_in(kte) + + IF (dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst(i,j) .ge. frp_threshold) ) then + WRITE(1000+mpiid,*) 'module_plumerise_before:xlat,xlong,curr_secs,ebu(kts),frp_inst',xlat(i,j), xlong(i,j), int(curr_secs),ebu(i,kts,j),frp_inst(i,j) + WRITE(1000+mpiid,*) 'module_plumerise_before:xlat,xlong,curr_secs,u(10),v(10),w(10),qv(10)',xlat(i,j), xlong(i,j),int(curr_secs), u_in(10),v_in(10),w_in(kte),qv_in(10) END IF ! RAR: the plume rise calculation step: @@ -127,7 +138,8 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & wind_eff_opt, & frp_inst(i,j), k_min(i,j), & k_max(i,j), dbg_opt, g, con_cp, & - con_rd, cpor, errmsg, errflg ) + con_rd, cpor, errmsg, errflg, & + icall, mpiid, xlat(i,j), xlong(i,j), curr_secs ) if(errflg/=0) return kp1= k_min(i,j) @@ -136,9 +148,13 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & ! SRB: Adding condition for overwriting plumerise levels uspdavg=SUM(uspd(kts:kpbl_thetav(i,j)))/kpbl_thetav(i,j) !Average wind speed within the boundary layer + +! SRB: Adding output + uspdavg2(i,j) = uspdavg + hpbl_thetav2(i,j) = z_lev(kpbl_thetav(i,j)) IF ((frp_inst(i,j) .gt. frp_threshold) .AND. (frp_inst(i,j) .le. frp_threshold500) .AND. & - (z_at_w(i,kpbl_thetav(i,j),j) .gt. zpbl_threshold) .AND. (wind_eff_opt .eq. 1)) THEN + (z_lev(kpbl_thetav(i,j)) .gt. zpbl_threshold) .AND. (wind_eff_opt .eq. 1)) THEN kp1=1 IF (uspdavg .ge. uspd_threshold) THEN ! Too windy kp2=kpbl_thetav(i,j)/3 @@ -157,11 +173,18 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & END IF ! SRB: End modification - IF ( dbg_opt ) then - WRITE(*,*) 'module_plumerise: i,j ',i,j - WRITE(*,*) 'module_plumerise: k_min(i,j), k_max(i,j) ',kp1, kp2 ! SRB: replaced k_min, k_max with kp1, kp2 + IF ( dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst(i,j) .ge. frp_threshold) ) then + WRITE(1000+mpiid,*) 'mod_plumerise_after:xlat,xlong,curr_secs,k_min(i,j), k_max(i,j) ',xlat(i,j),xlong(i,j),int(curr_secs),kp1,kp2 + WRITE(1000+mpiid,*) 'mod_plumerise_after:xlat,xlong,curr_secs,ebu(kts),frp_inst',xlat(i,j),xlong(i,j),int(curr_secs),ebu(i,kts,j),frp_inst(i,j) + WRITE(1000+mpiid,*) 'mod_plumerise_after:xlat,xlong,curr_secs,u(10),v(10),w(10),qv(10)',xlat(i,j),xlong(i,j),int(curr_secs),u_in(10),v_in(10),w_in(kte),qv_in(10) + WRITE(1000+mpiid,*) 'mod_plumerise_after:xlat,xlong,curr_secs,uspdavg,kpbl_thetav',xlat(i,j),xlong(i,j),int(curr_secs),uspdavg,kpbl_thetav(i,j) + IF ( frp_inst(i,j) .ge. 3.e+9 ) then + WRITE(1000+mpiid,*) 'mod_plumerise_after:High FRP at : xlat,xlong,curr_secs,frp_inst',xlat(i,j),xlong(i,j),int(curr_secs),frp_inst(i,j) + END IF + icall = icall + 1 END IF -! endif check_frp +! endif check_frp +! icall = icall + 1 enddo enddo diff --git a/physics/smoke_dust/module_smoke_plumerise.F90 b/physics/smoke_dust/module_smoke_plumerise.F90 index 0fca91de4..aa45890f4 100755 --- a/physics/smoke_dust/module_smoke_plumerise.F90 +++ b/physics/smoke_dust/module_smoke_plumerise.F90 @@ -14,10 +14,11 @@ module module_smoke_plumerise use machine , only : kind_phys - use plume_data_mod, only : num_frp_plume, p_frp_hr, p_frp_std + !use plume_data_mod, only : num_frp_plume, p_frp_hr, p_frp_std !tropical_forest, boreal_forest, savannah, grassland, & ! wind_eff USE module_zero_plumegen_coms + USE rrfs_smoke_config, only : n_dbg_lines !real(kind=kind_phys),parameter :: rgas=r_d !real(kind=kind_phys),parameter :: cpor=cp/r_d @@ -28,12 +29,13 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & up,vp,wp,theta,pp,dn0,rv,zt_rams,zm_rams, & wind_eff_opt, & frp_inst,k1,k2, dbg_opt, g, cp, rgas, & - cpor, errmsg, errflg ) + cpor, errmsg, errflg, icall, mpiid, lat, long, curr_secs ) implicit none LOGICAL, INTENT (IN) :: dbg_opt - INTEGER, INTENT (IN) :: wind_eff_opt + INTEGER, INTENT (IN) :: wind_eff_opt, mpiid + real(kind_phys), INTENT(IN) :: lat,long, curr_secs ! SRB ! INTEGER, PARAMETER :: ihr_frp=1, istd_frp=2!, imean_fsize=3, istd_fsize=4 ! RAR: @@ -70,7 +72,7 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & ! real(kind=kind_phys), dimension(nveg_agreg) :: firesize,mean_fct INTEGER :: wind_eff - + INTEGER, INTENT(IN) :: icall type(plumegen_coms), pointer :: coms ! Set wind effect from namelist @@ -162,19 +164,11 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & endif burnt_area= max(1.0e4,burnt_area) - IF (dbg_opt) THEN - WRITE(*,*) 'plumerise: m1 ', m1 - WRITE(*,*) 'plumerise: imm, FRP,burnt_area ', imm, FRP,burnt_area - ! WRITE(*,*) 'convert_smold_to_flam ',convert_smold_to_flam - WRITE(*,*) 'plumerise: zcon ', coms%zcon - WRITE(*,*) 'plumerise: zzcon ', coms%zzcon - END IF - - IF (dbg_opt) then - WRITE(*,*) 'plumerise: imm ', imm - WRITE(*,*) 'plumerise: burnt_area ',burnt_area - END IF - + IF ( dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst .ge. frp_threshold) ) THEN + WRITE(1000+mpiid,*) 'inside plumerise: xlat,xlong,curr_secs, m1 ', lat,long, int(curr_secs), m1 + WRITE(1000+mpiid,*) 'inside plumerise: xlat,xlong,curr_secs,imm,FRP,burnt_area ', lat, long, int(curr_secs), imm, FRP,burnt_area + END IF + !- get fire properties (burned area, plume radius, heating rates ...) call get_fire_properties(coms,imm,iveg_ag,burnt_area,FRP,errmsg,errflg) if(errflg/=0) return @@ -182,8 +176,8 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & !------ generates the plume rise ------ call makeplume (coms,kmt,ztopmax(imm),ixx,imm) - IF (dbg_opt) then - WRITE(*,*) 'plumerise after makeplume: imm,kmt,ztopmax(imm) ',imm,kmt,ztopmax(imm) + IF ( dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst .ge. frp_threshold) ) then + WRITE(1000+mpiid,*) 'inside plumerise after makeplume:xlat,xlong,curr_secs,imm,kmt,ztopmax(imm) ', lat, long, int(curr_secs), imm,kmt, ztopmax(imm) END IF enddo lp_minmax @@ -203,12 +197,12 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & ! enddo !enddo - IF (dbg_opt) then - WRITE(*,*) 'plumerise after set_flam_vert: nkp,k1,k2, ', nkp,k1,k2 - WRITE(*,*) 'plumerise after set_flam_vert: dzi ', dzi + !IF (dbg_opt) then + ! WRITE(*,*) 'plumerise after set_flam_vert: nkp,k1,k2, ', nkp,k1,k2 + ! WRITE(*,*) 'plumerise after set_flam_vert: dzi ', dzi !WRITE(*,*) 'plumerise after set_flam_vert: eburn_in(2) ', eburn_in(2) !WRITE(*,*) 'plumerise after set_flam_vert: eburn_out(:,2) ',eburn_out(:,2) - END IF + !END IF ! enddo lp_veg ! sub-grid vegetation, currently it's aggregated diff --git a/physics/smoke_dust/module_wetdep_ls.F90 b/physics/smoke_dust/module_wetdep_ls.F90 index 8ba8f67d9..2ef07e38c 100755 --- a/physics/smoke_dust/module_wetdep_ls.F90 +++ b/physics/smoke_dust/module_wetdep_ls.F90 @@ -31,11 +31,18 @@ subroutine wetdep_ls(dt,var,rain,moist, real(kind_phys) :: dvar,factor,clsum integer :: nv,i,j,k,km,kb,kbeg !real(kind_phys), parameter :: alpha = .5 ! scavenging factor + integer, save :: print_alpha = 0 wetdpr_smoke =0. wetdpr_dust =0. wetdpr_coarsepm=0. + !if ( print_alpha == 0 ) then + ! write(*,*) 'wetdep_ls, alpha = ',alpha + ! print_alpha = print_alpha + 1 + !endif + + do nv=1,nchem do i=its,ite do j=jts,jte @@ -76,11 +83,11 @@ subroutine wetdep_ls(dt,var,rain,moist, dvar=alpha*factor/(1+factor)*var(i,k,j,nv) ! Accumulate diags if (nv .eq. p_smoke ) then - wetdpr_smoke(i,j) = wetdpr_smoke(i,j) + dvar * rho(i,k,j) * dt * 1.E-6 + wetdpr_smoke(i,j) = wetdpr_smoke(i,j) + dvar * rho(i,k,j) / dt elseif (nv .eq. p_dust_1 ) then - wetdpr_dust(i,j) = wetdpr_dust(i,j) + dvar * rho(i,k,j) * dt * 1.E-6 + wetdpr_dust(i,j) = wetdpr_dust(i,j) + dvar * rho(i,k,j) / dt elseif (nv .eq. p_coarse_pm ) then - wetdpr_coarsepm(i,j) = wetdpr_coarsepm(i,j) + dvar * rho(i,k,j) * dt * 1.E-6 + wetdpr_coarsepm(i,j) = wetdpr_coarsepm(i,j) + dvar * rho(i,k,j) / dt endif var(i,k,j,nv)=max(1.e-16,var(i,k,j,nv)-dvar) endif diff --git a/physics/smoke_dust/plume_data_mod.F90 b/physics/smoke_dust/plume_data_mod.F90 deleted file mode 100755 index 3f0bcdecd..000000000 --- a/physics/smoke_dust/plume_data_mod.F90 +++ /dev/null @@ -1,51 +0,0 @@ -!>\file plume_data_mod.F90 -!! This file contains data for the fire plume rise module. - -module plume_data_mod - - use machine , only : kind_phys - - implicit none - - ! -- FRP parameters - integer, dimension(0:20), parameter :: & - catb = (/ & - 0, & - 2, 1, 2, 1, & !floresta tropical 2 and 4 / extra trop fores 1,3,5 - 2, 3, 3, 3, 3, & !cerrado/woody savanna :6 a 9 - 4, 4, 4, 4, 4, 0, 4, 0, 0, 0, 0 & !pastagem/lavouras: 10 ... - /) - - real(kind=kind_phys), dimension(0:4), parameter :: & - flaming = (/ & - 0.00, & ! - 0.45, & ! % biomass burned at flaming phase : tropical forest igbp 2 and 4 - 0.45, & ! % biomass burned at flaming phase : extratropical forest igbp 1 , 3 and 5 - 0.75, & ! % biomass burned at flaming phase : cerrado/woody savanna igbp 6 to 9 - 0.00 & ! % biomass burned at flaming phase : pastagem/lavoura: igbp 10 a 17 - /) - - real(kind=kind_phys), dimension(0:20), parameter :: & - msize= (/ & - 0.00021, & !0near water,1Evergreen needleleaf,2EvergreenBroadleaf,!3Deciduous Needleleaf,4Deciduous Broadleaf - 0.00021, 0.00021, 0.00021, 0.00021, & !5Mixed forest,6Closed shrublands,7Open shrublands,8Woody savannas,9Savannas, - 0.00023, 0.00022, 0.00022, 0.00022, 0.00029, &! 10Grassland,11Permanent wetlands,12cropland,13'Urban and Built-Up' - 0.00029, 0.00021, 0.00026, 0.00021, 0.00026, &!14cropland/natural vegetation mosaic,15Snow and ice,16Barren or sparsely vegetated - 0.00021, 0.00021, 0.00021, 0.00021, 0.00021, 0.00021 & !17Water,18Wooded Tundra,19Mixed Tundra,20Bare Ground Tundra - /) - - ! -- FRP buffer indices - integer, parameter :: p_frp_hr = 1 - integer, parameter :: p_frp_std = 2 - integer, parameter :: num_frp_plume = 2 - - ! -- plumerise parameters - integer, parameter :: tropical_forest = 1 - integer, parameter :: boreal_forest = 2 - integer, parameter :: savannah = 3 - integer, parameter :: grassland = 4 - integer, parameter :: nveg_agreg = 4 - - public - -end module plume_data_mod diff --git a/physics/smoke_dust/rrfs_smoke_config.F90 b/physics/smoke_dust/rrfs_smoke_config.F90 index d7478986b..dae4338bb 100755 --- a/physics/smoke_dust/rrfs_smoke_config.F90 +++ b/physics/smoke_dust/rrfs_smoke_config.F90 @@ -24,6 +24,7 @@ module rrfs_smoke_config integer :: addsmoke_flag = 1 integer :: smoke_forecast = 1 integer :: plumerisefire_frq=60 + integer :: n_dbg_lines = 3 integer :: wetdep_ls_opt = 1 integer :: drydep_opt = 1 integer :: pm_settling = 1 @@ -39,7 +40,7 @@ module rrfs_smoke_config ! -- integer, parameter :: CHEM_OPT_GOCART= 1 - integer, parameter :: num_moist=2, num_chem=20, num_emis_seas=5, num_emis_dust=5 + integer, parameter :: num_moist=3, num_chem=20, num_emis_seas=5, num_emis_dust=5 ! -- hydrometeors integer, parameter :: p_qv=1 diff --git a/physics/smoke_dust/rrfs_smoke_postpbl.meta b/physics/smoke_dust/rrfs_smoke_postpbl.meta index 50fbb4e03..8d7481ec4 100755 --- a/physics/smoke_dust/rrfs_smoke_postpbl.meta +++ b/physics/smoke_dust/rrfs_smoke_postpbl.meta @@ -1,8 +1,7 @@ [ccpp-table-properties] name = rrfs_smoke_postpbl type = scheme - dependencies = ../hooks/machine.F,dep_dry_simple_mod.F90,module_wetdep_ls.F90,dust_data_mod.F90,dust_fengsha_mod.F90,module_add_emiss_burn.F90,module_plumerise.F90,module_smoke_plumerise.F90,module_zero_plumegen_coms.F90,plume_data_mod.F90,rrfs_smoke_config.F90,seas_data_mod.F90,seas_mod.F90,seas_ngac_mod.F90,dep_dry_mod_emerson.F90,dep_data_mod.F90 - + dependencies = dep_dry_simple_mod.F90,module_wetdep_ls.F90,dust_data_mod.F90,dust_fengsha_mod.F90,module_add_emiss_burn.F90,module_plumerise.F90,module_smoke_plumerise.F90,module_zero_plumegen_coms.F90,rrfs_smoke_config.F90,seas_data_mod.F90,seas_mod.F90,seas_ngac_mod.F90,dep_dry_mod_emerson.F90,dep_data_mod.F90 ######################################################################## [ccpp-arg-table] name = rrfs_smoke_postpbl_run diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.F90 b/physics/smoke_dust/rrfs_smoke_wrapper.F90 index eb7f83af6..4daad7168 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.F90 +++ b/physics/smoke_dust/rrfs_smoke_wrapper.F90 @@ -12,7 +12,7 @@ module rrfs_smoke_wrapper ebb_dcycle, extended_sd_diags,add_fire_heat_flux, & num_moist, num_chem, num_emis_seas, num_emis_dust, & p_qv, p_atm_shum, p_atm_cldq, & - p_smoke, p_dust_1, p_coarse_pm, epsilc + p_smoke, p_dust_1, p_coarse_pm, epsilc, n_dbg_lines use dust_data_mod, only : dust_alpha, dust_gamma, dust_moist_opt, & dust_moist_correction, dust_drylimit_factor use seas_mod, only : gocart_seasalt_driver @@ -52,7 +52,7 @@ subroutine rrfs_smoke_wrapper_init( seas_opt_in, dust_moist_correction_in, dust_drylimit_factor_in, & ! Dust namelist aero_ind_fdb_in, & ! Feedback namelist extended_sd_diags_in,dbg_opt_in, & ! Other namelist - errmsg, errflg ) + errmsg, errflg, n_dbg_lines_in ) !>-- Namelist @@ -62,7 +62,7 @@ subroutine rrfs_smoke_wrapper_init( seas_opt_in, integer, intent(in) :: dust_opt_in,dust_moist_opt_in, wetdep_ls_opt_in, pm_settling_in, seas_opt_in integer, intent(in) :: drydep_opt_in logical, intent(in) :: aero_ind_fdb_in,dbg_opt_in, extended_sd_diags_in, add_fire_heat_flux_in - integer, intent(in) :: smoke_forecast_in, plume_wind_eff_in, plumerisefire_frq_in + integer, intent(in) :: smoke_forecast_in, plume_wind_eff_in, plumerisefire_frq_in, n_dbg_lines_in integer, intent(in) :: addsmoke_flag_in, ebb_dcycle_in logical, intent(in) :: do_plumerise_in, rrfs_sd character(len=*),intent(out):: errmsg @@ -100,6 +100,7 @@ subroutine rrfs_smoke_wrapper_init( seas_opt_in, !>-Other extended_sd_diags = extended_sd_diags_in dbg_opt = dbg_opt_in + n_dbg_lines = n_dbg_lines_in end subroutine rrfs_smoke_wrapper_init @@ -111,17 +112,19 @@ end subroutine rrfs_smoke_wrapper_init subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, & u10m, v10m, ustar, rlat, rlon, tskin, pb2d, t2m, dpt2m, & pr3d, ph3d,phl3d, prl3d, tk3d, us3d, vs3d, spechum, w, & - nsoil, smc, vegtype_dom, vegtype_frac, soiltyp, nlcat, & + nsoil, smc, tslb, vegtype_dom, vegtype_frac, soiltyp, nlcat, & dswsfc, zorl, snow, julian,recmol, & idat, rain_cpl, rainc_cpl, hf2d, g, pi, con_cp, con_rd, con_fv, & dust12m_in, emi_ant_in, smoke_RRFS, smoke2d_RRFS, & ntrac, qgrs, gq0, chem3d, tile_num, & ntsmoke, ntdust, ntcoarsepm, imp_physics, imp_physics_thompson, & nwfa, nifa, emanoc, emdust, emseas, drydep_flux_out, wetdpr, & - ebb_smoke_in, frp_output, coef_bb, ebu_smoke,fhist,min_fplume, & + ebb_smoke_in, frp_output, coef_bb, fire_type_out, & + ebu_smoke,fhist,min_fplume, & max_fplume, hwp, hwp_ave, wetness, ndvel, ddvel_inout,fire_in, & + peak_hr_out,lu_nofire_out,lu_qfire_out, & fire_heat_flux_out, frac_grid_burned_out, kpbl,oro, & - errmsg,errflg ) + uspdavg, hpbl_thetav, mpicomm, mpirank, mpiroot, errmsg,errflg ) implicit none @@ -135,7 +138,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, integer, parameter :: its=1,jts=1,jte=1, kts=1 integer, dimension(:), intent(in) :: land, vegtype_dom, soiltyp - real(kind_phys), dimension(:,:), intent(in) :: smc + real(kind_phys), dimension(:,:), intent(in) :: smc, tslb real(kind_phys), dimension(:,:,:), intent(in) :: dust12m_in real(kind_phys), dimension(:,:,:), intent(in) :: smoke_RRFS real(kind_phys), dimension(:,:), intent(in) :: smoke2d_RRFS @@ -153,14 +156,16 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, real(kind_phys), dimension(:,:), intent(inout) :: ebu_smoke real(kind_phys), dimension(:,:), intent(inout) :: fire_in real(kind_phys), dimension(:), intent(out) :: fire_heat_flux_out, frac_grid_burned_out - real(kind_phys), dimension(:), intent(inout) :: max_fplume, min_fplume - real(kind_phys), dimension(:), intent( out) :: hwp + real(kind_phys), dimension(:), intent(inout) :: max_fplume, min_fplume, uspdavg, hpbl_thetav + real(kind_phys), dimension(:), intent(inout) :: hwp, peak_hr_out real(kind_phys), dimension(:), intent(inout) :: hwp_ave real(kind_phys), dimension(:,:), intent(inout) :: nwfa, nifa real(kind_phys), dimension(:,:), intent(inout) :: ddvel_inout real(kind_phys), dimension(:,:), intent(inout) :: drydep_flux_out real(kind_phys), dimension(:,:), intent(inout) :: wetdpr real(kind_phys), dimension(:), intent(in) :: wetness + real(kind_phys), dimension(:), intent(out) :: lu_nofire_out,lu_qfire_out + integer, dimension(:), intent(out) :: fire_type_out integer, intent(in) :: imp_physics, imp_physics_thompson integer, dimension(:), intent(in) :: kpbl real(kind_phys), dimension(:), intent(in) :: oro @@ -187,16 +192,17 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, real(kind_phys), dimension(ims:im, jms:jme) :: ssm, rdrag, uthr, snowh ! fengsha dust real(kind_phys), dimension(ims:im, jms:jme) :: rmol, swdown, znt, clayf, sandf real(kind_phys), dimension(ims:im, nlcat, jms:jme) :: vegfrac - real(kind_phys), dimension(ims:im, nsoil, jms:jme) :: smois + real(kind_phys), dimension(ims:im, nsoil, jms:jme) :: smois, stemp real(kind_phys), dimension(ims:im, 1:1, jms:jme, 1:num_emis_dust) :: emis_dust integer, dimension(ims:im, jms:jme) :: isltyp, ivgtyp !>- plume variables ! -- buffers real(kind_phys), dimension(ims:im, jms:jme ) :: coef_bb_dc, flam_frac, frp_in, & fire_hist, peak_hr, lu_nofire, lu_qfire, ebu_in, & - fire_end_hr, hwp_day_avg, kpbl_thetav + fire_end_hr, hwp_day_avg, kpbl_thetav,& + uspdavg2, hpbl_thetav2 integer, dimension(ims:im, jms:jme ) :: min_fplume2, max_fplume2, fire_type - logical :: call_plume, reset_hwp_ave + logical :: call_plume, reset_hwp_ave, avg_hwp_ave !>- optical variables real(kind_phys), dimension(ims:im, jms:jme, ndvel) :: ddvel, settling_flux, drydep_flux_local real(kind_phys), dimension(ims:im, kms:kme, jms:jme, ndvel) :: vgrav @@ -218,6 +224,13 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, real(kind_phys) :: factor, factor2, factor3 integer :: nbegin, nv integer :: i, j, k, kp, n +! MPI variables + integer :: mpiid + integer, intent(in) :: mpicomm + integer, intent(in) :: mpirank + integer, intent(in) :: mpiroot + + mpiid = mpirank errmsg = '' errflg = 0 @@ -232,6 +245,8 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, min_fplume2 = 0 max_fplume2 = 0 + uspdavg2 = 0. + hpbl_thetav2 = 0. emis_seas = 0. emis_dust = 0. peak_hr = 0. @@ -260,12 +275,17 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, ! -- compute incremental convective and large-scale rainfall do i=its,ite rnav(i,1)=max((rain_cpl(i)-rainc_cpl(i))*1000., 0.) ! meter to mm +! coef_bb initializes as clear_val (from GFS_typedefs.F90) +! at ktau = 1, coef_bb_dc is set = 1.0 coef_bb_dc(i,1) = coef_bb(i) +! fhist initializes as 1. (from GFS_typedefs.F90) fire_hist (i,1) = fhist (i) + peak_hr (i,1) = peak_hr_out(i) enddo ! Is this a reset timestep (00:00 + dt)? reset_hwp_ave = mod(int(curr_secs-dt),3600) == 0 + avg_hwp_ave = mod(int(curr_secs),3600) == 0 ! plumerise frequency in minutes set up by the namelist input call_plume = (do_plumerise .and. (plumerisefire_frq > 0)) @@ -276,7 +296,8 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, ktau,current_month, current_hour, gmt, con_rd, con_fv, con_cp, & u10m,v10m,ustar,land,garea,rlat,rlon,tskin, & pr3d,ph3d,phl3d,tk3d,prl3d,us3d,vs3d,spechum,w, & - nsoil,smc,vegtype_dom,soiltyp,nlcat,vegtype_frac,dswsfc,zorl, & + nsoil,smc,tslb,vegtype_dom,soiltyp, & + nlcat,vegtype_frac,dswsfc,zorl, & snow,dust12m_in,emi_ant_in,smoke_RRFS,smoke2d_RRFS,coef_bb_dc, & hf2d, pb2d, g, pi, hour_int, peak_hr, & u10,v10,ust,tsk,xland,xlat,xlong,dxy, & @@ -287,20 +308,14 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, num_chem,num_moist, & ntsmoke, ntdust,ntcoarsepm, & moist,chem,ebu_in,kpbl_thetav,ebb_smoke_in, & - fhist,frp_in, hwp_day_avg, fire_end_hr, & - emis_anoc,smois,ivgtyp,isltyp,vegfrac,rmol,swdown,znt,hfx,pbl, & - snowh,clayf,rdrag,sandf,ssm,uthr,oro, hwp_local, & + fire_hist,frp_in, hwp_day_avg, fire_end_hr, & + emis_anoc,smois,stemp,ivgtyp,isltyp,vegfrac,rmol,swdown,znt, & + hfx,pbl,snowh,clayf,rdrag,sandf,ssm,uthr,oro, hwp_local, & t2m,dpt2m,wetness,kpbl, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) - do j=jts,jte - do i=its,ite - peak_hr(i,j)= fire_in(i,1) - enddo - enddo - IF (ktau==1) THEN ebu = 0. do j=jts,jte @@ -311,6 +326,13 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, enddo enddo enddo + ELSE + do k=kts,kte + do i=its,ite + ! ebu is divided by coef_bb_dc since it is applied in the output + ebu(i,k,1)=ebu_smoke(i,k) / coef_bb_dc(i,1) + enddo + enddo ENDIF !RAR: change this to the fractional LU type; fire_type: 0- no fires, 1- Ag @@ -320,6 +342,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, do i=its,ite if (ebu_in(i,j)<0.01) then fire_type(i,j) = 0 + lu_nofire(i,j) = 1.0 else ! Permanent wetlands, snow/ice, water, barren tundra lu_nofire(i,j) = vegfrac(i,11,j) + vegfrac(i,15,j) + vegfrac(i,17,j) + vegfrac(i,20,j) @@ -350,9 +373,10 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, endif !-- compute dust (opt=5) - if (dust_opt==5) then - call gocart_dust_fengsha_driver(dt,chem,rho_phy,smois,p8w,ssm, & - isltyp,snowh,xland,dxy,g,emis_dust,ust,znt, & + if (dust_opt==1) then + call gocart_dust_fengsha_driver(dt,chem,rho_phy, & + smois,stemp,p8w,ssm, & + isltyp,snowh,xland,dxy,g,emis_dust,ust,znt, & clayf,sandf,rdrag,uthr, & num_emis_dust,num_chem,nsoil, & ids,ide, jds,jde, kds,kde, & @@ -367,6 +391,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, ! Every hour (per namelist) the ebu_driver is called to calculate ebu, but ! the plumerise is controlled by the namelist option of plumerise_flag if (add_fire_heat_flux) then + WRITE(1000+mpiid,*) 'Entered add_fire_heat_flux at timestep:',ktau do i = its,ite if ( coef_bb_dc(i,1)*frp_in(i,1) .ge. 1.E7 ) then fire_heat_flux_out(i) = min(max(0.,0.88*coef_bb_dc(i,1)*frp_in(i,1) / & @@ -396,7 +421,8 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, kpbl_thetav, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte, errmsg, errflg ) + its,ite, jts,jte, kts,kte, errmsg, errflg, curr_secs, & + xlat, xlong, uspdavg2, hpbl_thetav2, mpiid ) if(errflg/=0) return end if @@ -409,7 +435,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, swdown,ebb_dcycle,ebu_in,ebu,fire_type, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte ) + its,ite, jts,jte, kts,kte , mpiid ) endif !>-- compute coarsepm setting if using simple dry dep option and @@ -431,7 +457,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, pm_settling,drydep_flux_local,settling_flux,dbg_opt, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte ) + its,ite, jts,jte, kts,kte, curr_secs, mpiid, xlat, xlong ) do nv=1,ndvel do i=its,ite ddvel_inout(i,nv)=ddvel(i,1,nv) @@ -470,10 +496,12 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, endif endif -! Smoke emisisons diagnostic +! Smoke emisisons diagnostic, RAR: let's multiply by coef_bb_dc before output +! Since ebu_smoke includes coef_bb_dc, we need to divide by coef_bb_dc when it +! comes back into the wrapper. do k=kts,kte do i=its,ite - ebu_smoke(i,k)=ebu(i,k,1) + ebu_smoke(i,k)=ebu(i,k,1) * coef_bb_dc(i,1) enddo enddo @@ -485,15 +513,21 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, do i=its,ite hwp(i)=hwp_local(i,1) hwp_ave(i) = hwp_ave(i) + hwp(i)*dt + if ( ktau == 1) then + hwp_ave(i) = hwp_ave(i) / dt + elseif ( avg_hwp_ave ) then + hwp_ave(i) = hwp_ave(i) / 3600._kind_phys + endif enddo - + + !---- diagnostic output of dry deposition & gravitational settling fluxes if ( drydep_opt == 1 .and. (extended_sd_diags .or. dbg_opt) ) then do nv = 1, ndvel do i=its,ite drydep_flux_out(i,nv) = drydep_flux_out(i,nv) + & - drydep_flux_local(i,1,nv) + & - settling_flux(i,1,nv) + drydep_flux_local(i,1,nv) !+ & + !settling_flux(i,1,nv) enddo enddo endif @@ -520,6 +554,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, !------------------------------------- !-- to output for diagnostics do i = 1, im +! RAR: let's remove the seas and ant. OC emseas (i) = emis_seas(i,1,1,1)*1.e+9 ! size bin 1 sea salt emission: ug/m2/s emanoc (i) = emis_anoc (i) ! anthropogenic organic carbon: ug/m2/s emdust (i) = emis_dust(i,1,1,1) + emis_dust(i,1,1,2) + & @@ -529,10 +564,13 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, fhist (i) = fire_hist (i,1) min_fplume (i) = real(min_fplume2(i,1)) max_fplume (i) = real(max_fplume2(i,1)) + fire_type_out(i)=fire_type(i,1) + lu_nofire_out(i)=lu_nofire(i,1) + lu_qfire_out (i)=lu_qfire(i,1) enddo do i = 1, im - fire_in(i,1) = peak_hr(i,1) + peak_hr_out(i) = peak_hr(i,1) enddo !-- to provide real aerosol emission for Thompson MP @@ -568,7 +606,7 @@ subroutine rrfs_smoke_prep( & ktau,current_month,current_hour,gmt,con_rd,con_fv,con_cp, & u10m,v10m,ustar,land,garea,rlat,rlon,ts2d, & pr3d,ph3d,phl3d,tk3d,prl3d,us3d,vs3d,spechum,w, & - nsoil,smc,vegtype_dom,soiltyp,nlcat,vegtype_frac,dswsfc,zorl, & + nsoil,smc,tslb,vegtype_dom,soiltyp,nlcat,vegtype_frac,dswsfc,zorl, & snow_cpl,dust12m_in,emi_ant_in,smoke_RRFS,smoke2d_RRFS,coef_bb_dc, & hf2d, pb2d, g, pi, hour_int, peak_hr, & u10,v10,ust,tsk,xland,xlat,xlong,dxy, & @@ -579,9 +617,9 @@ subroutine rrfs_smoke_prep( & num_chem, num_moist, & ntsmoke, ntdust, ntcoarsepm, & moist,chem,ebu_in,kpbl_thetav,ebb_smoke_in, & - fhist,frp_in, hwp_day_avg, fire_end_hr, & - emis_anoc,smois,ivgtyp,isltyp,vegfrac,rmol,swdown,znt,hfx,pbl, & - snowh,clayf,rdrag,sandf,ssm,uthr,oro,hwp_local, & + fire_hist,frp_in, hwp_day_avg, fire_end_hr, & + emis_anoc,smois,stemp,ivgtyp,isltyp,vegfrac,rmol,swdown, & + znt,hfx,pbl,snowh,clayf,rdrag,sandf,ssm,uthr,oro,hwp_local, & t2m,dpt2m,wetness,kpbl, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -599,7 +637,7 @@ subroutine rrfs_smoke_prep( & u10m, v10m, ustar, garea, rlat, rlon, ts2d, dswsfc, & zorl, snow_cpl, pb2d, hf2d, oro, t2m, dpt2m, wetness, recmol real(kind=kind_phys), dimension(ims:ime, nlcat), intent(in) :: vegtype_frac - real(kind=kind_phys), dimension(ims:ime, nsoil), intent(in) :: smc + real(kind=kind_phys), dimension(ims:ime, nsoil), intent(in) :: smc,tslb real(kind=kind_phys), dimension(ims:ime, 12, 5), intent(in) :: dust12m_in real(kind=kind_phys), dimension(ims:ime, 24, 2), intent(in) :: smoke_RRFS ! This is a place holder for ebb_dcycle == 2, currently set to hold a single @@ -633,8 +671,8 @@ subroutine rrfs_smoke_prep( & real(kind_phys), dimension(ims:ime, kms:kme, jms:jme, num_chem), intent(out) :: chem real(kind_phys), dimension(ims:ime, kms:kme, jms:jme), intent(out) :: z_at_w - real(kind_phys), dimension(ims:ime, nsoil, jms:jme), intent(out) :: smois - real(kind_phys), dimension(ims:ime,jms:jme), intent(inout) :: frp_in, fire_end_hr, fhist, coef_bb_dc + real(kind_phys), dimension(ims:ime, nsoil, jms:jme), intent(out) :: smois,stemp + real(kind_phys), dimension(ims:ime,jms:jme), intent(inout) :: frp_in, fire_end_hr, fire_hist, coef_bb_dc real(kind_phys), dimension(ims:ime,jms:jme), intent(inout) :: hwp_day_avg, peak_hr real(kind_phys), dimension(ims:ime), intent(inout) :: emis_anoc,ebb_smoke_in real(kind_phys), parameter :: conv_frp = 1.e+06_kind_phys ! FRP conversion factor, MW to W @@ -730,6 +768,7 @@ subroutine rrfs_smoke_prep( & do j=jts,jte do i=its,ite smois(i,k,j)=smc(i,k) + stemp(i,k,j)=tslb(i,k) enddo enddo enddo @@ -776,13 +815,14 @@ subroutine rrfs_smoke_prep( & rri(i,k,j)=1./rho_phy(i,k,j) vvel(i,k,j)=-w(i,kkp)*rri(i,k,j)/g moist(i,k,j,:)=0. - moist(i,k,j,1)=gq0(i,kkp,p_atm_shum) - if (t_phy(i,k,j) > 265.) then - moist(i,k,j,2)=gq0(i,kkp,p_atm_cldq) + moist(i,k,j,1)=gq0(i,kkp,1) + !if (t_phy(i,k,j) > 265.) then + moist(i,k,j,2)=gq0(i,kkp,2) if (moist(i,k,j,2) < 1.e-8) moist(i,k,j,2)=0. - else - moist(i,k,j,2)=0. - endif + !else + ! moist(i,k,j,2)=0. + if(moist(i,k,j,3) < 1.e-8)moist(i,k,j,3)=0. + !endif !-- zmid(i,k,j)=phl3d(i,kkp)/g enddo @@ -862,11 +902,11 @@ subroutine rrfs_smoke_prep( & if (hour_int .le. 24) then do j=jts,jte do i=its,ite - ebu_in (i,j) = smoke_RRFS(i,hour_int+1,1) ! smoke frp_in (i,j) = smoke_RRFS(i,hour_int+1,2)*conv_frp ! frp - fire_end_hr(i,j) = 0.0 - hwp_day_avg(i,j) = 0.0 + ! These 2 arrays aren't needed for this option + ! fire_end_hr(i,j) = 0.0 + ! hwp_day_avg(i,j) = 0.0 ebb_smoke_in (i) = ebu_in(i,j) enddo enddo @@ -890,7 +930,8 @@ subroutine rrfs_smoke_prep( & if (ktau==1) then do j=jts,jte do i=its,ite - fhist (i,j) = 1. + ! GFS_typedefs.F90 initializes this = 1, but should be OK to duplicate, RAR?? + fire_hist (i,j) = 1. coef_bb_dc (i,j) = 1. if (xlong(i,j)<230.) then peak_hr(i,j)= 0.0* 3600. ! peak at 24 UTC, fires in Alaska @@ -909,7 +950,7 @@ subroutine rrfs_smoke_prep( & enddo endif - ! We will add a namelist variable, real :: flam_frac_global + ! We will add a namelist variable, real :: flam_frac_global, RAR?? do k=kms,kte do i=ims,ime chem(i,k,jts,p_smoke )=max(epsilc,gq0(i,k,ntsmoke )) diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.meta b/physics/smoke_dust/rrfs_smoke_wrapper.meta index af61ac05e..fc3aa9fe6 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.meta +++ b/physics/smoke_dust/rrfs_smoke_wrapper.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = rrfs_smoke_wrapper type = scheme - dependencies = ../hooks/machine.F,dep_dry_simple_mod.F90,module_wetdep_ls.F90,dust_data_mod.F90,dust_fengsha_mod.F90,module_add_emiss_burn.F90,module_plumerise.F90,module_smoke_plumerise.F90,module_zero_plumegen_coms.F90,plume_data_mod.F90,rrfs_smoke_config.F90,seas_data_mod.F90,seas_mod.F90,seas_ngac_mod.F90,coarsepm_settling_mod.F90, dep_dry_mod_emerson.F90,dep_data_mod.F90 + dependencies = dep_dry_simple_mod.F90,module_wetdep_ls.F90,dust_data_mod.F90,dust_fengsha_mod.F90,module_add_emiss_burn.F90,module_plumerise.F90,module_smoke_plumerise.F90,module_zero_plumegen_coms.F90,rrfs_smoke_config.F90,seas_data_mod.F90,seas_mod.F90,seas_ngac_mod.F90,coarsepm_settling_mod.F90, dep_dry_mod_emerson.F90,dep_data_mod.F90 ######################################################################## [ccpp-arg-table] @@ -71,6 +71,13 @@ dimensions = () type = integer intent = in +[n_dbg_lines_in] + standard_name = smoke_debug_lines + long_name = rrfs smoke add smoke option + units = index + dimensions = () + type = integer + intent = in [plume_wind_eff_in] standard_name = option_for_wind_effects_on_smoke_plumerise long_name = wind effect plumerise option @@ -406,6 +413,14 @@ type = real kind = kind_phys intent = inout +[tslb] + standard_name = soil_temperature_for_land_surface_model + long_name = soil temperature for land surface model + units = K + dimensions = (horizontal_loop_extent,vertical_dimension_of_soil_internal_to_land_surface_scheme) + type = real + kind = kind_phys + intent = in [vegtype_dom] standard_name = vegetation_type_classification long_name = vegetation type at each grid cell @@ -715,6 +730,13 @@ type = real kind = kind_phys intent = inout +[fire_type_out] + standard_name = fire_type_out + long_name = type of fire + units = none + dimensions = (horizontal_loop_extent) + type = integer + intent = out [ebu_smoke] standard_name = ebu_smoke long_name = buffer of vertical fire emission @@ -747,6 +769,43 @@ type = real kind = kind_phys intent = inout +[mpicomm] + standard_name = mpi_communicator + long_name = MPI communicator + units = index + dimensions = () + type = integer + intent = in +[mpirank] + standard_name = mpi_rank + long_name = current MPI-rank + units = index + dimensions = () + type = integer + intent = in +[mpiroot] + standard_name = mpi_root + long_name = master MPI-rank + units = index + dimensions = () + type = integer + intent = in +[uspdavg] + standard_name = bl_averaged_wind_speed + long_name = average wind speed within the boundary layer + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[hpbl_thetav] + standard_name = pbl_height_thetav + long_name = pbl height based on modified parcel method + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout [drydep_flux_out] standard_name = dry_deposition_flux long_name = rrfs dry deposition flux @@ -810,6 +869,30 @@ type = real kind = kind_phys intent = inout +[peak_hr_out] + standard_name = peak_hr_fire + long_name = hour of peak fire emissions + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[lu_nofire_out] + standard_name = lu_nofire_out + long_name = land use of no fire pixels for type + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[lu_qfire_out] + standard_name = lu_qfire_out + long_name = land use of fire pixels for type + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [fire_heat_flux_out] standard_name = surface_fire_heat_flux long_name = heat flux of fire at the surface From c0544c218776d4a94169d45cb7dae102800594a1 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Tue, 23 Jan 2024 16:13:19 +0000 Subject: [PATCH 13/14] "update to address code reviewer's comments" --- physics/smoke_dust/dep_dry_mod_emerson.F90 | 14 ---------- physics/smoke_dust/module_add_emiss_burn.F90 | 27 +++++++++---------- physics/smoke_dust/rrfs_smoke_config.F90 | 2 +- physics/smoke_dust/rrfs_smoke_wrapper.F90 | 14 +++------- physics/smoke_dust/rrfs_smoke_wrapper.meta | 28 +++++++------------- 5 files changed, 27 insertions(+), 58 deletions(-) diff --git a/physics/smoke_dust/dep_dry_mod_emerson.F90 b/physics/smoke_dust/dep_dry_mod_emerson.F90 index 771801c44..e69d6bc3f 100755 --- a/physics/smoke_dust/dep_dry_mod_emerson.F90 +++ b/physics/smoke_dust/dep_dry_mod_emerson.F90 @@ -179,10 +179,6 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & eps0 = eps0_grs end if ! Set if snow greater than 1 cm -! if ( snowh(i,j) .gt. 0.01 ) then ! snow -! A = A_wat -! eps0 = eps0_wat -! endif ! Interception Ein = Cin * ( dp / A )**vv ! Surface resistance @@ -234,25 +230,15 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & IF (ndt_settl(nv) > 12) ndt_settl(nv) = 12 dt_settl(nv) = REAL(ntdt,kind=kind_phys) /REAL(ndt_settl(nv),kind=kind_phys) enddo - !do nv = 1, ndvel - ! chem_before(nv) = 0._kind_phys - ! do k = kts, kte - ! chem_before(nv) = chem_before(nv) + (cblk_col(k,nv) * rho_phy(i,k,j) * delz(i,k,j) ) ! ug/m2 - ! enddo - !enddo ! Perform gravitational settling if desired if ( settling_flag == 1 ) then call particle_settling(cblk_col,rho_col,delz_col,vg_col,dt_settl,ndt_settl,ndvel,kts,kte) endif ! Put cblk back into chem array do nv= 1, ndvel - !chem_after(nv) = 0._kind_phys - !settling_flux(i,j,nv) = 0._kind_phys do k = kts, kte chem(i,k,j,chem_pointers(nv)) = cblk_col(k,nv) - !chem_after(nv) = chem_after(nv) + (cblk_col(k,nv) * rho_phy(i,k,j) * delz(i,k,j) ) ! ug/m2 enddo ! k - !settling_flux(i,j,nv) = settling_flux(i,j,nv) + (chem_before(nv) - chem_after(nv)) ! ug/m2 enddo ! nv end do ! j end do ! i diff --git a/physics/smoke_dust/module_add_emiss_burn.F90 b/physics/smoke_dust/module_add_emiss_burn.F90 index 0a22fcfd7..f1bbaeee9 100755 --- a/physics/smoke_dust/module_add_emiss_burn.F90 +++ b/physics/smoke_dust/module_add_emiss_burn.F90 @@ -47,7 +47,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & INTEGER, PARAMETER :: kfire_max=51 ! max vertical level for BB plume rise - real(kind_phys) :: timeq, fire_age, age_hr, dt1,dt2,dtm ! For BB emis. diurnal cycle calculation + real(kind_phys) :: timeq, fire_age, age_hr, dt1,dt2,dtm, coef_con ! For BB emis. diurnal cycle calculation ! For Gaussian diurnal cycle real(kind_phys), PARAMETER :: sc_factor=1. ! to scale up the wildfire emissions, TBD later @@ -89,6 +89,8 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & if (ebb_dcycle==2) then ! Constants for the fire diurnal cycle calculation + coef_con = 1._kind_phys/((2._kind_phys*pi)**0.5_kind_phys * sigma_fire_dur(1) *fire_age) * & + exp(- ( log(fire_age) - avg_fire_dur(1))**2 /(2._kind_phys*sigma_fire_dur(1)**2 )) do j=jts,jte do i=its,ite fire_age= time_int + (fire_end_hr(i,j)-1._kind_phys)*3600._kind_phys !One hour delay is due to the latency of the RAVE files @@ -97,13 +99,12 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & SELECT CASE ( fire_type(i,j) ) !Ag, urban fires, bare land etc. CASE (1) ! these fires will have exponentially decreasing diurnal cycle, - coef_bb_dc(i,j) = 1._kind_phys/((2._kind_phys*pi)**0.5_kind_phys * sigma_fire_dur(1) *fire_age) * & - exp(- ( log(fire_age) - avg_fire_dur(1))**2 /(2._kind_phys*sigma_fire_dur(1)**2 )) + coef_bb_dc(i,j) = coef_con - ! IF ( dbg_opt .AND. time_int<5000.) then - ! WRITE(6,*) 'i,j,peak_hr(i,j) ',i,j,peak_hr(i,j) - ! WRITE(6,*) 'coef_bb_dc(i,j) ',coef_bb_dc(i,j) - ! END IF + IF ( dbg_opt .AND. time_int<5000.) then + WRITE(6,*) 'i,j,peak_hr(i,j) ',i,j,peak_hr(i,j) + WRITE(6,*) 'coef_bb_dc(i,j) ',coef_bb_dc(i,j) + END IF CASE (3) age_hr= fire_age/3600._kind_phys @@ -123,9 +124,6 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & dc_hwp= MAX(0._kind_phys,dc_hwp) dc_hwp= MIN(25._kind_phys,dc_hwp) - !coef_bb_dc(i,j)= sc_factor* fire_hist(i,j)* rate_ebb2(i,j)* (1. + log( - !hwp_(i,j)/ hwp_day_avg(i,j))) - ! RAR: Gaussian profile for wildfires dt1= abs(timeq - peak_hr(i,j)) dt2= timeq_max - peak_hr(i,j) + timeq ! peak hour is always <86400. @@ -134,13 +132,12 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & dc_gp = MAX(0._kind_phys,dc_gp) dc_fn = MIN(dc_hwp/dc_gp,3._kind_phys) - !coef_bb_dc(i,j) = fire_hist(i,j)* dc_fn coef_bb_dc(i,j) = fire_hist(i,j)* dc_hwp - ! IF ( dbg_opt .AND. time_int<5000.) then - ! WRITE(6,*) 'i,j,fire_hist(i,j),peak_hr(i,j) ', i,j,fire_hist(i,j),peak_hr(i,j) - ! WRITE(6,*) 'dc_gp,dc_hwp,dc_fn,coef_bb_dc(i,j) ',dc_gp,dc_hwp,dc_fn,coef_bb_dc(i,j) - ! END IF + IF ( dbg_opt .AND. time_int<5000.) then + WRITE(6,*) 'i,j,fire_hist(i,j),peak_hr(i,j) ', i,j,fire_hist(i,j),peak_hr(i,j) + WRITE(6,*) 'dc_gp,dc_hwp,dc_fn,coef_bb_dc(i,j) ',dc_gp,dc_hwp,dc_fn,coef_bb_dc(i,j) + END IF CASE DEFAULT END SELECT diff --git a/physics/smoke_dust/rrfs_smoke_config.F90 b/physics/smoke_dust/rrfs_smoke_config.F90 index dae4338bb..c20d6e2db 100755 --- a/physics/smoke_dust/rrfs_smoke_config.F90 +++ b/physics/smoke_dust/rrfs_smoke_config.F90 @@ -40,7 +40,7 @@ module rrfs_smoke_config ! -- integer, parameter :: CHEM_OPT_GOCART= 1 - integer, parameter :: num_moist=3, num_chem=20, num_emis_seas=5, num_emis_dust=5 + integer, parameter :: num_moist=2, num_chem=20, num_emis_seas=5, num_emis_dust=5 ! -- hydrometeors integer, parameter :: p_qv=1 diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.F90 b/physics/smoke_dust/rrfs_smoke_wrapper.F90 index 4daad7168..3842cba54 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.F90 +++ b/physics/smoke_dust/rrfs_smoke_wrapper.F90 @@ -121,7 +121,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, nwfa, nifa, emanoc, emdust, emseas, drydep_flux_out, wetdpr, & ebb_smoke_in, frp_output, coef_bb, fire_type_out, & ebu_smoke,fhist,min_fplume, & - max_fplume, hwp, hwp_ave, wetness, ndvel, ddvel_inout,fire_in, & + max_fplume, hwp, hwp_ave, wetness, ndvel, ddvel_inout, & peak_hr_out,lu_nofire_out,lu_qfire_out, & fire_heat_flux_out, frac_grid_burned_out, kpbl,oro, & uspdavg, hpbl_thetav, mpicomm, mpirank, mpiroot, errmsg,errflg ) @@ -154,8 +154,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, real(kind_phys), dimension(:), intent(inout) :: emdust, emseas, emanoc real(kind_phys), dimension(:), intent(inout) :: ebb_smoke_in,coef_bb, frp_output, fhist real(kind_phys), dimension(:,:), intent(inout) :: ebu_smoke - real(kind_phys), dimension(:,:), intent(inout) :: fire_in - real(kind_phys), dimension(:), intent(out) :: fire_heat_flux_out, frac_grid_burned_out + real(kind_phys), dimension(:), intent(out ) :: fire_heat_flux_out, frac_grid_burned_out real(kind_phys), dimension(:), intent(inout) :: max_fplume, min_fplume, uspdavg, hpbl_thetav real(kind_phys), dimension(:), intent(inout) :: hwp, peak_hr_out real(kind_phys), dimension(:), intent(inout) :: hwp_ave @@ -816,13 +815,8 @@ subroutine rrfs_smoke_prep( & vvel(i,k,j)=-w(i,kkp)*rri(i,k,j)/g moist(i,k,j,:)=0. moist(i,k,j,1)=gq0(i,kkp,1) - !if (t_phy(i,k,j) > 265.) then - moist(i,k,j,2)=gq0(i,kkp,2) - if (moist(i,k,j,2) < 1.e-8) moist(i,k,j,2)=0. - !else - ! moist(i,k,j,2)=0. - if(moist(i,k,j,3) < 1.e-8)moist(i,k,j,3)=0. - !endif + moist(i,k,j,2)=gq0(i,kkp,2) + if (moist(i,k,j,2) < 1.e-8) moist(i,k,j,2)=0. !-- zmid(i,k,j)=phl3d(i,kkp)/g enddo diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.meta b/physics/smoke_dust/rrfs_smoke_wrapper.meta index fc3aa9fe6..e00781ec1 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.meta +++ b/physics/smoke_dust/rrfs_smoke_wrapper.meta @@ -731,9 +731,9 @@ kind = kind_phys intent = inout [fire_type_out] - standard_name = fire_type_out + standard_name = fire_type long_name = type of fire - units = none + units = 1 dimensions = (horizontal_loop_extent) type = integer intent = out @@ -791,17 +791,17 @@ type = integer intent = in [uspdavg] - standard_name = bl_averaged_wind_speed + standard_name = mean_wind_speed_in_boundary_layer long_name = average wind speed within the boundary layer - units = none + units = m s-1 dimensions = (horizontal_loop_extent) type = real kind = kind_phys intent = inout [hpbl_thetav] - standard_name = pbl_height_thetav + standard_name = atmosphere_boundary_layer_thickness_from_modified_parcel long_name = pbl height based on modified parcel method - units = none + units = m dimensions = (horizontal_loop_extent) type = real kind = kind_phys @@ -861,24 +861,16 @@ type = real kind = kind_phys intent = inout -[fire_in] - standard_name = smoke_fire_auxiliary_input - long_name = smoke fire auxiliary input variables - units = various - dimensions = (horizontal_loop_extent,fire_auxiliary_data_extent) - type = real - kind = kind_phys - intent = inout [peak_hr_out] standard_name = peak_hr_fire - long_name = hour of peak fire emissions - units = none + long_name = time_of_peak_fire_emissions + units = s dimensions = (horizontal_loop_extent) type = real kind = kind_phys intent = out [lu_nofire_out] - standard_name = lu_nofire_out + standard_name = sum_of_land_use_fractions_for_no_fire_pixels long_name = land use of no fire pixels for type units = frac dimensions = (horizontal_loop_extent) @@ -886,7 +878,7 @@ kind = kind_phys intent = out [lu_qfire_out] - standard_name = lu_qfire_out + standard_name = sum_of_land_use_fractions_for_cropland_fire_pixels long_name = land use of fire pixels for type units = frac dimensions = (horizontal_loop_extent) From a0acaedeb7512f9c3cc062922b83466ffdfc2478 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Wed, 24 Jan 2024 16:35:59 +0000 Subject: [PATCH 14/14] "update to resolve code managers' comments" --- physics/smoke_dust/module_add_emiss_burn.F90 | 6 +++--- physics/smoke_dust/rrfs_smoke_wrapper.meta | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/physics/smoke_dust/module_add_emiss_burn.F90 b/physics/smoke_dust/module_add_emiss_burn.F90 index f1bbaeee9..80d91bb0e 100755 --- a/physics/smoke_dust/module_add_emiss_burn.F90 +++ b/physics/smoke_dust/module_add_emiss_burn.F90 @@ -89,8 +89,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & if (ebb_dcycle==2) then ! Constants for the fire diurnal cycle calculation - coef_con = 1._kind_phys/((2._kind_phys*pi)**0.5_kind_phys * sigma_fire_dur(1) *fire_age) * & - exp(- ( log(fire_age) - avg_fire_dur(1))**2 /(2._kind_phys*sigma_fire_dur(1)**2 )) + coef_con=1._kind_phys/((2._kind_phys*pi)**0.5_kind_phys) do j=jts,jte do i=its,ite fire_age= time_int + (fire_end_hr(i,j)-1._kind_phys)*3600._kind_phys !One hour delay is due to the latency of the RAVE files @@ -99,7 +98,8 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & SELECT CASE ( fire_type(i,j) ) !Ag, urban fires, bare land etc. CASE (1) ! these fires will have exponentially decreasing diurnal cycle, - coef_bb_dc(i,j) = coef_con + coef_bb_dc(i,j) = coef_con*1._kind_phys/(sigma_fire_dur(1) *fire_age) * & + exp(- ( log(fire_age) - avg_fire_dur(1))**2 /(2._kind_phys*sigma_fire_dur(1)**2 )) IF ( dbg_opt .AND. time_int<5000.) then WRITE(6,*) 'i,j,peak_hr(i,j) ',i,j,peak_hr(i,j) diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.meta b/physics/smoke_dust/rrfs_smoke_wrapper.meta index e00781ec1..271d2dd36 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.meta +++ b/physics/smoke_dust/rrfs_smoke_wrapper.meta @@ -872,7 +872,7 @@ [lu_nofire_out] standard_name = sum_of_land_use_fractions_for_no_fire_pixels long_name = land use of no fire pixels for type - units = frac + units = 1 dimensions = (horizontal_loop_extent) type = real kind = kind_phys @@ -880,7 +880,7 @@ [lu_qfire_out] standard_name = sum_of_land_use_fractions_for_cropland_fire_pixels long_name = land use of fire pixels for type - units = frac + units = 1 dimensions = (horizontal_loop_extent) type = real kind = kind_phys