diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index c6a613d449..b99be288e7 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -1764,6 +1764,7 @@ sub process_namelist_inline_logic { # namelist options for dust emissions ###################################### setup_logic_dust_emis($opts, $nl_flags, $definition, $defaults, $nl, $envxml_ref); + setup_logic_prigent_roughness($opts, $nl_flags, $definition, $defaults, $nl); ################################# # namelist group: megan_emis_nl # @@ -1874,10 +1875,15 @@ sub process_namelist_inline_logic { ######################################### setup_logic_initinterp($opts, $nl_flags, $definition, $defaults, $nl); - ############################### - # namelist group: exice_streams # - ############################### - setup_logic_exice($opts, $nl_flags, $definition, $defaults, $nl); + ################################# + # namelist group: exice_streams # + ################################# + setup_logic_exice($opts, $nl_flags, $definition, $defaults, $nl, $physv); + + ########################################## + # namelist group: clm_temperature_inparm # + ########################################## + setup_logic_coldstart_temp($opts,$nl_flags, $definition, $defaults, $nl); } #------------------------------------------------------------------------------- @@ -2381,7 +2387,6 @@ sub setup_logic_soilstate { add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'organic_frac_squared' ); add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_bedrock', 'use_fates'=>$nl_flags->{'use_fates'}, 'vichydro'=>$nl_flags->{'vichydro'} ); - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_excess_ice'); # excess ice flag should be read before stream vars my $var1 = "soil_layerstruct_predefined"; my $var2 = "soil_layerstruct_userdefined"; @@ -4071,9 +4076,6 @@ sub setup_logic_dust_emis { "$option is being set, need to change one or the other" ); } } - if ( $dust_emis_method eq "Leung_2023" ) { - $log->warning("dust_emis_method is Leung_2023 and that option has NOT been brought into CTSM yet"); - } } # Otherwise make sure dust settings are NOT being set in CLM } else { @@ -4930,13 +4932,31 @@ sub setup_logic_exice { # # excess ice streams # - my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; + my ($opts, $nl_flags, $definition, $defaults, $nl, $physv) = @_; + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_excess_ice', 'phys'=>$physv->as_string()); my $use_exice = $nl->get_value( 'use_excess_ice' ); my $use_exice_streams = $nl->get_value( 'use_excess_ice_streams' ); - # IF excess ice streams is on + my $finidat = $nl->get_value('finidat'); + # If coldstart and use_excess_ice is on: + if ( ( (not defined($use_exice_streams)) && value_is_true($use_exice) ) && string_is_undef_or_empty($finidat) ) { + $nl->set_variable_value('exice_streams', 'use_excess_ice_streams' , '.true.'); + $use_exice_streams = '.true.'; + # if excess ice is turned off + } elsif ( (not defined($use_exice_streams)) && (not value_is_true($use_exice)) ) { + $use_exice_streams = '.false.'; + # Checking for cold clm_start_type and not finidat here since finidat can be not set set in branch/hybrid runs and + # These cases are handled in the restart routines in the model + } elsif ( defined($use_exice_streams) && (not value_is_true($use_exice_streams)) && value_is_true($use_exice) && + ( $nl_flags->{'clm_start_type'} eq "'cold'" || $nl_flags->{'clm_start_type'} eq "'arb_ic'" )) { + $log->fatal_error("use_excess_ice_streams can NOT be FALSE when use_excess_ice is TRUE on the cold start" ); + } + + # Put use_exice_streams into nl_flags so can be referenced later + $nl_flags->{'use_excice_streams'} = $use_exice_streams; + # If excess ice streams is on if (defined($use_exice_streams) && value_is_true($use_exice_streams)) { # Can only be true if excess ice is also on, otherwise fail - if (defined($use_exice) && not value_is_true($use_exice)) { + if ( defined($use_exice) && (not value_is_true($use_exice)) ) { $log->fatal_error("use_excess_ice_streams can NOT be TRUE when use_excess_ice is FALSE" ); } # Otherwise if ice streams are off @@ -4973,6 +4993,46 @@ sub setup_logic_exice { } # end exice streams +sub setup_logic_coldstart_temp { + + my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; + + # set initial temperatures for excess ice gridcells: needs to be set whether excess ice is on or not + + my $use_exice = $nl->get_value( 'use_excess_ice' ); + my $finidat = $nl->get_value('finidat'); + + my @list = ( "excess_ice_coldstart_temp", "excess_ice_coldstart_depth" ); + + # Only needs to be set by the user if it's a coldstart + if ( ! string_is_undef_or_empty($finidat) ) { + foreach my $var ( @list ) { + my $val = $nl->get_value( $var ); + if ( defined($val) ) { + $log->warning("$var only needs to be set if this is a cold-start, although InitCold is always called"); + } + } + } + + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'excess_ice_coldstart_temp', + 'use_excess_ice'=>$use_exice); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'excess_ice_coldstart_depth', + 'use_excess_ice'=>$use_exice); + + my $use_exice_streams = $nl_flags->{'use_excice_streams'}; + my $exice_cs_temp = $nl->get_value( 'excess_ice_coldstart_temp' ); + my $exice_cs_depth = $nl->get_value( 'excess_ice_coldstart_depth' ); + + if (defined($use_exice_streams) && value_is_true($use_exice_streams)) { + if (defined($exice_cs_depth) && $exice_cs_depth <= 0.0 ) { + $log->fatal_error("excess_ice_coldstart_depth is <= 0.0" ); + } + if (defined($exice_cs_temp) && $exice_cs_temp >= 0.0 ) { + $log->fatal_error("excess_ice_coldstart_temp is >= 0.0, no excess ice will be present in this run" ); + } + } +} + #------------------------------------------------------------------------------- sub setup_logic_z0param { @@ -5019,6 +5079,33 @@ sub setup_logic_misc { #------------------------------------------------------------------------------- +sub setup_logic_prigent_roughness { + # + # The Prigent roughness stream data set read in if needed + # + my ($opts, $nl_flags, $definition, $defaults, $nl) = @_; + my $var = "use_prigent_roughness"; + my $dust_emis_method = remove_leading_and_trailing_quotes( $nl->get_value('dust_emis_method') ); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var, + 'dust_emis_method'=>$dust_emis_method ); + my $use_prigent = $nl->get_value($var); + if ( &value_is_true($use_prigent) ) { + if ( $dust_emis_method ne "Leung_2023" ) { + # The Prigent dataset could be used for other purposes + # (such as roughness as in https://github.com/ESCOMP/CTSM/issues/2349) + $log->warning( "$var does NOT need to on without dust_emis_method being Leung_2023, it simply won't be used" ); + } + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldfilename_prigentroughness' ); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_meshfile_prigentroughness' ); + } elsif ( $dust_emis_method eq "Leung_2023" ) { + # In the future we WILL allow it to be turned off for testing and Paleo work + # see: https://github.com/ESCOMP/CTSM/issues/2381) + $log->fatal_error("variable \"$var\" MUST be true when Leung_2023 dust emission method is being used" ); + } +} + +#------------------------------------------------------------------------------- + sub write_output_files { my ($opts, $nl_flags, $defaults, $nl) = @_; @@ -5068,8 +5155,10 @@ sub write_output_files { push @groups, "lifire_inparm"; push @groups, "ch4finundated"; push @groups, "exice_streams"; + push @groups, "clm_temperature_inparm"; push @groups, "soilbgc_decomp"; push @groups, "clm_canopy_inparm"; + push @groups, "prigentroughness"; push @groups, "zendersoilerod"; if (remove_leading_and_trailing_quotes($nl->get_value('snow_cover_fraction_method')) eq 'SwensonLawrence2012') { push @groups, "scf_swenson_lawrence_2012_inparm"; diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 3594fe5c38..c4c056195d 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -2565,6 +2565,17 @@ lnd/clm2/surfdata_esmf/NEON/surfdata_1x1_NEON_TOOL_hist_78pfts_CMIP6_simyr2000_c lnd/clm2/dustemisdata/dust_2x2_ESMFmesh_cdf5_c230730.nc + + + + +.true. +.false. +lnd/clm2/dustemisdata/Prigent_2005_roughness_0.25x0.25_cdf5_c240127.nc +lnd/clm2/dustemisdata/dust_0.25x0.25_ESMFmesh_cdf5_c240222.nc + @@ -2635,7 +2646,10 @@ lnd/clm2/surfdata_esmf/NEON/surfdata_1x1_NEON_TOOL_hist_78pfts_CMIP6_simyr2000_c .false. - +-1.0 +0.5 +-3.15 +0.5 lnd/clm2/paramdata/exice_init_0.125x0.125_c20220516.nc lnd/clm2/paramdata/exice_init_0.125x0.125_ESMFmesh_cdf5_c20220802.nc bilinear diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 9cbcc5e692..4c80f78fd3 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -1694,6 +1694,25 @@ Mapping method from Nitrogen deposition input file to the model resolution copy = copy using the same indices + + + + + +If TRUE use the Prigent roughness dataset + + + +Filename of input stream data for aeolian roughness length (from Prigent's roughness dataset) + + + +mesh filename of input stream data for aeolian roughness length (from Prigent's roughness dataset) + + @@ -1783,6 +1802,7 @@ If true, will ignore the prescribed soilm data for that point and let the model prescribed data. + @@ -2946,6 +2966,18 @@ use case.) If TRUE turn on the excess ice physics, (Lee et al., 2014; Cai et al., 2020) + +Initial soil temperature to use for gridcells with excess ice present during a run starting with coldstart (deg C). Value only applys if use_excess_ice is true. + + + +Soil depth below which initial excess ice concentration will be applied during a run starting with coldstart (m). Value only applys if use_excess_ice is true. +If this is set below depth of the soil depth, only the last soil layer will get excess ice. + + + If TRUE and use_excess_ice is TRUE, use the excess ice stream to determine the initial values of the excess ice field diff --git a/bld/unit_testers/build-namelist_test.pl b/bld/unit_testers/build-namelist_test.pl index 5bccb3b77c..c51b0b0cb3 100755 --- a/bld/unit_testers/build-namelist_test.pl +++ b/bld/unit_testers/build-namelist_test.pl @@ -163,7 +163,7 @@ sub cat_and_create_namelistinfile { # # Figure out number of tests that will run # -my $ntests = 3329; +my $ntests = 3339; if ( defined($opts{'compare'}) ) { $ntests += 1999; @@ -322,7 +322,7 @@ sub cat_and_create_namelistinfile { "-res 0.9x1.25 -namelist '&a irrigate=.true./'", "-res 0.9x1.25 -verbose", "-res 0.9x1.25 -ssp_rcp SSP2-4.5", "-res 0.9x1.25 -test", "-res 0.9x1.25 -sim_year 1850", "-res 0.9x1.25 -namelist '&a use_lai_streams=.true.,use_soil_moisture_streams=.true./'", "-res 0.9x1.25 -namelist '&a use_excess_ice=.true. use_excess_ice_streams=.true./'", - "-res 0.9x1.25 -namelist '&a use_excess_ice=.true. use_excess_ice_streams=.false./'", + "-res 0.9x1.25 --clm_start_type cold -namelist '&a use_excess_ice=.true. use_excess_ice_streams=.true./'", "-res 0.9x1.25 -use_case 1850_control", "-res 1x1pt_US-UMB -clm_usr_name 1x1pt_US-UMB -namelist '&a fsurdat=\"/dev/null\"/'", "-res 1x1_brazil", @@ -517,6 +517,10 @@ sub cat_and_create_namelistinfile { namelst=>"use_crop=.true.", phys=>"clm4_5", }, + "LeungDust_WO_Prigent" =>{ options=>" -envxml_dir . -bgc sp", + namelst=>"use_prigent_roughness=.false.", + phys=>"clm5_1", + }, "soilm_stream off w file" =>{ options=>"-res 0.9x1.25 -envxml_dir .", namelst=>"use_soil_moisture_streams = .false.,stream_fldfilename_soilm='file_provided_when_off'", phys=>"clm5_0", @@ -537,6 +541,18 @@ sub cat_and_create_namelistinfile { namelst=>"use_excess_ice=.true., use_excess_ice_streams = .false.,stream_mapalgo_exice='bilinear'", phys=>"clm5_0", }, + "coldstart exice on wo stream"=>{ options=>"-res 0.9x1.25 -envxml_dir . --clm_start_type cold", + namelst=>"use_excess_ice=.true., use_excess_ice_streams = .false.", + phys=>"clm6_0", + }, + "coldstart exice on bad temp" =>{ options=>"-res 0.9x1.25 -envxml_dir . --clm_start_type cold", + namelst=>"use_excess_ice=.true., use_excess_ice_streams = .true., excess_ice_coldstart_temp=0.0", + phys=>"clm6_0", + }, + "coldstart exice on bad depth" =>{ options=>"-res 0.9x1.25 -envxml_dir . --clm_start_type cold", + namelst=>"use_excess_ice=.true., use_excess_ice_streams = .true., excess_ice_coldstart_depth=0.0", + phys=>"clm6_0", + }, "clm50CNDVwtransient" =>{ options=>" -envxml_dir . -use_case 20thC_transient -dynamic_vegetation -res 10x15 -ignore_warnings", namelst=>"", phys=>"clm5_0", @@ -1220,10 +1236,6 @@ sub cat_and_create_namelistinfile { my %warntest = ( # Warnings without the -ignore_warnings option given - "dustemisLeung" =>{ options=>"-envxml_dir .", - namelst=>"dust_emis_method = 'Leung_2023'", - phys=>"clm5_1", - }, "coldwfinidat" =>{ options=>"-envxml_dir . -clm_start_type cold", namelst=>"finidat = 'testfile.nc'", phys=>"clm5_0", @@ -1256,6 +1268,18 @@ sub cat_and_create_namelistinfile { namelst=>"use_fun=.true.,use_flexiblecn=.false.", phys=>"clm6_0", }, + "Set coldtemp wo coldstart" =>{ options=>"-envxml_dir . --clm_start_type startup", + namelst=>"use_excess_ice=.true.,excess_ice_coldstart_temp=-10.", + phys=>"clm6_0", + }, + "Set colddepth wo coldstart" =>{ options=>"-envxml_dir . --clm_start_type startup", + namelst=>"use_excess_ice=.true.,excess_ice_coldstart_depth=0.5", + phys=>"clm6_0", + }, + "PrigentOnWOLeung" =>{ options=>"-envxml_dir . -bgc sp", + namelst=>"use_prigent_roughness=.true.,dust_emis_method='Zender_2003'", + phys=>"clm6_0", + }, "NotNEONbutNEONlightres" =>{ options=>"--res CLM_USRDAT --clm_usr_name regional --envxml_dir . --bgc bgc --light_res 106x174", namelst=>"fsurdat='build-namelist_test.pl'", phys=>"clm6_0", diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index 951e0d8d5b..cd4c41f23c 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -100,7 +100,7 @@ - + @@ -109,7 +109,7 @@ - + @@ -129,7 +129,7 @@ - + @@ -138,15 +138,6 @@ - - - - - - - - - @@ -156,7 +147,7 @@ - + @@ -247,7 +238,7 @@ - + diff --git a/cime_config/testdefs/testmods_dirs/clm/clm45cam4LndTuningModeZDustSoilErod/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/clm45cam4LndTuningModeZDustSoilErod/user_nl_clm index 93b7ee2e48..d97fbbac57 100644 --- a/cime_config/testdefs/testmods_dirs/clm/clm45cam4LndTuningModeZDustSoilErod/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/clm/clm45cam4LndTuningModeZDustSoilErod/user_nl_clm @@ -1,3 +1,4 @@ ! Turn on using the soil eroditability file in CTSM dust_emis_method = 'Zender_2003' zender_soil_erod_source = 'lnd' +hist_fincl1 += 'FV' diff --git a/cime_config/testdefs/testmods_dirs/clm/clm50cam5LndTuningModeZDustSoilErod/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/clm50cam5LndTuningModeZDustSoilErod/user_nl_clm index 93b7ee2e48..d97fbbac57 100644 --- a/cime_config/testdefs/testmods_dirs/clm/clm50cam5LndTuningModeZDustSoilErod/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/clm/clm50cam5LndTuningModeZDustSoilErod/user_nl_clm @@ -1,3 +1,4 @@ ! Turn on using the soil eroditability file in CTSM dust_emis_method = 'Zender_2003' zender_soil_erod_source = 'lnd' +hist_fincl1 += 'FV' diff --git a/cime_config/testdefs/testmods_dirs/clm/clm51cam6LndTuningModeZDustSoilErod/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/clm51cam6LndTuningModeZDustSoilErod/user_nl_clm index 93b7ee2e48..d97fbbac57 100644 --- a/cime_config/testdefs/testmods_dirs/clm/clm51cam6LndTuningModeZDustSoilErod/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/clm/clm51cam6LndTuningModeZDustSoilErod/user_nl_clm @@ -1,3 +1,4 @@ ! Turn on using the soil eroditability file in CTSM dust_emis_method = 'Zender_2003' zender_soil_erod_source = 'lnd' +hist_fincl1 += 'FV' diff --git a/cime_config/testdefs/testmods_dirs/clm/clm60cam7LndTuningModeLDust/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/clm60cam7LndTuningModeLDust/include_user_mods new file mode 100644 index 0000000000..ef8619d930 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/clm60cam7LndTuningModeLDust/include_user_mods @@ -0,0 +1 @@ +../clm60cam7LndTuningMode diff --git a/cime_config/testdefs/testmods_dirs/clm/clm60cam7LndTuningModeLDust/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/clm60cam7LndTuningModeLDust/user_nl_clm new file mode 100644 index 0000000000..09887c81b8 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/clm60cam7LndTuningModeLDust/user_nl_clm @@ -0,0 +1,5 @@ +dust_emis_method = 'Leung_2023' +! Add all of the Leung optional history fields +hist_fincl1 += 'FV', 'DUST_EMIS_COEFF', 'WND_FRC_FT', 'WND_FRC_FT_DRY', 'WND_FRC_IT', 'WND_FRC_SOIL', 'LND_FRC_MBLE', 'GWC', 'LIQ_FRAC', + 'U_S_MEAN', 'U_S_MEAN', 'ZETAOBU', 'U_FT', 'U_IT', 'ALPHA_TC_RATE', 'P_IT', 'ETA', 'SSR', 'VAI_OKIN', + 'FRC_THR_RGHN_FCT', 'WND_FRC_FT_STD', 'DPFCT_ROCK' diff --git a/doc/ChangeLog b/doc/ChangeLog index 03a3c8999b..c1f5dfaa5a 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,12 +1,11 @@ =============================================================== -Tag name: ctsm5.2.0?? +Tag name: ctsm5.2.020 Originator(s): HuiWangWanderInGitHub (Hui Wang,huiw16@uci.edu) -Date: Thu 01 Aug 2024 05:09:46 PM MDT +Date: Mon 12 Aug 2024 10:51:04 AM MDT One-line Summary: MEGAN updates (MEGAN-CLM6) Purpose and description of changes ---------------------------------- - We add new features to MEGANv2.1 for simulating isoprene emissions based on three recent studies conducted at the BAI lab in the University of California, Irvine. The first one is about the drought impact on isoprene (Wang et al., 2022). The second study investigates the effect of temperature history on the emission factors of boreal broadleaf deciduous shrubs (Wang et al., 2024a). The third study explores a different temperature response curve for C3 Arctic grass (Wang et al., 2024b, under press). These changes improved the model's representation of isoprene emissions during drought and in high-latitude ecosystems. This work relates to issue #1323 and is based on the old MEGANv2.1 framework, but incorporates new scientific insights. @@ -32,15 +31,10 @@ Does this tag change answers significantly for any of the following physics conf Notes of particular relevance for users --------------------------------------- - Changes to documentation: Hui Wang has offered to provide the modifications to the documentation that we need. The PR https://github.com/ESCOMP/CTSM/pull/2588 includes a list of references. -Substantial timing or memory changes: -[e.g., check PFS test in the test suite and look at timings, if you -expect possible significant timing changes] - Testing summary: ---------------- @@ -58,14 +52,241 @@ Answer changes Changes answers relative to baseline: Yes Summarize any changes to answers, i.e., - - what code configurations: All that calculate BVOCs (so not FATES) + - what code configurations: All that calculate BVOCs - what platforms/compilers: All - nature of change: Larger than roundoff, only in BVOC output and specifically limited to MEG_ and _voc variables. Other details ------------- Pull Requests that document the changes (include PR ids): - https://github.com/ESCOMP/CTSM/pull/2588 + https://github.com/ESCOMP/CTSM/pull/2588 New updates on MEGANv2.1 + +=============================================================== +=============================================================== +Tag name: ctsm5.2.019 +Originator(s): erik (Erik Kluzek,UCAR/TSS,303-497-1326) +Date: Sun 11 Aug 2024 09:00:43 AM MDT +One-line Summary: Add in an additional dust emission method Leung_2023, by default off + +Purpose and description of changes +---------------------------------- + +Author: Danny Leung +This is a scheme that builds upon a switch of CESM2's default dust emission scheme (Zender et al., 2003) to +@jfkok's more physical and less empirical one (Kok et al., 2014). This brings in additional modifications and add new aeolian +physics to the Kok's scheme, most notably, by adding the roughness effect (or called drag partition effect) which discounts surface +soil erosion by winds due to the presence of local-scale land-surface roughness elements (mostly plants and rocks). We use a hybrid +approach to account for both roughness from rocks (with a 2-D time-invariant dataset provided by Prigent et al., 2005; 2012) and +roughness from plants (time-varying, as a function of CLM's LAI). We further include the dust emission intermittency effects due to +boundary-layer turbulence. + +Significant changes to scientifically-supported configurations +-------------------------------------------------------------- + +Does this tag change answers significantly for any of the following physics configurations? +(Details of any changes will be given in the "Answer changes" section below.) + + [Put an [X] in the box for any configuration with significant answer changes.] + +[ ] clm6_0 + +[ ] clm5_1 + +[ ] clm5_0 + +[ ] ctsm5_0-nwp + +[ ] clm4_5 + + +Bugs fixed +---------- + +List of CTSM issues fixed (include CTSM Issue # and description) [one per line]: + #1604 -- A new physically based dust emission scheme with more aeolian physics (updated) + +Notes of particular relevance for users +--------------------------------------- + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): New namelist items + dust_emis_method can now be set to Leung_2023 + + use_prigent_roughness ---------------- Logical to use the Prigent dataset (default on when Leung_2023 is used) + stream_fldfilename_prigentroughness -- Prigent streams dataset + stream_meshile_prigentroughness ------ Mesh file for the dataset + + History field FV can now be output for any compset, but it's default off + +Changes to the datasets (e.g., parameter, surface or initial files): + New Prigent streams dataset with aeolian roughness length + +Notes of particular relevance for developers: +--------------------------------------------- + +Caveats for developers (e.g., code that is duplicated that requires double maintenance): + Some issues will be worked on later + + #2568 -- Turn on Leung by default for clm6_0 physics + + #2381 -- Allow Prigent streams to be off in Leung_2023 dust emissions for Paleo work + (Fixing this will allow a unit-test module to be removed) + + #2680 -- prigentroughnessmapalgo is hardcoded to bilinear add it to the namelist so it can change (change name to snake-case) + + #2681 -- Add global dust variables to the DustEmisBase class that apply to both Zender and Leung dust emission methods bfb enhancement next + + #2682 -- Small function and saved constants for convective velocity code cleanup enhancement priority + +Changes to tests or testing: + New test mod that explicitly turns on Leung_2023 dust emissions: clm60cam7LndTuningModeLDust + Some of the tests with clm60cam7LndTuningMode were converted to this test mod + PF unit tester was added for Leung_2023 + More work could be done in the unit-testers to exercise more of the code + +Testing summary: regular +---------------- + + [PASS means all tests PASS; OK means tests PASS other than expected fails.] + + build-namelist tests (if CLMBuildNamelist.pm has changed): + + derecho - PASS (1002 are different) + + python testing (if python code has changed; see instructions in python/README.md; document testing done): + + derecho - PASS + + regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing): + + derecho ----- OK + izumi ------- OK + +Answer changes +-------------- + +Changes answers relative to baseline: no bit-for-bit (unless you turn the new dust scheme on) + Tests with BGC change answers because fieldlists are different (FV is no longer on history files by default) + +Other details +------------- + +Pull Requests that document the changes (include PR ids): +(https://github.com/ESCOMP/ctsm/pull) + #1897 -- A new physically based dust emission scheme with more aeolian physics (updated) + +=============================================================== +=============================================================== +Tag name: ctsm5.2.018 +Originator(s): mvdebolskiy +Date: Fri 02 Aug 2024 09:26:33 AM MDT +One-line Summary: Fix/excess ice cold start + +Purpose and description of changes +---------------------------------- + +Changed the way soil temperature is initialized when excess ice is on and the model starts from cold. + +Specific notes + +TemperatureType now has a private ReadNL subroutine that reads two new namelist options: excess_ice_coldstart_depth and excess_ice_coldstart_temp which control top depth (higher layers get default initial temperature) and soil temperature for soil layers ONLY in columns where excess ice is present. Other columns get their default soil temperature. + +New namelist options belong to a new group "clm_temperature_inparm" with its own logic routine. This is done so in the future hardcoded cold start temperatures can be moved to the namelists. Two namelist tests have been added to check for invalid values logic as well as a test for use_excess_ice_streams logic. use_excess_ice_streams has no default value (due to options for restart) and default value set in the CLMBuildNamelist.pm. + +excessicestream_type has been taken out of waterstate_type and its routines are called directly in clm_inst%Init. Checks for UseExcessIceStreams() are still in place in WaterStateType.F90 for double-checking. clm_inst%Init now has a new local variable with excess ice concentration that are read from streams (or zeros if use_excess_ice is false). Arguments for temperature_type%Init and water_type%Init (and children types) have been changed to include this new variable. Fortran unit tests are also updated to account for these new arguments. + + + + +Significant changes to scientifically-supported configurations +-------------------------------------------------------------- + +Does this tag change answers significantly for any of the following physics configurations? +(Details of any changes will be given in the "Answer changes" section below.) + + [Put an [X] in the box for any configuration with significant answer changes.] + +[ ] clm6_0 + +[ ] clm5_1 + +[ ] clm5_0 + +[ ] ctsm5_0-nwp + +[ ] clm4_5 + + +Bugs fixed +---------- + +List of CTSM issues fixed (include CTSM Issue # and description) [one per line]: + Fixes #2384 -- Cold start temperature init when excess ice is on + Fixes #2373 -- SMS_Lm3_D_Mmpi-serial.1x1_brazil.I2000Clm50FatesCruRsGs.izumi_intel.clm-FatesColdHydro fails + +Notes of particular relevance for users +--------------------------------------- + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): + Added: excess_ice_coldstart_temp and excess_ice_coldstart_depth + +Changes made to namelist defaults (e.g., changed parameter values): + Set new namelist items differently when excess ice is on or off + +Notes of particular relevance for developers: +--------------------------------------------- + +Caveats for developers (e.g., code that is duplicated that requires double maintenance): + Note that the coldstart variables are always used even without excess ice or with an finidat file. + InitCold is always called so the variables are always set. + +Changes to tests or testing: + New tests for build-namelist unit tester + +Testing summary: regular +---------------- + + [PASS means all tests PASS; OK means tests PASS other than expected fails.] + + build-namelist tests (if CLMBuildNamelist.pm has changed): + + derecho - PASS (998 are different from baseline) + + python testing (if python code has changed; see instructions in python/README.md; document testing done): + + derecho - PASS + + regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing): + + derecho ----- OK + izumi ------- OK + + any other testing (give details below): + + SMS_Lm12.f09_f09_mg17.I1850Clm60Sp.derecho_intel.clm-ExcessIceStartup_output_sp_exice + derecho ---- PASS + +Answer changes +-------------- + +Changes answers relative to baseline: + + Summarize any changes to answers, i.e., + - what code configurations: excess ice on, with coldstart or with finidat file interpolation + - what platforms/compilers: all + - nature of change: different at startup + Different for cold-start and can be different for a few points where the cold-start + values are still used on interpoaltion of an existing finidat file + + Tests that compare different to baseline: + ERS_D.f10_f10_mg37.I1850Clm60Sp.izumi_nag.clm-ExcessIceStreams + failed baseline comparison (though answers for default physics configurations have not changed). + +Other details +------------- + +Pull Requests that document the changes (include PR ids): +(https://github.com/ESCOMP/ctsm/pull) + #2465 -- fix excess ice cold starts =============================================================== =============================================================== diff --git a/doc/ChangeSum b/doc/ChangeSum index 315a5cb470..d035fe4861 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,6 +1,8 @@ Tag Who Date Summary ============================================================================================================================ - ctsm5.2.0?? slevis 08/01/2024 MEGAN updates (MEGAN-CLM6) + ctsm5.2.020 slevis 08/12/2024 MEGAN updates (MEGAN-CLM6) + ctsm5.2.019 erik 08/11/2024 Add in an additional dust emission method Leung_2023, by default off + ctsm5.2.018 mvdebols 08/02/2024 Fix/excess ice cold start ctsm5.2.017 erik 07/30/2024 Dust emissions control moved to cmeps ctsm5.2.016 samrabin 07/27/2024 Enable new crop calendars for clm60 compsets ctsm5.2.015 multiple 07/22/2024 Update submodule tags to pass runoff from cism to rof diff --git a/src/biogeochem/CMakeLists.txt b/src/biogeochem/CMakeLists.txt index 0763be096d..270e85838b 100644 --- a/src/biogeochem/CMakeLists.txt +++ b/src/biogeochem/CMakeLists.txt @@ -7,6 +7,7 @@ list(APPEND clm_sources CNSpeciesMod.F90 CNDVType.F90 DustEmisBase.F90 + DustEmisLeung2023.F90 DustEmisZender2003.F90 DustEmisFactory.F90 CropReprPoolsMod.F90 diff --git a/src/biogeochem/DustEmisBase.F90 b/src/biogeochem/DustEmisBase.F90 index 47f2f32688..84467c701d 100644 --- a/src/biogeochem/DustEmisBase.F90 +++ b/src/biogeochem/DustEmisBase.F90 @@ -16,7 +16,8 @@ module DustEmisBase use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) - use clm_varpar , only : dst_src_nbr, ndst, sz_nbr + use clm_varpar , only : dst_src_nbr, ndst, sz_nbr, & + natpft_lb, natpft_ub, natpft_size use clm_varcon , only : grav, spval use landunit_varcon , only : istcrop, istsoil use clm_varctl , only : iulog @@ -68,6 +69,8 @@ module DustEmisBase procedure , public :: GetConstVars ! Get important constant variables procedure , public :: CleanBase ! Base object deallocation (allows extension) procedure , public :: Clean => CleanBase ! Deallocation used by callers + procedure , public :: SetDragPartitionBase ! Base SetDragPartition method that just aborts + procedure , public :: SetDragPartition => SetDragPartitionBase ! SetDrgPartiotion used by callers procedure , private :: InitAllocateBase procedure , private :: InitHistoryBase procedure , private :: InitDustVars ! Initialize variables used in DustEmission method @@ -236,6 +239,23 @@ subroutine InitHistoryBase(this, bounds) end subroutine InitHistoryBase + !------------------------------------------------------------------------ + + subroutine SetDragPartitionBase(this, bounds, drag_partition) + ! + ! !DESCRIPTION: + ! Set the drag partition for testing -- only aborts as only used by the Leung instance + ! + ! !USES: + ! !ARGUMENTS: + class(dust_emis_base_type) :: this + type(bounds_type), intent(in) :: bounds + real(r8), intent(in) :: drag_partition + + call endrun(msg="SetDragPartition is NOT allowed for this dust emission class type") + + end subroutine SetDragPartitionBase + !----------------------------------------------------------------------- subroutine WritePatchToLog(this, p) @@ -271,6 +291,7 @@ subroutine GetPatchVars(this, p, flx_mss_vrt_dst, flx_mss_vrt_dst_tot, vlc_trb, real(r8), optional, intent(out) :: vlc_trb_3 real(r8), optional, intent(out) :: vlc_trb_4 + if ( present(flx_mss_vrt_dst ) ) flx_mss_vrt_dst = this%flx_mss_vrt_dst_patch(p,:) if ( present(flx_mss_vrt_dst_tot) ) flx_mss_vrt_dst_tot = this%flx_mss_vrt_dst_tot_patch(p) if ( present(vlc_trb ) ) vlc_trb = this%vlc_trb_patch(p,:) @@ -787,6 +808,6 @@ subroutine InitDustVars(this, bounds) end subroutine InitDustVars - !------------------------------------------------------------------------ + !============================================================================== end module DustEmisBase diff --git a/src/biogeochem/DustEmisFactory.F90 b/src/biogeochem/DustEmisFactory.F90 index 344596cdb0..cdae8bc20a 100644 --- a/src/biogeochem/DustEmisFactory.F90 +++ b/src/biogeochem/DustEmisFactory.F90 @@ -38,6 +38,7 @@ function create_dust_emissions(bounds, NLFilename) result(dust_emis) !--------------------------------------------------------------------------- use DustEmisBase , only : dust_emis_base_type use DustEmisZender2003, only : dust_emis_zender2003_type + use DustEmisLeung2023 , only : dust_emis_leung2023_type use decompMod , only : bounds_type use shr_kind_mod , only : CL => shr_kind_cl use shr_dust_emis_mod , only : is_dust_emis_zender, is_dust_emis_leung @@ -46,14 +47,13 @@ function create_dust_emissions(bounds, NLFilename) result(dust_emis) class(dust_emis_base_type), allocatable :: dust_emis type(bounds_type), intent(in) :: bounds character(len=*), intent(in) :: NLFilename - ! Local variables if ( is_dust_emis_zender() )then allocate(dust_emis, source=dust_emis_zender2003_type() ) - - ! This will be added when the Leung2023 comes in - !else if ( is_dust_emis_leung() ) - ! allocate(dust_emis, source=dust_emis_zender2003_type() ) + + else if ( is_dust_emis_leung() )then + allocate(dust_emis, source=dust_emis_leung2023_type() ) + else write(iulog,*) 'ERROR: unknown dust_emis_method: ', & errMsg(sourcefile, __LINE__) diff --git a/src/biogeochem/DustEmisLeung2023.F90 b/src/biogeochem/DustEmisLeung2023.F90 new file mode 100644 index 0000000000..869c18c0dc --- /dev/null +++ b/src/biogeochem/DustEmisLeung2023.F90 @@ -0,0 +1,889 @@ +module DustEmisLeung2023 + + !----------------------------------------------------------------------- + ! !DESCRIPTION: + ! Routines in this module calculate Dust mobilization and dry deposition for dust. + ! Simulates dust mobilization due to wind from the surface into the + ! lowest atmospheric layer. On output flx_mss_vrt_dst(ndst) is the surface dust + ! emission (kg/m**2/s) [ + = to atm]. + ! Calculates the turbulent component of dust dry deposition, (the turbulent deposition + ! velocity through the lowest atmospheric layer). CAM will calculate the settling + ! velocity through the whole atmospheric column. The two calculations will determine + ! the dust dry deposition flux to the surface. + ! + ! !USES: + use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_log_mod , only : errMsg => shr_log_errMsg + use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=), shr_infnan_isnan + use clm_varpar , only : dst_src_nbr, ndst + use clm_varcon , only : grav, spval + use landunit_varcon , only : istcrop, istsoil + use clm_varctl , only : iulog + use abortutils , only : endrun + use decompMod , only : bounds_type, subgrid_level_landunit + use atm2lndType , only : atm2lnd_type + use SoilStateType , only : soilstate_type + use CanopyStateType , only : canopystate_type + use WaterStateBulkType , only : waterstatebulk_type + use WaterDiagnosticBulkType, only : waterdiagnosticbulk_type + use FrictionVelocityMod , only : frictionvel_type + use LandunitType , only : lun + use PatchType , only : patch + use DustEmisBase , only : dust_emis_base_type + use pftconMod , only : noveg + use PrigentRoughnessStreamType, only : prigent_roughness_stream_type + ! + ! !PUBLIC TYPES + implicit none + private + ! + ! !PRIVATE DATA: + ! + ! + ! !PUBLIC DATA TYPES: + ! + type, public, extends(dust_emis_base_type) :: dust_emis_leung2023_type + + real(r8), pointer, private :: dst_emiss_coeff_patch (:) ! dust emission coefficient (unitless) + real(r8), pointer, private :: wnd_frc_thr_patch (:) ! wet fluid threshold (m/s) + real(r8), pointer, private :: wnd_frc_thr_dry_patch (:) ! dry fluid threshold (m/s) + real(r8), pointer, private :: wnd_frc_thr_it_patch (:) ! impact threshold (m/s) + real(r8), pointer, private :: lnd_frc_mble_patch (:) ! land mobile fraction + real(r8), pointer, private :: liq_frac_patch (:) ! liquid fraction of total water + real(r8), pointer, private :: wnd_frc_soil_patch (:) ! soil wind friction velocity (m/s) + real(r8), pointer, private :: gwc_patch (:) ! gravimetric water content (kg/kg) + real(r8), pointer, private :: intrmtncy_fct_patch (:) ! intermittency factor, accounting for turbulence shutting down dust emissions (unitless) + real(r8), pointer, private :: stblty_patch (:) ! stability parameter for checking stability condition (stblty < 0 is unstable atmosphere) + real(r8), pointer, private :: u_mean_slt_patch (:) ! wind speed 0.1 m level of dust saltation (m/s) + real(r8), pointer, private :: u_sd_slt_patch (:) ! sd of wind speed 0.1 m level of dust saltation (m/s) + real(r8), pointer, private :: u_fld_thr_patch (:) ! fluid threshold wind speed 0.1 m level of dust saltation (m/s) + real(r8), pointer, private :: u_impct_thr_patch (:) ! impact threshold wind speed at 0.1 m level of dust saltation (m/s) + real(r8), pointer, private :: thr_crs_rate_patch (:) ! threshold crossing rate (unitless) + real(r8), pointer, private :: prb_crs_fld_thr_patch (:) ! probability of wind speed crossing fluid threshold + real(r8), pointer, private :: prb_crs_impct_thr_patch (:) ! probability of wind speed crossing impact threshold + real(r8), pointer, private :: ssr_patch (:) ! [dimless] integrated shear stress ratiio, defined by Okin (2008) and then integrated by Caroline Pierre et al. (2014) + real(r8), pointer, private :: vai_Okin_patch (:) ! [m2 leaf /m2 land] LAI+SAI for calculating Okin drag partition + real(r8), pointer, private :: frc_thr_rghn_fct_patch (:) ! [dimless] hybrid drag partition (or called roughness) factor + real(r8), pointer, private :: wnd_frc_thr_std_patch (:) ! standardized fluid threshold friction velocity (m/s) + type(prigent_roughness_stream_type), private :: prigent_roughness_stream ! Prigent roughness stream data + real(r8), pointer, private :: dpfct_rock_patch (:) ! [fraction] rock drag partition factor, time-constant + + contains + + procedure , public :: Init => InitLeung2023 + procedure , public :: DustEmission ! Dust mobilization + procedure , public :: Clean => CleanLeung2023 + procedure , private :: InitAllocate ! Allocate data + procedure , private :: InitHistory ! History initialization + procedure , private :: InitCold + procedure , private :: CalcDragPartition ! Calculate drag partitioning based on Prigent roughness stream + ! Public for unit testing + procedure , public :: SetDragPartition ! Set drag partitioning for testing + + end type dust_emis_leung2023_type + + interface dust_emis_leung2023_type + ! initialize a new dust emission object + module procedure constructor + end interface dust_emis_leung2023_type + !------------------------------------------------------------------------ + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + +contains + + !----------------------------------------------------------------------- + type(dust_emis_leung2023_type) function constructor() + ! + ! Creates a dust emission object for Leung-2023 type + ! For now this is just a placeholder + !----------------------------------------------------------------------- + + end function constructor + + !------------------------------------------------------------------------ + + subroutine InitLeung2023(this, bounds, NLFilename) + + ! Initialization for this extended class, calling base level initiation and adding to it + class(dust_emis_leung2023_type) :: this + type(bounds_type), intent(in) :: bounds + character(len=*), intent(in) :: NLFilename + + call this%InitBase(bounds, NLFilename) + call this%prigent_roughness_stream%Init( bounds, NLFilename ) + call this%InitAllocate (bounds) + call this%InitHistory (bounds) + call this%InitCold (bounds) + + end subroutine InitLeung2023 + + !------------------------------------------------------------------------ + + subroutine InitAllocate(this, bounds) + ! + ! !ARGUMENTS: + class (dust_emis_leung2023_type) :: this + type(bounds_type), intent(in) :: bounds + ! + ! !LOCAL VARIABLES: + integer :: begp,endp + !------------------------------------------------------------------------ + + begp = bounds%begp ; endp = bounds%endp + allocate(this%dst_emiss_coeff_patch (begp:endp)) ; this%dst_emiss_coeff_patch (:) = nan + allocate(this%wnd_frc_thr_patch (begp:endp)) ; this%wnd_frc_thr_patch (:) = nan + allocate(this%wnd_frc_thr_dry_patch (begp:endp)) ; this%wnd_frc_thr_dry_patch (:) = nan + allocate(this%wnd_frc_thr_it_patch (begp:endp)) ; this%wnd_frc_thr_it_patch (:) = nan + allocate(this%lnd_frc_mble_patch (begp:endp)) ; this%lnd_frc_mble_patch (:) = nan + allocate(this%wnd_frc_soil_patch (begp:endp)) ; this%wnd_frc_soil_patch (:) = nan + allocate(this%gwc_patch (begp:endp)) ; this%gwc_patch (:) = nan + allocate(this%liq_frac_patch (begp:endp)) ; this%liq_frac_patch (:) = nan + allocate(this%intrmtncy_fct_patch (begp:endp)) ; this%intrmtncy_fct_patch (:) = nan + allocate(this%stblty_patch (begp:endp)) ; this%stblty_patch (:) = nan + allocate(this%u_mean_slt_patch (begp:endp)) ; this%u_mean_slt_patch (:) = nan + allocate(this%u_sd_slt_patch (begp:endp)) ; this%u_sd_slt_patch (:) = nan + allocate(this%u_fld_thr_patch (begp:endp)) ; this%u_fld_thr_patch (:) = nan + allocate(this%u_impct_thr_patch (begp:endp)) ; this%u_impct_thr_patch (:) = nan + allocate(this%thr_crs_rate_patch (begp:endp)) ; this%thr_crs_rate_patch (:) = nan + allocate(this%prb_crs_fld_thr_patch (begp:endp)) ; this%prb_crs_fld_thr_patch (:) = nan + allocate(this%prb_crs_impct_thr_patch (begp:endp)) ; this%prb_crs_impct_thr_patch (:) = nan + allocate(this%ssr_patch (begp:endp)) ; this%ssr_patch (:) = nan + allocate(this%vai_Okin_patch (begp:endp)) ; this%vai_Okin_patch (:) = nan + allocate(this%frc_thr_rghn_fct_patch (begp:endp)) ; this%frc_thr_rghn_fct_patch (:) = nan + allocate(this%wnd_frc_thr_std_patch (begp:endp)) ; this%wnd_frc_thr_std_patch (:) = nan + allocate(this%dpfct_rock_patch (begp:endp)) ; this%dpfct_rock_patch (:) = nan + + end subroutine InitAllocate + + !------------------------------------------------------------------------ + + subroutine CleanLeung2023(this) + ! + ! Deallocation for this extended class, calling base level deallocation and adding to it + ! !ARGUMENTS: + class (dust_emis_leung2023_type) :: this + ! + ! !LOCAL VARIABLES: + !------------------------------------------------------------------------ + + call this%CleanBase() + call this%prigent_roughness_stream%Clean( ) + deallocate(this%dst_emiss_coeff_patch ) + deallocate(this%wnd_frc_thr_patch ) + deallocate(this%wnd_frc_thr_dry_patch ) + deallocate(this%wnd_frc_thr_it_patch ) + deallocate(this%lnd_frc_mble_patch ) + deallocate(this%wnd_frc_soil_patch ) + deallocate(this%gwc_patch ) + deallocate(this%liq_frac_patch ) + deallocate(this%intrmtncy_fct_patch ) + deallocate(this%stblty_patch ) + deallocate(this%u_mean_slt_patch ) + deallocate(this%u_sd_slt_patch ) + deallocate(this%u_fld_thr_patch ) + deallocate(this%u_impct_thr_patch ) + deallocate(this%thr_crs_rate_patch ) + deallocate(this%prb_crs_fld_thr_patch ) + deallocate(this%prb_crs_impct_thr_patch ) + deallocate(this%ssr_patch ) + deallocate(this%vai_Okin_patch ) + deallocate(this%frc_thr_rghn_fct_patch ) + deallocate(this%wnd_frc_thr_std_patch ) + deallocate(this%dpfct_rock_patch ) + + end subroutine CleanLeung2023 + + !------------------------------------------------------------------------ + + subroutine InitHistory(this, bounds) + ! + ! !USES: + use histFileMod, only : hist_addfld1d + ! + ! + ! !ARGUMENTS: + class (dust_emis_leung2023_type) :: this + type(bounds_type), intent(in) :: bounds + ! + ! !LOCAL VARIABLES: + integer :: begp,endp + !------------------------------------------------------------------------ + + begp = bounds%begp; endp = bounds%endp + this%dst_emiss_coeff_patch(begp:endp) = spval + call hist_addfld1d (fname='DUST_EMIS_COEFF', units='dimensionless', & + avgflag='A', long_name='soil erodibility or dust emission coefficient for Kok emission scheme', & + ptr_patch=this%dst_emiss_coeff_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%wnd_frc_thr_patch(begp:endp) = spval + call hist_addfld1d (fname='WND_FRC_FT', units='m/s', & + avgflag='A', long_name='fluid threshold friction velocity', & + ptr_patch=this%wnd_frc_thr_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%wnd_frc_thr_dry_patch(begp:endp) = spval + call hist_addfld1d (fname='WND_FRC_FT_DRY', units='m/s', & + avgflag='A', long_name='dry fluid threshold friction velocity', & + ptr_patch=this%wnd_frc_thr_dry_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%wnd_frc_thr_it_patch(begp:endp) = spval + call hist_addfld1d (fname='WND_FRC_IT', units='m/s', & + avgflag='A', long_name='impact threshold friction velocity', & + ptr_patch=this%wnd_frc_thr_it_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%wnd_frc_soil_patch(begp:endp) = spval + call hist_addfld1d (fname='WND_FRC_SOIL', units='m/s', & + avgflag='A', long_name='soil surface wind friction velocity', & + ptr_patch=this%wnd_frc_soil_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%lnd_frc_mble_patch(begp:endp) = spval + call hist_addfld1d (fname='LND_FRC_MBLE', units='dimensionless', & + avgflag='A', long_name='land mobile fraction', & + ptr_patch=this%lnd_frc_mble_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%gwc_patch(begp:endp) = spval + call hist_addfld1d (fname='GWC', units='kg/kg', & + avgflag='A', long_name='gravimetric soil moisture at the topmost soil layer', & + ptr_patch=this%gwc_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%liq_frac_patch(begp:endp) = spval + call hist_addfld1d (fname='LIQ_FRAC', units='dimensionless', & + avgflag='A', long_name='fraction of total water that is liquid', & + ptr_patch=this%liq_frac_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%u_mean_slt_patch(begp:endp) = spval + call hist_addfld1d (fname='U_S_MEAN', units='m/s', & + avgflag='A', long_name='mean wind velocity at saltation level', & + ptr_patch=this%u_mean_slt_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%u_sd_slt_patch(begp:endp) = spval + call hist_addfld1d (fname='U_S_SIGMA', units='m/s', & + avgflag='A', long_name='sd of wind velocity at saltation level', & + ptr_patch=this%u_sd_slt_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%stblty_patch(begp:endp) = spval + call hist_addfld1d (fname='ZETAOBU', units='', & + avgflag='A', long_name='stability parameter', & + ptr_patch=this%stblty_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%u_fld_thr_patch(begp:endp) = spval + call hist_addfld1d (fname='U_FT', units='m/s', & + avgflag='A', long_name='fluid threshold velocity at saltation level', & + ptr_patch=this%u_fld_thr_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%u_impct_thr_patch(begp:endp) = spval + call hist_addfld1d (fname='U_IT', units='m/s', & + avgflag='A', long_name='impact threshold velocity at saltation level', & + ptr_patch=this%u_impct_thr_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%thr_crs_rate_patch(begp:endp) = spval + call hist_addfld1d (fname='ALPHA_TC_RATE', units='', & + avgflag='A', long_name='threshold crossing rate', & + ptr_patch=this%thr_crs_rate_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%prb_crs_fld_thr_patch(begp:endp) = spval + call hist_addfld1d (fname='P_FT', units='', & + avgflag='A', long_name='probability of winds crossing fluid threshold', & + ptr_patch=this%prb_crs_fld_thr_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%prb_crs_impct_thr_patch(begp:endp) = spval + call hist_addfld1d (fname='P_IT', units='', & + avgflag='A', long_name='probability of winds crossing impact threshold', & + ptr_patch=this%prb_crs_impct_thr_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%intrmtncy_fct_patch(begp:endp) = spval + call hist_addfld1d (fname='ETA', units='', & + avgflag='A', long_name='intermittency factor', & + ptr_patch=this%intrmtncy_fct_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%ssr_patch(begp:endp) = spval + call hist_addfld1d (fname='SSR', units='m/s', & + avgflag='A', long_name='Okin-Pierre vegetation shear stress ratio (drag partition factor)', & + ptr_patch=this%ssr_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%vai_Okin_patch(begp:endp) = spval + call hist_addfld1d (fname='VAI_OKIN', units='m/s', & + avgflag='A', long_name='vegetation area index used in the Okin-Pierre plant drag partition scheme', & + ptr_patch=this%vai_Okin_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%frc_thr_rghn_fct_patch(begp:endp) = spval + call hist_addfld1d (fname='FRC_THR_RGHN_FCT', units='dimensionless', & + avgflag='A', long_name='hybrid drag partition (or roughness) factor', & + ptr_patch=this%frc_thr_rghn_fct_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%wnd_frc_thr_std_patch(begp:endp) = spval + call hist_addfld1d (fname='WND_FRC_FT_STD', units='m/s', & + avgflag='A', long_name='standardized fluid threshold friction velocity', & + ptr_patch=this%wnd_frc_thr_std_patch, set_lake=0.0_r8, set_urb=0.0_r8) + this%dpfct_rock_patch(begp:endp) = spval + call hist_addfld1d (fname='DPFCT_ROCK', units='m/s', & + avgflag='A', long_name='rock drag partition factor', & + ptr_patch=this%dpfct_rock_patch) + + end subroutine InitHistory + + !----------------------------------------------------------------------- + + subroutine InitCold(this, bounds) + ! + ! Initialize values from a cold start + ! !ARGUMENTS: + class (dust_emis_leung2023_type) :: this + type(bounds_type), intent(in) :: bounds + ! + ! !LOCAL VARIABLES: + !----------------------------------------------------------------------- + ! Calculate Drag Partition factor (Marticorena and Bergametti 1995 formulation) + if ( this%prigent_roughness_stream%useStreams() )then !if usestreams == true, and it should be always true + call this%CalcDragPartition( bounds ) + else + + call endrun( "ERROR:: dus_emis_Leung_2023 requires requires a streams file of aeolian roughness length to calculate drag partitioning" ) + + end if + + end subroutine InitCold + + !------------------------------------------------------------------------ + + subroutine DustEmission (this, bounds, & + num_nolakep, filter_nolakep, & + atm2lnd_inst, soilstate_inst, canopystate_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, & + frictionvel_inst) + ! + ! !DESCRIPTION: + ! Dust mobilization. This code simulates dust mobilization due to wind + ! from the surface into the lowest atmospheric layer + ! On output flx_mss_vrt_dst(ndst) is the surface dust emission + ! (kg/m**2/s) [ + = to atm] + ! + ! !USES + use shr_const_mod, only : SHR_CONST_RHOFW + use subgridaveMod, only : p2g + ! + ! !ARGUMENTS: + class (dust_emis_leung2023_type) :: this + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_nolakep ! number of column non-lake points in patch filter + integer , intent(in) :: filter_nolakep(num_nolakep) ! patch filter for non-lake points + type(atm2lnd_type) , intent(in) :: atm2lnd_inst + type(soilstate_type) , intent(in) :: soilstate_inst + type(canopystate_type) , intent(in) :: canopystate_inst + type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst + type(waterdiagnosticbulk_type) , intent(in) :: waterdiagnosticbulk_inst + type(frictionvel_type) , intent(in) :: frictionvel_inst + + ! + ! !LOCAL VARIABLES + integer :: fp,p,c,l,g,m,n ! indices + real(r8) :: liqfrac ! fraction of total water that is liquid + real(r8) :: wnd_frc_rat ! [frc] Wind friction threshold over wind friction + real(r8) :: wnd_frc_slt_dlt ! [m s-1] Friction velocity increase from saltatn + real(r8) :: wnd_rfr_dlt ! [m s-1] Reference windspeed excess over threshld + real(r8) :: dst_slt_flx_rat_ttl + real(r8) :: flx_mss_hrz_slt_ttl + real(r8) :: flx_mss_vrt_dst_ttl(bounds%begp:bounds%endp) + real(r8) :: frc_thr_wet_fct + real(r8) :: frc_thr_rgh_fct + real(r8) :: wnd_frc_slt + real(r8) :: lnd_frc_mbl(bounds%begp:bounds%endp) + real(r8) :: bd + real(r8) :: gwc_sfc + real(r8) :: ttlai(bounds%begp:bounds%endp) + real(r8) :: tlai_lu(bounds%begl:bounds%endl) + real(r8) :: sumwt(bounds%begl:bounds%endl) ! sum of weights + logical :: found ! temporary for error check + integer :: index + + real(r8) :: tmp2 ! calculates the dry fluid threshold using Shao and Lu (2000) scheme; + ! replace the tmp1 (Iversen and White, 1982) that was passed from Dustini to DustEmission; + ! tmp2 will be calculated here + real(r8) :: wnd_frc_thr_slt_std ! [m/s] The soil threshold friction speed at standard air density (1.2250 kg/m3) + real(r8) :: frag_expt ! fragmentation exponent + real(r8) :: wnd_frc_thr_slt_it ! [m/s] created for impact threshold friction velocity + real(r8) :: wnd_frc_thr_slt ! [m/s] used for wet fluid threshold friction velocity + real(r8) :: K_length ! [dimless] normalized mean interobstacle distance, or called gap length (Okin, 2008) + real(r8) :: bare_frc ! LUH2 bare soil land cover fraction + real(r8) :: veg_frc ! LUH2 natural vegetation + crop land cover fraction + ! + ! constants + ! + real(r8), parameter :: cst_slt = 2.61_r8 ! [frc] Saltation constant + real(r8), parameter :: flx_mss_fdg_fct = 5.0e-4_r8 ! [frc] Empir. mass flx tuning eflx_lh_vegt + real(r8), parameter :: vai_mbl_thr = 1.0_r8 ! [m2 m-2] new VAI threshold; Danny M. Leung suggests 1, and Zender's scheme uses 0.3 + + real(r8), parameter :: Cd0 = 4.4e-5_r8 ! [dimless] proportionality constant in calculation of dust emission coefficient + real(r8), parameter :: Ca = 2.7_r8 ! [dimless] proportionality constant in scaling of dust emission exponent + real(r8), parameter :: Ce = 2.0_r8 ! [dimless] proportionality constant scaling exponential dependence of dust emission coefficient on standardized soil threshold friction speed + real(r8), parameter :: C_tune = 0.05_r8 ! [dimless] global tuning constant for vertical dust flux; set to produce ~same global dust flux in control sim (I_2000) as old parameterization + real(r8), parameter :: wnd_frc_thr_slt_std_min = 0.16_r8 ! [m/s] minimum standardized soil threshold friction speed + real(r8), parameter :: forc_rho_std = 1.2250_r8 ! [kg/m3] density of air at standard pressure (101325) and temperature (293 K) + real(r8), parameter :: dns_slt = 2650.0_r8 ! [kg m-3] Density of optimal saltation particles + real(r8), parameter :: B_it = 0.82_r8 ! [dimless] ratio = u_star_it / u_star_ft0 + real(r8), parameter :: k = 0.4_r8 ! [dimless] von Karman constant + real(r8), parameter :: f_0 = 0.32_r8 ! [dimless] SSR in the immediate lee of a plant + real(r8), parameter :: c_e = 4.8_r8 ! [dimless] e-folding distance velocity recovery + real(r8), parameter :: D_p = 130e-6_r8 ! [m] Medium soil particle diameter, assuming a global constant of ~130 um following Leung et al. (2023). dmleung 16 Feb 2024 + real(r8), parameter :: gamma_Shao = 1.65e-4_r8 ! [kg s-2] interparticle cohesion: fitting parameter in Shao and Lu (2000) (S&L00). dmleung 16 Feb 2024 + real(r8), parameter :: A_Shao = 0.0123_r8 ! [dimless] coefficient for aerodynamic force: fitting parameter in Shao and Lu (2000). dmleung 16 Feb 2024 + real(r8), parameter :: frag_expt_thr = 5.0_r8 ! [dimless] threshold for fragmentation exponent defined in Leung et al. (2023), somewhere within 3 to 5. It is used to prevent a local AOD blowup (over Patagonia, Argentina), but one can test larger values and relax the threshold if wanted. dmleung 16 Feb 2024 + real(r8), parameter :: z0a_glob = 1e-4_r8 ! [m] assumed globally constant aeolian roughness length value in Leung et al. (2023), for the log law of the wall for Comola et al. (2019) intermittency scheme. dmleung 20 Feb 2024 + real(r8), parameter :: hgt_sal = 0.1_r8 ! [m] saltation height used by Comola et al. (2019) intermittency scheme for the log law of the wall. dmleung 20 Feb 2024 + real(r8), parameter :: vai0_Okin = 0.1_r8 ! [m2/m2] minimum VAI needed for Okin-Pierre's vegetation drag partition equation. lai=0 in the equation will lead to infinity, so a small value is added into this lai dmleung defined. + real(r8), parameter :: zii = 1000.0_r8 ! [m] convective boundary layer height added by dmleung 20 Feb 2024, following other CTSM modules (e.g., CanopyFluxesMod). Should we transfer PBL height (PBLH) from CAM? + real(r8) :: numer ! Numerator term for threshold crossing rate + real(r8) :: denom ! Denominator term for threshold crossing rate + !------------------------------------------------------------------------ + + associate( & + forc_rho => atm2lnd_inst%forc_rho_downscaled_col , & ! Input: [real(r8) (:) ] downscaled density (kg/m**3) + + gwc_thr => soilstate_inst%gwc_thr_col , & ! Input: [real(r8) (:) ] threshold gravimetric soil moisture based on clay content + mss_frc_cly_vld => soilstate_inst%mss_frc_cly_vld_col , & ! Input: [real(r8) (:) ] [frc] Mass fraction clay limited to 0.20 + watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] saturated volumetric soil water + + tlai => canopystate_inst%tlai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index, no burying by snow + tsai => canopystate_inst%tsai_patch , & ! Input: [real(r8) (:) ] one-sided stem area index, no burying by snow + + frac_sno => waterdiagnosticbulk_inst%frac_sno_col, & ! Input: [real(r8) (:) ] fraction of ground covered by snow (0 to 1) + h2osoi_vol => waterstatebulk_inst%h2osoi_vol_col , & ! Input: [real(r8) (:,:) ] volumetric soil water (0<=h2osoi_vol<=watsat) + h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid soil water (kg/m2) + h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] frozen soil water (kg/m2) + + fv => frictionvel_inst%fv_patch , & ! Input: [real(r8) (:) ] friction velocity (m/s) (for dust model) + obu => frictionvel_inst%obu_patch , & ! Input: [real(r8) (:) ] Monin-Obukhov length from the friction Velocity module obu => frictionvel_inst%obu_patch & ! Input: [real(r8) (:) ] Monin-Obukhov length from the friction Velocity module + + dpfct_rock => this%dpfct_rock_patch , & ! Input: rock drag partition factor defined in Marticorena and Bergametti 1995. A fraction between 0 and 1. + + flx_mss_vrt_dst => this%flx_mss_vrt_dst_patch , & ! Output: [real(r8) (:,:) ] surface dust emission (kg/m**2/s) + flx_mss_vrt_dst_tot => this%flx_mss_vrt_dst_tot_patch , & ! Output: [real(r8) (:) ] total dust flux back to atmosphere (pft) + ! below variables are defined in Kok et al. (2014) or (mostly) Leung et al. (2023) dust emission scheme. dmleung 16 Feb 2024 + dst_emiss_coeff => this%dst_emiss_coeff_patch , & ! Output: dust emission coefficient + wnd_frc_thr => this%wnd_frc_thr_patch , & ! Output: fluid threshold + wnd_frc_thr_dry => this%wnd_frc_thr_dry_patch , & ! Output: dry fluid threshold + wnd_frc_thr_it => this%wnd_frc_thr_it_patch , & ! Output: impact threshold + lnd_frc_mble => this%lnd_frc_mble_patch , & ! Output: bare land fraction + wnd_frc_soil => this%wnd_frc_soil_patch , & ! Output: soil friction velocity u_*s = (u_*)(f_eff) + gwc => this%gwc_patch , & ! output: gravimetric water content + liq_frac => this%liq_frac_patch , & ! Output: fraction of liquid moisture + intrmtncy_fct => this%intrmtncy_fct_patch , & ! Output: intermittency factor eta (fraction of time that dust emission is active within a timestep) + stblty => this%stblty_patch , & ! Output: stability in similarity theory (no need to output) + u_mean_slt => this%u_mean_slt_patch , & ! Output: mean wind speed at 0.1 m height translated from friction velocity using the log law of the wall, assuming neutral condition + u_sd_slt => this%u_sd_slt_patch , & ! Output: standard deviation of wind speed from similarity theory + u_fld_thr => this%u_fld_thr_patch , & ! Output: fluid threshold wind speed at 0.1 m height translated from the log law of the wall + u_impct_thr => this%u_impct_thr_patch , & ! Output: impact threshold wind speed at 0.1 m height translated from the log law of the wall + thr_crs_rate => this%thr_crs_rate_patch , & ! Output: threshold crossing rate in Comola 2019 intermittency parameterization + prb_crs_fld_thr => this%prb_crs_fld_thr_patch , & ! Output: probability of instantaneous wind crossing fluid threshold in Comola 2019 intermittency parameterization + prb_crs_impct_thr => this%prb_crs_impct_thr_patch , & ! Output: probability of instantaneous wind crossing impact threshold in Comola 2019 intermittency parameterization + ssr => this%ssr_patch , & ! Output: vegetation drag partition factor in Okin 2008 vegetation roughness effect (called shear stress ratio, SSR in Okin 2008) + vai_Okin => this%vai_Okin_patch , & ! Output: vegetation area index for calculating Okin-Pierre vegetation drag partitioning. vai=0 in the ssr equation will lead to infinity, so a small value is added into this vai dmleung defined. (no need to output) 16 Feb 2024 + frc_thr_rghn_fct => this%frc_thr_rghn_fct_patch , & ! Output: hybrid/total drag partition factor considering both rock and vegetation drag partition factors. + wnd_frc_thr_std => this%wnd_frc_thr_std_patch & ! Output: standardized dust emission threshold friction velocity defined in Jasper Kok et al. (2014). + ) + + ttlai(bounds%begp : bounds%endp) = 0.0_r8 + ! make lai average at landunit level + do fp = 1,num_nolakep + p = filter_nolakep(fp) + ttlai(p) = tlai(p)+tsai(p) + enddo + + tlai_lu(bounds%begl : bounds%endl) = spval + sumwt(bounds%begl : bounds%endl) = 0.0_r8 + do p = bounds%begp,bounds%endp + if (ttlai(p) /= spval .and. patch%active(p) .and. patch%wtlunit(p) /= 0.0_r8) then + c = patch%column(p) + l = patch%landunit(p) + if (sumwt(l) == 0.0_r8) tlai_lu(l) = 0.0_r8 + tlai_lu(l) = tlai_lu(l) + ttlai(p) * patch%wtlunit(p) + sumwt(l) = sumwt(l) + patch%wtlunit(p) + end if + end do + found = .false. + do l = bounds%begl,bounds%endl + if (sumwt(l) > 1.0_r8 + 1.e-6_r8) then + found = .true. + index = l + exit + else if (sumwt(l) /= 0.0_r8) then + tlai_lu(l) = tlai_lu(l)/sumwt(l) + end if + end do + if (found) then + write(iulog,*) 'error: sumwt is greater than 1.0 at l= ',index + call endrun(subgrid_index=index, subgrid_level=subgrid_level_landunit, msg=errMsg(sourcefile, __LINE__)) + return + end if + + ! Loop through patches + + ! initialize variables which get passed to the atmosphere + flx_mss_vrt_dst(bounds%begp:bounds%endp,:)=0.0_r8 + + do fp = 1,num_nolakep + p = filter_nolakep(fp) + c = patch%column(p) + l = patch%landunit(p) + + ! the following code from subr. lnd_frc_mbl_get was adapted for lsm use + ! purpose: return fraction of each gridcell suitable for dust mobilization + + ! the "bare ground" fraction of the current sub-gridscale cell decreases + ! linearly from 1 to 0 as VAI(=tlai+tsai) increases from 0 to vai_mbl_thr + ! if ice sheet, wetland, or lake, no dust allowed + + if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then + if (tlai_lu(l) < vai_mbl_thr) then + lnd_frc_mbl(p) = 1.0_r8 - (tlai_lu(l))/vai_mbl_thr + else + lnd_frc_mbl(p) = 0.0_r8 + endif + lnd_frc_mbl(p) = lnd_frc_mbl(p) * (1.0_r8 - frac_sno(c)) + else + lnd_frc_mbl(p) = 0.0_r8 + end if + end do + + do fp = 1,num_nolakep + p = filter_nolakep(fp) + if (lnd_frc_mbl(p)>1.0_r8 .or. lnd_frc_mbl(p)<0.0_r8) then + write(iulog,*)'Error dstmbl: pft= ',p,' lnd_frc_mbl(p)= ',lnd_frc_mbl(p), & + errMsg(sourcefile, __LINE__) + call endrun("Bad value for dust mobilization fraction") + return + end if + end do + + ! dmleung add output for bare_frc and veg_frc here if wanted !!!!---------------------- + + ! reset history output variables before next if-statement to avoid output = inf + + do fp = 1,num_nolakep + p = filter_nolakep(fp) + flx_mss_vrt_dst_tot(p) = 0.0_r8 + dst_emiss_coeff(p) = 0.0_r8 + wnd_frc_thr(p) = 0.0_r8 + wnd_frc_thr_dry(p) = 0.0_r8 + lnd_frc_mble(p) = 0.0_r8 + wnd_frc_soil(p) = 0.0_r8 + gwc(p) = 0.0_r8 + liq_frac(p) = 0.0_r8 + u_mean_slt(p) = 0.0_r8 + u_sd_slt(p) = 0.0_r8 + stblty(p) = 0.0_r8 + u_fld_thr(p) = 0.0_r8 + u_impct_thr(p) = 0.0_r8 + thr_crs_rate(p) = 0.0_r8 + prb_crs_fld_thr(p) = 0.0_r8 + prb_crs_impct_thr(p) = 0.0_r8 + intrmtncy_fct(p) = 0.0_r8 + ssr(p) = 0.0_r8 + vai_Okin(p) = 0.0_r8 + frc_thr_rghn_fct(p) = 0.0_r8 + wnd_frc_thr_std(p) = 0.0_r8 + end do + do n = 1, ndst + do fp = 1,num_nolakep + p = filter_nolakep(fp) + flx_mss_vrt_dst(p,n) = 0.0_r8 + end do + end do + + do fp = 1,num_nolakep + p = filter_nolakep(fp) + c = patch%column(p) + l = patch%landunit(p) + g = patch%gridcell(p) + + !-------------------------------------------------------------------------------------------------- + ! put dust emission calculation here to output threshold friction velocity for the whole globe, + ! not just when lnd_frc_mbl = 0. Danny M. Leung 27 Nov 2021 + + !#################################################################################################### + ! calculate soil moisture effect for dust emission threshold + ! following Fecan, Marticorena et al. (1999) + ! also see Zender et al. (2003) for DUST emission scheme and Kok et al. (2014b) for K14 emission scheme in CESM + bd = (1.0_r8-watsat(c,1))*dns_slt ![kg m-3] Bulk density of dry surface soil (dmleung changed from 2700 to dns_slt, soil particle density, on 16 Feb 2024. Note that dn s_slt=2650 kg m-3 so the value is changed by a tiny bit from 2700 to 2650. dns_slt has been here for many years so dns_slt should be used here instead of explicitly typing the value out. dmleung 16 Feb 2024) + + ! use emission threshold to calculate standardized threshold and dust emission coefficient + + ! Here convert h2osoi_vol (H2OSOI) at the topmost CTSM soil layer from volumetric (m3 water / m3 soil) to gravimetric soil moisture (kg water / kg soil) + gwc_sfc = h2osoi_vol(c,1)*SHR_CONST_RHOFW/bd ![kg kg-1] Gravimetric H2O cont + if (gwc_sfc > gwc_thr(c)) then + frc_thr_wet_fct = sqrt(1.0_r8 + 1.21_r8 * (100.0_r8*(gwc_sfc - gwc_thr(c)))**0.68_r8) ! dmleung's comment: this is an empirical equation by Fecan, Marticorena et al. (1999) on relating the soil moisture factor on enhancing dust emission threshold to gravimetric soil moisture. 1.21 and 0.68 are fitting parameters in the regression done by Fecan; 100 is to convert gracimetric soil moisture from fraction (kg water / kg soil) to percentage. Note that gwc_thr was defined in SoilStateInitConst.F90 as a fraction. 1.0_r8 means there is no soil moisture effect on enhancing dust emission threhsold. dmleung 16 Feb 2024. + else + frc_thr_wet_fct = 1.0_r8 + end if + + ! output moisture variables + gwc(p) = gwc_sfc ! output surface gravimetric water content + + ! slevis: adding liqfrac here, because related to effects from soil water + liqfrac = max( 0.0_r8, min( 1.0_r8, h2osoi_liq(c,1) / (h2osoi_ice(c,1)+h2osoi_liq(c,1)+1.0e-6_r8) ) ) + ! dmleung: output liquid fraction + liq_frac(p) = liqfrac + + !####################################################################################################### + ! calculate Shao & Lu (2000) dust emission threshold scheme here + ! use tmp1 from DUSTini for Iversen and White I&W (1982) (~75 um is optimal); use tmp2 for S&L (2000) (~80 um is optimal) + ! see Danny M. Leung et al. (2023) + !####################################################################################################### + tmp2 = sqrt(A_Shao * (dns_slt*grav*D_p + gamma_Shao/D_p)) ! calculate S&L (2000) scheme here for threshold; gamma = 1.65e-4 following S&L00, D_p = 127 um ~ 130 um following Leung et al. (2023). dmleung use defined parameters instead of typing numerical values 16 Feb 2024 + wnd_frc_thr_dry(p) = tmp2 / sqrt(forc_rho(c)) ! dry fluid threshold + wnd_frc_thr_slt = tmp2 / sqrt(forc_rho(c)) * frc_thr_wet_fct !* frc_thr_rgh_fct ! fluid threshold. dmleung commented out frc_thr_rgh_fct since it is used to modify the wind, not the wind threshold. + wnd_frc_thr_slt_it = B_it * tmp2 / sqrt(forc_rho(c)) ! define impact threshold + + ! the above formula is true for Iversen and White (1982) and Shao and Lu (2000) scheme + wnd_frc_thr(p) = wnd_frc_thr_slt ! output fluid threshold + wnd_frc_thr_it(p) = wnd_frc_thr_slt_it ! output impact threshold + + !############################################################################################## + ! dmleung: here, calculate quantities relevant to the fluid threshold + ! standardized fluid threshold + wnd_frc_thr_slt_std = wnd_frc_thr_slt * sqrt(forc_rho(c) / forc_rho_std) ! standardized soil threshold friction speed (defined using fluid threshold) + wnd_frc_thr_std(p) = wnd_frc_thr_slt_std ! output standardized fluid threshold + ! dust emission coefficient or soil erodibility coefficient (this is analogous to the soil erodibility map or prefenertial source filter in Zender; see zendersoilerodstream) + dst_emiss_coeff(p) = Cd0 * exp(-Ce * (wnd_frc_thr_slt_std - wnd_frc_thr_slt_std_min) / wnd_frc_thr_slt_std_min) ! save dust emission coefficient here for all grids + + ! framentation exponent (dependent on fluid threshold) + frag_expt = (Ca * (wnd_frc_thr_slt_std - wnd_frc_thr_slt_std_min) / wnd_frc_thr_slt_std_min) ! fragmentation exponent, defined in Kok et al. (2014a) + if (frag_expt > frag_expt_thr) then ! set fragmentation exponent to be 3 or 5 at maximum, to avoid local AOD blowup + frag_expt = frag_expt_thr + end if + + !############################################################################################## + !################ drag partition effect, and soil-surface friction velocity ################### + ! subsection on computing vegetation drag partition and hybrid drag partition factors + ! in Leung et al. (2023), drag partition effect is applied on the wind instead of the threshold + !############################################################################################## + ! the following comes from subr. frc_thr_rgh_fct_get + ! purpose: compute factor by which surface roughness increases threshold + ! friction velocity (currently a constant) + + if (lnd_frc_mbl(p) > 0.0_r8 .AND. tlai_lu(l)<= vai_mbl_thr) then + + vai_Okin(p) = tlai_lu(l)+vai0_Okin ! LAI+SAI averaged to landunit level; the equation is undefined at lai=0, and LAI in CTSM has some zeros over deserts, so we add in a small number. + if (vai_Okin(p) > 1.0_r8) then + vai_Okin(p) = 1.0_r8 ! setting LAI = 1 to be a max value (since K_length goes to negative when LAI>1) + end if + + + ! calculate Okin's shear stress ratio (SSR, which is vegetation drag partition factor) using Pierre's equation + K_length = 2.0_r8 * (1.0_r8/vai_Okin(p) - 1.0_r8) ! Here LAI has to be non-zero to avoid blowup, and < 1 to avoid -ve K_length. See this equation in Leung et al. (2023). This line is Okin's formulation + ssr(p) = (K_length+f_0*c_e)/(K_length+c_e) ! see this equation in Caroline Pierre et al. (2014) or Leung et al. (2023). This line is Pierre's formulation. + + ! calculation of the hybrid/total drag partition effect considering both rock and vegetation drag partitioning using LUH2 bare and veg fractions within a grid + if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then + if (patch%itype(p) == noveg) then ! if bare, uses rock drag partition factor + if (shr_infnan_isnan(dpfct_rock(p)) ) then ! dmleung added 24 May 2024: dpfct_rock(p) could be NaN; CLM could run when DEBUG=FALSE in env_build.xml but dies when DEBUG=TRUE (usually when checking if wnd_frc_slt > wnd_frc_thr_slt_it and if numer/denom < 30 below) + frc_thr_rgh_fct = 0.001_r8 ! Set drag partition effect to be a very small value (or zero) such that there is no emission whenever dpfct_rock(p) = NaN; dmleung 24 May 2024 + else + frc_thr_rgh_fct = dpfct_rock(p) + end if + else ! if vegetation, uses vegetation drag partition factor + frc_thr_rgh_fct = ssr(p) + end if + else + frc_thr_rgh_fct = 1.0_r8 + end if + + wnd_frc_slt = fv(p) * frc_thr_rgh_fct ! wnd_frc_slt is the drag-parition-modified wind speed and will be used in the dust emission equation below + + frc_thr_rghn_fct(p) = frc_thr_rgh_fct ! save and output hybrid drag partition factor + + else ! for lnd_frc_mbl=0, do not change friction velocity and assume drag partition factor = 0 + wnd_frc_slt = fv(p) * frc_thr_rghn_fct(p) ! The value here is not important since once lnd_frc_mbl(p) <= 0.0_r8 there will be no emission. dmleung added 5 Jul 2024 + frc_thr_rghn_fct(p) = 0.0_r8 ! When LAI > vai_mbl_thr, the drag partition effect is zero. dmleung 16 Feb 2024. + end if + + !########## end of drag partition effect ####################################################### + + ! save soil friction velocity and roughness effect before the if-statement + wnd_frc_soil(p) = wnd_frc_slt ! save soil friction velocity for CLM output, which has drag partition and Owen effect + ! 20 Feb 2024: dmleung notes that Leung does not consider the Owen's effect. This is Jasper Kok's decision. The Owen's effect should be in Zender's DUST emission scheme. + + ! save land mobile fraction + lnd_frc_mble(p) = lnd_frc_mbl(p) ! save land mobile fraction first, before the if-statement + ! only perform the following calculations if lnd_frc_mbl is non-zero + + if (lnd_frc_mbl(p) > 0.0_r8) then ! if bare land fraction is larger than 0 then calculate the dust emission equation + + ! reset these variables which will be updated in the following if-block + + flx_mss_hrz_slt_ttl = 0.0_r8 + flx_mss_vrt_dst_ttl(p) = 0.0_r8 + + ! the following comes from subr. flx_mss_hrz_slt_ttl_Whi79_get + ! purpose: compute vertically integrated streamwise mass flux of particles + + if (wnd_frc_slt > wnd_frc_thr_slt_it) then! if using Leung's scheme, use impact threshold for dust emission equation; if Zender, uses fluid threshold (wnd_frc_thr_slt) for dust emission equation + + !################### for Leung et al. (2023) ################################################ + ! dmleung: instead of using mss_frc_cly_vld(c) with a range of [0,0.2] , which makes dust too sensitive to input clay surface dataset), for now use 0.1 + mss_frc_cly_vld(c) * 0.1 / 0.20 with a range of [0.1,0.2]. So, instead of scaling dust emission to 1/20 times for El Djouf (western Sahara) because of its 1 % clay fraction, scale its emission to 1/2 times. This reduces the sensitivity of dust emission to the clay input dataset. In particular, because dust emission is a small-scale process and the grid-averaged clay from the new WISE surface dataset over El Djouf is 1 %, much lower than the former FAO estimation of 5 %, dmleung is not sure if the clay value of 1 % suppresses too much of El Djouf's small-scale dust emission process. Similar to Charlie Zender's feeling suspicious about the soil texture datasets (or the dust emission schemes' sandblasting process), dmleung feels the same and for now decides to still allow dust emission to weakly scale with clay fraction, however limiting the scaling factor to 0.1-0.2. See modification in SoilStateInitTimeConst.F90. dmleung 5 Jul 2024 + flx_mss_vrt_dst_ttl(p) = dst_emiss_coeff(p) * mss_frc_cly_vld(c) * forc_rho(c) * ((wnd_frc_slt**2.0_r8 - wnd_frc_thr_slt_it**2.0_r8) / wnd_frc_thr_slt_it) * (wnd_frc_slt / wnd_frc_thr_slt_it)**frag_expt ! Leung et al. (2022) uses Kok et al. (2014) dust emission euqation for emission flux + + ! account for bare soil fraction, frozen soil fraction, and apply global tuning parameter (Kok et al. 2014) + flx_mss_vrt_dst_ttl(p) = flx_mss_vrt_dst_ttl(p) * lnd_frc_mbl(p) * C_tune * liqfrac + !######################################################################################## + end if + + !############## Danny M. Leung added the intermittency calculation ################################# + ! subsection for intermittency factor calculation (only used by Leung's scheme, not Zender's scheme) + ! Leung et al. (2023) uses the Comola et al. (2019) intermittency scheme for the calculation of intermittent dust emissions. + ! This part takes care of the sub-timestep, high-frequency (< 1 minute period) turblent wind fluctuations occuring at the planetary boundary layer (PBL) near surface. Subtimestep wind gusts and episodes are important for generating emissions in marginal dust source regions, such as semiarid areas and high-latitude polar deserts. + ! 2 Dec 2021: assume no buoyancy contribution to the wind fluctuation (u_sd_slt), so no obul(p) is needed. It is shown to be important for the wind fluctuations contribute little to the intermittency factor. We might add this back in the future revisions. + ! 20 Feb 2024: dmleung notes that dmleung may revise Comola's scheme in the future to improve Comola's formulation of the statistical parameterization. + + ! mean lowpass-filtered wind speed at hgt_sal = 0.1 m saltation height (assuming aerodynamic roughness length z0a_glob = 1e-4 m globally for ease; also assuming neutral condition) + u_mean_slt(p) = (wnd_frc_slt/k) * log(hgt_sal / z0a_glob) ! translating from ustar (velocity scale) to actual wind + + if ( obu(p) == 0.0_r8 )then + call endrun(msg='Input obu is zero and can NOT be' ) + return + end if + stblty(p) = zii / obu(p) ! -dmleung 20 Feb 2024: use obu from CTSM and PBL height = zii (= 1000_r8) which is default in CTSM. Should we transfer PBL height from CAM? + if ((12.0_r8 - 0.5_r8 * stblty(p)) .GE. 0.001_r8) then ! should have used 0 theoretically; used 0.001 here to avoid undefined values + u_sd_slt(p) = wnd_frc_slt * (12.0_r8 - 0.5_r8 * stblty(p))**0.333_r8 + else + u_sd_slt(p) = wnd_frc_slt * (0.001_r8)**0.333_r8 ! should have used 0 theoretically; used 0.001 here to avoid undefined values + end if + + ! threshold velocities + ! Here wnd_frc_thr_slt is the fluid threshold; wnd_frc_thr_dry(p) is the dry fluid threshold; B_it*wnd_frc_thr_dry(p) is the impact threshold + ! fluid threshold wind at 0.1 m saltation height + u_fld_thr(p) = (wnd_frc_thr_slt/k) * log(hgt_sal / z0a_glob) ! assume a globally constant z0a value for the log law of the wall, but it can be z0m from CLM or, better, z0a from Prigent's roughness dataset. Danny M. Leung et al. (2023) chose to assume a global constant z0a = 1e-4 m. dmleung 20 Feb 2024 + ! impact threshold wind at 0.1 m saltation height + u_impct_thr(p) = (wnd_frc_thr_slt_it/k) * log(hgt_sal / z0a_glob) + + ! threshold crossing rate + numer = (u_fld_thr(p)**2.0_r8 - u_impct_thr(p)**2.0_r8 - 2.0_r8 * u_mean_slt(p) * (u_fld_thr(p) - u_impct_thr(p))) + denom = (2.0_r8 * u_sd_slt(p)**2.0_r8) ! note that u_sd_slt should be always positive + ! Truncate to zero if the expression inside exp is becoming too large + if ( numer/denom < 30.0_r8 ) then ! set numer/denom to be < 30 given exp(30) below is already very large; also denom itself should be non-zero and non-negative given the standard deviation (u_sd_slt) of the subtimestep wind fluctuation is non-negative. dmleung 28 May 2024 + thr_crs_rate(p) = (exp((u_fld_thr(p)**2.0_r8 - u_impct_thr(p)**2.0_r8 - 2.0_r8 * u_mean_slt(p) * (u_fld_thr(p) - u_impct_thr(p))) / (2.0_r8 * u_sd_slt(p)**2.0_r8)) + 1.0_r8)**(-1.0_r8) + else + thr_crs_rate(p) = 0.0_r8 + end if + + ! probability that lowpass-filtered wind speed does not exceed u_ft + prb_crs_fld_thr(p) = 0.5_r8 * (1.0_r8 + erf((u_fld_thr(p) - u_mean_slt(p)) / ( sqrt(2.0_r8) * u_sd_slt(p)))) + ! probability that lowpass-filtered wind speed does not exceed u_it + prb_crs_impct_thr(p) = 0.5_r8 * (1.0_r8 + erf((u_impct_thr(p) - u_mean_slt(p)) / ( sqrt(2.0_r8) * u_sd_slt(p)))) + + ! intermittency factor (eta; ranging from 0 to 1) + intrmtncy_fct(p) = 1.0_r8 - prb_crs_fld_thr(p) + thr_crs_rate(p) * (prb_crs_fld_thr(p) - prb_crs_impct_thr(p)) + + ! multiply dust emission flux by intermittency factor + if ( shr_infnan_isnan(intrmtncy_fct(p)) ) then ! if intrmtncy_fct(p) is not NaN then multiply by intermittency factor; this statement is needed because dust emission flx_mss_vrt_dst_ttl(p) has to be non NaN (at least zero) to be outputted + flx_mss_vrt_dst_ttl(p) = flx_mss_vrt_dst_ttl(p) + else + flx_mss_vrt_dst_ttl(p) = flx_mss_vrt_dst_ttl(p) * intrmtncy_fct(p) ! multiply dust flux by intermittency + end if + + !############ end the intermittency subsection here; only use for Leung's scheme ########################## + + end if ! lnd_frc_mbl > 0.0 + + end do + + ! the following comes from subr. flx_mss_vrt_dst_prt in C. Zender's code + ! purpose: partition total vertical mass flux of dust into transport bins + + do n = 1, ndst + do m = 1, dst_src_nbr + do fp = 1,num_nolakep + p = filter_nolakep(fp) + if (lnd_frc_mbl(p) > 0.0_r8) then + flx_mss_vrt_dst(p,n) = flx_mss_vrt_dst(p,n) + this%ovr_src_snk_mss(m,n) * flx_mss_vrt_dst_ttl(p) + end if + end do + end do + end do + + do n = 1, ndst + do fp = 1,num_nolakep + p = filter_nolakep(fp) + if (lnd_frc_mbl(p) > 0.0_r8) then + flx_mss_vrt_dst_tot(p) = flx_mss_vrt_dst_tot(p) + flx_mss_vrt_dst(p,n) + end if + end do + end do + + end associate + + end subroutine DustEmission + + !------------------------------------------------------------------------ + subroutine CalcDragPartition(this, bounds) + ! + ! !DESCRIPTION: + ! Commented below by Danny M. Leung 31 Dec 2022 + ! Calculate the drag partition effect of friction velocity due to surface roughness following + ! Leung et al. (2023). This module is used in the dust emission module DUSTMod.F90 for + ! calculating drag partitioning. The drag partition equation comes from Marticorena and + ! Bergametti (1995) with constants modified by Darmenova et al. (2009). Here it is assumed + ! that this equation is used only over arid/desertic regions, such that Catherine Prigent's + ! roughness measurements represents mostly rocks. For more vegetated areas, the vegetation + ! roughness and drag partitioning are calculated in the DustEmission subroutine. This + ! subroutine is used in the InitCold subroutine of DUSTMod.F90. + ! + ! !USES: + use PatchType , only : patch + use landunit_varcon , only : istdlak + use LandunitType , only : lun + ! + ! !ARGUMENTS: + implicit none + class (dust_emis_leung2023_type) :: this + type(bounds_type), intent(in) :: bounds + ! + ! !LOCAL VARIABLES: + integer :: g, p, fp, l ! Indices + real(r8) :: z0s ! smooth roughness length (m) + + ! constants + real(r8), parameter :: D_p = 130e-6_r8 ! [m] Medium soil particle diameter, assuming a global constant + ! of ~130 um following Leung et al. (2023) + real(r8), parameter :: X = 10_r8 ! [m] distance downwind of the roughness element (rock). Assume + ! estimating roughness effect at a distance of 10 m following Leung et al. (2023) + real(r8), parameter :: b1 = 0.7_r8 ! [dimless] first fitting coefficient for the drag partition equation by Marticorena and Bergametti (1995), later modified by Darmenova et al. (2009). + real(r8), parameter :: b2 = 0.8_r8 ! [dimless] second fitting coefficient for the drag partition equation by Marticorena and Bergametti (1995), later modified by Darmenova et al. (2009). + !--------------------------------------------------------------------- + + ! Make sure we've initialized the Prigent roughness streams + if ( .not. this%prigent_roughness_stream%IsStreamInit() )then + write(iulog,*)'Error : Prigent roughness stream is NOT on: ', errMsg(sourcefile, __LINE__) + call endrun(msg=' ERROR: Streams have not been initialized, make sure Init is called first' & + //', and streams are on') + return + end if + + ! dmleung: this loop calculates the drag partition effect (or roughness effect) of rocks. + ! We save the drag partition factor as a patch level quantity. + ! TODO: EBK 02/13/2024: Several magic numbers here that should become parameters so the meaning is preserved + z0s = 2.0_r8 * D_p / 30.0_r8 ! equation for smooth roughness length for soil grain. See Danny M. Leung et al. (2023) and Martina Klose et al. (2021) for instance. 1/15 is a coefficient that relates roughness to soil particle diameter D_p. + ! Here we assume soil medium size is a global constant, and so is smooth roughness length. + do p = bounds%begp,bounds%endp + g = patch%gridcell(p) + l = patch%landunit(p) + if (lun%itype(l) /= istdlak) then + ! Calculating rock drag partition factor using the Marticorena and Bergametti (1995) formulation. + ! 0.01 is used to convert Prigent's roughness length dataset from centimeter to meter. + this%dpfct_rock_patch(p) = 1.0_r8 - ( log(this%prigent_roughness_stream%prigent_rghn(g)*0.01_r8/z0s) & + / log(b1 * (X/z0s)**b2 ) ) + end if + end do + + end subroutine CalcDragPartition + + + !------------------------------------------------------------------------ + + subroutine SetDragPartition(this, bounds, drag_partition) + ! + ! !DESCRIPTION: + ! Set the drag partition for testing + ! + ! !USES: + ! + ! !ARGUMENTS: + implicit none + class (dust_emis_leung2023_type) :: this + type(bounds_type), intent(in) :: bounds + real(r8), intent(in) :: drag_partition + ! + ! !LOCAL VARIABLES: + integer :: p ! Indices + + !--------------------------------------------------------------------- + + do p = bounds%begp,bounds%endp + this%dpfct_rock_patch(p) = drag_partition + end do + + end subroutine SetDragPartition + + !============================================================================== + +end module DustEmisLeung2023 \ No newline at end of file diff --git a/src/biogeochem/test/DustEmis_test/CMakeLists.txt b/src/biogeochem/test/DustEmis_test/CMakeLists.txt index 9239d16e22..c705d6c2e3 100644 --- a/src/biogeochem/test/DustEmis_test/CMakeLists.txt +++ b/src/biogeochem/test/DustEmis_test/CMakeLists.txt @@ -1,5 +1,7 @@ set (pfunit_sources - test_DustEmisZender2003.pf) + test_DustEmisZender2003.pf + test_DustEmisLeung2023.pf +) add_pfunit_ctest(DustEmis TEST_SOURCES "${pfunit_sources}" diff --git a/src/biogeochem/test/DustEmis_test/test_DustEmisLeung2023.pf b/src/biogeochem/test/DustEmis_test/test_DustEmisLeung2023.pf new file mode 100644 index 0000000000..77e0cba4c0 --- /dev/null +++ b/src/biogeochem/test/DustEmis_test/test_DustEmisLeung2023.pf @@ -0,0 +1,389 @@ +module test_DustEmisLeung2023 + + ! Tests of DustEmisLeung2023 + + use funit + use unittestDustEmisInputs, only : unittest_dust_emis_input_type + use unittestSubgridMod, only : bounds + use DustEmisBase + use DustEmisLeung2023 + use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) + use DustEmisFactory, only : create_dust_emissions + use shr_dust_emis_mod, only : dust_emis_set_options + + implicit none + + real(r8), parameter :: tol = 1.e-16_r8 + + @TestCase + type, extends(TestCase) :: TestDustEmisLeung2023 + class(dust_emis_base_type), allocatable :: dust_emis + type(unittest_dust_emis_input_type) :: input + contains + procedure :: setUp + procedure :: tearDown + procedure :: print_values + procedure :: validate_patch + end type TestDustEmisLeung2023 + +contains + + !----------------------------------------------------------------------- + + subroutine setUp(this) + class(TestDustEmisLeung2023), intent(inout) :: this + ! Allocate and initialize the test object for input objects dust-emission needs + character(len=5) :: NLFilename = 'none' + + call dust_emis_set_options( 'Leung_2023', 'none') + call this%input%setUp( ) + + ! Create the dust emission object last + allocate(this%dust_emis, source = create_dust_emissions(bounds, NLFilename)) + end subroutine setUp + + !----------------------------------------------------------------------- + + subroutine tearDown(this) + class(TestDustEmisLeung2023), intent(inout) :: this + + call this%dust_emis%Clean() + call this%input%tearDown() + end subroutine tearDown + + !----------------------------------------------------------------------- + + subroutine print_values(this) + ! For debugging + use PatchType, only : patch + class(TestDustEmisLeung2023), intent(inout) :: this + real(r8) :: SaltationFactor + integer :: p, c + + call this%input%print_values() + call this%dust_emis%GetConstVars( SaltationFactor ) + do c = bounds%begc, bounds%endc + print *, 'saltation per rho = ', (SaltationFactor / this%input%atm2lnd_inst%forc_rho_downscaled_col(c)) + end do + do p = bounds%begp, bounds%endp + c = patch%column(p) + print *, 'Wind threshold fraction = ', (SaltationFactor / this%input%atm2lnd_inst%forc_rho_downscaled_col(c)) & + / this%input%frictionvel_inst%fv_patch(p) + call this%dust_emis%WritePatchToLog( p ) + end do + end subroutine print_values + + !----------------------------------------------------------------------- + + subroutine validate_patch(this, p) + class(TestDustEmisLeung2023), intent(inout) :: this + integer, intent(in) :: p + + call this%dust_emis%CheckDustEmisIsValid( p ) + end subroutine validate_patch + + !----------------------------------------------------------------------- + + @Test + subroutine check_dust_emis(this) + ! Check dust emissions for default values + class(TestDustEmisLeung2023), intent(inout) :: this + integer :: p + real(r8) :: flx_mss_vrt_dst_tot + real(r8) :: vlc_trb_1 + real(r8) :: vlc_trb_2 + real(r8) :: vlc_trb_3 + real(r8) :: vlc_trb_4 + + call this%input%create_atm2lnd() + call this%input%create_fv() + call this%dust_emis%DustEmission(bounds, this%input%num_nolakep, this%input%filter_nolakep, this%input%atm2lnd_inst, & + this%input%soilstate_inst, this%input%canopystate_inst, this%input%water_inst%waterstatebulk_inst, & + this%input%water_inst%waterdiagnosticbulk_inst, this%input%frictionvel_inst) + call this%dust_emis%DustDryDep(bounds, this%input%atm2lnd_inst, this%input%frictionvel_inst) + call this%print_values() ! Call print subroutine just to make sure it can be used for debugging + do p = bounds%begp, bounds%endp + call this%validate_patch(p) + call this%dust_emis%GetPatchVars( p, flx_mss_vrt_dst_tot=flx_mss_vrt_dst_tot, & + vlc_trb_1=vlc_trb_1, vlc_trb_2=vlc_trb_2, vlc_trb_3=vlc_trb_3, & + vlc_trb_4=vlc_trb_4) + @assertEqual( flx_mss_vrt_dst_tot, 1.305340725852736d-006, tolerance=tol ) + @assertEqual( vlc_trb_1, 3.407721147709135d-003, tolerance=tol ) + @assertEqual( vlc_trb_2, 4.961153753164878d-003, tolerance=tol ) + @assertEqual( vlc_trb_3, 4.980100969983446d-003, tolerance=tol ) + @assertEqual( vlc_trb_4, 4.977071672163210d-003, tolerance=tol ) + end do + + end subroutine check_dust_emis + + !----------------------------------------------------------------------- + + @Test + subroutine dust_zero_for_fixed_ratio(this) + ! Check dust emissions are zero for a fixed ratio + class(TestDustEmisLeung2023), intent(inout) :: this + integer :: p + real(r8) :: flx_mss_vrt_dst_tot + real(r8) :: fv + real(r8) :: SaltationFactor + + call this%input%create_atm2lnd() + call this%dust_emis%GetConstVars( SaltationFactor ) + ! Figure out what fv needs to be so that the wind threshold will result in zero dust emission + fv = ( SaltationFactor / sqrt( this%input%atm2lnd_inst%forc_rho_downscaled_col(bounds%begc)) ) - 1.d-15 + call this%input%create_fv( fv=fv ) + call this%dust_emis%DustEmission(bounds, this%input%num_nolakep, this%input%filter_nolakep, this%input%atm2lnd_inst, & + this%input%soilstate_inst, this%input%canopystate_inst, this%input%water_inst%waterstatebulk_inst, & + this%input%water_inst%waterdiagnosticbulk_inst, this%input%frictionvel_inst) + call this%dust_emis%DustDryDep(bounds, this%input%atm2lnd_inst, this%input%frictionvel_inst) + do p = bounds%begp, bounds%endp + call this%validate_patch(p) + call this%dust_emis%GetPatchVars( p, flx_mss_vrt_dst_tot=flx_mss_vrt_dst_tot ) + @assertEqual( flx_mss_vrt_dst_tot, 0.0_r8 ) + end do + + end subroutine dust_zero_for_fixed_ratio + + !----------------------------------------------------------------------- + + @Test + subroutine dust_zero_when_fsno_one(this) + ! Check dust emissions are zero when snow fraction is identically 1 + class(TestDustEmisLeung2023), intent(inout) :: this + integer :: p + real(r8) :: flx_mss_vrt_dst_tot + + call this%input%create_atm2lnd() + this%input%water_inst%waterdiagnosticbulk_inst%frac_sno_col(:) = 1.0_r8 + call this%input%create_fv( ) + call this%dust_emis%DustEmission(bounds, this%input%num_nolakep, this%input%filter_nolakep, this%input%atm2lnd_inst, & + this%input%soilstate_inst, this%input%canopystate_inst, this%input%water_inst%waterstatebulk_inst, & + this%input%water_inst%waterdiagnosticbulk_inst, this%input%frictionvel_inst) + call this%dust_emis%DustDryDep(bounds, this%input%atm2lnd_inst, this%input%frictionvel_inst) + do p = bounds%begp, bounds%endp + call this%dust_emis%GetPatchVars( p, flx_mss_vrt_dst_tot=flx_mss_vrt_dst_tot ) + @assertEqual( flx_mss_vrt_dst_tot, 0.0_r8 ) + end do + + end subroutine dust_zero_when_fsno_one + + !----------------------------------------------------------------------- + + @Test + subroutine dust_zero_non_veg_lu(this) + ! Check dust emissions are zero for non-veg landunits + use landunit_varcon, only: istcrop, max_lunit + use LandunitType, only : lun + class(TestDustEmisLeung2023), intent(inout) :: this + integer :: p, l + real(r8) :: flx_mss_vrt_dst_tot + + call this%input%create_atm2lnd() + call this%input%create_fv( ) + ! Set the lanunit type for + do l = istcrop+1, max_lunit + lun%itype(bounds%begl:bounds%endl) = l + call this%dust_emis%DustEmission(bounds, this%input%num_nolakep, this%input%filter_nolakep, this%input%atm2lnd_inst, & + this%input%soilstate_inst, this%input%canopystate_inst, this%input%water_inst%waterstatebulk_inst, & + this%input%water_inst%waterdiagnosticbulk_inst, this%input%frictionvel_inst) + call this%dust_emis%DustDryDep(bounds, this%input%atm2lnd_inst, this%input%frictionvel_inst) + do p = bounds%begp, bounds%endp + call this%dust_emis%GetPatchVars( p, flx_mss_vrt_dst_tot=flx_mss_vrt_dst_tot ) + @assertEqual( flx_mss_vrt_dst_tot, 0.0_r8 ) + end do + end do + + end subroutine dust_zero_non_veg_lu + + !----------------------------------------------------------------------- + + @Test + subroutine aborts_on_bad_dust_mobility(this) + ! Check dust abort when dust mobility is bad + class(TestDustEmisLeung2023), intent(inout) :: this + real(r8) :: flx_mss_vrt_dst_tot + character(100) :: expected_msg + + call this%input%create_atm2lnd() + call this%input%create_fv( ) + ! Dust should abort with an error on dust mobility when snow fraction greater than 1 + this%input%water_inst%waterdiagnosticbulk_inst%frac_sno_col(:) = 1.0_r8 + 1.e-15_r8 + call this%dust_emis%DustEmission(bounds, this%input%num_nolakep, this%input%filter_nolakep, this%input%atm2lnd_inst, & + this%input%soilstate_inst, this%input%canopystate_inst, this%input%water_inst%waterstatebulk_inst, & + this%input%water_inst%waterdiagnosticbulk_inst, this%input%frictionvel_inst) + expected_msg = "ABORTED: Bad value for dust mobilization fraction" + @assertExceptionRaised(expected_msg) + + end subroutine aborts_on_bad_dust_mobility + + !----------------------------------------------------------------------- + + @Test + subroutine dust_zero_when_tlai_high(this) + use PatchType, only : patch + ! Check dust emissions are zero when LAI is high enough + class(TestDustEmisLeung2023), intent(inout) :: this + integer :: p + real(r8) :: flx_mss_vrt_dst_tot + + ! Explicitly set the patch type to a hard-coded 1 (so NOT bare-soil) + ! pft indices can't be used without reading them from the parameter file + ! + ! To do this fully the subgrid setup in unittestDustEmisInputs to baresoil + ! should really be run again. But, just doing this is likely sufficient for testing + patch%itype(bounds%begp:bounds%endp) = 1 + call this%input%create_atm2lnd() + call this%input%create_fv( ) + this%input%canopystate_inst%tlai_patch(:) = 1.0_r8 + call this%dust_emis%DustEmission(bounds, this%input%num_nolakep, this%input%filter_nolakep, this%input%atm2lnd_inst, & + this%input%soilstate_inst, this%input%canopystate_inst, this%input%water_inst%waterstatebulk_inst, & + this%input%water_inst%waterdiagnosticbulk_inst, this%input%frictionvel_inst) + call this%dust_emis%DustDryDep(bounds, this%input%atm2lnd_inst, this%input%frictionvel_inst) + do p = bounds%begp, bounds%endp + call this%validate_patch(p) + call this%dust_emis%GetPatchVars( p, flx_mss_vrt_dst_tot=flx_mss_vrt_dst_tot ) + @assertEqual( flx_mss_vrt_dst_tot, 0.0_r8 ) + end do + + end subroutine dust_zero_when_tlai_high + + !----------------------------------------------------------------------- + @Test + subroutine dust_handles_dpfct_rock_nan(this) + use PatchType, only : patch + ! Check dust emissions can handle dpftc_rock being set to NaN + class(TestDustEmisLeung2023), intent(inout) :: this + integer :: p, i + real(r8) :: flx_mss_vrt_dst_tot, drag_partition, previous_dst_tot + + call this%input%create_atm2lnd() + call this%input%create_fv( obu = -500.0_r8, fv=15.0_r8 ) + ! Loop over twice first with nan and then with the tiny value that nan should be changed to + do i = 1, 2 + if ( i == 1 )then + drag_partition = nan + else if ( i == 2 ) then + drag_partition = 0.001 + else + @assertEqual( i, 1, "Iteration too high, should just be up to 2" ) + end if + call this%dust_emis%SetDragPartition(bounds, drag_partition) + call this%dust_emis%DustEmission(bounds, this%input%num_nolakep, this%input%filter_nolakep, this%input%atm2lnd_inst, & + this%input%soilstate_inst, this%input%canopystate_inst, this%input%water_inst%waterstatebulk_inst, & + this%input%water_inst%waterdiagnosticbulk_inst, this%input%frictionvel_inst) + do p = bounds%begp, bounds%endp + call this%dust_emis%GetPatchVars( p, flx_mss_vrt_dst_tot=flx_mss_vrt_dst_tot ) + ! TODO: To have a more robust test the input should be adjusted so that there is dust emissions rather than zero + @assertEqual( flx_mss_vrt_dst_tot, 0.0_r8 ) + if ( i == 1 )then + previous_dst_tot = flx_mss_vrt_dst_tot + else + ! Verify that the result with nan is equal to the result with the tiny value as expected + @assertEqual( flx_mss_vrt_dst_tot, previous_dst_tot ) + end if + end do + end do + + end subroutine dust_handles_dpfct_rock_nan + + !----------------------------------------------------------------------- + + @Test + subroutine dust_constant_low_stability(this) + use PatchType, only : patch + ! Check that dust emissions stay constant until stability high enough + class(TestDustEmisLeung2023), intent(inout) :: this + integer :: p, i + real(r8) :: flx_mss_vrt_dst_tot, obu, previous_dst_tot + + call this%input%create_atm2lnd() + do i = 1, 3 + if ( i == 1 )then + obu = 100._r8 + else if ( i == 2 ) then + obu = 50._r8 + else if ( i == 3 ) then + obu = 41.67013917826486_r8 + else + @assertEqual( i, 1, "Iteration too high, should just be up to 3" ) + end if + call this%input%create_fv( obu=obu ) + call this%dust_emis%DustEmission(bounds, this%input%num_nolakep, this%input%filter_nolakep, this%input%atm2lnd_inst, & + this%input%soilstate_inst, this%input%canopystate_inst, this%input%water_inst%waterstatebulk_inst, & + this%input%water_inst%waterdiagnosticbulk_inst, this%input%frictionvel_inst) + do p = bounds%begp, bounds%endp + call this%dust_emis%GetPatchVars( p, flx_mss_vrt_dst_tot=flx_mss_vrt_dst_tot ) + print *, ' obu = ', obu, ' dust_tot = ', flx_mss_vrt_dst_tot + @assertEqual( flx_mss_vrt_dst_tot, 1.305341766366198d-006, tolerance=tol ) + if ( i == 1 )then + previous_dst_tot = flx_mss_vrt_dst_tot + else + ! Verify that the previous result is identical to the others + @assertEqual( flx_mss_vrt_dst_tot, previous_dst_tot, tolerance=tol ) + end if + end do + end do + + end subroutine dust_constant_low_stability + + !----------------------------------------------------------------------- + + @Test + subroutine dust_aborts_on_zero_obu(this) + use PatchType, only : patch + ! Check dust emissions can handle dpftc_rock being set to NaN + class(TestDustEmisLeung2023), intent(inout) :: this + character(100) :: expected_msg + + call this%input%create_atm2lnd() + call this%input%create_fv( obu=0.0_r8 ) + call this%dust_emis%DustEmission(bounds, this%input%num_nolakep, this%input%filter_nolakep, this%input%atm2lnd_inst, & + this%input%soilstate_inst, this%input%canopystate_inst, this%input%water_inst%waterstatebulk_inst, & + this%input%water_inst%waterdiagnosticbulk_inst, this%input%frictionvel_inst) + expected_msg = "ABORTED: Input obu is zero and can NOT be" + @assertExceptionRaised(expected_msg) + + end subroutine dust_aborts_on_zero_obu + !----------------------------------------------------------------------- + @Test + subroutine check_dust_emis_increasing_fv(this) + ! Check dust emissions with increasing friction velocity + class(TestDustEmisLeung2023), intent(inout) :: this + integer :: p, c + real(r8) :: flx_mss_vrt_dst_tot + real(r8) :: fv = 4.0_r8 + real(r8) :: total_dust0, total_dust_higher + + ! Run baseline fv + call this%input%create_atm2lnd() + call this%input%create_fv( fv=fv ) + call this%dust_emis%DustEmission(bounds, this%input%num_nolakep, this%input%filter_nolakep, this%input%atm2lnd_inst, & + this%input%soilstate_inst, this%input%canopystate_inst, this%input%water_inst%waterstatebulk_inst, & + this%input%water_inst%waterdiagnosticbulk_inst, this%input%frictionvel_inst) + call this%dust_emis%DustDryDep(bounds, this%input%atm2lnd_inst, this%input%frictionvel_inst) + do p = bounds%begp, bounds%endp + call this%validate_patch(p) + call this%dust_emis%GetPatchVars( p, flx_mss_vrt_dst_tot=flx_mss_vrt_dst_tot ) + total_dust0 = flx_mss_vrt_dst_tot + @assertEqual( flx_mss_vrt_dst_tot, 1.064145664977026d-5, tolerance=tol ) + end do + ! Double fv and show result is higher + call this%input%create_fv( fv=fv*2.0_r8) + call this%dust_emis%DustEmission(bounds, this%input%num_nolakep, this%input%filter_nolakep, this%input%atm2lnd_inst, & + this%input%soilstate_inst, this%input%canopystate_inst, this%input%water_inst%waterstatebulk_inst, & + this%input%water_inst%waterdiagnosticbulk_inst, this%input%frictionvel_inst) + call this%dust_emis%DustDryDep(bounds, this%input%atm2lnd_inst, this%input%frictionvel_inst) + do p = bounds%begp, bounds%endp + call this%validate_patch(p) + call this%dust_emis%GetPatchVars( p, flx_mss_vrt_dst_tot=flx_mss_vrt_dst_tot ) + total_dust_higher = flx_mss_vrt_dst_tot + @assertEqual( flx_mss_vrt_dst_tot, 8.309603071671919d-5, tolerance=tol ) + end do + @assertGreaterThan( total_dust_higher, total_dust0 ) + + end subroutine check_dust_emis_increasing_fv + + !----------------------------------------------------------------------- + +end module test_DustEmisLeung2023 diff --git a/src/biogeochem/test/DustEmis_test/test_DustEmisZender2003.pf b/src/biogeochem/test/DustEmis_test/test_DustEmisZender2003.pf index 0289cadabc..6c883d0dfe 100644 --- a/src/biogeochem/test/DustEmis_test/test_DustEmisZender2003.pf +++ b/src/biogeochem/test/DustEmis_test/test_DustEmisZender2003.pf @@ -35,10 +35,10 @@ contains ! Allocate and initialize the test object for input objects dust-emission needs character(len=5) :: NLFilename = 'none' + call dust_emis_set_options( 'Zender_2003', 'atm') call this%input%setUp() ! Create the dust emission object last - call dust_emis_set_options( 'Zender_2003', 'atm') allocate(this%dust_emis, source = create_dust_emissions(bounds, NLFilename)) end subroutine setUp @@ -290,4 +290,4 @@ contains !----------------------------------------------------------------------- -end module test_DustEmisZender2003 \ No newline at end of file +end module test_DustEmisZender2003 diff --git a/src/biogeophys/FrictionVelocityMod.F90 b/src/biogeophys/FrictionVelocityMod.F90 index 042415f545..44e1bf6294 100644 --- a/src/biogeophys/FrictionVelocityMod.F90 +++ b/src/biogeophys/FrictionVelocityMod.F90 @@ -262,12 +262,10 @@ subroutine InitHistory(this, bounds) ptr_patch=this%ram1_patch, default='inactive') end if - if (use_cn) then this%fv_patch(begp:endp) = spval call hist_addfld1d (fname='FV', units='m/s', & avgflag='A', long_name='friction velocity', & - ptr_patch=this%fv_patch) - end if + ptr_patch=this%fv_patch, default='inactive') call hist_addfld1d (fname='RAH1', units='s/m', & avgflag='A', long_name='aerodynamical resistance ', & @@ -455,6 +453,11 @@ subroutine Restart(this, bounds, ncid, flag) long_name='ground momentum roughness length', units='m', & interpinic_flag='interp', readvar=readvar, data=this%z0mg_col) + call restartvar(ncid=ncid, flag=flag, varname='OBU', xtype=ncd_double, & + dim1name='pft', & + long_name='Monin-Obukhov length', units='m', & + interpinic_flag='interp', readvar=readvar, data=this%obu_patch) + if(use_luna)then call restartvar(ncid=ncid, flag=flag, varname='rb10', xtype=ncd_double, & dim1name='pft', long_name='10-day mean boundary layer resistance at the pacth', units='s/m', & diff --git a/src/biogeophys/SoilStateInitTimeConstMod.F90 b/src/biogeophys/SoilStateInitTimeConstMod.F90 index e47577c9b9..e6fcca0f27 100644 --- a/src/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/biogeophys/SoilStateInitTimeConstMod.F90 @@ -10,6 +10,8 @@ module SoilStateInitTimeConstMod use LandunitType , only : lun use ColumnType , only : col use PatchType , only : patch + use abortUtils , only : endrun + use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) ! implicit none private @@ -22,6 +24,7 @@ module SoilStateInitTimeConstMod public :: ThresholdSoilMoistZender2003 public :: ThresholdSoilMoistKok2014 public :: MassFracClay + public :: MassFracClayLeung2023 ! ! !PRIVATE MEMBER FUNCTIONS: private :: ReadNL @@ -183,7 +186,8 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) use organicFileMod , only : organicrd use FuncPedotransferMod , only : pedotransf, get_ipedof use RootBiophysMod , only : init_vegrootfr - use GridcellType , only : grc + use GridcellType , only : grc + use shr_dust_emis_mod , only : is_dust_emis_zender, is_dust_emis_leung ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds @@ -701,14 +705,21 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) end do ! -------------------------------------------------------------------- - ! Initialize threshold soil moisture and mass fracion of clay limited to 0.20 + ! Initialize threshold soil moisture, and mass fraction of clay as + ! scaling coefficient of dust emission flux (kg/m2/s) in each DustEmisType + ! module. See the comments in each function. ! -------------------------------------------------------------------- do c = begc,endc g = col%gridcell(c) soilstate_inst%gwc_thr_col(c) = ThresholdSoilMoistZender2003( clay3d(g,1) ) - soilstate_inst%mss_frc_cly_vld_col(c) = MassFracClay( clay3d(g,1) ) + if ( is_dust_emis_leung() )then + soilstate_inst%mss_frc_cly_vld_col(c) = MassFracClayLeung2023( clay3d(g,1) ) + else + soilstate_inst%mss_frc_cly_vld_col(c) = MassFracClay( clay3d(g,1) ) + end if + end do ! -------------------------------------------------------------------- @@ -725,13 +736,15 @@ end subroutine SoilStateInitTimeConst real(r8) function ThresholdSoilMoistZender2003( clay ) !------------------------------------------------------------------------------ ! - ! Calculate the threshold soil moisture needed for dust emission, based on clay content + ! Calculate the threshold gravimetric water content needed for dust emission, based on clay content ! This was the original equation with a = 1 / (%clay) being the tuning factor for soil - ! moisture effect in Zender's 2003 dust emission scheme. + ! moisture effect in Zender's 2003 dust emission scheme (only for top layer). ! ! 0.17 and 0.14 are fitting coefficients in Fecan et al. (1999), and 0.01 is used to ! convert surface clay fraction from percentage to fraction. + ! The equation comes from Eq. 14 of Fecan et al. (1999; https://doi.org/10.1007/s00585-999-0149-7). ! + ! NOTE: dmleung 19 Feb 2024. !------------------------------------------------------------------------------ ! For future developments Danny M. Leung decided (Dec, 2023) that the Leung et al. (2023) o ! dust emission scheme in the CESM will use Zender's tuning as well, which overall @@ -743,11 +756,13 @@ real(r8) function ThresholdSoilMoistZender2003( clay ) ! ! Also see the notes below for ThresholdSoilMoistKok2014. ! - ! Notes from: dmleung 19 Feb 2024. + ! NOTE on Leung dust emissions: dmleung Jul 5 2024: + ! + ! dmleung followed Zender (2003) DUST scheme to again avoid dust flux being too sensitive to the choice + ! of clay dataset. This is different from what Leung et al. (2023) intended to do. + ! NOTE: This might need to be adjusted for tuning in the future. ! !------------------------------------------------------------------------------ - use abortUtils , only : endrun - use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) real(r8), intent(IN) :: clay ! Fraction of clay in the soil (%) if ( clay < 0.0_r8 .or. clay > 100.0_r8 )then @@ -764,6 +779,8 @@ real(r8) function ThresholdSoilMoistKok2014( clay ) !------------------------------------------------------------------------------ ! Calculate the threshold soil moisture needed for dust emission, based on clay content ! + ! NOTE: dmleung 24 May 2024. + ! ! The below calculates the threshold gravimetric water content for the dust emission ! calculation in DustEmis. The equation comes from Eq. 14 of Fecan et al. ! (1999; https://doi.org/10.1007/s00585-999-0149-7). @@ -775,7 +792,6 @@ real(r8) function ThresholdSoilMoistKok2014( clay ) ! Charlie Zender (2003a) chose: a = 1/clay3d, which gives the ThresholdSoilMoistZender2003 ! function above. ! - ! Notes from: dmleung 24 May 2024. !------------------------------------------------------------------------------ real(r8), intent(IN) :: clay ! Fraction of clay in the soil (%) @@ -793,4 +809,17 @@ end function MassFracClay !------------------------------------------------------------------------------ + real(r8) function MassFracClayLeung2023( clay ) + ! Calculate the mass fraction of clay needed for dust emission, based on clay content + ! Based on the base Zender_2003 version, with a slight modification for Leung_2023 + ! dmleung modified 5 Jul 2024, reducing sensitivity of dust emission + ! flux to clay fraction. + ! NOTE: This might need to be adjusted for tuning in the future. + real(r8), intent(IN) :: clay ! Fraction of lay in the soil (%) + + MassFracClayLeung2023 = 0.1_r8 + MassFracClay( clay ) * 0.1_r8 / 0.20_r8 ! dmleung added this line to reduce the sensitivity of dust emission flux to clay fraction in DUSTMod. 5 Jul 2024 + end function MassFracClayLeung2023 + + !------------------------------------------------------------------------------ + end module SoilStateInitTimeConstMod diff --git a/src/biogeophys/SoilStateType.F90 b/src/biogeophys/SoilStateType.F90 index 2e7225f962..4b9ced3466 100644 --- a/src/biogeophys/SoilStateType.F90 +++ b/src/biogeophys/SoilStateType.F90 @@ -24,7 +24,7 @@ module SoilStateType ! sand/ clay/ organic matter real(r8), pointer :: sandfrac_patch (:) ! patch sand fraction - real(r8), pointer :: clayfrac_patch (:) ! patch clay fraction + real(r8), pointer :: clayfrac_patch (:) ! patch clay fraction real(r8), pointer :: mss_frc_cly_vld_col (:) ! col mass fraction clay limited to 0.20 real(r8), pointer :: cellorg_col (:,:) ! col organic matter for gridcell containing column (1:nlevsoi) real(r8), pointer :: cellsand_col (:,:) ! sand value for gridcell containing column (1:nlevsoi) diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index 707218cc27..899da6b882 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/src/biogeophys/TemperatureType.F90 @@ -9,7 +9,7 @@ module TemperatureType use abortutils , only : endrun use clm_varctl , only : use_cndv, iulog, use_luna, use_crop, use_biomass_heat_storage use clm_varctl , only : flush_gdd20 - use clm_varpar , only : nlevsno, nlevgrnd, nlevlak, nlevurb, nlevmaxurbgrnd + use clm_varpar , only : nlevsno, nlevgrnd, nlevlak, nlevurb, nlevmaxurbgrnd, nlevsoi use clm_varcon , only : spval, ispval use GridcellType , only : grc use LandunitType , only : lun @@ -119,6 +119,10 @@ module TemperatureType real(r8), pointer :: xmf_h2osfc_col (:) ! latent heat of phase change of surface water real(r8), pointer :: fact_col (:,:) ! used in computing tridiagonal matrix real(r8), pointer :: c_h2osfc_col (:) ! heat capacity of surface water + + ! Namelist parameters for initialization + real(r8) :: excess_ice_coldstart_depth ! depth below which excess ice will be present + real(r8) :: excess_ice_coldstart_temp ! coldstart temperature of layers with excess ice present contains @@ -132,6 +136,8 @@ module TemperatureType procedure, public :: UpdateAccVars procedure, private :: UpdateAccVars_CropGDDs + procedure, private :: ReadNL + end type temperature_type character(len=*), parameter, private :: sourcefile = & @@ -143,7 +149,7 @@ module TemperatureType !------------------------------------------------------------------------ subroutine Init(this, bounds, & em_roof_lun, em_wall_lun, em_improad_lun, em_perroad_lun, & - is_simple_buildtemp, is_prog_buildtemp) + is_simple_buildtemp, is_prog_buildtemp, exice_init_conc_col, NLFileName) ! ! !DESCRIPTION: ! @@ -158,7 +164,16 @@ subroutine Init(this, bounds, & real(r8) , intent(in) :: em_perroad_lun(bounds%begl:) logical , intent(in) :: is_simple_buildtemp ! Simple building temp is being used logical , intent(in) :: is_prog_buildtemp ! Prognostic building temp is being used + real(r8) , intent(in) :: exice_init_conc_col(bounds%begc:) ! initial coldstart excess ice concentration (from the stream file) + character(len=*) , intent(in), optional :: NLFilename ! Namelist filename + + if ( present(NLFilename) ) then + call this%ReadNL(NLFilename) + else ! this is needed for testing + this%excess_ice_coldstart_depth = 0.5_r8 + this%excess_ice_coldstart_temp = -5.0_r8 + endif call this%InitAllocate ( bounds ) call this%InitHistory ( bounds, is_simple_buildtemp, is_prog_buildtemp ) call this%InitCold ( bounds, & @@ -166,7 +181,8 @@ subroutine Init(this, bounds, & em_wall_lun(bounds%begl:bounds%endl), & em_improad_lun(bounds%begl:bounds%endl), & em_perroad_lun(bounds%begl:bounds%endl), & - is_simple_buildtemp, is_prog_buildtemp) + is_simple_buildtemp, is_prog_buildtemp, & + exice_init_conc_col(bounds%begc:bounds%endc) ) end subroutine Init @@ -654,7 +670,7 @@ end subroutine InitHistory !----------------------------------------------------------------------- subroutine InitCold(this, bounds, & em_roof_lun, em_wall_lun, em_improad_lun, em_perroad_lun, & - is_simple_buildtemp, is_prog_buildtemp) + is_simple_buildtemp, is_prog_buildtemp, exice_init_conc_col) ! ! !DESCRIPTION: ! Initialize cold start conditions for module variables @@ -662,11 +678,13 @@ subroutine InitCold(this, bounds, & ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use shr_const_mod , only : SHR_CONST_TKFRZ - use clm_varcon , only : denice, denh2o + use clm_varcon , only : denice, denh2o, zisoi use landunit_varcon, only : istwet, istsoil, istdlak, istice, istcrop use column_varcon , only : icol_road_imperv, icol_roof, icol_sunwall use column_varcon , only : icol_shadewall, icol_road_perv - use clm_varctl , only : iulog, use_vancouver, use_mexicocity, use_excess_ice + use clm_varctl , only : iulog, use_vancouver, use_mexicocity + use clm_varctl , only : use_excess_ice + use initVerticalMod , only : find_soil_layer_containing_depth ! ! !ARGUMENTS: class(temperature_type) :: this @@ -677,6 +695,7 @@ subroutine InitCold(this, bounds, & real(r8) , intent(in) :: em_perroad_lun(bounds%begl:) logical , intent(in) :: is_simple_buildtemp ! Simple building temp is being used logical , intent(in) :: is_prog_buildtemp ! Prognostic building temp is being used + real(r8) , intent(in) :: exice_init_conc_col(bounds%begc:) ! initial coldstart excess ice concentration (from the stream file) ! ! !LOCAL VARIABLES: integer :: j,l,c,p ! indices @@ -684,6 +703,7 @@ subroutine InitCold(this, bounds, & real(r8) :: snowbd ! temporary calculation of snow bulk density (kg/m3) real(r8) :: fmelt ! snowbd/100 integer :: lev + integer :: nexice_start ! layer number containing 0.5 meter depth !----------------------------------------------------------------------- SHR_ASSERT_ALL_FL((ubound(em_roof_lun) == (/bounds%endl/)), sourcefile, __LINE__) @@ -757,8 +777,23 @@ subroutine InitCold(this, bounds, & end if else this%t_soisno_col(c,1:nlevgrnd) = 272._r8 - if (use_excess_ice .and. (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop)) then - this%t_soisno_col(c,1:nlevgrnd) = SHR_CONST_TKFRZ - 5.0_r8 !needs to be below freezing to properly initiate excess ice + if (use_excess_ice .and. exice_init_conc_col(c) > 0.0_r8) then + nexice_start = nlevsoi - 1 + if (this%excess_ice_coldstart_depth <= 0.0_r8) then + ! we double check this here, and when building namelists + call endrun(msg="ERROR excess_ice_coldstart_depth <= 0.0. Set a positive value in the namelist"//errmsg(sourcefile, __LINE__)) + endif + if (zisoi(nlevsoi) >= this%excess_ice_coldstart_depth) then + call find_soil_layer_containing_depth(this%excess_ice_coldstart_depth,nexice_start) + else + nexice_start=nlevsoi-1 + endif + if (this%excess_ice_coldstart_temp >= 0.0_r8) then + ! this is here, since we care about excess_ice_coldstart_temp only when there is excess ice in the gridcell + ! which happens only when the streams are read. + call endrun(msg="ERROR excess_ice_coldstart_temp is not below freezing point"//errmsg(sourcefile, __LINE__)) + endif + this%t_soisno_col(c,nexice_start:nlevgrnd) = SHR_CONST_TKFRZ + this%excess_ice_coldstart_temp end if endif endif @@ -1718,4 +1753,67 @@ subroutine Clean(this) end subroutine Clean + !----------------------------------------------------------------------- + subroutine ReadNL( this, NLFilename ) + ! + ! !DESCRIPTION: + ! Read namelist for Temperature type + ! right now (17.07.2024) it only reads variables related to excess ice coldstart initialization + ! but can be extended to replace hardocded values in InitCold by namelist variables + ! + ! !USES: + use shr_mpi_mod , only : shr_mpi_bcast + use shr_log_mod , only : errMsg => shr_log_errMsg + use spmdMod , only : masterproc, mpicom + use fileutils , only : getavu, relavu, opnfil + use clm_nlUtilsMod , only : find_nlgroup_name + use clm_varctl , only : iulog + use abortutils , only : endrun + ! + ! !ARGUMENTS: + class(temperature_type) :: this + character(len=*), intent(in) :: NLFilename ! Namelist filename + ! + ! !LOCAL VARIABLES: + integer :: ierr ! error code + integer :: unitn ! unit for namelist file + real(r8) :: excess_ice_coldstart_temp = spval ! coldstart temperature of layers with excess ice present (deg C) + real(r8) :: excess_ice_coldstart_depth = spval ! depth below which excess ice will be present (m) + character(len=32) :: subname = 'Temperature_readnl' ! subroutine name + !----------------------------------------------------------------------- + + namelist / clm_temperature_inparm / excess_ice_coldstart_depth, excess_ice_coldstart_temp + + if ( masterproc )then + + unitn = getavu() + write(iulog,*) 'Read in clm_temperature_inparm namelist' + call opnfil (NLFilename, unitn, 'F') + call find_nlgroup_name(unitn, 'clm_temperature_inparm', status=ierr) + if (ierr == 0) then + read(unitn, nml=clm_temperature_inparm, iostat=ierr) + if (ierr /= 0) then + call endrun(msg="ERROR reading clm_temperature_inparm namelist"//errmsg(sourcefile, __LINE__)) + end if + else + call endrun(msg="ERROR finding clm_temperature_inparm namelist"//errmsg(sourcefile, __LINE__)) + end if + call relavu( unitn ) + ! namelist might be read but the values not properly set + if ( excess_ice_coldstart_depth == spval ) then + call endrun(msg="ERROR exice_coldstart_depth namelist value is not properly set"//errmsg(sourcefile, __LINE__)) + endif + if ( excess_ice_coldstart_temp == spval ) then + call endrun(msg="ERROR exice_coldstart_temp namelist value is not properly set"//errmsg(sourcefile, __LINE__)) + endif + end if + + call shr_mpi_bcast(excess_ice_coldstart_depth, mpicom) + call shr_mpi_bcast(excess_ice_coldstart_temp, mpicom) + + this%excess_ice_coldstart_depth = excess_ice_coldstart_depth + this%excess_ice_coldstart_temp = excess_ice_coldstart_temp + + end subroutine ReadNL + end module TemperatureType diff --git a/src/biogeophys/WaterStateBulkType.F90 b/src/biogeophys/WaterStateBulkType.F90 index ba3f0513c5..4cd425c976 100644 --- a/src/biogeophys/WaterStateBulkType.F90 +++ b/src/biogeophys/WaterStateBulkType.F90 @@ -47,7 +47,7 @@ module WaterStateBulkType !------------------------------------------------------------------------ subroutine InitBulk(this, bounds, info, vars, & - h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer, NLFilename) + h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer, exice_coldstart_depth, exice_init_conc_col) class(waterstatebulk_type), intent(inout) :: this type(bounds_type) , intent(in) :: bounds @@ -57,7 +57,8 @@ subroutine InitBulk(this, bounds, info, vars, & real(r8) , intent(in) :: watsat_col(bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) real(r8) , intent(in) :: t_soisno_col(bounds%begc:, -nlevsno+1:) ! col soil temperature (Kelvin) logical , intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run - character(len=*) , intent(in) :: NLFilename ! Namelist filename + real(r8) , intent(in) :: exice_coldstart_depth ! depth below which excess ice will be present + real(r8) , intent(in) :: exice_init_conc_col(bounds%begc:) ! initial coldstart excess ice concentration (from the stream file) call this%Init(bounds = bounds, & info = info, & @@ -66,7 +67,7 @@ subroutine InitBulk(this, bounds, info, vars, & watsat_col = watsat_col, & t_soisno_col = t_soisno_col, & use_aquifer_layer = use_aquifer_layer, & - NLFilename = NLFilename) + exice_coldstart_depth = exice_coldstart_depth, exice_init_conc_col = exice_init_conc_col(bounds%begc:bounds%endc)) call this%InitBulkAllocate(bounds) diff --git a/src/biogeophys/WaterStateType.F90 b/src/biogeophys/WaterStateType.F90 index 390e9e8691..35441d65d9 100644 --- a/src/biogeophys/WaterStateType.F90 +++ b/src/biogeophys/WaterStateType.F90 @@ -56,8 +56,6 @@ module WaterStateType real(r8), pointer :: excess_ice_col (:,:) ! col excess ice (kg/m2) (new) (-nlevsno+1:nlevgrnd) real(r8), pointer :: exice_bulk_init (:) ! inital value for excess ice (new) (unitless) - type(excessicestream_type), private :: exicestream ! stream type for excess ice initialization NUOPC only - ! Hillslope stream variables real(r8), pointer :: stream_water_volume_lun(:) ! landunit volume of water in the streams (m3) @@ -82,7 +80,7 @@ module WaterStateType !------------------------------------------------------------------------ subroutine Init(this, bounds, info, tracer_vars, & - h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer, NLFilename) + h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer, exice_coldstart_depth, exice_init_conc_col) class(waterstate_type), intent(inout) :: this type(bounds_type) , intent(in) :: bounds @@ -91,8 +89,9 @@ subroutine Init(this, bounds, info, tracer_vars, & real(r8) , intent(in) :: h2osno_input_col(bounds%begc:) real(r8) , intent(in) :: watsat_col(bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) real(r8) , intent(in) :: t_soisno_col(bounds%begc:, -nlevsno+1:) ! col soil temperature (Kelvin) - logical , intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run - character(len=*) , intent(in) :: NLFilename ! Namelist filename + logical , intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run + real(r8) , intent(in) :: exice_coldstart_depth ! depth below which excess ice will be present + real(r8) , intent(in) :: exice_init_conc_col(bounds%begc:bounds%endc) ! initial coldstart excess ice concentration (from the stream file) this%info => info @@ -104,7 +103,7 @@ subroutine Init(this, bounds, info, tracer_vars, & watsat_col = watsat_col, & t_soisno_col = t_soisno_col, & use_aquifer_layer = use_aquifer_layer, & - NLFilename = NLFilename) + exice_coldstart_depth = exice_coldstart_depth , exice_init_conc_col = exice_init_conc_col) end subroutine Init @@ -323,7 +322,7 @@ end subroutine InitHistory !----------------------------------------------------------------------- subroutine InitCold(this, bounds, & - h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer, NLFilename) + h2osno_input_col, watsat_col, t_soisno_col, use_aquifer_layer, exice_coldstart_depth, exice_init_conc_col) ! ! !DESCRIPTION: ! Initialize time constant variables and cold start conditions @@ -343,11 +342,12 @@ subroutine InitCold(this, bounds, & real(r8) , intent(in) :: watsat_col(bounds%begc:, 1:) ! volumetric soil water at saturation (porosity) real(r8) , intent(in) :: t_soisno_col(bounds%begc:, -nlevsno+1:) ! col soil temperature (Kelvin) logical , intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run - character(len=*) , intent(in) :: NLFilename ! Namelist filename + real(r8) , intent(in) :: exice_coldstart_depth ! depth below which excess ice will be present + real(r8) , intent(in) :: exice_init_conc_col(bounds%begc:bounds%endc) ! initial coldstart excess ice concentration (from the stream file) ! ! !LOCAL VARIABLES: integer :: c,j,l,nlevs,g - integer :: nbedrock, n05m ! layer containing 0.5 m + integer :: nbedrock, nexice ! layer containing 0.5 m real(r8) :: ratio !----------------------------------------------------------------------- @@ -550,52 +550,40 @@ subroutine InitCold(this, bounds, & this%dynbal_baseline_ice_col(bounds%begc:bounds%endc) = 0._r8 !Initialize excess ice - if (use_excess_ice .and. NLFilename /= '') then - ! enforce initialization with 0 for everything - this%excess_ice_col(bounds%begc:bounds%endc,-nlevsno+1:nlevmaxurbgrnd)=0.0_r8 - this%exice_bulk_init(bounds%begc:bounds%endc)=0.0_r8 - call this%exicestream%Init(bounds, NLFilename) ! get initial fraction of excess ice per column - ! - ! If excess ice is being read from streams, use the streams to - ! initialize - ! - if ( UseExcessIceStreams() )then - call this%exicestream%CalcExcessIce(bounds, this%exice_bulk_init) - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - l = col%landunit(c) - if (.not. lun%lakpoi(l)) then !not lake - if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then - if (zisoi(nlevsoi) >= 0.5_r8) then - call find_soil_layer_containing_depth(0.5_r8,n05m) - else - n05m=nlevsoi-1 - endif - if (use_bedrock .and. col%nbedrock(c) <=nlevsoi) then - nbedrock = col%nbedrock(c) - else - nbedrock = nlevsoi - endif - do j = 2, nlevmaxurbgrnd ! ignore first layer - if (n05m= n05m .and. j= exice_coldstart_depth) then + call find_soil_layer_containing_depth(exice_coldstart_depth,nexice) + else + nexice=nlevsoi-1 + endif + if (use_bedrock .and. col%nbedrock(c) <=nlevsoi) then + nbedrock = col%nbedrock(c) + else + nbedrock = nlevsoi endif - else ! just in case zeros for lakes and other columns - this%excess_ice_col(c,-nlevsno+1:nlevmaxurbgrnd) = 0.0_r8 + do j = 2, nlevmaxurbgrnd ! ignore first layer + if (nexice= nexice .and. j shr_kind_r8, CL => shr_kind_cl + use shr_log_mod , only : errMsg => shr_log_errMsg + use spmdMod , only : mpicom, masterproc + use clm_varctl , only : iulog + use abortutils , only : endrun + use decompMod , only : bounds_type + + ! !PUBLIC TYPES: + implicit none + private + + type, public :: prigent_roughness_stream_type + real(r8), pointer, public :: prigent_rghn (:) ! Prigent et al. (1997) roughness length (m) + contains + + ! !PUBLIC MEMBER FUNCTIONS: + procedure, public :: Init ! Initialize and read data in + procedure, public :: UseStreams ! If Prigent rougness streams will be used + procedure, public :: IsStreamInit ! If the streams have been initialized and read in, so data can be used + procedure, public :: Clean ! Clean and deallocate the object + + ! !PRIVATE MEMBER FUNCTIONS: + procedure, private :: InitAllocate ! Allocate data + + end type prigent_roughness_stream_type + + ! ! PRIVATE DATA: + type, private :: streamcontrol_type + character(len=CL) :: stream_fldFileName_prigentroughness ! data Filename + character(len=CL) :: stream_meshfile_prigentroughness ! mesh Filename + character(len=CL) :: prigentroughnessmapalgo ! map algo + contains + procedure, private :: ReadNML ! Read in control namelist + end type streamcontrol_type + + type(streamcontrol_type), private :: control ! Stream control data + logical , private :: NMLRead = .false. ! If namelist has been read + logical , private :: InitDone = .false. ! If initialization of streams has been done + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + +!============================================================================== +contains +!============================================================================== + + subroutine Init(this, bounds, NLFilename) + ! + ! Initialize the prigent roughness stream object + ! + ! Uses: + use spmdMod , only : iam + use lnd_comp_shr , only : mesh, model_clock + use dshr_strdata_mod , only : shr_strdata_init_from_inline, shr_strdata_print + use dshr_strdata_mod , only : shr_strdata_advance + use dshr_methods_mod , only : dshr_fldbun_getfldptr + ! + ! arguments + implicit none + class(prigent_roughness_stream_type) :: this + type(bounds_type), intent(in) :: bounds + character(len=*), intent(in) :: NLFilename ! Namelist filename + ! + ! local variables + integer :: ig, g, n ! Indices + integer :: year ! year (0, ...) for nstep+1 + integer :: mon ! month (1, ..., 12) for nstep+1 + integer :: day ! day of month (1, ..., 31) for nstep+1 + integer :: sec ! seconds into current date for nstep+1 + integer :: mcdate ! Current model date (yyyymmdd) + type(shr_strdata_type) :: sdat_rghn ! input data stream + character(len=16), allocatable :: stream_varnames(:) ! array of stream field names + integer :: rc ! error code + real(r8), pointer :: dataptr1d(:) ! temporary pointer + character(len=*), parameter :: stream_name = 'prigent_roughness' + character(len=*), parameter :: subname = 'PrigentRoughnessStream::Init' + !----------------------------------------------------------------------- + + call control%ReadNML( bounds, NLFileName ) + + call this%InitAllocate( bounds ) + + if ( this%useStreams() )then + + + allocate(stream_varnames(1)) + stream_varnames = (/"Z0a"/) ! varname in cdf5_Z0a_Prigent-Globe-025x025-09262022.nc, in centimeter + + if (masterproc) then + write(iulog,*) ' stream_varnames = ',stream_varnames + end if + + ! Initialize the cdeps data type sdat_rghn + call shr_strdata_init_from_inline(sdat_rghn, & + my_task = iam, & + logunit = iulog, & + compname = 'LND', & + model_clock = model_clock, & + model_mesh = mesh, & + stream_meshfile = control%stream_meshfile_prigentroughness, & + stream_lev_dimname = 'null', & + stream_mapalgo = control%prigentroughnessmapalgo, & + stream_filenames = (/trim(control%stream_fldFileName_prigentroughness)/), & + stream_fldlistFile = stream_varnames, & + stream_fldListModel = stream_varnames, & + stream_yearFirst = 1997, & + stream_yearLast = 1997, & + stream_yearAlign = 1, & + stream_offset = 0, & + stream_taxmode = 'extend', & + stream_dtlimit = 1.0e30_r8, & + stream_tintalgo = 'linear', & + stream_name = 'Prigent roughness', & + rc = rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + ! Explicitly set current date to a hardcoded constant value. Otherwise + ! using the real date can cause roundoff differences that are + ! detrected as issues with exact restart. EBK M05/20/2017 + ! call get_curr_date(year, mon, day, sec) + year = 1997 + mon = 12 + day = 31 + sec = 0 + mcdate = year*10000 + mon*100 + day + + call shr_strdata_advance(sdat_rghn, ymd=mcdate, tod=sec, logunit=iulog, istr='prigentrghn', rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + + ! Get pointer for stream data that is time and spatially interpolate to model time and grid + do n = 1,size(stream_varnames) + call dshr_fldbun_getFldPtr(sdat_rghn%pstrm(1)%fldbun_model, stream_varnames(n), fldptr1=dataptr1d, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then + call ESMF_Finalize(endflag=ESMF_END_ABORT) + end if + if (trim(stream_varnames(n)) == 'Z0a') then + ig = 0 + do g = bounds%begg,bounds%endg + ig = ig+1 + this%prigent_rghn(g) = dataptr1d(ig) + end do + + end if + + end do + deallocate(stream_varnames) + InitDone = .true. + end if + + end subroutine Init + + !============================================================================== + logical function UseStreams(this) + ! + ! !DESCRIPTION: + ! Return true if the Prigent Roughness stream is being used + ! + ! !USES: + ! + ! !ARGUMENTS: + implicit none + class(prigent_roughness_stream_type) :: this + ! + character(len=*), parameter :: subname = 'PrigentRoughnessStream::UseStreams' + ! + if ( this%IsStreamInit() )then + if ( trim(control%stream_fldFileName_prigentroughness) == '' )then + UseStreams = .false. ! Prigent streams are off without a filename given + else + UseStreams = .true. + end if + else + UseStreams = .true. + end if + end function UseStreams + + !============================================================================== + logical function IsStreamInit(this) + ! + ! !DESCRIPTION: + ! Return true if the streams have been initialized + ! + ! !USES: + ! + ! !ARGUMENTS: + implicit none + class(prigent_roughness_stream_type) :: this + ! + character(len=*), parameter :: subname = 'PrigentRoughnessStream::IsStreamInit' + ! + if ( .not. NMLRead )then + call endrun(msg=subname//' ERROR Namelist has NOT been read first, call Init before this') + end if + if ( InitDone )then + IsStreamInit = .true. + else + IsStreamInit = .false. + end if + end function IsStreamInit + + !============================================================================== + subroutine Clean(this) + ! + ! Deallocate and clean the object + ! + ! Uses: + ! + ! arguments + implicit none + class(prigent_roughness_stream_type) :: this + ! + ! local variables + !----------------------------------------------------------------------- + deallocate(this%prigent_rghn) + this%prigent_rghn => NULL() + InitDone = .false. + + end subroutine Clean + + !============================================================================== + subroutine InitAllocate(this, bounds) + ! + ! !DESCRIPTION: + ! Allocate module variables and data structures + ! + ! !USES: + use shr_infnan_mod, only: nan => shr_infnan_nan, assignment(=) + ! + ! !ARGUMENTS: + implicit none + class(prigent_roughness_stream_type) :: this + type(bounds_type), intent(in) :: bounds + ! + ! !LOCAL VARIABLES: + integer :: begg, endg + !--------------------------------------------------------------------- + + begg = bounds%begg; endg = bounds%endg + + if ( this%useStreams() )then + allocate(this%prigent_rghn(begg:endg)) + else + allocate(this%prigent_rghn(0)) + end if + this%prigent_rghn(:) = nan + + end subroutine InitAllocate + + !============================================================================== + subroutine ReadNML(this, bounds, NLFilename) + ! + ! Read the namelist data stream information. + ! + ! Uses: + use shr_nl_mod , only : shr_nl_find_group_name + use shr_log_mod , only : errMsg => shr_log_errMsg + use shr_mpi_mod , only : shr_mpi_bcast + ! + ! arguments + implicit none + class(streamcontrol_type) :: this + type(bounds_type), intent(in) :: bounds + character(len=*), intent(in) :: NLFilename ! Namelist filename + ! + ! local variables + integer :: nu_nml ! unit for namelist file + integer :: nml_error ! namelist i/o error flag + logical :: use_prigent_roughness = .true. + character(len=CL) :: stream_fldFileName_prigentroughness = ' ' + character(len=CL) :: stream_meshfile_prigentroughness = ' ' + character(len=CL) :: prigentroughnessmapalgo = 'bilinear' + character(len=*), parameter :: namelist_name = 'prigentroughness' ! MUST agree with group name in namelist definition to read. + character(len=*), parameter :: subName = "('prigentroughness::ReadNML')" + !----------------------------------------------------------------------- + + namelist /prigentroughness/ & ! MUST agree with namelist_name above + prigentroughnessmapalgo, stream_fldFileName_prigentroughness, stream_meshfile_prigentroughness, & + use_prigent_roughness + + ! Default values for namelist + + ! Read prigentroughness namelist + if (masterproc) then + open( newunit=nu_nml, file=trim(NLFilename), status='old', iostat=nml_error ) + call shr_nl_find_group_name(nu_nml, namelist_name, status=nml_error) + if (nml_error == 0) then + read(nu_nml, nml=prigentroughness,iostat=nml_error) ! MUST agree with namelist_name above + if (nml_error /= 0) then + call endrun(msg=' ERROR reading '//namelist_name//' namelist'//errMsg(sourcefile, __LINE__)) + end if + else + call endrun(msg=' ERROR finding '//namelist_name//' namelist'//errMsg(sourcefile, __LINE__)) + end if + close(nu_nml) + endif + + call shr_mpi_bcast(use_prigent_roughness , mpicom) + call shr_mpi_bcast(prigentroughnessmapalgo , mpicom) + call shr_mpi_bcast(stream_fldFileName_prigentroughness , mpicom) + call shr_mpi_bcast(stream_meshfile_prigentroughness , mpicom) + + ! Error checking + if ( .not. use_prigent_roughness )then + if ( len_trim(stream_fldFileName_prigentroughness) /= 0 )then + call endrun(msg=' ERROR stream_fldFileName_prigentroughness is set, but use_prigent_roughness is FALSE' & + //errMsg(sourcefile, __LINE__)) + end if + if ( len_trim(stream_meshfile_prigentroughness) /= 0 )then + call endrun(msg=' ERROR stream_meshfile_prigentroughness is set, but use_prigent_roughness is FALSE' & + //errMsg(sourcefile, __LINE__)) + end if + else + if ( len_trim(stream_fldFileName_prigentroughness) == 0 )then + call endrun(msg=' ERROR stream_fldFileName_prigentroughness is NOT set, but use_prigent_roughness is TRUE' & + //errMsg(sourcefile, __LINE__)) + end if + if ( len_trim(stream_meshfile_prigentroughness) == 0 )then + call endrun(msg=' ERROR stream_meshfile_prigentroughness is NOT set, but use_prigent_roughness is TRUE' & + //errMsg(sourcefile, __LINE__)) + end if + end if + + if (masterproc) then + write(iulog,*) ' ' + write(iulog,*) namelist_name, ' stream settings:' + write(iulog,*) ' use_prigent_roughness = ',use_prigent_roughness + if ( use_prigent_roughness )then + write(iulog,*) ' stream_fldFileName_prigentroughness = ',trim(stream_fldFileName_prigentroughness) + write(iulog,*) ' stream_meshfile_prigentroughness = ',trim(stream_meshfile_prigentroughness) + write(iulog,*) ' prigentroughnessmapalgo = ',trim(prigentroughnessmapalgo) + end if + endif + this%stream_fldFileName_prigentroughness = stream_fldFileName_prigentroughness + this%stream_meshfile_prigentroughness = stream_meshfile_prigentroughness + this%prigentroughnessmapalgo = prigentroughnessmapalgo + + ! Mark namelist read as having been done + NMLRead = .true. + + end subroutine ReadNML + +end module PrigentRoughnessStreamType diff --git a/src/main/clm_instMod.F90 b/src/main/clm_instMod.F90 index a134265e02..ae83a1b31f 100644 --- a/src/main/clm_instMod.F90 +++ b/src/main/clm_instMod.F90 @@ -11,6 +11,7 @@ module clm_instMod use clm_varctl , only : use_cn, use_c13, use_c14, use_lch4, use_cndv, use_fates, use_fates_bgc use clm_varctl , only : iulog use clm_varctl , only : use_crop, snow_cover_fraction_method, paramfile + use clm_varctl , only : use_excess_ice use SoilBiogeochemDecompCascadeConType , only : mimics_decomp, no_soil_decomp, century_decomp, decomp_method use clm_varcon , only : bdsno, c13ratio, c14ratio use landunit_varcon , only : istice, istsoil @@ -193,6 +194,8 @@ subroutine clm_instInit(bounds) use SoilBiogeochemDecompCascadeBGCMod , only : init_decompcascade_bgc use SoilBiogeochemDecompCascadeContype , only : init_decomp_cascade_constants use SoilBiogeochemCompetitionMod , only : SoilBiogeochemCompetitionInit + use clm_varctl , only : use_excess_ice + use ExcessIceStreamType , only : excessicestream_type, UseExcessIceStreams use initVerticalMod , only : initVertical use SnowHydrologyMod , only : InitSnowLayers @@ -220,6 +223,8 @@ subroutine clm_instInit(bounds) type(file_desc_t) :: params_ncid ! pio netCDF file id for parameter file real(r8), allocatable :: h2osno_col(:) real(r8), allocatable :: snow_depth_col(:) + real(r8), allocatable :: exice_init_conc_col(:) ! initial coldstart excess ice concentration (from the stream file or 0.0) (-) + type(excessicestream_type) :: exice_stream integer :: dummy_to_make_pgi_happy !---------------------------------------------------------------------- @@ -297,12 +302,23 @@ subroutine clm_instInit(bounds) ! Initialization of public data types + ! If excess ice is read from the stream, it has to be read before we coldstart the temperature + allocate(exice_init_conc_col(bounds%begc:bounds%endc)) + exice_init_conc_col(bounds%begc:bounds%endc) = 0.0_r8 + if (use_excess_ice) then + call exice_stream%Init(bounds, NLFilename) + if (UseExcessIceStreams()) then + call exice_stream%CalcExcessIce(bounds, exice_init_conc_col(bounds%begc:bounds%endc)) + endif + endif + call temperature_inst%Init(bounds, & em_roof_lun=urbanparams_inst%em_roof(begl:endl), & em_wall_lun=urbanparams_inst%em_wall(begl:endl), & em_improad_lun=urbanparams_inst%em_improad(begl:endl), & em_perroad_lun=urbanparams_inst%em_perroad(begl:endl), & - is_simple_buildtemp=IsSimpleBuildTemp(), is_prog_buildtemp=IsProgBuildTemp() ) + is_simple_buildtemp=IsSimpleBuildTemp(), is_prog_buildtemp=IsProgBuildTemp(), & + exice_init_conc_col=exice_init_conc_col(bounds%begc:bounds%endc) , NLFileName=NLFilename) call active_layer_inst%Init(bounds) @@ -316,7 +332,9 @@ subroutine clm_instInit(bounds) snow_depth_col = snow_depth_col(begc:endc), & watsat_col = soilstate_inst%watsat_col(begc:endc, 1:), & t_soisno_col = temperature_inst%t_soisno_col(begc:endc, -nlevsno+1:), & - use_aquifer_layer = use_aquifer_layer()) + use_aquifer_layer = use_aquifer_layer(), & + exice_coldstart_depth = temperature_inst%excess_ice_coldstart_depth, & + exice_init_conc_col = exice_init_conc_col(begc:endc)) call glacier_smb_inst%Init(bounds) @@ -454,6 +472,7 @@ subroutine clm_instInit(bounds) deallocate (h2osno_col) deallocate (snow_depth_col) + deallocate (exice_init_conc_col) ! ------------------------------------------------------------------------ ! Initialize accumulated fields diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index 73d528fbc5..cb7e2e3931 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -431,7 +431,7 @@ module clm_varctl !---------------------------------------------------------- - ! excess ice physics switch + ! excess ice physics switch and params !---------------------------------------------------------- logical, public :: use_excess_ice = .false. ! true. => use excess ice physics diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index 787b827605..0b4366a0cb 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -259,7 +259,8 @@ subroutine control_init(dtime) namelist /clm_inparm/ use_soil_moisture_streams - namelist /clm_inparm/ use_excess_ice + ! excess ice flag + namelist /clm_inparm/ use_excess_ice namelist /clm_inparm/ use_lai_streams diff --git a/src/unit_test_shr/unittestDustEmisInputs.F90 b/src/unit_test_shr/unittestDustEmisInputs.F90 index 6b15680ef7..771eb410f7 100644 --- a/src/unit_test_shr/unittestDustEmisInputs.F90 +++ b/src/unit_test_shr/unittestDustEmisInputs.F90 @@ -16,6 +16,8 @@ module unittestDustEmisInputs use FrictionVelocityMod, only : frictionvel_type use unittestWaterTypeFactory, only : unittest_water_type_factory_type use SoilStateInitTimeConstMod, only : ThresholdSoilMoistZender2003, MassFracClay + use SoilStateInitTimeConstMod, only : MassFracClayLeung2023 + use abortutils, only : endrun implicit none private @@ -50,11 +52,13 @@ subroutine setUp(this) character(len=5) :: NLFilename = 'none' real(r8), allocatable :: urb_em(:) + real(r8), allocatable :: exice_init_conc_col(:) integer :: begl, endl, begc, endc integer :: c type(atm2lnd_params_type) :: atm2lnd_params integer, parameter :: snl = 3 + !----------------------------------------------------------------------- ! Settings needed for clm_varpar soil_layerstruct_predefined = '20SL_8.5m' create_crop_landunit = .true. @@ -80,6 +84,7 @@ subroutine setUp(this) glcmec_downscale_longwave = .false., & lapse_rate = 0.01_r8 & ! arbitrary (this is unused for these tests) ) + call this%atm2lnd_inst%InitForTesting(bounds, atm2lnd_params) ! Water and soil state -- after the subgrid setup call this%water_factory%setup_after_subgrid(snl = snl) @@ -91,12 +96,15 @@ subroutine setUp(this) call this%frictionvel_inst%InitForTesting(bounds) allocate( urb_em(begl:endl) ) urb_em(begl:endl) = 0.99_r8 ! Arbitrary won't matter here + allocate( exice_init_conc_col(begc:endc) ) + exice_init_conc_col(begc:endc) = 0.0_r8 ! zero, so it doesn't affect anything. call this%temperature_inst%Init(bounds, & em_roof_lun=urb_em(begl:endl), & em_wall_lun=urb_em(begl:endl), & em_improad_lun=urb_em(begl:endl), & em_perroad_lun=urb_em(begl:endl), & - is_simple_buildtemp=.true., is_prog_buildtemp=.false.) + is_simple_buildtemp=.true., is_prog_buildtemp=.false., & + exice_init_conc_col = exice_init_conc_col(begc:endc)) deallocate( urb_em ) end subroutine setUp @@ -124,6 +132,7 @@ subroutine setupSoilState(this) ! use ColumnType, only : col use GridcellType, only : grc + use shr_dust_emis_mod , only : is_dust_emis_zender, is_dust_emis_leung class(unittest_dust_emis_input_type), intent(in) :: this ! integer :: c,j @@ -145,7 +154,13 @@ subroutine setupSoilState(this) ! These are needed for dust emissions initialization do c = bounds%begc, bounds%endc this%soilstate_inst%gwc_thr_col(c) = ThresholdSoilMoistZender2003( clay ) - this%soilstate_inst%mss_frc_cly_vld_col(c) = MassFracClay( clay ) + if ( is_dust_emis_zender() )then + this%soilstate_inst%mss_frc_cly_vld_col(c) = MassFracClay( clay ) + else if ( is_dust_emis_leung() )then + this%soilstate_inst%mss_frc_cly_vld_col(c) = MassFracClayLeung2023( clay ) + else + call endrun("ERROR: do NOT know about this dust_emis_method") + end if end do end subroutine setupSoilState @@ -187,17 +202,19 @@ end subroutine create_atm2lnd !----------------------------------------------------------------------- - subroutine create_fv(this, fv, u10, ram1) + subroutine create_fv(this, fv, u10, ram1, obu) ! Initializes some fields needed for dust emissions in this%frictionvel_inst, and sets ! fields based on inputs. Excluded inputs are given a default value class(unittest_dust_emis_input_type), intent(inout) :: this real(r8), intent(in), optional :: fv real(r8), intent(in), optional :: u10 real(r8), intent(in), optional :: ram1 + real(r8), intent(in), optional :: obu real(r8), parameter :: fv_default = 2.0_r8 real(r8), parameter :: u10_default = 4._r8 real(r8), parameter :: ram1_default = 200._r8 + real(r8), parameter :: obu_default = -100._r8 ! ------------------------------------------------------------------------ if (present(fv)) then @@ -218,6 +235,12 @@ subroutine create_fv(this, fv, u10, ram1) this%frictionvel_inst%ram1_patch(bounds%begp:bounds%endp) = ram1_default end if + if (present(obu)) then + this%frictionvel_inst%obu_patch(bounds%begp:bounds%endp) = obu + else + this%frictionvel_inst%obu_patch(bounds%begp:bounds%endp) = obu_default + end if + end subroutine create_fv !----------------------------------------------------------------------- diff --git a/src/unit_test_shr/unittestWaterTypeFactory.F90 b/src/unit_test_shr/unittestWaterTypeFactory.F90 index f5f417f9ce..2d83fd0058 100644 --- a/src/unit_test_shr/unittestWaterTypeFactory.F90 +++ b/src/unit_test_shr/unittestWaterTypeFactory.F90 @@ -91,6 +91,7 @@ end subroutine setup_after_subgrid subroutine create_water_type(this, water_inst, & t_soisno_col, watsat_col, & + exice_coldstart_depth, exice_init_conc_col, & enable_consistency_checks, enable_isotopes) ! Initialize water_inst ! @@ -113,6 +114,13 @@ subroutine create_water_type(this, water_inst, & ! 1..nlevgrnd. If not provided, we use arbitrary values for this variable. real(r8), intent(in), optional :: watsat_col( bounds%begc: , 1: ) + ! Depth below which excess ice is present on coldstart. + ! If not provided, we use an arbitrary value. + real(r8), intent(in), optional :: exice_coldstart_depth + ! Excess initial concentration for each column. + ! If not provided, is set to 0. + real(r8), intent(in), optional :: exice_init_conc_col(bounds%begc:) + logical, intent(in), optional :: enable_consistency_checks logical, intent(in), optional :: enable_isotopes @@ -121,6 +129,8 @@ subroutine create_water_type(this, water_inst, & type(water_params_type) :: params real(r8) :: l_watsat_col(bounds%begc:bounds%endc, nlevgrnd) real(r8) :: l_t_soisno_col(bounds%begc:bounds%endc, -nlevsno+1:nlevgrnd) + real(r8) :: l_exice_coldstart_depth + real(r8) :: l_exice_init_conc_col(bounds%begc:bounds%endc) if (present(t_soisno_col)) then SHR_ASSERT_ALL_FL((ubound(t_soisno_col) == [bounds%endc, nlevgrnd]), sourcefile, __LINE__) @@ -129,6 +139,10 @@ subroutine create_water_type(this, water_inst, & SHR_ASSERT_ALL_FL((ubound(watsat_col) == [bounds%endc, nlevgrnd]), sourcefile, __LINE__) end if + if (present(exice_init_conc_col)) then + SHR_ASSERT_ALL_FL((ubound(exice_init_conc_col,1) == bounds%endc), sourcefile, __LINE__) + endif + if (present(enable_consistency_checks)) then l_enable_consistency_checks = enable_consistency_checks else @@ -157,12 +171,24 @@ subroutine create_water_type(this, water_inst, & l_t_soisno_col(:,:) = 275._r8 end if + if (present(exice_coldstart_depth)) then + l_exice_coldstart_depth = exice_coldstart_depth + else + l_exice_coldstart_depth = 0.5_r8 + endif + if (present(exice_init_conc_col)) then + l_exice_init_conc_col(:) = exice_init_conc_col + else + l_exice_init_conc_col(:) = 0.0_r8 + endif + call water_inst%InitForTesting(bounds, params, & h2osno_col = col_array(0._r8), & snow_depth_col = col_array(0._r8), & watsat_col = l_watsat_col, & t_soisno_col = l_t_soisno_col, & - NLFilename = ' ') + exice_coldstart_depth = l_exice_coldstart_depth, & + exice_init_conc_col = l_exice_init_conc_col) end subroutine create_water_type subroutine teardown(this, water_inst) diff --git a/src/unit_test_stubs/share_esmf/CMakeLists.txt b/src/unit_test_stubs/share_esmf/CMakeLists.txt index 5ffce0ba94..1d767543ea 100644 --- a/src/unit_test_stubs/share_esmf/CMakeLists.txt +++ b/src/unit_test_stubs/share_esmf/CMakeLists.txt @@ -1,5 +1,6 @@ list(APPEND clm_sources ExcessIceStreamType.F90 + PrigentRoughnessStreamType.F90 ZenderSoilErodStreamType.F90 ) diff --git a/src/unit_test_stubs/share_esmf/PrigentRoughnessStreamType.F90 b/src/unit_test_stubs/share_esmf/PrigentRoughnessStreamType.F90 new file mode 100644 index 0000000000..d9ded09c30 --- /dev/null +++ b/src/unit_test_stubs/share_esmf/PrigentRoughnessStreamType.F90 @@ -0,0 +1,125 @@ +module PrigentRoughnessStreamType + + !----------------------------------------------------------------------- + ! !DESCRIPTION: + ! UNIT-TEST STUB for Prigent Roughtness Stream + ! This just allows the Leung dust emission code to be tested without + ! reading in the streams data, by faking it and setting it to a + ! constant value. + ! NOTE: THIS SHOULD BE REMOVED WHEN THE MAIN VERSION ALLOWS IT TO BE OFF. + ! https://github.com/ESCOMP/CTSM/issues/2381 + ! Removing this version will remove testing code duplication. + ! !USES + use shr_kind_mod , only : r8 => shr_kind_r8, CL => shr_kind_cl + use shr_log_mod , only : errMsg => shr_log_errMsg + use clm_varctl , only : iulog + use abortutils , only : endrun + use decompMod , only : bounds_type + + ! !PUBLIC TYPES: + implicit none + private + + type, public :: prigent_roughness_stream_type + real(r8), pointer, public :: prigent_rghn (:) ! Prigent et al. (1997) roughness length (m) + contains + + ! !PUBLIC MEMBER FUNCTIONS: + procedure, public :: Init ! Initialize and read data in + procedure, public :: UseStreams ! If Prigent rougness streams will be used + procedure, public :: IsStreamInit ! If the streams have been initialized and read in, so data can be used + procedure, public :: Clean ! Clean and deallocate + + end type prigent_roughness_stream_type + + logical, private :: InitDone = .false. ! If initialization of streams has been done + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + +!============================================================================== +contains +!============================================================================== + + subroutine Init(this, bounds, NLFilename) + ! + ! Initialize the prigent roughness stream object + ! + ! Uses: + ! + ! arguments + implicit none + class(prigent_roughness_stream_type) :: this + type(bounds_type), intent(in) :: bounds + character(len=*), intent(in) :: NLFilename ! Namelist filename + ! + ! local variables + !----------------------------------------------------------------------- + allocate(this%prigent_rghn(bounds%begg:bounds%endg)) + this%prigent_rghn(:) = 1.0_r8 ! This should likely be a smaller value like 0.1 based on the Prigent dataset average values + InitDone = .true. + + end subroutine Init + + !============================================================================== + logical function UseStreams(this) + ! + ! !DESCRIPTION: + ! Return true if the Prigent Roughness stream is being used + ! + ! !USES: + ! + ! !ARGUMENTS: + implicit none + class(prigent_roughness_stream_type) :: this + ! + ! + if ( .not. this%IsStreamInit() )then + UseStreams = .false. + call endrun( "Not initialized" ) + else + UseStreams = .true. + end if + end function UseStreams + + !============================================================================== + logical function IsStreamInit(this) + ! + ! !DESCRIPTION: + ! Return true if the streams have been initialized + ! + ! !USES: + ! + ! !ARGUMENTS: + implicit none + class(prigent_roughness_stream_type) :: this + ! + character(len=*), parameter :: subname = 'PrigentRoughnessStream::IsStreamInit' + ! + if ( InitDone )then + IsStreamInit = .true. + else + IsStreamInit = .false. + end if + end function IsStreamInit + + !============================================================================== + subroutine Clean(this) + ! + ! Deallocate and clean the object + ! + ! Uses: + ! + ! arguments + implicit none + class(prigent_roughness_stream_type) :: this + ! + ! local variables + !----------------------------------------------------------------------- + deallocate(this%prigent_rghn) + this%prigent_rghn => NULL() + InitDone = .false. + + end subroutine Clean + +end module PrigentRoughnessStreamType