diff --git a/.github/workflows/ci_ctests.yml b/.github/workflows/ci_ctests.yml new file mode 100644 index 000000000..3c557ecf4 --- /dev/null +++ b/.github/workflows/ci_ctests.yml @@ -0,0 +1,26 @@ +name: Physics Scheme Unit Tests (CTest) + +on: [push, pull_request] + +jobs: + ctests: + if: github.repository == 'NCAR/ccpp-physics' || github.repository == 'ufs-community/ccpp-physics' || github.repository == 'bbakernoaa/ccpp-physics' + runs-on: ubuntu-latest + + steps: + - name: Checkout code + uses: actions/checkout@v4 + + - name: Install dependencies + run: | + sudo apt-get update + sudo apt-get install -y gfortran cmake + + - name: Build and run smoke and dust tests + run: | + cd physics/smoke_dust/tests + mkdir build + cd build + cmake .. + make + ctest --output-on-failure diff --git a/physics/smoke_dust/dust_fengsha_mod.F90 b/physics/smoke_dust/dust_fengsha_mod.F90 index 6ec8f8d4a..c02ccaed1 100755 --- a/physics/smoke_dust/dust_fengsha_mod.F90 +++ b/physics/smoke_dust/dust_fengsha_mod.F90 @@ -8,9 +8,14 @@ module dust_fengsha_mod ! ! 07/16/2019 - Adapted for NUOPC/GOCART, R. Montuoro ! 02/01/2020 - Adapted for FV3/CCPP, Haiqin Li +! Refactored for Fortran 2008 compliance and best practices. use machine , only : kind_phys - use dust_data_mod + use dust_data_mod, only : ndust, reff_dust, lo_dust, up_dust, & + dust_calcdrag, dust_moist_opt, dust_alpha, dust_gamma, & + dust_moist_correction, dust_drylimit_factor, & + p_dust_1, p_dust_2, p_dust_3, p_dust_4, p_dust_5, & + p_edust1, p_edust2, p_edust3, p_edust4, p_edust5 implicit none @@ -18,8 +23,75 @@ module dust_fengsha_mod public :: gocart_dust_fengsha_driver + ! -- unified densities for FENGSHA + real(kind_phys), parameter :: rho_soil = 2650.0_kind_phys + real(kind_phys), parameter :: rho_water = 1000.0_kind_phys + + ! -- FENGSHA parameters + type :: fengsha_params_type + real(kind_phys) :: mmd_dust = 3.4e-6_kind_phys !< median mass diameter (m) + real(kind_phys) :: gsd_dust = 3.0_kind_phys !< geom. std deviation + real(kind_phys) :: lambda = 12.0e-6_kind_phys !< crack propagation length (m) + real(kind_phys) :: cv = 12.62e-6_kind_phys !< normalization constant + real(kind_phys) :: z0s = 1.0e-4_kind_phys !< Surface roughness for ideal bare surface (m) + real(kind_phys) :: clay_thresh = 0.2_kind_phys !< clay fraction threshold + real(kind_phys) :: cmb = 1.0_kind_phys !< constant of proportionality + real(kind_phys) :: kvhmax = 2.0e-4_kind_phys !< max vertical to horizontal flux ratio + real(kind_phys) :: frozen_soil_thresh = 268.0_kind_phys !< frozen soil threshold (K) + real(kind_phys) :: znt_limit = 0.2_kind_phys !< roughness length limit (m) + real(kind_phys) :: snow_limit = 0.0_kind_phys !< snow depth limit (m) + end type fengsha_params_type + + !> Parameter set for FENGSHA scheme + type(fengsha_params_type), parameter :: fparams = fengsha_params_type() + contains + !> \brief Driver for the GOCART FENGSHA dust emission scheme + !! \section arg_table_gocart_dust_fengsha_driver Arguments + !! \htmlinclude gocart_dust_fengsha_driver.html + !! \param[in] dt physics time step (s) + !! \param[inout] chem constituent mixing ratio (kg kg-1) + !! \param[in] rho_phy air density (kg m-3) + !! \param[in] smois volumetric soil moisture (m3 m-3) + !! \param[in] stemp soil temperature (K) + !! \param[in] p8w air pressure at interfaces (Pa) + !! \param[in] ssm sediment supply map (none) + !! \param[in] isltyp dominant soil type (index) + !! \param[in] snowh snow depth (m) + !! \param[in] xland land mask (1 for land, 2 for water) + !! \param[in] area grid cell area (m2) + !! \param[in] g gravitational acceleration (m s-2) + !! \param[inout] emis_dust optional dust emission flux (kg m-2 s-1) + !! \param[in] ust friction velocity (m s-1) + !! \param[in] znt surface roughness length (m) + !! \param[in] clay clay fraction (none) + !! \param[in] sand sand fraction (none) + !! \param[in] rdrag drag partition correction (none) + !! \param[in] uthr dry threshold velocity (m s-1) + !! \param[in] num_emis_dust number of dust bins + !! \param[in] num_chem number of chemistry tracers + !! \param[in] num_soil_layers number of soil layers + !! \param[in] ids horizontal dimension start index + !! \param[in] ide horizontal dimension end index + !! \param[in] jds horizontal dimension 2 start index + !! \param[in] jde horizontal dimension 2 end index + !! \param[in] kds vertical dimension start index + !! \param[in] kde vertical dimension end index + !! \param[in] ims horizontal dimension ims + !! \param[in] ime horizontal dimension ime + !! \param[in] jms horizontal dimension 2 jms + !! \param[in] jme horizontal dimension 2 jme + !! \param[in] kms vertical dimension kms + !! \param[in] kme vertical dimension kme + !! \param[in] its horizontal dimension its + !! \param[in] ite horizontal dimension ite + !! \param[in] jts horizontal dimension 2 jts + !! \param[in] jte horizontal dimension 2 jte + !! \param[in] kts vertical dimension kts + !! \param[in] kte vertical dimension kte + !! \param[out] errmsg ccpp error message + !! \param[out] errflg ccpp error code subroutine gocart_dust_fengsha_driver(dt, & chem,rho_phy,smois,stemp,p8w,ssm, & isltyp,snowh,xland,area,g,emis_dust, & @@ -27,600 +99,413 @@ subroutine gocart_dust_fengsha_driver(dt, & num_emis_dust,num_chem,num_soil_layers, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte) - IMPLICIT NONE - INTEGER, INTENT(IN ) :: & + its,ite, jts,jte, kts,kte, & + errmsg, errflg) + + integer, intent(in) :: & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & num_emis_dust,num_chem,num_soil_layers ! 2d input variables - REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: ssm ! Sediment supply map - REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: snowh ! snow height (m) - REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: xland ! dominant land use type - REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: area ! area of grid cell - REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: ust ! friction velocity - REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: znt ! Surface Roughness length (m) - REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: clay ! Clay Fraction (-) - REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: sand ! Sand Fraction (-) - REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: rdrag ! Drag Partition (-) - REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: uthr ! Dry Threshold Velocity (m/s) - - INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: isltyp ! soil type + real(kind_phys), dimension( ims:ime , jms:jme ), intent(in) :: ssm ! Sediment supply map + real(kind_phys), dimension( ims:ime , jms:jme ), intent(in) :: snowh ! snow height (m) + real(kind_phys), dimension( ims:ime , jms:jme ), intent(in) :: xland ! dominant land use type + real(kind_phys), dimension( ims:ime , jms:jme ), intent(in) :: area ! area of grid cell + real(kind_phys), dimension( ims:ime , jms:jme ), intent(in) :: ust ! friction velocity + real(kind_phys), dimension( ims:ime , jms:jme ), intent(in) :: znt ! Surface Roughness length (m) + real(kind_phys), dimension( ims:ime , jms:jme ), intent(in) :: clay ! Clay Fraction (-) + real(kind_phys), dimension( ims:ime , jms:jme ), intent(in) :: sand ! Sand Fraction (-) + real(kind_phys), dimension( ims:ime , jms:jme ), intent(in) :: rdrag ! Drag Partition (-) + real(kind_phys), dimension( ims:ime , jms:jme ), intent(in) :: uthr ! Dry Threshold Velocity (m/s) + + integer, dimension( ims:ime , jms:jme ), intent(in) :: isltyp ! soil type ! 3d input variables - REAL(kind_phys), DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN) :: p8w - 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, stemp + real(kind_phys), dimension( ims:ime , kms:kme , jms:jme ), intent(in) :: p8w + 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, stemp !0d input variables - REAL(kind_phys), INTENT(IN) :: dt ! time step - REAL(kind_phys), INTENT(IN) :: g ! gravity (m/s**2) - + real(kind_phys), intent(in) :: dt ! time step + real(kind_phys), intent(in) :: g ! gravity (m/s**2) + ! CCPP error handling + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg ! Local variables - integer :: nmx,i,j,k,imx,jmx,lmx - integer :: ilwi - real(kind_phys) :: airden ! air density - REAL(kind_phys) :: airmas ! dry air mass - real(kind_phys) :: dxy - real(kind_phys) :: conver,converi ! conversion values - real(kind_phys) :: R ! local drag partition - real(kind_phys) :: ustar - real(kind_phys), DIMENSION (num_emis_dust) :: tc - real(kind_phys), DIMENSION (num_emis_dust) :: bems - real(kind_phys), DIMENSION (num_emis_dust) :: distribution - real(kind_phys), dimension (3) :: massfrac - real(kind_phys) :: erodtot - real(kind_phys) :: moist_volumetric - - ! conversion values - conver=1.e-9 - converi=1.e9 - - ! Number of dust bins - - imx=1 - jmx=1 - lmx=1 - nmx=ndust - - k=kts - do j=jts,jte - do i=its,ite - - ! Don't do dust over water!!! - - ilwi=0 - if(xland(i,j).lt.1.5)then - ilwi=1 - - ! Total concentration at lowest model level. This is still hardcoded for 5 bins. - - ! if(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) then - ! tc(:)=1.e-16*conver - ! else - tc(1)=chem(i,kts,j,p_dust_1)*conver - tc(2)=chem(i,kts,j,p_dust_2)*conver - tc(3)=chem(i,kts,j,p_dust_3)*conver - tc(4)=chem(i,kts,j,p_dust_4)*conver - tc(5)=chem(i,kts,j,p_dust_5)*conver - ! endif - - ! Air mass and density at lowest model level. - - airmas=-(p8w(i,kts+1,j)-p8w(i,kts,j))*area(i,j)/g - airden=rho_phy(i,kts,j) - ustar=ust(i,j) - dxy=area(i,j) - - ! Mass fractions of clay, silt, and sand. - massfrac(1)=clay(i,j) - massfrac(2)=1-(clay(i,j)+sand(i,j)) - massfrac(3)=sand(i,j) - - - ! Total erodibility. - - erodtot = ssm(i,j) ! SUM(erod(i,j,:)) - - ! Don't allow roughness lengths greater than 20 cm to be lofted. - ! This kludge accounts for land use types like urban areas and - ! forests which would otherwise show up as high dust emitters. - ! This is a placeholder for a more widely accepted kludge - ! factor in the literature, which reduces lofting for rough areas. - ! Forthcoming... - - IF (znt(i,j) .gt. 0.2) then - ilwi=0 - endif + integer :: i, j + real(kind_phys), parameter :: conver = 1.0e-9_kind_phys + real(kind_phys), parameter :: converi = 1.0e9_kind_phys + + ! Initialize error handling + errmsg = '' + errflg = 0 + + ! OpenMP parallelization over the grid + !$omp parallel do collapse(2) default(shared) private(i, j) + do j = jts, jte + do i = its, ite + block + integer :: ilwi + integer :: nmx + real(kind_phys) :: airden ! air density + real(kind_phys) :: airmas ! dry air mass + real(kind_phys) :: dxy + real(kind_phys) :: R ! local drag partition + real(kind_phys) :: ustar + real(kind_phys) :: tc(num_emis_dust) + real(kind_phys) :: bems(num_emis_dust) + real(kind_phys) :: massfrac(3) + real(kind_phys) :: erodtot + real(kind_phys) :: moist_volumetric + + nmx = ndust + + ! Don't do dust over water!!! + ilwi = 0 + if (xland(i,j) < 1.5_kind_phys) then + ilwi = 1 + + ! Total concentration at lowest model level. + tc(1) = chem(i,kts,j,p_dust_1) * conver + tc(2) = chem(i,kts,j,p_dust_2) * conver + tc(3) = chem(i,kts,j,p_dust_3) * conver + tc(4) = chem(i,kts,j,p_dust_4) * conver + tc(5) = chem(i,kts,j,p_dust_5) * conver + + ! Air mass and density at lowest model level. + airmas = -(p8w(i,kts+1,j) - p8w(i,kts,j)) * area(i,j) / g + airden = rho_phy(i,kts,j) + ustar = ust(i,j) + dxy = area(i,j) + + ! Mass fractions of clay, silt, and sand. + massfrac(1) = clay(i,j) + massfrac(2) = 1.0_kind_phys - (clay(i,j) + sand(i,j)) + massfrac(3) = sand(i,j) + + ! Total erodibility. + erodtot = ssm(i,j) + + ! Don't allow roughness lengths greater than limit to be lofted. + if (znt(i,j) > fparams%znt_limit) then + ilwi = 0 + endif - ! limit where there is lots of vegetation + ! limit where there is snow on the ground + if (snowh(i,j) > fparams%snow_limit) then + ilwi = 0 + endif - ! limit where there is snow on the ground - if (snowh(i,j) .gt. 0) then - ilwi = 0 - endif + ! Don't emit over frozen soil + if (stemp(i,1,j) < fparams%frozen_soil_thresh) then + 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) == 15 .or. isltyp(i,j) == 16 .or. & + isltyp(i,j) == 18 .or. isltyp(i,j) == 0) then + ilwi = 0 + endif - ! Do not allow areas with bedrock, lava, or land-ice to loft + if (ilwi /= 0) then + ! get drag partition + if (dust_calcdrag /= 1) then + call fengsha_drag(znt(i,j), R) + else + ! use the precalculated version + if (rdrag(i,j) > 0.0_kind_phys) then + R = rdrag(i,j) + else + ilwi = 0 + endif + endif + endif - IF (isltyp(i,j) .eq. 15 .or. isltyp(i,j) .eq. 16. .or. & - isltyp(i,j) .eq. 18) then - ilwi=0 - ENDIF - IF (isltyp(i,j) .eq. 0)then - ilwi=0 - endif - if(ilwi == 0 ) cycle - - ! get drag partition - ! FENGSHA uses the drag partition correction of MacKinnon et al 2004 - ! doi:10.1016/j.geomorph.2004.03.009 - if (dust_calcdrag .ne. 1) then - call fengsha_drag(znt(i,j),R) - else - ! use the precalculated version derived from ASCAT; Prigent et al. (2012,2015) - ! doi:10.1109/TGRS.2014.2338913 & doi:10.5194/amt-5-2703-2012 - ! pick only valid values - if (rdrag(i,j) > 0.) then - R = real(rdrag(i,j), kind=kind_phys) - else - cycle + if (ilwi /= 0) then + ! soil moisture correction factor + moist_volumetric = dust_moist_correction * smois(i,2,j) + + ! Call dust emission routine. + call source_dust(nmx, dt, tc, ustar, massfrac, & + erodtot, dxy, moist_volumetric, airden, airmas, bems, g, dust_alpha, dust_gamma, & + R, uthr(i,j)) + + ! convert back to concentration + chem(i,kts,j,p_dust_1) = tc(1) * converi + chem(i,kts,j,p_dust_2) = tc(2) * converi + chem(i,kts,j,p_dust_3) = tc(3) * converi + chem(i,kts,j,p_dust_4) = tc(4) * converi + chem(i,kts,j,p_dust_5) = tc(5) * converi + + ! For output diagnostics + if (present(emis_dust)) then + emis_dust(i,1,j,p_edust1) = bems(1) + emis_dust(i,1,j,p_edust2) = bems(2) + emis_dust(i,1,j,p_edust3) = bems(3) + emis_dust(i,1,j,p_edust4) = bems(4) + emis_dust(i,1,j,p_edust5) = bems(5) + endif endif endif - - ! soil moisture correction factor - moist_volumetric = dust_moist_correction * smois(i,2,j) - - ! Call dust emission routine. - - call source_dust(imx,jmx, lmx, nmx, dt, tc, ustar, massfrac, & - erodtot, dxy, moist_volumetric, airden, airmas, bems, g, dust_alpha, dust_gamma, & - R, uthr(i,j)) - - ! convert back to concentration - - chem(i,kts,j,p_dust_1)=tc(1)*converi - chem(i,kts,j,p_dust_2)=tc(2)*converi - chem(i,kts,j,p_dust_3)=tc(3)*converi - chem(i,kts,j,p_dust_4)=tc(4)*converi - chem(i,kts,j,p_dust_5)=tc(5)*converi - - ! For output diagnostics - - emis_dust(i,1,j,p_edust1)=bems(1) - emis_dust(i,1,j,p_edust2)=bems(2) - emis_dust(i,1,j,p_edust3)=bems(3) - emis_dust(i,1,j,p_edust4)=bems(4) - emis_dust(i,1,j,p_edust5)=bems(5) - endif + end block enddo enddo - ! end subroutine gocart_dust_fengsha_driver - subroutine source_dust(imx, jmx, lmx, nmx, dt1, tc, ustar, massfrac, & + !> \brief Evaluates the source of each dust particle size bin + !! \param[in] nmx number of dust bins + !! \param[in] dt1 time step (s) + !! \param[inout] tc total concentration of dust (kg kg-1) + !! \param[in] ustar friction velocity (m s-1) + !! \param[in] massfrac fraction of mass in each of 3 soil classes + !! \param[in] erod fraction of erodible grid cell + !! \param[in] dxy grid cell area (m2) + !! \param[in] smois volumetric soil moisture (m3 m-3) + !! \param[in] airden density of air (kg m-3) + !! \param[in] airmas mass of air for each grid box (kg) + !! \param[out] bems source of each dust type (ug m-2 s-1) + !! \param[in] g0 gravitational acceleration (m s-2) + !! \param[in] alpha scaling factor + !! \param[in] gamma scaling factor + !! \param[in] R drag partition + !! \param[in] uthres dry threshold velocity (m s-1) + subroutine source_dust(nmx, dt1, tc, ustar, massfrac, & erod, dxy, smois, airden, airmas, bems, g0, alpha, gamma, & R, uthres) - ! **************************************************************************** - ! * Evaluate the source of each dust particles size bin by soil emission - ! * - ! * Input: - ! * EROD Fraction of erodible grid cell (-) - ! * smois Volumetric soil moisture (m3/m3) - ! * ALPHA Constant to fudge the total emission of dust (1/m) - ! * GAMMA Tuning constant for erodibility (-) - ! * DXY Surface of each grid cell (m2) - ! * AIRMAS Mass of air for each grid box (kg) - ! * AIRDEN Density of air for each grid box (kg/m3) - ! * USTAR Friction velocity (m/s) - ! * DT1 Time step (s) - ! * NMX Number of dust bins (-) - ! * IMX Number of I points (-) - ! * JMX Number of J points (-) - ! * LMX Number of L points (-) - ! * R Drag Partition (-) - ! * UTHRES FENGSHA Dry Threshold Velocities (m/s) - ! * - ! * Data: - ! * MASSFRAC Fraction of mass in each of 3 soil classes (-) (clay silt sand) - ! * DEN_DUST Dust density (kg/m3) - ! * DEN_SALT Saltation particle density (kg/m3) - ! * REFF_SALT Reference saltation particle diameter (m) - ! * REFF_DUST Reference dust particle diameter (m) - ! * LO_DUST Lower diameter limits for dust bins (m) - ! * UP_DUST Upper diameter limits for dust bins (m) - ! * FRAC_SALT Soil class mass fraction for saltation bins (-) - ! * - ! * Parameters: - ! * CMB Constant of proportionality (-) - ! * MMD_DUST Mass median diameter of dust (m) - ! * GSD_DUST Geometric standard deviation of dust (-) - ! * LAMBDA Side crack propagation length (m) - ! * CV Normalization constant (-) - ! * G0 Gravitational acceleration (m/s2) - ! * - ! * Working: - ! * RHOA Density of air in cgs (g/cm3) - ! * DS_REL Saltation surface area distribution (-) - ! * DLNDP Dust bin width (-) - ! * EMIT Total vertical mass flux (kg/m2/s) - ! * EMIT_VOL Total vertical volume flux (m/s) - ! * DSRC Mass of emitted dust (kg/timestep/cell) - ! * - ! * Output: - ! * TC Total concentration of dust (kg/kg/timestep/cell) - ! * BEMS Source of each dust type (kg/timestep/cell) - ! * - ! **************************************************************************** - implicit none - - ! Input - INTEGER, INTENT(IN) :: imx,jmx,lmx,nmx - REAL(kind_phys), INTENT(IN) :: dt1 - REAL(kind_phys), INTENT(IN) :: ustar - REAL(kind_phys), INTENT(IN) :: massfrac(3) - REAL(kind_phys), INTENT(IN) :: erod - REAL(kind_phys), INTENT(IN) :: dxy - REAL(kind_phys), INTENT(IN) :: smois - REAL(kind_phys), INTENT(IN) :: airden - REAL(kind_phys), INTENT(IN) :: airmas - REAL(kind_phys), INTENT(IN) :: g0 - REAL(kind_phys), INTENT(IN) :: alpha - REAL(kind_phys), INTENT(IN) :: gamma - REAL(kind_phys), INTENT(IN) :: R - REAL(kind_phys), INTENT(IN) :: uthres + integer, intent(in) :: nmx + real(kind_phys), intent(in) :: dt1 + real(kind_phys), intent(in) :: ustar + real(kind_phys), intent(in) :: massfrac(3) + real(kind_phys), intent(in) :: erod + real(kind_phys), intent(in) :: dxy + real(kind_phys), intent(in) :: smois + real(kind_phys), intent(in) :: airden + real(kind_phys), intent(in) :: airmas + real(kind_phys), intent(in) :: g0 + real(kind_phys), intent(in) :: alpha + real(kind_phys), intent(in) :: gamma + real(kind_phys), intent(in) :: R + real(kind_phys), intent(in) :: uthres ! Output - REAL(kind_phys), INTENT(INOUT) :: tc(nmx) + real(kind_phys), intent(inout) :: tc(nmx) + real(kind_phys), intent(out) :: bems(nmx) ! Local Variables - REAL(kind_phys), INTENT(OUT) :: bems(nmx) - - REAL(kind_phys) :: dvol(nmx) - REAL(kind_phys) :: distr_dust(nmx) - REAL(kind_phys) :: dlndp(nmx) - REAL(kind_phys) :: dsrc - REAL(kind_phys) :: dvol_tot - REAL(kind_phys) :: emit - REAL(kind_phys) :: emit_vol - REAL(kind_phys) :: rhoa - INTEGER :: i, j, n - - ! Constant of proportionality from Marticorena et al, 1997 (unitless) - ! Arguably more ~consistent~ fudge than alpha, which has many walnuts - ! sprinkled throughout the literature. - GC - - REAL(kind_phys), PARAMETER :: cmb=1.0 - REAL(kind_phys), PARAMETER :: kvhmax=2.0e-4 - - ! Parameters used in Kok distribution function. Advise not to play with - ! these without the expressed written consent of someone who knows what - ! they're doing. - GC - - REAL(kind_phys), PARAMETER :: mmd_dust=3.4D-6 ! median mass diameter (m) - REAL(kind_phys), PARAMETER :: gsd_dust=3.0 ! geom. std deviation - REAL(kind_phys), PARAMETER :: lambda=12.0D-6 ! crack propagation length (m) - REAL(kind_phys), PARAMETER :: cv=12.62D-6 ! normalization constant - REAL(kind_phys), PARAMETER :: RHOSOIL=2650. - + real(kind_phys) :: dvol(nmx) + real(kind_phys) :: distr_dust(nmx) + real(kind_phys) :: dlndp(nmx) + real(kind_phys) :: dsrc + real(kind_phys) :: dvol_tot + real(kind_phys) :: emit + integer :: n ! calculate the total vertical dust flux + emit = 0.0_kind_phys - emit = 0.0 - - call DustEmissionFENGSHA(smois,massfrac(1),massfrac(3), massfrac(2), & - erod, R, airden, ustar, uthres, alpha, gamma, kvhmax, & - g0, RHOSOIL, emit) + call DustEmissionFENGSHA(smois, massfrac(1), massfrac(3), massfrac(2), & + erod, R, airden, ustar, uthres, alpha, gamma, & + g0, emit) ! Now that we have the total dust emission, distribute into dust bins using - ! lognormal distribution (Dr. Jasper Kok, in press), and - ! calculate total mass emitted over the grid box over the timestep. - ! - ! In calculating the Kok distribution, we assume upper and lower limits to each bin. - ! For reff_dust=(/0.73D-6,1.4D-6,2.4D-6,4.5D-6,8.0D-6/) (default), - ! lower limits were ASSUMED at lo_dust=(/0.1D-6,1.0D-6,1.8D-6,3.0D-6,6.0D-6/) - ! upper limits were ASSUMED at up_dust=(/1.0D-6,1.8D-6,3.0D-6,6.0D-6,10.0D-6/) - ! These may be changed within module_data_gocart_dust.F, but make sure it is - ! consistent with reff_dust values. These values were taken from the original - ! GOCART bin configuration. We use them here to calculate dust bin width, dlndp. - ! dVol is the volume distribution. You know...if you were wondering. GC - - dvol_tot=0. - DO n=1,nmx - dlndp(n)=LOG(up_dust(n)/lo_dust(n)) - dvol(n)=(2.0*reff_dust(n)/cv)*(1.+ERF(LOG(2.0*reff_dust(n)/mmd_dust)/(SQRT(2.)*LOG(gsd_dust))))*& - EXP(-(2.0*reff_dust(n)/lambda)**3.0)*dlndp(n) - dvol_tot=dvol_tot+dvol(n) - ! Convert mass flux to volume flux - !emit_vol=emit/den_dust(n) ! (m s^-1) - END DO - DO n=1,nmx - distr_dust(n)=dvol(n)/dvol_tot - !print *,"distr_dust(",n,")=",distr_dust(n) - END DO + ! lognormal distribution (Dr. Jasper Kok) + dvol_tot = 0.0_kind_phys + do n = 1, nmx + dlndp(n) = log(up_dust(n) / lo_dust(n)) + dvol(n) = (2.0_kind_phys * reff_dust(n) / fparams%cv) * & + (1.0_kind_phys + erf(log(2.0_kind_phys * reff_dust(n) / fparams%mmd_dust) / & + (sqrt(2.0_kind_phys) * log(fparams%gsd_dust)))) * & + exp(-(2.0_kind_phys * reff_dust(n) / fparams%lambda)**3.0_kind_phys) * dlndp(n) + dvol_tot = dvol_tot + dvol(n) + end do + + do n = 1, nmx + distr_dust(n) = dvol(n) / dvol_tot + end do ! Now distribute total vertical emission into dust bins and update concentration. - - DO n=1,nmx + do n = 1, nmx ! Calculate total mass emitted - dsrc = emit*distr_dust(n)*dxy*dt1 ! (kg) - IF (dsrc < 0.0) dsrc = 0.0 + dsrc = emit * distr_dust(n) * dxy * dt1 ! (kg) + if (dsrc < 0.0_kind_phys) dsrc = 0.0_kind_phys ! Update dust mixing ratio at first model level. tc(n) = tc(n) + dsrc / airmas ! (kg/kg) - ! bems(i,j,n) = dsrc ! diagnostic - !bems(i,j,n) = 1000.*dsrc/(dxy(j)*dt1) ! diagnostic (g/m2/s) - bems(n) = 1.e+9*dsrc/(dxy*dt1) ! diagnostic (ug/m2/s) !lzhang - - END DO - tc(1)=tc(1)+0.286*tc(2) ! This is just for RRFS-SD. DO NOT use in other models!!! - tc(5)=0.714*tc(2)+tc(3)+tc(4) ! This is just for RRFS-SD. DO NOT use in other models!!! + bems(n) = 1.0e9_kind_phys * dsrc / (dxy * dt1) ! diagnostic (ug/m2/s) + end do - END SUBROUTINE source_dust + ! RRFS-SD Kludge + tc(1) = tc(1) + 0.286_kind_phys * tc(2) + tc(5) = 0.714_kind_phys * tc(2) + tc(3) + tc(4) + end subroutine source_dust - subroutine fengsha_drag(z0,R) - implicit none - real(kind_phys), intent(in) :: z0 + !> \brief Calculates the MacKinnon et al. 2004 Drag Partition Correction + !! \param[in] z0 surface roughness length (m) + !! \param[out] R drag partition correction + subroutine fengsha_drag(z0, R) + real(kind_phys), intent(in) :: z0 real(kind_phys), intent(out) :: R - real(kind_phys), parameter :: z0s = 1.0e-04 !Surface roughness for ideal bare surface [m] - ! ------------------------------------------------------------------------ - ! Function: Calculates the MacKinnon et al. 2004 Drag Partition Correction - ! - ! R = 1.0 - log(z0 / z0s) / log( 0.7 * (12255./z0s) ** 0.8) - ! - !-------------------------------------------------------------------------- - ! Drag partition correction. See MacKinnon et al. (2004), - ! doi:10.1016/j.geomorph.2004.03.009 - R = 1.0 - log(z0 / z0s) / log( 0.7 * (12255./z0s) ** 0.8) - ! Drag partition correction. See Marticorena et al. (1997), - ! doi:10.1029/96JD02964 - !R = 1.0 - log(z0 / z0s) / log( 0.7 * (10./z0s) ** 0.8) + ! Drag partition correction. See MacKinnon et al. (2004), + ! doi:10.1016/j.geomorph.2004.03.009 + R = 1.0_kind_phys - log(z0 / fparams%z0s) / log(0.7_kind_phys * (12255.0_kind_phys / fparams%z0s)**0.8_kind_phys) - return end subroutine fengsha_drag + !> \brief Computes dust emissions using NOAA/ARL FENGSHA model + !! \param[in] slc volumetric soil moisture fraction + !! \param[in] clay fractional clay content + !! \param[in] sand fractional sand content + !! \param[in] silt fractional silt content + !! \param[in] ssm erosion map + !! \param[in] rdrag drag partition + !! \param[in] airdens air density at lowest level (kg m-3) + !! \param[in] ustar friction velocity (m s-1) + !! \param[in] uthrs threshold velocity (m s-1) + !! \param[in] alpha scaling factor + !! \param[in] gamma scaling factor + !! \param[in] grav gravitational acceleration (m s-2) + !! \param[inout] emissions total surface emissions (kg m-2 s-1) subroutine DustEmissionFENGSHA(slc, clay, sand, silt, & ssm, rdrag, airdens, ustar, uthrs, alpha, gamma, & - kvhmax, grav, rhop, emissions) - - ! !USES: - implicit NONE - -! !INPUT PARAMETERS: - REAL(kind_phys), intent(in) :: slc ! liquid water content of soil layer, volumetric fraction [1] - REAL(kind_phys), intent(in) :: clay ! fractional clay content [1] - REAL(kind_phys), intent(in) :: sand ! fractional sand content [1] - REAL(kind_phys), intent(in) :: silt ! fractional silt content [1] - REAL(kind_phys), intent(in) :: ssm ! erosion map [1] - REAL(kind_phys), intent(in) :: rdrag ! drag partition [1/m] - REAL(kind_phys), intent(in) :: airdens ! air density at lowest level [kg/m^3] - REAL(kind_phys), intent(in) :: ustar ! friction velocity [m/sec] - REAL(kind_phys), intent(in) :: uthrs ! threshold velocity [m/2] - REAL(kind_phys), intent(in) :: alpha ! scaling factor [1] - REAL(kind_phys), intent(in) :: gamma ! scaling factor [1] - REAL(kind_phys), intent(in) :: kvhmax ! max. vertical to horizontal mass flux ratio [1] - REAL(kind_phys), intent(in) :: grav ! gravity [m/sec^2] - REAL(kind_phys), intent(in) :: rhop ! soil class density [kg/m^3] + grav, emissions) - ! !OUTPUT PARAMETERS: - REAL(kind_phys), intent(inout) :: emissions ! binned surface emissions [kg/(m^2 sec)] + real(kind_phys), intent(in) :: slc ! liquid water content of soil layer, volumetric fraction [1] + real(kind_phys), intent(in) :: clay ! fractional clay content [1] + real(kind_phys), intent(in) :: sand ! fractional sand content [1] + real(kind_phys), intent(in) :: silt ! fractional silt content [1] + real(kind_phys), intent(in) :: ssm ! erosion map [1] + real(kind_phys), intent(in) :: rdrag ! drag partition [1/m] + real(kind_phys), intent(in) :: airdens ! air density at lowest level [kg/m^3] + real(kind_phys), intent(in) :: ustar ! friction velocity [m/sec] + real(kind_phys), intent(in) :: uthrs ! threshold velocity [m/s] + real(kind_phys), intent(in) :: alpha ! scaling factor [1] + real(kind_phys), intent(in) :: gamma ! scaling factor [1] + real(kind_phys), intent(in) :: grav ! gravity [m/sec^2] - ! !DESCRIPTION: Compute dust emissions using NOAA/ARL FENGSHA model - ! - ! !REVISION HISTORY: - ! - ! 22Feb2020 B.Baker/NOAA - Original implementation - ! 29Mar2021 R.Montuoro/NOAA - Refactored for process library - ! 09Aug2022 B.Baker/NOAA - Adapted for CCPP-Physics + real(kind_phys), intent(inout) :: emissions ! surface emissions [kg/(m^2 sec)] - ! !Local Variables - real(kind_phys) :: alpha_grav - real(kind_phys) :: h - real(kind_phys) :: kvh - real(kind_phys) :: q - real(kind_phys) :: rustar - real(kind_phys) :: total_emissions - real(kind_phys) :: u_sum, u_thresh + ! Local Variables + real(kind_phys) :: alpha_grav + real(kind_phys) :: h + real(kind_phys) :: kvh + real(kind_phys) :: q + real(kind_phys) :: rustar + real(kind_phys) :: u_sum, u_thresh -!EOP -!------------------------------------------------------------------------- -! Begin - -! Initialize emissions -! -------------------- - emissions = 0. - -! Prepare scaling factor -! ---------------------- - alpha_grav = alpha / grav - - ! Compute vertical-to-horizontal mass flux ratio - ! ---------------------------------------------- - kvh = DustFluxV2HRatioMB95(clay, kvhmax) - - ! Compute total emissions - ! ----------------------- - emissions = alpha_grav * (ssm ** gamma) * airdens * kvh - - ! Compute threshold wind friction velocity using drag partition - ! ------------------------------------------------------------- - rustar = rdrag * ustar - - ! Now compute size-dependent total emission flux - ! ---------------------------------------------- - - if (dust_moist_opt .eq. 1) then - - ! Fecan moisture correction - ! ------------------------- - h = moistureCorrectionFecan(slc, sand, clay) - else - ! shao soil moisture correction - h = moistureCorrectionShao(slc) - end if - ! Adjust threshold - ! ---------------- - u_thresh = uthrs * h - - u_sum = rustar + u_thresh - - ! Compute Horizontal Saltation Flux according to Eq (9) in Webb et al. (2020) - ! --------------------------------------------------------------------------- - q = max(0., rustar - u_thresh) * u_sum * u_sum - - ! Distribute emissions to bins and convert to mass flux (kg s-1) - ! -------------------------------------------------------------- - emissions = emissions * q - - - end subroutine DustEmissionFENGSHA -!----------------------------------------------------------------- - real function soilMoistureConvertVol2Grav(vsoil, sandfrac) - -! !USES: - implicit NONE + emissions = 0.0_kind_phys + alpha_grav = alpha / grav -! !INPUT PARAMETERS: - REAL(kind_phys), intent(in) :: vsoil ! volumetric soil moisture fraction [1] - REAL(kind_phys), intent(in) :: sandfrac ! fractional sand content [1] + ! Compute vertical-to-horizontal mass flux ratio + kvh = DustFluxV2HRatioMB95(clay) -! !DESCRIPTION: Convert soil moisture fraction from volumetric to gravimetric. -! -! !REVISION HISTORY: -! -! 02Apr2020, B.Baker/NOAA - Original implementation -! 01Apr2020, R.Montuoro/NOAA - Adapted for GOCART process library + ! Compute total emissions base + emissions = alpha_grav * (ssm ** gamma) * airdens * kvh -! !Local Variables - real :: vsat + ! Compute threshold wind friction velocity using drag partition + rustar = rdrag * ustar -! !CONSTANTS: - REAL(kind_phys), parameter :: rhow = 1000. ! density of water [kg m-3] - REAL(kind_phys), parameter :: rhop = 1700. ! density of dry soil -!EOP -!------------------------------------------------------------------------- -! Begin... + if (dust_moist_opt == 1) then + ! Fecan moisture correction + h = moistureCorrectionFecan(slc, sand, clay) + else + ! shao soil moisture correction + h = moistureCorrectionShao(slc) + end if -! Saturated volumetric water content (sand-dependent) ! [m3 m-3] - vsat = 0.489 - 0.126 * sandfrac + ! Adjust threshold + u_thresh = uthrs * h + u_sum = rustar + u_thresh + + ! Compute Horizontal Saltation Flux according to Eq (9) in Webb et al. (2020) + q = max(0.0_kind_phys, rustar - u_thresh) * u_sum * u_sum + + ! Distribute emissions and convert to mass flux (kg m-2 s-1) + emissions = emissions * q + + end subroutine DustEmissionFENGSHA + + !> \brief Convert soil moisture fraction from volumetric to gravimetric + !! \param[in] vsoil volumetric soil moisture fraction + !! \param[in] sandfrac fractional sand content + !! \return gravimetric soil moisture fraction + function soilMoistureConvertVol2Grav(vsoil, sandfrac) result(res) + real(kind_phys), intent(in) :: vsoil ! volumetric soil moisture fraction [1] + real(kind_phys), intent(in) :: sandfrac ! fractional sand content [1] + real(kind_phys) :: res + real(kind_phys) :: vsat -! Gravimetric soil content - soilMoistureConvertVol2Grav = 100.0 * (vsoil * rhow / rhop / ( 1. - vsat)) - - end function soilMoistureConvertVol2Grav -!---------------------------------------------------------------- - real function moistureCorrectionFecan(slc, sand, clay) - -! !USES: - implicit NONE - -! !INPUT PARAMETERS: - REAL(kind_phys), intent(in) :: slc ! liquid water content of top soil layer, volumetric fraction [1] - REAL(kind_phys), intent(in) :: sand ! fractional sand content [1] - REAL(kind_phys), intent(in) :: clay ! fractional clay content [1] - -! !DESCRIPTION: Compute correction factor to account for Fecal soil moisture -! -! !REVISION HISTORY: -! -! 02Apr2020, B.Baker/NOAA - Original implementation -! 01Apr2020, R.Montuoro/NOAA - Adapted for GOCART process library + ! Saturated volumetric water content (sand-dependent) [m3 m-3] + vsat = 0.489_kind_phys - 0.126_kind_phys * sandfrac -! !Local Variables - real :: grvsoilm - real :: drylimit + ! Gravimetric soil content + res = 100.0_kind_phys * (vsoil * rho_water / rho_soil / (1.0_kind_phys - vsat)) -!EOP -!--------------------------------------------------------------- -! Begin... + end function soilMoistureConvertVol2Grav -! Convert soil moisture from volumetric to gravimetric + !> \brief Compute correction factor to account for Fecan soil moisture + !! \param[in] slc liquid water content of top soil layer + !! \param[in] sand fractional sand content + !! \param[in] clay fractional clay content + !! \return correction factor + function moistureCorrectionFecan(slc, sand, clay) result(res) + real(kind_phys), intent(in) :: slc ! liquid water content of top soil layer + real(kind_phys), intent(in) :: sand ! fractional sand content + real(kind_phys), intent(in) :: clay ! fractional clay content + real(kind_phys) :: res + + real(kind_phys) :: grvsoilm + real(kind_phys) :: drylimit + + ! Convert soil moisture from volumetric to gravimetric grvsoilm = soilMoistureConvertVol2Grav(slc, sand) -! Compute fecan dry limit - drylimit = dust_drylimit_factor * clay * (14.0 * clay + 17.0) + ! Compute fecan dry limit + drylimit = dust_drylimit_factor * clay * (14.0_kind_phys * clay + 17.0_kind_phys) -! Compute soil moisture correction - moistureCorrectionFecan = sqrt(1.0 + 1.21 * max(0., grvsoilm - drylimit)**0.68) + ! Compute soil moisture correction + res = sqrt(1.0_kind_phys + 1.21_kind_phys * max(0.0_kind_phys, grvsoilm - drylimit)**0.68_kind_phys) end function moistureCorrectionFecan -!---------------------------------------------------------------- - real function moistureCorrectionShao(slc) - -! !USES: - implicit NONE - -! !INPUT PARAMETERS: - REAL(kind_phys), intent(in) :: slc ! liquid water content of top soil layer, volumetric fraction [1] - -! !DESCRIPTION: Compute correction factor to account for Fecal soil moisture -! -! !REVISION HISTORY: -! -! 02Apr2020, B.Baker/NOAA - Original implementation -! 01Apr2020, R.Montuoro/NOAA - Adapted for GOCART process library -! !Local Variables - real :: grvsoilm - real :: drylimit + !> \brief Compute correction factor to account for Shao soil moisture + !! \param[in] slc liquid water content of top soil layer + !! \return correction factor + function moistureCorrectionShao(slc) result(res) + real(kind_phys), intent(in) :: slc + real(kind_phys) :: res -!EOP -!--------------------------------------------------------------- -! Begin... - - if (slc < 0.03) then - moistureCorrectionShao = exp(22.7 * slc) + if (slc < 0.03_kind_phys) then + res = exp(22.7_kind_phys * slc) else - moistureCorrectionShao = exp(95.3 * slc - 2.029) + res = exp(95.3_kind_phys * slc - 2.029_kind_phys) end if end function moistureCorrectionShao -!--------------------------------------------------------------- - real function DustFluxV2HRatioMB95(clay, kvhmax) - -! !USES: - implicit NONE - -! !INPUT PARAMETERS: - REAL(kind_phys), intent(in) :: clay ! fractional clay content [1] - REAL(kind_phys), intent(in) :: kvhmax ! maximum flux ratio [1] - -! !CONSTANTS: - REAL(kind_phys), parameter :: clay_thresh = 0.2 ! clay fraction above which the maximum flux ratio is returned -! !DESCRIPTION: Computes the vertical-to-horizontal dust flux ratio according to -! B.Marticorena, G.Bergametti, J.Geophys.Res., 100(D8), 164! doi:10.1029/95JD00690 -! -! !REVISION HISTORY: -! -! 22Feb2020 B.Baker/NOAA - Original implementation -! 01Apr2021 R.Montuoro/NOAA - Adapted for GOCART process library -! -!EOP -!------------------------------------------------------------------------- -! Begin... + !> \brief Computes the vertical-to-horizontal dust flux ratio + !! \param[in] clay fractional clay content + !! \return flux ratio + function DustFluxV2HRatioMB95(clay) result(res) + real(kind_phys), intent(in) :: clay ! fractional clay content + real(kind_phys) :: res - if (clay > clay_thresh) then - DustFluxV2HRatioMB95 = kvhmax + if (clay > fparams%clay_thresh) then + res = fparams%kvhmax else - DustFluxV2HRatioMB95 = 10.0**(13.4*clay-6.0) + res = 10.0_kind_phys**(13.4_kind_phys * clay - 6.0_kind_phys) end if end function DustFluxV2HRatioMB95 diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.F90 b/physics/smoke_dust/rrfs_smoke_wrapper.F90 index 7254304b8..afd68fc60 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.F90 +++ b/physics/smoke_dust/rrfs_smoke_wrapper.F90 @@ -421,7 +421,9 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land num_emis_dust,num_chem,nsoil, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte) + its,ite, jts,jte, kts,kte, & + errmsg, errflg) + if (errflg /= 0) return end if ! compute wild-fire plumes diff --git a/physics/smoke_dust/tests/CMakeLists.txt b/physics/smoke_dust/tests/CMakeLists.txt new file mode 100644 index 000000000..9eeeb737c --- /dev/null +++ b/physics/smoke_dust/tests/CMakeLists.txt @@ -0,0 +1,53 @@ +cmake_minimum_required(VERSION 3.10) +project(smoke_dust_tests LANGUAGES Fortran) + +# Find OpenMP +find_package(OpenMP) + +# Paths to source files +set(HOOKS_DIR ${CMAKE_CURRENT_SOURCE_DIR}/../../hooks) +set(SMOKE_DUST_DIR ${CMAKE_CURRENT_SOURCE_DIR}/..) + +# FENGSHA Test +set(FENGSHA_SOURCES + ${HOOKS_DIR}/machine.F + ${SMOKE_DUST_DIR}/rrfs_smoke_config.F90 + ${SMOKE_DUST_DIR}/dust_data_mod.F90 + ${SMOKE_DUST_DIR}/dust_fengsha_mod.F90 + test_dust_fengsha.F90 +) +add_executable(test_dust_fengsha ${FENGSHA_SOURCES}) + +# GOCART Sea Salt Test +set(SEAS_SOURCES + ${HOOKS_DIR}/machine.F + ${SMOKE_DUST_DIR}/seas_data_mod.F90 + ${SMOKE_DUST_DIR}/seas_ngac_mod.F90 + ${SMOKE_DUST_DIR}/seas_mod.F90 + test_seas_gocart.F90 +) +add_executable(test_seas_gocart ${SEAS_SOURCES}) + +# Large-scale Wet Deposition Test +set(WETDEP_SOURCES + ${HOOKS_DIR}/machine.F + ${SMOKE_DUST_DIR}/rrfs_smoke_config.F90 + ${SMOKE_DUST_DIR}/module_wetdep_ls.F90 + test_wetdep_ls.F90 +) +add_executable(test_wetdep_ls ${WETDEP_SOURCES}) + +# Compile flags and OpenMP linking +foreach(target test_dust_fengsha test_seas_gocart test_wetdep_ls) + set_target_properties(${target} PROPERTIES + Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/modules + ) + if(OpenMP_Fortran_FOUND) + target_link_libraries(${target} PUBLIC OpenMP::OpenMP_Fortran) + endif() +endforeach() + +enable_testing() +add_test(NAME run_fengsha_test COMMAND test_dust_fengsha) +add_test(NAME run_seas_test COMMAND test_seas_gocart) +add_test(NAME run_wetdep_test COMMAND test_wetdep_ls) diff --git a/physics/smoke_dust/tests/test_dust_fengsha.F90 b/physics/smoke_dust/tests/test_dust_fengsha.F90 new file mode 100644 index 000000000..cee8be431 --- /dev/null +++ b/physics/smoke_dust/tests/test_dust_fengsha.F90 @@ -0,0 +1,146 @@ +!> \file test_dust_fengsha.F90 +!! Comprehensive unit test for the FENGSHA dust emission scheme. + +program test_dust_fengsha + use machine, only : kind_phys + use dust_fengsha_mod, only : gocart_dust_fengsha_driver + implicit none + + ! Parameters + integer, parameter :: im = 1, jm = 1, km = 1 + integer, parameter :: num_chem = 20 + integer, parameter :: num_emis_dust = 5 + integer, parameter :: num_soil_layers = 2 + integer, parameter :: p_dust_1 = 10 + integer, parameter :: p_dust_2 = 11 + integer, parameter :: p_dust_3 = 12 + integer, parameter :: p_dust_4 = 13 + integer, parameter :: p_dust_5 = 14 + + ! Global data for tests + real(kind_phys) :: dt, g + real(kind_phys), dimension(im, km, jm, num_chem) :: chem + real(kind_phys), dimension(im, km, jm) :: rho_phy + real(kind_phys), dimension(im, num_soil_layers, jm) :: smois, stemp + real(kind_phys), dimension(im, km+1, jm) :: p8w + real(kind_phys), dimension(im, jm) :: ssm, snowh, xland, area, ust, znt, clay, sand, rdrag, uthr + integer, dimension(im, jm) :: isltyp + real(kind_phys), dimension(im, 1, jm, num_emis_dust) :: emis_dust + + character(len=128) :: errmsg + integer :: errflg + integer :: total_passed, total_failed + + total_passed = 0 + total_failed = 0 + + dt = 600.0_kind_phys + g = 9.80665_kind_phys + + print *, "Starting FENGSHA comprehensive unit test suite..." + + ! Case 1: Standard land emission (should be non-zero) + call reset_inputs() + call run_case("Standard land emission", .true.) + + ! Case 2: Water (xland=2.0) (should be zero) + call reset_inputs() + xland = 2.0_kind_phys + call run_case("Water (xland=2.0)", .false.) + + ! Case 3: Frozen soil (stemp < 268) (should be zero) + call reset_inputs() + stemp = 260.0_kind_phys + call run_case("Frozen soil", .false.) + + ! Case 4: Snow cover (snowh > 0) (should be zero) + call reset_inputs() + snowh = 0.1_kind_phys + call run_case("Snow cover", .false.) + + ! Case 5: Roughness too high (znt > 0.2) (should be zero) + call reset_inputs() + znt = 0.3_kind_phys + call run_case("High roughness length", .false.) + + ! Case 6: Invalid soil type (isltyp=15) (should be zero) + call reset_inputs() + isltyp = 15 + call run_case("Invalid soil type (15)", .false.) + + ! Case 7: Wind below threshold (should be zero) + call reset_inputs() + ust = 0.01_kind_phys + call run_case("Wind below threshold", .false.) + + ! Case 8: Very dry soil (should have higher emission than moist) + call reset_inputs() + smois = 0.01_kind_phys + call run_case("Very dry soil", .true.) + + print *, "---------------------------------------" + print *, "Test Summary: ", total_passed, " PASSED, ", total_failed, " FAILED" + print *, "---------------------------------------" + + if (total_failed > 0) stop 1 + +contains + + subroutine reset_inputs() + chem = 1.0e-12_kind_phys + rho_phy = 1.2_kind_phys + smois = 0.1_kind_phys + stemp = 300.0_kind_phys + p8w(:,1,:) = 101325.0_kind_phys + p8w(:,2,:) = 100000.0_kind_phys + ssm = 0.5_kind_phys + isltyp = 1 + snowh = 0.0_kind_phys + xland = 1.0_kind_phys ! land + area = 1000000.0_kind_phys + ust = 1.0_kind_phys + znt = 0.01_kind_phys + clay = 0.2_kind_phys + sand = 0.4_kind_phys + rdrag = 0.8_kind_phys + uthr = 0.3_kind_phys + emis_dust = 0.0_kind_phys + errmsg = '' + errflg = 0 + end subroutine reset_inputs + + subroutine run_case(name, expect_non_zero) + character(len=*), intent(in) :: name + logical, intent(in) :: expect_non_zero + logical :: is_non_zero + real(kind_phys) :: total_emis + + call gocart_dust_fengsha_driver(dt, & + 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, & + 1, im, 1, jm, 1, km, & + 1, im, 1, jm, 1, km, & + 1, im, 1, jm, 1, km, & + errmsg, errflg) + + if (errflg /= 0) then + print *, "FAIL: ", name, " - Driver returned error: ", trim(errmsg) + total_failed = total_failed + 1 + return + endif + + total_emis = sum(emis_dust(1,1,1,:)) + is_non_zero = total_emis > 1.0e-20_kind_phys + + if (is_non_zero .eqv. expect_non_zero) then + print *, "PASS: ", name, " (Emis: ", total_emis, ")" + total_passed = total_passed + 1 + else + print *, "FAIL: ", name, " (Emis: ", total_emis, ", Expected non-zero: ", expect_non_zero, ")" + total_failed = total_failed + 1 + endif + end subroutine run_case + +end program test_dust_fengsha diff --git a/physics/smoke_dust/tests/test_seas_gocart.F90 b/physics/smoke_dust/tests/test_seas_gocart.F90 new file mode 100644 index 000000000..dbfae08d4 --- /dev/null +++ b/physics/smoke_dust/tests/test_seas_gocart.F90 @@ -0,0 +1,110 @@ +!> \file test_seas_gocart.F90 +!! Comprehensive unit test for the GOCART sea salt emission scheme. + +program test_seas_gocart + use machine, only : kind_phys + use seas_mod, only : gocart_seasalt_driver + implicit none + + ! Parameters + integer, parameter :: im = 1, jm = 1, km = 1 + integer, parameter :: num_chem = 30 + integer, parameter :: num_emis_seas = 5 + + ! Global data for tests + real(kind_phys) :: dt, g, pi + real(kind_phys), dimension(im, km, jm, num_chem) :: chem + real(kind_phys), dimension(im, km, jm) :: rho_phy, alt, t_phy, dz8w, u_phy, v_phy + real(kind_phys), dimension(im, km+1, jm) :: p8w + real(kind_phys), dimension(im, jm) :: u10, v10, ustar, tsk, xland, xlat, xlong, area, seashelp + real(kind_phys), dimension(im, 1, jm, num_emis_seas) :: emis_seas + + integer :: total_passed, total_failed + + total_passed = 0 + total_failed = 0 + + dt = 600.0_kind_phys + g = 9.80665_kind_phys + pi = acos(-1.0_kind_phys) + + print *, "Starting GOCART sea salt comprehensive unit test suite..." + + ! Case 1: Standard ocean emission (xland=0.0) (should be non-zero) + call reset_inputs() + xland = 0.0_kind_phys + call run_case_opt("Standard ocean emission (GOCART)", .true., 1) + + ! Case 2: Land (xland=1.0) (should be zero) + call reset_inputs() + xland = 1.0_kind_phys + call run_case_opt("Land (xland=1.0)", .false., 1) + + ! Case 3: High wind speed + call reset_inputs() + xland = 0.0_kind_phys + u10 = 20.0_kind_phys + call run_case_opt("High wind speed", .true., 1) + + ! Case 4: seas_opt = 2 (NGAC scheme) + call reset_inputs() + xland = 0.0_kind_phys + call run_case_opt("NGAC sea salt scheme", .true., 2) + + print *, "---------------------------------------" + print *, "Test Summary: ", total_passed, " PASSED, ", total_failed, " FAILED" + print *, "---------------------------------------" + + if (total_failed > 0) stop 1 + +contains + + subroutine reset_inputs() + chem = 1.0e-12_kind_phys + rho_phy = 1.2_kind_phys + alt = 1.0_kind_phys / rho_phy + t_phy = 300.0_kind_phys + u_phy = 10.0_kind_phys + v_phy = 0.0_kind_phys + dz8w = 50.0_kind_phys + p8w(:,1,:) = 101325.0_kind_phys + p8w(:,2,:) = 100000.0_kind_phys + u10 = 10.0_kind_phys + v10 = 0.0_kind_phys + ustar = 0.5_kind_phys + tsk = 290.0_kind_phys + xland = 0.0_kind_phys ! ocean + xlat = 45.0_kind_phys + xlong = -30.0_kind_phys + area = 1000000.0_kind_phys + emis_seas = 0.0_kind_phys + end subroutine reset_inputs + + subroutine run_case_opt(name, expect_non_zero, opt) + character(len=*), intent(in) :: name + logical, intent(in) :: expect_non_zero + integer, intent(in) :: opt + logical :: is_non_zero + real(kind_phys) :: total_emis + + call gocart_seasalt_driver(dt,alt,t_phy,u_phy, & + v_phy,chem,rho_phy,dz8w,u10,v10,ustar,p8w,tsk, & + xland,xlat,xlong,area,g,emis_seas,pi, & + seashelp,num_emis_seas,num_chem,opt, & + 1, im, 1, jm, 1, km, & + 1, im, 1, jm, 1, km, & + 1, im, 1, jm, 1, km ) + + total_emis = sum(emis_seas(1,1,1,:)) + is_non_zero = total_emis > 1.0e-20_kind_phys + + if (is_non_zero .eqv. expect_non_zero) then + print *, "PASS: ", name, " (Emis: ", total_emis, ")" + total_passed = total_passed + 1 + else + print *, "FAIL: ", name, " (Emis: ", total_emis, ", Expected non-zero: ", expect_non_zero, ")" + total_failed = total_failed + 1 + endif + end subroutine run_case_opt + +end program test_seas_gocart diff --git a/physics/smoke_dust/tests/test_wetdep_ls.F90 b/physics/smoke_dust/tests/test_wetdep_ls.F90 new file mode 100644 index 000000000..e2a0d5fb5 --- /dev/null +++ b/physics/smoke_dust/tests/test_wetdep_ls.F90 @@ -0,0 +1,95 @@ +!> \file test_wetdep_ls.F90 +!! Comprehensive unit test for the large-scale wet deposition scheme. + +program test_wetdep_ls + use machine, only : kind_phys + use module_wetdep_ls, only : wetdep_ls + use rrfs_smoke_config, only : p_smoke, p_dust_1, p_coarse_pm, p_qc + implicit none + + ! Parameters + integer, parameter :: im = 1, jm = 1, km = 5 + integer, parameter :: num_chem = 20 + integer, parameter :: num_moist = 5 + integer, parameter :: ndvel = 1 + + ! Global data for tests + real(kind_phys) :: dt + real(kind_phys), dimension(im, km, jm, num_moist) :: moist + real(kind_phys), dimension(im, km, jm) :: rho, dz8w, vvel + real(kind_phys), dimension(im, km, jm, num_chem) :: var + real(kind_phys), dimension(im, jm) :: rain, wetdpr_smoke, wetdpr_dust, wetdpr_coarsepm + + integer :: total_passed, total_failed + + total_passed = 0 + total_failed = 0 + + dt = 600.0_kind_phys + + print *, "Starting Large-scale Wet Deposition comprehensive unit test suite..." + + ! Case 1: Standard rain (should have wet deposition) + call reset_inputs() + rain = 1.0e-3_kind_phys ! Some rain + moist(:, :, :, p_qc) = 1.0e-4_kind_phys ! some cloud water + call run_case("Standard rain wet deposition", .true.) + + ! Case 2: No rain (should have no wet deposition) + call reset_inputs() + rain = 0.0_kind_phys + moist(:, :, :, p_qc) = 1.0e-4_kind_phys + call run_case("No rain", .false.) + + ! Case 3: No cloud water (should have no wet deposition) + call reset_inputs() + rain = 1.0e-3_kind_phys + moist(:, :, :, p_qc) = 0.0_kind_phys + call run_case("No cloud water", .false.) + + print *, "---------------------------------------" + print *, "Test Summary: ", total_passed, " PASSED, ", total_failed, " FAILED" + print *, "---------------------------------------" + + if (total_failed > 0) stop 1 + +contains + + subroutine reset_inputs() + var = 1.0e-6_kind_phys + moist = 0.0_kind_phys + rho = 1.0_kind_phys + dz8w = 100.0_kind_phys + vvel = 1.0_kind_phys ! vertical velocity + rain = 0.0_kind_phys + wetdpr_smoke = 0.0_kind_phys + wetdpr_dust = 0.0_kind_phys + wetdpr_coarsepm = 0.0_kind_phys + end subroutine reset_inputs + + subroutine run_case(name, expect_non_zero) + character(len=*), intent(in) :: name + logical, intent(in) :: expect_non_zero + logical :: is_non_zero + real(kind_phys) :: total_dep + + call wetdep_ls(dt, var, rain, moist, & + rho, num_chem, num_moist, ndvel, dz8w, vvel,& + wetdpr_smoke, wetdpr_dust, wetdpr_coarsepm, & + 1, im, 1, jm, 1, km, & + 1, im, 1, jm, 1, km, & + 1, im, 1, jm, 1, km ) + + total_dep = sum(wetdpr_smoke) + sum(wetdpr_dust) + sum(wetdpr_coarsepm) + is_non_zero = total_dep > 1.0e-20_kind_phys + + if (is_non_zero .eqv. expect_non_zero) then + print *, "PASS: ", name, " (Dep: ", total_dep, ")" + total_passed = total_passed + 1 + else + print *, "FAIL: ", name, " (Dep: ", total_dep, ", Expected non-zero: ", expect_non_zero, ")" + total_failed = total_failed + 1 + endif + end subroutine run_case + +end program test_wetdep_ls