diff --git a/bld/namelist_files/namelist_defaults_clm4_5.xml b/bld/namelist_files/namelist_defaults_clm4_5.xml index 8177668942..4301e9c4fb 100644 --- a/bld/namelist_files/namelist_defaults_clm4_5.xml +++ b/bld/namelist_files/namelist_defaults_clm4_5.xml @@ -112,8 +112,8 @@ attributes from the config_cache.xml file (with keys converted to upper-case). 1 0 -.true. -.false. +FAST +NONE .false. diff --git a/bld/namelist_files/namelist_defaults_clm4_5_tools.xml b/bld/namelist_files/namelist_defaults_clm4_5_tools.xml index d6cece9e92..579a83f3c0 100644 --- a/bld/namelist_files/namelist_defaults_clm4_5_tools.xml +++ b/bld/namelist_files/namelist_defaults_clm4_5_tools.xml @@ -624,10 +624,7 @@ attributes from the config_cache.xml file (with keys converted to upper-case). -atm/cam/ggas/ghg_hist_1765-2005_c091218.nc -atm/cam/ggas/ghg_rcp26_1765-2500_c100405.nc -atm/cam/ggas/ghg_rcp45_1765-2500_c100405.nc -atm/cam/ggas/ghg_rcp60_1765-2500_c100901.nc -atm/cam/ggas/ghg_rcp85_1765-2500_c100203.nc +atm/waccm/lb/LBC_17500116-20150116_CMIP6_0p5degLat_c180905.nc +/gpfs/fs1/p/acom/acom-climate/cesm2/inputdata/atm/waccm/lb/LBC_1750-2015_CMIP6_GlobAnnAvg_c180905.nc diff --git a/bld/namelist_files/namelist_definition_clm4_5.xml b/bld/namelist_files/namelist_definition_clm4_5.xml index 0a4bbd9a28..fdf2f5e9e2 100644 --- a/bld/namelist_files/namelist_definition_clm4_5.xml +++ b/bld/namelist_files/namelist_definition_clm4_5.xml @@ -868,9 +868,13 @@ If TRUE, urban traffic flux will be activated (Currently NOT implemented). 1 = prognostic calculation of interior building temp (clm5_0) - -If TRUE, human stress indices will be calculated + +Human heat stress indices: + ALL = All indices will be calculated + FAST = A subset of indices will be calculated (will not include the computationally + expensive wet bulb calculation and associated indices) + NONE = No indices will be calculated + +Toggle to turn on on diagnostic Snow Radiative Effect + + + diff --git a/bld/namelist_files/use_cases/1850_control.xml b/bld/namelist_files/use_cases/1850_control.xml index 1ab17dfaf2..18c5959c78 100644 --- a/bld/namelist_files/use_cases/1850_control.xml +++ b/bld/namelist_files/use_cases/1850_control.xml @@ -35,5 +35,9 @@ lnd/clm2/ndepdata/fndep_clm_WACCM6_CMIP6piControl001_y21-50avg_1850monthly_0.95x1.25_c180802.nc + +cycle + +NDEP_month diff --git a/cime_config/testdefs/ExpectedTestFails.xml b/cime_config/testdefs/ExpectedTestFails.xml index da8d8e1080..fe7373dc38 100644 --- a/cime_config/testdefs/ExpectedTestFails.xml +++ b/cime_config/testdefs/ExpectedTestFails.xml @@ -2,13 +2,12 @@ - FAIL ERP_D_Ld10_P36x2.f10_f10_musgs.IHistClm50BgcCrop.cheyenne_intel.clm-ciso_decStart RUN - FAIL ERS_Lm20_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropGs.cheyenne_gnu.clm-monthly RUN - FAIL ERS_Lm20_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropGs.cheyenne_intel.clm-monthly RUN - FAIL ERS_D_Ld3.f10_f10_musgs.I1850Clm50BgcCrop.cheyenne_gnu.clm-default COMPARE_base_rest - FAIL ERP_D_Ld5.f09_g17.I2000Clm50Vic.cheyenne_intel.clm-vrtlay RUN - FAIL SMS.f10_f10_musgs.I2000Clm50BgcCrop.hobart_pgi.clm-crop RUN - FAIL SMS_D.f10_f10_musgs.I2000Clm50BgcCrop.hobart_pgi.clm-crop RUN + FAIL ERP_D_Ld10_P36x2.f10_f10_musgs.IHistClm50BgcCrop.cheyenne_intel.clm-ciso_decStart RUN + FAIL ERS_Lm20_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropGs.cheyenne_gnu.clm-monthly RUN + FAIL ERS_Lm20_Mmpi-serial.1x1_smallvilleIA.I2000Clm50BgcCropGs.cheyenne_intel.clm-monthly RUN + FAIL SMS.f10_f10_musgs.I2000Clm50BgcCrop.hobart_pgi.clm-crop RUN + FAIL SMS_D.f10_f10_musgs.I2000Clm50BgcCrop.hobart_pgi.clm-crop RUN + FAIL ERS_D_Ln9_P480x3.f19_g16.I2000Clm50SpGs.cheyenne_intel.clm-waccmx_offline MODEL_BUILD FAIL ERS_Ld60.f45_f45_mg37.I2000Clm45Fates.hobart_nag.clm-Fates COMPARE_base_rest diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index e98b0b8753..96d93654ae 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -107,6 +107,14 @@ + + + + + + + + @@ -952,7 +960,7 @@ - + diff --git a/cime_config/testdefs/testmods_dirs/clm/clm50KSinkMOut/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/clm50KSinkMOut/user_nl_clm index b2a51bd5d5..3c3d59e9c3 100644 --- a/cime_config/testdefs/testmods_dirs/clm/clm50KSinkMOut/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/clm/clm50KSinkMOut/user_nl_clm @@ -1,2 +1,3 @@ hist_nhtfrq = 0,-240 hist_mfilt = 1,1 + calc_human_stress_indices = 'ALL' diff --git a/cime_config/testdefs/testmods_dirs/clm/clm50KitchenSink/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/clm50KitchenSink/user_nl_clm index b99b882838..f587667b57 100644 --- a/cime_config/testdefs/testmods_dirs/clm/clm50KitchenSink/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/clm/clm50KitchenSink/user_nl_clm @@ -4,3 +4,4 @@ use_bedrock = .true. use_luna = .true. use_flexibleCN = .true. use_fun = .true. +calc_human_stress_indices = 'ALL' diff --git a/cime_config/testdefs/testmods_dirs/clm/default/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/default/user_nl_clm index 6b0ba531b7..e113c389ce 100644 --- a/cime_config/testdefs/testmods_dirs/clm/default/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/clm/default/user_nl_clm @@ -24,4 +24,4 @@ 'SoilAlpha_U', 'SWup', 'LWup', 'URBAN_AC', 'URBAN_HEAT' - + use_ssre = .true. diff --git a/cime_config/testdefs/testmods_dirs/clm/reduceOutput/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/reduceOutput/user_nl_clm index 3d94564c47..8c98af2f4e 100644 --- a/cime_config/testdefs/testmods_dirs/clm/reduceOutput/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/clm/reduceOutput/user_nl_clm @@ -10,3 +10,5 @@ hist_fincl1 = 'SNOWLIQ','SNOWICE' , 'TG', 'RH2M_U', 'RH2M_R', 'QRUNOFF' +calc_human_stress_indices = 'NONE' +use_ssre = .false. diff --git a/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/README b/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/README new file mode 100644 index 0000000000..a1c390600d --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/README @@ -0,0 +1,4 @@ +The purpose of this test is to run with the changes needed for the WACCM-X configuration +but with CTSM standalone. Since, WACCM-X is unlikely to have initial conditions that are +in equalibrium purposely use a file that will require the model to flush out a shock to the system +at startup. diff --git a/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/include_user_mods new file mode 100644 index 0000000000..fe0e18cf88 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/include_user_mods @@ -0,0 +1 @@ +../default diff --git a/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/shell_commands b/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/shell_commands new file mode 100755 index 0000000000..09dc866921 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/shell_commands @@ -0,0 +1,2 @@ +./xmlchange USE_ESMF_LIB=TRUE,ATM_NCPL=288,CALENDAR=GREGORIAN + diff --git a/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/user_nl_clm new file mode 100644 index 0000000000..2822111c11 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/waccmx_offline/user_nl_clm @@ -0,0 +1,5 @@ +! Add in an initial condition file that isn't in equilibrium and will need to be brought into equilibrium. +! We purposely want to add a shock to ensure the model can properly respond to it, skipping balance checks +! at the beginning and then coming into balance after enough time. + finidat = '$DIN_LOC_ROOT/lnd/clm2/initdata_map/clmi.I1850Clm45BgcGs.0901-01-01.0.9x1.25_gx1v6_simyr1850_c180204.nc' + use_init_interp = .true. diff --git a/doc/ChangeLog b/doc/ChangeLog index 9dac3ba0df..5bf6f59553 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,494 @@ =============================================================== +Tag name: ctsm1.0.dev012 +Originator(s): erik (Erik Kluzek) +Date: Sat Sep 29 11:49:35 MDT 2018 +One-line Summary: Add snow-free fields for snowmip, fix several issues + +Purpose of changes +------------------ + +Bring in new diagnostic fields added by Justin Perket, Sean Swenson and Mark Flanner +for Snow-MIP. Most of those are fields that represent "Snow Free" data. + +Also bring in fixes for a list of issues. Add handling of the new CO2 which includes +both latitude-band and global average versions. Add some changes to make it easier +for input data processing including NOT doing the slow 1km map file creation. Have +the number of steps that are skipped at startup dependent on the time-step size. Add +a test for some requirements of WACCMX (time-step and use of ESMF). Calculations of +local time are now done in a global subroutine, that can handle negative longitudes. +Fix how FFIX_TO_SMINN is handled for history output. The namelist logical "calc_human_stress_indices" +changed from logical to a character string of three values: FAST, NONE, ALL. FAST +is the default so the less expensive ones are output, NONE turns them all off, and ] +ALL does all of them including the expensive ones. + + +Bugs fixed or introduced +------------------------ + +Issues fixed (include CTSM Issue #): #428, #474, #475, #476, #450, #482, #481, #491 + Fix #428 -- Update getco2_historical.ncl to handle latitude varying CO2 + Fix #474 -- Add ability to send GRIDFILE to regridbatch.sh script + Fix #475 -- Have number of steps to skip balance-check based on time + Fix #476 -- Add a test for WACCMX standalone + Fix #450 -- Add option to use global average of terrain standard deviation on surfdata files + (partial fix with simplest option) + Fix #482 -- Add extra field on CO2 streams file for global/time-averaged data + Fix #481 -- FFIX_TO_SMINN needs to be output when FUN is on + Fix #491 -- Calculations of local noon assume that longitude is 0 to 360 rather than -180 to 180 + +Known bugs introduced in this tag (include github issue ID): cime#2801 + cime#2801 -- Problem building with ESMF_LIB + +Known bugs found since the previous tag (include github issue ID): #507, #505 + #507 -- Albedo's are bad at night with negative longitudes + #505 -- CTSM input data-set tools can only work on 0-360 grids, and require monotone increasing longitude + + +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.] + +[ ] clm5_0 + +[ ] clm4_5 + +[ ] clm4_0 + +Notes of particular relevance for users +--------------------------------------- + +Caveats for users (e.g., need to interpolate initial conditions): None + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): + + New namelist: + use_ssre -- Turn on show free fields needed for SnowMIP + + Changed namelist: + calc_human_stress_indices changed from logical to character with options: ALL, FAST, NONE + + New history fields: + Mostly added Snow Free (SF) fields + ALBDSF + ALBISF + FSRSF + FSRSFND + FSRSFNDLN + FSRSFNI + FSRSFVD + FSRSFVDLN + FSRSFVI + SSRE_FSR + SSRE_FSRND + SSRE_FSRNDLN + SSRE_FSRNI + SSRE_FSRVD + SSRE_FSRVDLN + SSRE_FSRVI + +Changes made to namelist defaults (e.g., changed parameter values): + calc_human_stress_indices = 'FAST' is now the default + +Changes to the datasets (e.g., parameter, surface or initial files): + mkghg_bndtvghg -- Update with new CO2 files, both monthly, lat-bands and yearly, global + +Substantial timing or memory changes: None + +Notes of particular relevance for developers: (including Code reviews and testing) +--------------------------------------------- +NOTE: Be sure to review the steps in .CTSMTrunkChecklist as well as the coding style in the Developers Guide + +Caveats for developers (e.g., code that is duplicated that requires double maintenance): + I was able to reduce the duplication in SurfaceAlbedoMod where the original implementation added a cut + and paste copy of code. But, there is still a lot of duplication in this file that could be improved, by + making smaller functions/subroutines to do sections of code that are essentially repeated many times. + There's a bit of an increase in complexity to reduce the duplication, but reducing the duplication was worth it. + +Changes to tests or testing: + Add a new waccmx_offline test mods and test with it + New test expected fail because of cime issue: ERS_D_Ln9_P480x3.f19_g16.I2000Clm50SpGs.cheyenne_intel.clm-waccmx_offline + Turn use_ssre on for most tests, off for reducedOutput + And set calc_human_stress_indices=NONE for reducedOutput, FAST for most, and ALL for KitchenSink and KSMOut tests + +Code reviewed by: self, @olyson, @billsacks + + +CTSM testing: regular + + [PASS means all tests PASS and OK means tests PASS other than expected fails.] + + build-namelist tests: + + cheyenne - PASS + + unit-tests (components/clm/src): + + cheyenne - PASS + + regular tests (aux_clm): + + cheyenne_intel ---- OK + cheyenne_gnu ------ OK + hobart_nag -------- OK + hobart_pgi -------- OK + hobart_intel ------ OK + +CTSM tag used for the baseline comparisons: ctsm1.0.dev011 + + +Answer changes +-------------- + +Changes answers relative to baseline: no bit-for-bit + +Detailed list of changes +------------------------ + +List any externals directories updated (cime, rtm, mosart, cism, fates, etc.): No + +Pull Requests that document the changes (include PR ids): #462 #449 +(https://github.com/ESCOMP/ctsm/pull) + + #462 -- Add namelist item to calculate all heat stress indices only if requested; to speed up model + #449 -- Bring in snowmip diagnostic fields + +=============================================================== +=============================================================== +Tag name: ctsm1.0.dev011 +Originator(s): sacks (Bill Sacks), mvr (Mathew Rothstein) +Date: Wed Sep 12 10:50:31 MDT 2018 +One-line Summary: Add water tracer consistency checks, and other water tracer work + +Purpose of changes +------------------ + +1. Add water tracer consistency checks + +2. Add infrastructure for looping over all water tracers - currently + just used for the tracer consistency checks + +3. Breakout of atm2lnd and lnd2atm water variables, needed for water tracers + +4. Add some namelist control over the addition of water tracers + +5. Add a system test that exercises the water tracer consistency checks + +6. Add a 'ratio' variable for each water tracer + +7. Add some unit tests of the new water tracer infrastructure + +Bugs fixed or introduced +------------------------ + +Issues fixed (include CTSM Issue #): +- Partially addresses #357 +- Resolves #479 +- Resolves #492 + +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.] + +[ ] clm5_0 + +[ ] clm4_5 + +[ ] clm4_0 + +Notes of particular relevance for users +--------------------------------------- + +Caveats for users (e.g., need to interpolate initial conditions): +- Don't be fooled by the new namelist variable, enable_water_isotopes: + This is just a place-holder for now, not implying that water isotopes + are actually working. + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): +- New namelist variables: enable_water_tracer_consistency_checks and + enable_water_isotopes. The latter is just a place-holder for now. + +Changes made to namelist defaults (e.g., changed parameter values): none + +Changes to the datasets (e.g., parameter, surface or initial files): none + +Substantial timing or memory changes: none + +Notes of particular relevance for developers: (including Code reviews and testing) +--------------------------------------------- +NOTE: Be sure to review the steps in .CTSMTrunkChecklist as well as the coding style in the Developers Guide + +Caveats for developers (e.g., code that is duplicated that requires double maintenance): none + +Changes to tests or testing: New test that runs the water tracer consistency check + I ran this test on cheyenne_gnu and cheyenne_intel along with the + out-of-the-box hobart_nag version + +Code reviewed by: Mat Rothstein and I have worked together on many of +these changes, but not all code has been reviewed by the other. + +CTSM testing: + + [PASS means all tests PASS and OK means tests PASS other than expected fails.] + + build-namelist tests: + + cheyenne - ok (tests pass, clm4_5 and clm5_0 namelists differ from + baseline as expected) + + unit-tests (components/clm/src): + + cheyenne - pass + + tools-tests (components/clm/test/tools): + + cheyenne - not run + + PTCLM testing (components/clm/tools/shared/PTCLM/test): + + cheyenne - not run + + regular tests (aux_clm): + + cheyenne_intel ---- pass + cheyenne_gnu ------ pass + hobart_nag -------- pass + hobart_pgi -------- pass + hobart_intel ------ pass + +CTSM tag used for the baseline comparisons: ctsm1.0.dev010 + + +Answer changes +-------------- + +Changes answers relative to baseline: NO + +Detailed list of changes +------------------------ + +List any externals directories updated (cime, rtm, mosart, cism, fates, etc.): none + +Pull Requests that document the changes (include PR ids): +- https://github.com/ESCOMP/ctsm/pull/497 + +=============================================================== +=============================================================== +Tag name: ctsm1.0.dev009 +Originator(s): sacks (Bill Sacks) +Date: Wed Aug 22 20:32:36 MDT 2018 +One-line Summary: Fix initialization of AnnET in InitAccVars + +Purpose of changes +------------------ + +InitAccVars was mistakenly setting qflx_evap_tot_col rather than +AnnET. This fix allows us to remove now-redundant cold start and restart +code for AnnET. + +Bugs fixed or introduced +------------------------ + +Issues fixed (include CTSM Issue #): +- Fixes #480 +- Partially addresses #285 + + +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.] + +[ ] clm5_0 + +[ ] clm4_5 + +[ ] clm4_0 + +Notes of particular relevance for users +--------------------------------------- + +Caveats for users (e.g., need to interpolate initial conditions): none + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): none + +Changes made to namelist defaults (e.g., changed parameter values): none + +Changes to the datasets (e.g., parameter, surface or initial files): none + +Substantial timing or memory changes: none + +Notes of particular relevance for developers: (including Code reviews and testing) +--------------------------------------------- +NOTE: Be sure to review the steps in .CTSMTrunkChecklist as well as the coding style in the Developers Guide + +Caveats for developers (e.g., code that is duplicated that requires double maintenance): none + +Changes to tests or testing: none + +Code reviewed by: basic proposed changes reviewed by Erik Kluzek + + +CTSM testing: + + [PASS means all tests PASS and OK means tests PASS other than expected fails.] + + build-namelist tests: + + cheyenne - not run + + unit-tests (components/clm/src): + + cheyenne - pass + + tools-tests (components/clm/test/tools): + + cheyenne - not run + + PTCLM testing (components/clm/tools/shared/PTCLM/test): + + cheyenne - not run + + regular tests (aux_clm): + + cheyenne_intel ---- pass + cheyenne_gnu ------ pass + hobart_nag -------- pass + hobart_pgi -------- pass + hobart_intel ------ pass + +CTSM tag used for the baseline comparisons: ctsm1.0.dev008 + + +Answer changes +-------------- + +Changes answers relative to baseline: NO + +Detailed list of changes +------------------------ + +List any externals directories updated (cime, rtm, mosart, cism, fates, etc.): none + +Pull Requests that document the changes (include PR ids): none + +=============================================================== +=============================================================== +Tag name: ctsm1.0.dev008 +Originator(s): erik (Erik Kluzek) +Date: Tue Aug 14 10:25:12 MDT 2018 +One-line Summary: Update 1850 ndep file and last year for streams for Historical transient cases + +Purpose of changes +------------------ + +Bring in changes from release-clm5.0.05. Update to latest Nitrogen Deposition file from simulations with WACCM for 1850. +Also fix an issue with the last year for historical transient cases. + + +Bugs fixed or introduced +------------------------ + +Issues fixed (include CTSM Issue #): 461 + #461 -- increase last year in streams for transient + +Known bugs found since the previous tag (include github issue ID): [If none, remove this line] + #478 -- Bare soil g1 should be missing value or zero + + +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.] + +[ ] clm5_0 + +[ ] clm4_5 + +[ ] clm4_0 + +Notes of particular relevance for users +--------------------------------------- + +Caveats for users (e.g., need to interpolate initial conditions): None + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): None + +Changes made to namelist defaults (e.g., changed parameter values): Last year extended for transient datasets + +Changes to the datasets (e.g., parameter, surface or initial files): New ndep dataset for 1850 + +Substantial timing or memory changes: + +Notes of particular relevance for developers: (including Code reviews and testing) +--------------------------------------------- +NOTE: Be sure to review the steps in .CTSMTrunkChecklist as well as the coding style in the Developers Guide + +Caveats for developers (e.g., code that is duplicated that requires double maintenance): None + +Changes to tests or testing: Lengthen some tests + +Code reviewed by: self + + +CTSM testing: regular + + [PASS means all tests PASS and OK means tests PASS other than expected fails.] + + build-namelist tests: + + cheyenne - PASS (11 show differences for 1850_control and 20thC_transient) + + unit-tests (components/clm/src): + + cheyenne - PASS + + regular tests (aux_clm): + + cheyenne_intel ---- OK + cheyenne_gnu ------ OK + hobart_nag -------- OK + hobart_pgi -------- OK + hobart_intel ------ OK + +CTSM tag used for the baseline comparisons: ctsm1.0.dev007 + + +Answer changes +-------------- + +Changes answers relative to baseline: Yes! + + Summarize any changes to answers, i.e., + - what code configurations: 1850_control or 20thC_transient for Clm50 + - what platforms/compilers: all + - nature of change: similar climate + +Detailed list of changes +------------------------ + +List any externals directories updated (cime, rtm, mosart, cism, fates, etc.): None + +Pull Requests that document the changes (include PR ids): +(https://github.com/ESCOMP/ctsm/pull) + + #477 -- Move changes from release-clm5.0.05 onto master + +=============================================================== +=============================================================== Tag name: ctsm1.0.dev007 Originator(s): sacks (Bill Sacks) Date: Sun Aug 5 21:03:28 MDT 2018 diff --git a/doc/ChangeSum b/doc/ChangeSum index 1fdca6b3cd..5df2922ee9 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,15 +1,23 @@ Tag Who Date Summary ============================================================================================================================ +release-clm5.0.09 erik 10/23/2018 Bring in bit-for-bit changes from master up to ctsm1.0.dev012: AnnEt init, snowmip fields + ctsm1.0.dev012 erik 09/29/2018 Add snow-free fields for snowmip, fix several issues + ctsm1.0.dev011 sacks 09/12/2018 SKIPPED ON RELEASE BRANCH -- water tracer consistency checks + ctsm1.0.dev010 sacks 08/30/2018 SKIPPED ON RELEASE BRANCH -- cime update + ctsm1.0.dev009 sacks 08/22/2018 Fix initialization of AnnET in InitAccVars release-clm5.0.08 erik 09/28/2018 Updated CMIP6 ndep file for historical transient Bgc cases, 1850_control same as before release-clm5.0.07 erik 08/08/2018 Bring in some simple fixes from ctsm1.0.dev006 and avoid glacier adjustment at startup from ctsm1.0.dev007 release-clm5.0.06 erik 08/07/2018 Bring in some simple fixes from ctsm1.0.dev006 and avoid glacier adjustment at startup from ctsm1.0.dev007 + ctsm1.0.dev008 erik 08/14/2018 Update 1850 ndep file and last year for streams for Historical transient cases (same as release-clm5.0.05) release-clm5.0.05 erik 08/05/2018 Update 1850 ndep file, and last year for transient streams ctsm1.0.dev007 sacks 08/05/2018 Avoid glacier dynamic landunit adjustments in first time step ctsm1.0.dev006 sacks 08/04/2018 Minor bug fixes, cleanup, documentation and enhancements release-clm5.0.04 erik 07/18/2018 Fix some NFIX variables, update cime/cism for upgraded hobart new glade, new diagnostic fields, update cmip6 output + ctsm1.0.dev005 sacks 08/03/2018 SKIPPED ON RELEASE BRANCH -- water types refactor ctsm1.0.dev004 erik 07/18/2018 Add some new diagnostic fields, fix a few issues, update cmip6 output ctsm1.0.dev003 erik 07/15/2018 Update cime/cism to work on upgraded hobart and with glade changes on cheyenne ctsm1.0.dev002 erik 07/06/2018 Fix NFIX flux variables so special land-units are zeroed out, tools update, add some *_MAX fields on mksurfdata_map for transient cases + ctsm1.0.dev001 sacks 06/22/2018 SKIPPED ON RELEASE BRANCH -- SoilHydrology release-clm5.0.03 erik 06/12/2018 Second release branch tag for CESM2.0 release, fixing DA and tools and README files, identical to clm5.0.dev013 release-clm5.0.02 erik 06/12/2018 Mistake tag identical to previous version clm5.0.dev013 erik 06/12/2018 cleanup and update cime and cism diff --git a/doc/release-clm5.0.ChangeLog b/doc/release-clm5.0.ChangeLog index be078c167a..c0fd701b3c 100644 --- a/doc/release-clm5.0.ChangeLog +++ b/doc/release-clm5.0.ChangeLog @@ -1,4 +1,141 @@ =============================================================== +Tag name: release-clm5.0.09 +Originator(s): erik (Erik Kluzek) +Date: Tue Oct 23 00:00:50 MDT 2018 +One-line Summary: Bring in bit-for-bit changes from master up to ctsm1.0.dev012: AnnEt init, snowmip fields + +Purpose of this version: +------------------------ + +Bring in new diagnostic fields added by Justin Perket, Sean Swenson and Mark Flanner +for Snow-MIP. Most of those are fields that represent "Snow Free" data. + +Also bring in fixes for a list of issues. Add handling of the new CO2 which includes +both latitude-band and global average versions. Add some changes to make it easier +for input data processing including NOT doing the slow 1km map file creation. Have +the number of steps that are skipped at startup dependent on the time-step size. Add +a test for some requirements of WACCMX (time-step and use of ESMF). Calculations of +local time are now done in a global subroutine, that can handle negative longitudes. +Fix how FFIX_TO_SMINN is handled for history output. The namelist logical "calc_human_stress_indices" +changed from logical to a character string of three values: FAST, NONE, ALL. FAST +is the default so the less expensive ones are output, NONE turns them all off, and ] +ALL does all of them including the expensive ones. + +InitAccVars was mistakenly setting qflx_evap_tot_col rather than +AnnET. This fix allows us to remove now-redundant cold start and restart +code for AnnET. + +CTSM Master Tag This Corresponds To: ctsm1.0.dev012 + +Summary of changes: +------------------- + +Issues fixed (include CTSM Issue #): + Fix #428 -- Update getco2_historical.ncl to handle latitude varying CO2 + Fix #474 -- Add ability to send GRIDFILE to regridbatch.sh script + Fix #475 -- Have number of steps to skip balance-check based on time + Fix #476 -- Add a test for WACCMX standalone + Fix #450 -- Add option to use global average of terrain standard deviation on surfdata files + (partial fix with simplest option) + Fix #482 -- Add extra field on CO2 streams file for global/time-averaged data + Fix #481 -- FFIX_TO_SMINN needs to be output when FUN is on + Fix #491 -- Calculations of local noon assume that longitude is 0 to 360 rather than -180 to 180 + Fix #480 InitAccVars for AnnET initializing the wrong variable + Fix #285 Remove an un-needed restart variable (partially) + +Science changes since: release-clm5.0.08 + + * None + +Software changes since: release-clm5.0.08 + + New fields and fix a list of issues + +Changes to User Interface since: release-clm5.0.08 + + New namelist: + use_ssre -- Turn on show free fields needed for SnowMIP + + Changed namelist: + calc_human_stress_indices changed from logical to character with options: ALL, FAST, NONE + + New history fields: + Mostly added Snow Free (SF) fields + ALBDSF + ALBISF + FSRSF + FSRSFND + FSRSFNDLN + FSRSFNI + FSRSFVD + FSRSFVDLN + FSRSFVI + SSRE_FSR + SSRE_FSRND + SSRE_FSRNDLN + SSRE_FSRNI + SSRE_FSRVD + SSRE_FSRVDLN + SSRE_FSRVI + +Testing: +-------- + + [PASS means all tests PASS and OK means tests PASS other than expected fails.] + + build-namelist tests: + + cheyenne - PASS + + unit-tests (components/clm/src): + + cheyenne - PASS + hobart ---PASS + + tools-tests (components/clm/test/tools): + + cheyenne - OK (PTCLM tests fail) + + PTCLM testing (components/clm/tools/shared/PTCLM/test): + + cheyenne - PASS + + regular tests (aux_clm): + + cheyenne_intel ---- OK + cheyenne_gnu ------ OK + hobart_nag -------- OK + hobart_pgi -------- OK + hobart_intel ------ OK + +Summary of Answer changes: +------------------------- + +Baseline version for comparison: release-clm5.0.08 + +Changes answers relative to baseline: No bit-for-bit + +Detailed list of changes: +------------------------ + +Externals being used: + + cism: release-cesm2.0.04 + rtm: release-cesm2.0.00 + mosart: release-cesm2.0.00 + cime: cime_cesm2_0_rel_05 + FATES: fates_s1.8.1_a3.0.0 + PTCLM: PTCLM2_180611 + +CTSM Tag versions pulled over from master development branch: + ctsm1.0.dev008, ctsm1.0.dev009, ctsm1.0.dev012 + +Pull Requests that document the changes (include PR ids): +(https://github.com/ESCOMP/ctsm/pull) + #543 -- Update release branch to ctsm1.0.dev012 + +=============================================================== +=============================================================== Tag name: release-clm5.0.08 Originator(s): erik (Erik Kluzek,UCAR/TSS,303-497-1326) Date: Fri Sep 28 14:17:52 MDT 2018 diff --git a/src/biogeochem/CNVegetationFacade.F90 b/src/biogeochem/CNVegetationFacade.F90 index 9f555f6b14..1e24a98193 100644 --- a/src/biogeochem/CNVegetationFacade.F90 +++ b/src/biogeochem/CNVegetationFacade.F90 @@ -179,13 +179,16 @@ module CNVegetationFacade procedure, private :: CNReadNML ! Read in the CN general namelist end type cn_vegetation_type + ! !PRIVATE DATA MEMBERS: + + integer, private :: skip_steps ! Number of steps to skip at startup character(len=*), parameter, private :: sourcefile = & __FILE__ contains !----------------------------------------------------------------------- - subroutine Init(this, bounds, NLFilename) + subroutine Init(this, bounds, NLFilename, nskip_steps) ! ! !DESCRIPTION: ! Initialize a CNVeg object. @@ -199,7 +202,8 @@ subroutine Init(this, bounds, NLFilename) ! !ARGUMENTS: class(cn_vegetation_type), intent(inout) :: this type(bounds_type), intent(in) :: bounds - character(len=*) , intent(in) :: NLFilename ! namelist filename + character(len=*) , intent(in) :: NLFilename ! namelist filename + integer , intent(in) :: nskip_steps ! Number of steps to skip at startup ! ! !LOCAL VARIABLES: integer :: begp, endp @@ -213,6 +217,8 @@ subroutine Init(this, bounds, NLFilename) ! Note - always initialize the memory for cnveg_state_inst (used in biogeophys/) call this%cnveg_state_inst%Init(bounds) + skip_steps = nskip_steps + if (use_cn) then ! Read in the general CN namelist @@ -1002,7 +1008,7 @@ subroutine BalanceCheck(this, bounds, num_soilc, filter_soilc, & !----------------------------------------------------------------------- DA_nstep = get_nstep_since_startup_or_lastDA_restart_or_pause() - if (DA_nstep < 2 )then + if (DA_nstep <= skip_steps )then if (masterproc) then write(iulog,*) '--WARNING-- skipping CN balance check for first timesteps after startup or data assimilation' end if diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index 91fdd30ad3..1db44262e7 100644 --- a/src/biogeophys/BalanceCheckMod.F90 +++ b/src/biogeophys/BalanceCheckMod.F90 @@ -35,11 +35,19 @@ module BalanceCheckMod ! !PUBLIC TYPES: implicit none save + private ! ! !PUBLIC MEMBER FUNCTIONS: + public :: BalanceCheckInit ! Initialization of Water and energy balance check public :: BeginWaterBalance ! Initialize water balance check public :: BalanceCheck ! Water and energy balance check + public :: GetBalanceCheckSkipSteps ! Get the number of steps to skip for the balance check + public :: BalanceCheckClean ! Clean up for BalanceCheck + + ! !PRIVATE MEMBER DATA: + real(r8), private, parameter :: skip_size = 3600.0_r8 ! Time steps to skip the balance check at startup (sec) + integer, private :: skip_steps = -999 ! Number of time steps to skip the balance check for character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -47,6 +55,63 @@ module BalanceCheckMod contains + !----------------------------------------------------------------------- + subroutine BalanceCheckInit( ) + !----------------------------------------------------------------------- + ! + ! !DESCRIPTION: + ! Initialize balance check + ! + ! !USES: + use spmdMod , only : masterproc + use clm_time_manager , only : get_step_size + ! !ARGUMENTS: + ! + ! !LOCAL VARIABLES: + real(r8) :: dtime ! land model time step (sec) + !----------------------------------------------------------------------- + dtime = get_step_size() + ! Skip a minimum of two time steps, but otherwise skip the number of time-steps in the skip_size rounded to the nearest integer + skip_steps = max(2, nint( (skip_size / dtime) ) ) + + if ( masterproc ) write(iulog,*) ' Skip balance checking for the first ', skip_steps, ' time steps' + + end subroutine BalanceCheckInit + + !----------------------------------------------------------------------- + subroutine BalanceCheckClean( ) + !----------------------------------------------------------------------- + ! + ! !DESCRIPTION: + ! Clean up BalanceCheck + ! + ! !USES: + ! !ARGUMENTS: + ! + ! !LOCAL VARIABLES: + !----------------------------------------------------------------------- + skip_steps = -999 + end subroutine BalanceCheckClean + + !----------------------------------------------------------------------- + function GetBalanceCheckSkipSteps( ) result( get_skip ) + !----------------------------------------------------------------------- + ! + ! !DESCRIPTION: + ! Get the number of steps to skip for the balance check + ! + ! !ARGUMENTS: + integer :: get_skip ! function result + ! !LOCAL VARIABLES: + if ( skip_steps > 0 )then + get_skip = skip_steps + else + get_skip = -999 + call endrun('ERROR:: GetBalanceCheckSkipSteps called before BalanceCheckInit') + end if + end function GetBalanceCheckSkipSteps + !----------------------------------------------------------------------- + !----------------------------------------------------------------------- subroutine BeginWaterBalance(bounds, & num_nolakec, filter_nolakec, num_lakec, filter_lakec, & @@ -97,7 +162,7 @@ end subroutine BeginWaterBalance !----------------------------------------------------------------------- subroutine BalanceCheck( bounds, & atm2lnd_inst, solarabs_inst, waterflux_inst, waterstate_inst, & - irrigation_inst, glacier_smb_inst, energyflux_inst, canopystate_inst) + irrigation_inst, glacier_smb_inst, surfalb_inst, energyflux_inst, canopystate_inst) ! ! !DESCRIPTION: ! This subroutine accumulates the numerical truncation errors of the water @@ -117,8 +182,8 @@ subroutine BalanceCheck( bounds, & use clm_varcon , only : spval use clm_time_manager , only : get_step_size, get_nstep use clm_time_manager , only : get_nstep_since_startup_or_lastDA_restart_or_pause - use clm_instMod , only : surfalb_inst use CanopyStateType , only : canopystate_type + use SurfaceAlbedoType , only : surfalb_type use subgridAveMod ! ! !ARGUMENTS: @@ -129,6 +194,7 @@ subroutine BalanceCheck( bounds, & type(waterstate_type) , intent(inout) :: waterstate_inst type(irrigation_type) , intent(in) :: irrigation_inst type(glacier_smb_type), intent(in) :: glacier_smb_inst + type(surfalb_type) , intent(inout) :: surfalb_inst type(energyflux_type) , intent(inout) :: energyflux_inst type(canopystate_type), intent(inout) :: canopystate_inst ! @@ -308,7 +374,7 @@ subroutine BalanceCheck( bounds, & if ((col%itype(indexc) == icol_roof .or. & col%itype(indexc) == icol_road_imperv .or. & col%itype(indexc) == icol_road_perv) .and. & - abs(errh2o(indexc)) > 1.e-5_r8 .and. (DAnstep > 2) ) then + abs(errh2o(indexc)) > 1.e-5_r8 .and. (DAnstep > skip_steps) ) then write(iulog,*)'clm urban model is stopping - error is greater than 1e-5 (mm)' write(iulog,*)'nstep = ',nstep @@ -337,7 +403,7 @@ subroutine BalanceCheck( bounds, & write(iulog,*)'clm model is stopping' call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__)) - else if (abs(errh2o(indexc)) > 1.e-5_r8 .and. (DAnstep > 2) ) then + else if (abs(errh2o(indexc)) > 1.e-5_r8 .and. (DAnstep > skip_steps) ) then write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)' write(iulog,*)'nstep = ',nstep @@ -436,7 +502,7 @@ subroutine BalanceCheck( bounds, & ' lun%itype= ',lun%itype(col%landunit(indexc)), & ' errh2osno= ',errh2osno(indexc) - if (abs(errh2osno(indexc)) > 1.e-5_r8 .and. (DAnstep > 2) ) then + if (abs(errh2osno(indexc)) > 1.e-5_r8 .and. (DAnstep > skip_steps) ) then write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)' write(iulog,*)'nstep = ',nstep write(iulog,*)'errh2osno = ',errh2osno(indexc) @@ -528,7 +594,7 @@ subroutine BalanceCheck( bounds, & end if end if end do - if ( found .and. (DAnstep > 2) ) then + if ( found .and. (DAnstep > skip_steps) ) then write(iulog,*)'WARNING:: BalanceCheck, solar radiation balance error (W/m2)' write(iulog,*)'nstep = ',nstep write(iulog,*)'errsol = ',errsol(indexp) @@ -558,7 +624,7 @@ subroutine BalanceCheck( bounds, & end if end if end do - if ( found .and. (DAnstep > 2) ) then + if ( found .and. (DAnstep > skip_steps) ) then write(iulog,*)'WARNING: BalanceCheck: longwave energy balance error (W/m2)' write(iulog,*)'nstep = ',nstep write(iulog,*)'errlon = ',errlon(indexp) @@ -581,7 +647,7 @@ subroutine BalanceCheck( bounds, & end if end if end do - if ( found .and. (DAnstep > 2) ) then + if ( found .and. (DAnstep > skip_steps) ) then write(iulog,*)'WARNING: BalanceCheck: surface flux energy balance error (W/m2)' write(iulog,*)'nstep = ' ,nstep write(iulog,*)'errseb = ' ,errseb(indexp) @@ -624,7 +690,7 @@ subroutine BalanceCheck( bounds, & write(iulog,*)'WARNING: BalanceCheck: soil balance error (W/m2)' write(iulog,*)'nstep = ',nstep write(iulog,*)'errsoi_col = ',errsoi_col(indexc) - if (abs(errsoi_col(indexc)) > 1.e-4_r8 .and. (DAnstep > 2) ) then + if (abs(errsoi_col(indexc)) > 1.e-4_r8 .and. (DAnstep > skip_steps) ) then write(iulog,*)'clm model is stopping' call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__)) end if diff --git a/src/biogeophys/BareGroundFluxesMod.F90 b/src/biogeophys/BareGroundFluxesMod.F90 index cfbc99b740..5b914a3a07 100644 --- a/src/biogeophys/BareGroundFluxesMod.F90 +++ b/src/biogeophys/BareGroundFluxesMod.F90 @@ -51,7 +51,8 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, & use FrictionVelocityMod , only : FrictionVelocity, MoninObukIni use QSatMod , only : QSat use SurfaceResistanceMod , only : do_soilevap_beta,do_soil_resistance_sl14 - use HumanIndexMod , only : calc_human_stress_indices, Wet_Bulb, Wet_BulbS, HeatIndex, AppTemp, & + use HumanIndexMod , only : all_human_stress_indices, fast_human_stress_indices, & + Wet_Bulb, Wet_BulbS, HeatIndex, AppTemp, & swbgt, hmdex, dis_coi, dis_coiS, THIndex, & SwampCoolEff, KtoC, VaporPres use FrictionVelocityMod , only : frictionvel_parms_inst @@ -390,37 +391,42 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, & end if ! Human Heat Stress - if ( calc_human_stress_indices )then + if ( all_human_stress_indices .or. fast_human_stress_indices ) then call KtoC(t_ref2m(p), tc_ref2m(p)) call VaporPres(rh_ref2m(p), e_ref2m, vap_ref2m(p)) - call Wet_Bulb(t_ref2m(p), vap_ref2m(p), forc_pbot(c), rh_ref2m(p), q_ref2m(p), & - teq_ref2m(p), ept_ref2m(p), wb_ref2m(p)) call Wet_BulbS(tc_ref2m(p),rh_ref2m(p), wbt_ref2m(p)) call HeatIndex(tc_ref2m(p), rh_ref2m(p), nws_hi_ref2m(p)) call AppTemp(tc_ref2m(p), vap_ref2m(p), u10_clm(p), appar_temp_ref2m(p)) call swbgt(tc_ref2m(p), vap_ref2m(p), swbgt_ref2m(p)) call hmdex(tc_ref2m(p), vap_ref2m(p), humidex_ref2m(p)) - call dis_coi(tc_ref2m(p), wb_ref2m(p), discomf_index_ref2m(p)) call dis_coiS(tc_ref2m(p), rh_ref2m(p), wbt_ref2m(p), discomf_index_ref2mS(p)) - call THIndex(tc_ref2m(p), wb_ref2m(p), thic_ref2m(p), thip_ref2m(p)) - call SwampCoolEff(tc_ref2m(p), wb_ref2m(p), swmp80_ref2m(p), swmp65_ref2m(p)) + if ( all_human_stress_indices ) then + call Wet_Bulb(t_ref2m(p), vap_ref2m(p), forc_pbot(c), rh_ref2m(p), q_ref2m(p), & + teq_ref2m(p), ept_ref2m(p), wb_ref2m(p)) + call dis_coi(tc_ref2m(p), wb_ref2m(p), discomf_index_ref2m(p)) + call THIndex(tc_ref2m(p), wb_ref2m(p), thic_ref2m(p), thip_ref2m(p)) + call SwampCoolEff(tc_ref2m(p), wb_ref2m(p), swmp80_ref2m(p), swmp65_ref2m(p)) + end if if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then - teq_ref2m_r(p) = teq_ref2m(p) - ept_ref2m_r(p) = ept_ref2m(p) - wb_ref2m_r(p) = wb_ref2m(p) wbt_ref2m_r(p) = wbt_ref2m(p) nws_hi_ref2m_r(p) = nws_hi_ref2m(p) appar_temp_ref2m_r(p) = appar_temp_ref2m(p) swbgt_ref2m_r(p) = swbgt_ref2m(p) humidex_ref2m_r(p) = humidex_ref2m(p) - discomf_index_ref2m_r(p) = discomf_index_ref2m(p) discomf_index_ref2mS_r(p) = discomf_index_ref2mS(p) - thic_ref2m_r(p) = thic_ref2m(p) - thip_ref2m_r(p) = thip_ref2m(p) - swmp80_ref2m_r(p) = swmp80_ref2m(p) - swmp65_ref2m_r(p) = swmp65_ref2m(p) + if ( all_human_stress_indices ) then + teq_ref2m_r(p) = teq_ref2m(p) + ept_ref2m_r(p) = ept_ref2m(p) + wb_ref2m_r(p) = wb_ref2m(p) + discomf_index_ref2m_r(p) = discomf_index_ref2m(p) + thic_ref2m_r(p) = thic_ref2m(p) + thip_ref2m_r(p) = thip_ref2m(p) + swmp80_ref2m_r(p) = swmp80_ref2m(p) + swmp65_ref2m_r(p) = swmp65_ref2m(p) + end if end if + end if end do diff --git a/src/biogeophys/CMakeLists.txt b/src/biogeophys/CMakeLists.txt index 55ad7e1c65..d91cc1cdf1 100644 --- a/src/biogeophys/CMakeLists.txt +++ b/src/biogeophys/CMakeLists.txt @@ -4,6 +4,9 @@ list(APPEND clm_sources AerosolMod.F90 DaylengthMod.F90 + BalanceCheckMod.F90 + CanopyStateType.F90 + EnergyFluxType.F90 GlacierSurfaceMassBalanceMod.F90 HumanIndexMod.F90 IrrigationMod.F90 @@ -15,6 +18,8 @@ list(APPEND clm_sources SoilHydrologyType.F90 SoilStateType.F90 SoilWaterRetentionCurveMod.F90 + SolarAbsorbedType.F90 + SurfaceAlbedoType.F90 TemperatureType.F90 TotalWaterAndHeatMod.F90 UrbanParamsType.F90 diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index 9b662bb372..8efbb2e6ed 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -176,7 +176,8 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, use QSatMod , only : QSat use CLMFatesInterfaceMod, only : hlm_fates_interface_type use FrictionVelocityMod, only : FrictionVelocity, MoninObukIni, frictionvel_parms_inst - use HumanIndexMod , only : calc_human_stress_indices, Wet_Bulb, Wet_BulbS, HeatIndex, AppTemp, & + use HumanIndexMod , only : all_human_stress_indices, fast_human_stress_indices, & + Wet_Bulb, Wet_BulbS, HeatIndex, AppTemp, & swbgt, hmdex, dis_coi, dis_coiS, THIndex, & SwampCoolEff, KtoC, VaporPres use SoilWaterRetentionCurveMod, only : soil_water_retention_curve_type @@ -1193,36 +1194,39 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, rh_ref2m_r(p) = rh_ref2m(p) ! Human Heat Stress - if ( calc_human_stress_indices )then - + if ( all_human_stress_indices .or. fast_human_stress_indices ) then call KtoC(t_ref2m(p), tc_ref2m(p)) call VaporPres(rh_ref2m(p), e_ref2m, vap_ref2m(p)) - call Wet_Bulb(t_ref2m(p), vap_ref2m(p), forc_pbot(c), rh_ref2m(p), q_ref2m(p), & - teq_ref2m(p), ept_ref2m(p), wb_ref2m(p)) call Wet_BulbS(tc_ref2m(p),rh_ref2m(p), wbt_ref2m(p)) call HeatIndex(tc_ref2m(p), rh_ref2m(p), nws_hi_ref2m(p)) call AppTemp(tc_ref2m(p), vap_ref2m(p), u10_clm(p), appar_temp_ref2m(p)) call swbgt(tc_ref2m(p), vap_ref2m(p), swbgt_ref2m(p)) call hmdex(tc_ref2m(p), vap_ref2m(p), humidex_ref2m(p)) - call dis_coi(tc_ref2m(p), wb_ref2m(p), discomf_index_ref2m(p)) call dis_coiS(tc_ref2m(p), rh_ref2m(p), wbt_ref2m(p), discomf_index_ref2mS(p)) - call THIndex(tc_ref2m(p), wb_ref2m(p), thic_ref2m(p), thip_ref2m(p)) - call SwampCoolEff(tc_ref2m(p), wb_ref2m(p), swmp80_ref2m(p), swmp65_ref2m(p)) - - teq_ref2m_r(p) = teq_ref2m(p) - ept_ref2m_r(p) = ept_ref2m(p) - wb_ref2m_r(p) = wb_ref2m(p) + if ( all_human_stress_indices ) then + call Wet_Bulb(t_ref2m(p), vap_ref2m(p), forc_pbot(c), rh_ref2m(p), q_ref2m(p), & + teq_ref2m(p), ept_ref2m(p), wb_ref2m(p)) + call dis_coi(tc_ref2m(p), wb_ref2m(p), discomf_index_ref2m(p)) + call THIndex(tc_ref2m(p), wb_ref2m(p), thic_ref2m(p), thip_ref2m(p)) + call SwampCoolEff(tc_ref2m(p), wb_ref2m(p), swmp80_ref2m(p), swmp65_ref2m(p)) + end if wbt_ref2m_r(p) = wbt_ref2m(p) nws_hi_ref2m_r(p) = nws_hi_ref2m(p) appar_temp_ref2m_r(p) = appar_temp_ref2m(p) swbgt_ref2m_r(p) = swbgt_ref2m(p) humidex_ref2m_r(p) = humidex_ref2m(p) - discomf_index_ref2m_r(p) = discomf_index_ref2m(p) discomf_index_ref2mS_r(p) = discomf_index_ref2mS(p) - thic_ref2m_r(p) = thic_ref2m(p) - thip_ref2m_r(p) = thip_ref2m(p) - swmp80_ref2m_r(p) = swmp80_ref2m(p) - swmp65_ref2m_r(p) = swmp65_ref2m(p) + if ( all_human_stress_indices ) then + teq_ref2m_r(p) = teq_ref2m(p) + ept_ref2m_r(p) = ept_ref2m(p) + wb_ref2m_r(p) = wb_ref2m(p) + discomf_index_ref2m_r(p) = discomf_index_ref2m(p) + thic_ref2m_r(p) = thic_ref2m(p) + thip_ref2m_r(p) = thip_ref2m(p) + swmp80_ref2m_r(p) = swmp80_ref2m(p) + swmp65_ref2m_r(p) = swmp65_ref2m(p) + end if + end if ! Downward longwave radiation below the canopy diff --git a/src/biogeophys/CanopyHydrologyMod.F90 b/src/biogeophys/CanopyHydrologyMod.F90 index ecb35092c7..1ed2776aa8 100644 --- a/src/biogeophys/CanopyHydrologyMod.F90 +++ b/src/biogeophys/CanopyHydrologyMod.F90 @@ -81,7 +81,6 @@ subroutine CanopyHydrology_readnl( NLFilename ) integer :: unitn ! unit for namelist file character(len=32) :: subname = 'CanopyHydrology_readnl' ! subroutine name !----------------------------------------------------------------------- - namelist /clm_canopyhydrology_inparm/ & oldfflag, & interception_fraction, & @@ -920,4 +919,5 @@ logical function IsSnowvegFlagOnRad( ) end function IsSnowvegFlagOnRad + end module CanopyHydrologyMod diff --git a/src/biogeophys/HumanIndexMod.F90 b/src/biogeophys/HumanIndexMod.F90 index b2f3d8abff..a231c0919d 100644 --- a/src/biogeophys/HumanIndexMod.F90 +++ b/src/biogeophys/HumanIndexMod.F90 @@ -39,7 +39,13 @@ module HumanIndexMod ! ! !PUBLIC MEMBER DATA: - logical, public :: calc_human_stress_indices = .true. ! If should calculate the set of human stress indices + character(len= *), parameter, public :: calc_human_stress_indices_all = 'ALL' + character(len= *), parameter, public :: calc_human_stress_indices_fast = 'FAST' + character(len= *), parameter, public :: calc_human_stress_indices_none = 'NONE' + character(len= 16), public :: calc_human_stress_indices = calc_human_stress_indices_fast + logical, public :: all_human_stress_indices = .false. ! If should calculate the full set of human stress indices + logical, public :: fast_human_stress_indices = .true. ! If should calculate the fast (limited) set of human + ! stress indices type, public :: humanindex_type real(r8), pointer :: tc_ref2m_patch (:) ! Patch 2 m height surface air temperature (C) real(r8), pointer :: vap_ref2m_patch (:) ! Patch 2 m height vapor pressure (Pa) @@ -151,15 +157,27 @@ subroutine Init(this, bounds ) type(bounds_type) :: bounds_tmp !EOP !----------------------------------------------------------------------- - if ( calc_human_stress_indices ) then - call this%InitAllocate ( bounds ) - call this%InitHistory ( bounds ) - else + if (trim(calc_human_stress_indices) == calc_human_stress_indices_all) then + all_human_stress_indices = .true. + fast_human_stress_indices = .false. + else if (trim(calc_human_stress_indices) == calc_human_stress_indices_fast) then + all_human_stress_indices = .false. + fast_human_stress_indices = .true. + else if (trim(calc_human_stress_indices) == calc_human_stress_indices_none) then + all_human_stress_indices = .false. + fast_human_stress_indices = .false. + end if + + ! Allocation always needs to be done... + if (trim(calc_human_stress_indices) == calc_human_stress_indices_none) then ! Associate statements need humanindex_inst to be allocated ! So allocate with size 1 when not being used bounds_tmp%begp = 1 bounds_tmp%endp = 1 + call this%InitAllocate ( bounds_tmp ) + else call this%InitAllocate ( bounds ) + call this%InitHistory ( bounds ) end if end subroutine Init @@ -183,50 +201,52 @@ subroutine InitAllocate(this, bounds) begp = bounds%begp; endp= bounds%endp - allocate(this%vap_ref2m_patch (begp:endp)) ; this%vap_ref2m_patch (:) = nan - allocate(this%humidex_ref2m_u_patch (begp:endp)) ; this%humidex_ref2m_u_patch (:) = nan - allocate(this%humidex_ref2m_patch (begp:endp)) ; this%humidex_ref2m_patch (:) = nan - allocate(this%humidex_ref2m_r_patch (begp:endp)) ; this%humidex_ref2m_r_patch (:) = nan - allocate(this%nws_hi_ref2m_patch (begp:endp)) ; this%nws_hi_ref2m_patch (:) = nan - allocate(this%nws_hi_ref2m_r_patch (begp:endp)) ; this%nws_hi_ref2m_r_patch (:) = nan - allocate(this%thip_ref2m_patch (begp:endp)) ; this%thip_ref2m_patch (:) = nan - allocate(this%thip_ref2m_r_patch (begp:endp)) ; this%thip_ref2m_r_patch (:) = nan - allocate(this%thic_ref2m_patch (begp:endp)) ; this%thic_ref2m_patch (:) = nan - allocate(this%thic_ref2m_r_patch (begp:endp)) ; this%thic_ref2m_r_patch (:) = nan - allocate(this%nws_hi_ref2m_u_patch (begp:endp)) ; this%nws_hi_ref2m_u_patch (:) = nan - allocate(this%thip_ref2m_u_patch (begp:endp)) ; this%thip_ref2m_u_patch (:) = nan - allocate(this%thic_ref2m_u_patch (begp:endp)) ; this%thic_ref2m_u_patch (:) = nan - allocate(this%tc_ref2m_patch (begp:endp)) ; this%tc_ref2m_patch (:) = nan - allocate(this%appar_temp_ref2m_patch (begp:endp)) ; this%appar_temp_ref2m_patch (:) = nan - allocate(this%appar_temp_ref2m_r_patch (begp:endp)) ; this%appar_temp_ref2m_r_patch (:) = nan - allocate(this%swbgt_ref2m_patch (begp:endp)) ; this%swbgt_ref2m_patch (:) = nan - allocate(this%swbgt_ref2m_r_patch (begp:endp)) ; this%swbgt_ref2m_r_patch (:) = nan - allocate(this%wbt_ref2m_patch (begp:endp)) ; this%wbt_ref2m_patch (:) = nan - allocate(this%wbt_ref2m_r_patch (begp:endp)) ; this%wbt_ref2m_r_patch (:) = nan - allocate(this%wb_ref2m_patch (begp:endp)) ; this%wb_ref2m_patch (:) = nan - allocate(this%wb_ref2m_r_patch (begp:endp)) ; this%wb_ref2m_r_patch (:) = nan - allocate(this%teq_ref2m_patch (begp:endp)) ; this%teq_ref2m_patch (:) = nan - allocate(this%teq_ref2m_r_patch (begp:endp)) ; this%teq_ref2m_r_patch (:) = nan - allocate(this%ept_ref2m_patch (begp:endp)) ; this%ept_ref2m_patch (:) = nan - allocate(this%ept_ref2m_r_patch (begp:endp)) ; this%ept_ref2m_r_patch (:) = nan - allocate(this%discomf_index_ref2m_patch (begp:endp)) ; this%discomf_index_ref2m_patch (:) = nan - allocate(this%discomf_index_ref2m_r_patch(begp:endp)) ; this%discomf_index_ref2m_r_patch(:) = nan - allocate(this%discomf_index_ref2mS_patch(begp:endp)) ; this%discomf_index_ref2mS_patch (:) = nan - allocate(this%discomf_index_ref2mS_r_patch(begp:endp)) ; this%discomf_index_ref2mS_r_patch(:) = nan - allocate(this%discomf_index_ref2mS_u_patch(begp:endp)) ; this%discomf_index_ref2mS_u_patch(:) = nan - allocate(this%swmp65_ref2m_patch (begp:endp)) ; this%swmp65_ref2m_patch (:) = nan - allocate(this%swmp65_ref2m_r_patch (begp:endp)) ; this%swmp65_ref2m_r_patch (:) = nan - allocate(this%swmp80_ref2m_patch (begp:endp)) ; this%swmp80_ref2m_patch (:) = nan - allocate(this%swmp80_ref2m_r_patch (begp:endp)) ; this%swmp80_ref2m_r_patch (:) = nan - allocate(this%swmp80_ref2m_u_patch (begp:endp)) ; this%swmp80_ref2m_u_patch (:) = nan - allocate(this%appar_temp_ref2m_u_patch (begp:endp)) ; this%appar_temp_ref2m_u_patch (:) = nan - allocate(this%swbgt_ref2m_u_patch (begp:endp)) ; this%swbgt_ref2m_u_patch (:) = nan - allocate(this%wbt_ref2m_u_patch (begp:endp)) ; this%wbt_ref2m_u_patch (:) = nan - allocate(this%wb_ref2m_u_patch (begp:endp)) ; this%wb_ref2m_u_patch (:) = nan - allocate(this%teq_ref2m_u_patch (begp:endp)) ; this%teq_ref2m_u_patch (:) = nan - allocate(this%ept_ref2m_u_patch (begp:endp)) ; this%ept_ref2m_u_patch (:) = nan - allocate(this%discomf_index_ref2m_u_patch(begp:endp)) ; this%discomf_index_ref2m_u_patch(:) = nan - allocate(this%swmp65_ref2m_u_patch (begp:endp)) ; this%swmp65_ref2m_u_patch (:) = nan + allocate(this%vap_ref2m_patch (begp:endp)) ; this%vap_ref2m_patch (:) = nan + allocate(this%tc_ref2m_patch (begp:endp)) ; this%tc_ref2m_patch (:) = nan + allocate(this%humidex_ref2m_patch (begp:endp)) ; this%humidex_ref2m_patch (:) = nan + allocate(this%humidex_ref2m_u_patch (begp:endp)) ; this%humidex_ref2m_u_patch (:) = nan + allocate(this%humidex_ref2m_r_patch (begp:endp)) ; this%humidex_ref2m_r_patch (:) = nan + allocate(this%nws_hi_ref2m_patch (begp:endp)) ; this%nws_hi_ref2m_patch (:) = nan + allocate(this%nws_hi_ref2m_r_patch (begp:endp)) ; this%nws_hi_ref2m_r_patch (:) = nan + allocate(this%nws_hi_ref2m_u_patch (begp:endp)) ; this%nws_hi_ref2m_u_patch (:) = nan + allocate(this%appar_temp_ref2m_patch (begp:endp)) ; this%appar_temp_ref2m_patch (:) = nan + allocate(this%appar_temp_ref2m_r_patch (begp:endp)) ; this%appar_temp_ref2m_r_patch (:) = nan + allocate(this%appar_temp_ref2m_u_patch (begp:endp)) ; this%appar_temp_ref2m_u_patch (:) = nan + allocate(this%swbgt_ref2m_patch (begp:endp)) ; this%swbgt_ref2m_patch (:) = nan + allocate(this%swbgt_ref2m_r_patch (begp:endp)) ; this%swbgt_ref2m_r_patch (:) = nan + allocate(this%swbgt_ref2m_u_patch (begp:endp)) ; this%swbgt_ref2m_u_patch (:) = nan + allocate(this%wbt_ref2m_patch (begp:endp)) ; this%wbt_ref2m_patch (:) = nan + allocate(this%wbt_ref2m_r_patch (begp:endp)) ; this%wbt_ref2m_r_patch (:) = nan + allocate(this%wbt_ref2m_u_patch (begp:endp)) ; this%wbt_ref2m_u_patch (:) = nan + allocate(this%discomf_index_ref2mS_patch (begp:endp)) ; this%discomf_index_ref2mS_patch (:) = nan + allocate(this%discomf_index_ref2mS_r_patch (begp:endp)) ; this%discomf_index_ref2mS_r_patch(:) = nan + allocate(this%discomf_index_ref2mS_u_patch (begp:endp)) ; this%discomf_index_ref2mS_u_patch(:) = nan + + allocate(this%wb_ref2m_patch (begp:endp)) ; this%wb_ref2m_patch (:) = nan + allocate(this%wb_ref2m_r_patch (begp:endp)) ; this%wb_ref2m_r_patch (:) = nan + allocate(this%wb_ref2m_u_patch (begp:endp)) ; this%wb_ref2m_u_patch (:) = nan + allocate(this%teq_ref2m_patch (begp:endp)) ; this%teq_ref2m_patch (:) = nan + allocate(this%teq_ref2m_r_patch (begp:endp)) ; this%teq_ref2m_r_patch (:) = nan + allocate(this%teq_ref2m_u_patch (begp:endp)) ; this%teq_ref2m_u_patch (:) = nan + allocate(this%ept_ref2m_patch (begp:endp)) ; this%ept_ref2m_patch (:) = nan + allocate(this%ept_ref2m_r_patch (begp:endp)) ; this%ept_ref2m_r_patch (:) = nan + allocate(this%ept_ref2m_u_patch (begp:endp)) ; this%ept_ref2m_u_patch (:) = nan + allocate(this%discomf_index_ref2m_patch (begp:endp)) ; this%discomf_index_ref2m_patch (:) = nan + allocate(this%discomf_index_ref2m_r_patch (begp:endp)) ; this%discomf_index_ref2m_r_patch (:) = nan + allocate(this%discomf_index_ref2m_u_patch (begp:endp)) ; this%discomf_index_ref2m_u_patch (:) = nan + allocate(this%thip_ref2m_patch (begp:endp)) ; this%thip_ref2m_patch (:) = nan + allocate(this%thip_ref2m_r_patch (begp:endp)) ; this%thip_ref2m_r_patch (:) = nan + allocate(this%thip_ref2m_u_patch (begp:endp)) ; this%thip_ref2m_u_patch (:) = nan + allocate(this%thic_ref2m_patch (begp:endp)) ; this%thic_ref2m_patch (:) = nan + allocate(this%thic_ref2m_r_patch (begp:endp)) ; this%thic_ref2m_r_patch (:) = nan + allocate(this%thic_ref2m_u_patch (begp:endp)) ; this%thic_ref2m_u_patch (:) = nan + allocate(this%swmp65_ref2m_patch (begp:endp)) ; this%swmp65_ref2m_patch (:) = nan + allocate(this%swmp65_ref2m_r_patch (begp:endp)) ; this%swmp65_ref2m_r_patch (:) = nan + allocate(this%swmp65_ref2m_u_patch (begp:endp)) ; this%swmp65_ref2m_u_patch (:) = nan + allocate(this%swmp80_ref2m_patch (begp:endp)) ; this%swmp80_ref2m_patch (:) = nan + allocate(this%swmp80_ref2m_r_patch (begp:endp)) ; this%swmp80_ref2m_r_patch (:) = nan + allocate(this%swmp80_ref2m_u_patch (begp:endp)) ; this%swmp80_ref2m_u_patch (:) = nan + end subroutine InitAllocate @@ -250,21 +270,6 @@ subroutine InitHistory(this, bounds) begp = bounds%begp; endp= bounds%endp - this%appar_temp_ref2m_patch(begp:endp) = spval - call hist_addfld1d (fname='APPAR_TEMP', units='C', & - avgflag='A', long_name='2 m apparent temperature', & - ptr_patch=this%appar_temp_ref2m_patch, default='inactive') - - this%appar_temp_ref2m_u_patch(begp:endp) = spval - call hist_addfld1d (fname='APPAR_TEMP_U', units='C', & - avgflag='A', long_name='Urban 2 m apparent temperature', & - ptr_patch=this%appar_temp_ref2m_u_patch, set_nourb=spval, default='inactive') - - this%appar_temp_ref2m_r_patch(begp:endp) = spval - call hist_addfld1d (fname='APPAR_TEMP_R', units='C', & - avgflag='A', long_name='Rural 2 m apparent temperature', & - ptr_patch=this%appar_temp_ref2m_r_patch, set_spec=spval, default='inactive') - this%swbgt_ref2m_patch(begp:endp) = spval call hist_addfld1d (fname='SWBGT', units='C', & avgflag='A', long_name='2 m Simplified Wetbulb Globe Temp', & @@ -310,81 +315,6 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='Rural 2 m Stull Wet Bulb', & ptr_patch=this%wbt_ref2m_r_patch, set_spec=spval) - this%wb_ref2m_patch(begp:endp) = spval - call hist_addfld1d (fname='WBA', units='C', & - avgflag='A', long_name='2 m Wet Bulb', & - ptr_patch=this%wb_ref2m_patch, default='inactive') - - this%wb_ref2m_u_patch(begp:endp) = spval - call hist_addfld1d (fname='WBA_U', units='C', & - avgflag='A', long_name='Urban 2 m Wet Bulb', & - ptr_patch=this%wb_ref2m_u_patch, set_nourb=spval, default='inactive') - - this%wb_ref2m_r_patch(begp:endp) = spval - call hist_addfld1d (fname='WBA_R', units='C', & - avgflag='A', long_name='Rural 2 m Wet Bulb', & - ptr_patch=this%wb_ref2m_r_patch, set_spec=spval, default='inactive') - - this%teq_ref2m_patch(begp:endp) = spval - call hist_addfld1d (fname='TEQ', units='K', & - avgflag='A', long_name='2 m Equiv Temp', & - ptr_patch=this%teq_ref2m_patch, default='inactive') - - this%teq_ref2m_u_patch(begp:endp) = spval - call hist_addfld1d (fname='TEQ_U', units='K', & - avgflag='A', long_name='Urban 2 m Equiv Temp', & - ptr_patch=this%teq_ref2m_u_patch, set_nourb=spval, default='inactive') - - this%teq_ref2m_r_patch(begp:endp) = spval - call hist_addfld1d (fname='TEQ_R', units='K', & - avgflag='A', long_name='Rural 2 m Equiv Temp', & - ptr_patch=this%teq_ref2m_r_patch, set_spec=spval, default='inactive') - - this%ept_ref2m_patch(begp:endp) = spval - call hist_addfld1d (fname='EPT', units='K', & - avgflag='A', long_name='2 m Equiv Pot Temp', & - ptr_patch=this%ept_ref2m_patch, default='inactive') - - this%ept_ref2m_u_patch(begp:endp) = spval - call hist_addfld1d (fname='EPT_U', units='K', & - avgflag='A', long_name='Urban 2 m Equiv Pot Temp', & - ptr_patch=this%ept_ref2m_u_patch, set_nourb=spval, default='inactive') - - this%ept_ref2m_r_patch(begp:endp) = spval - call hist_addfld1d (fname='EPT_R', units='K', & - avgflag='A', long_name='Rural 2 m Equiv Pot Temp', & - ptr_patch=this%ept_ref2m_r_patch, set_spec=spval, default='inactive') - - this%discomf_index_ref2m_patch(begp:endp) = spval - call hist_addfld1d (fname='DISCOI', units='C', & - avgflag='A', long_name='2 m Discomfort Index', & - ptr_patch=this%discomf_index_ref2m_patch, default='inactive') - - this%discomf_index_ref2m_u_patch(begp:endp) = spval - call hist_addfld1d (fname='DISCOI_U', units='C', & - avgflag='A', long_name='Urban 2 m Discomfort Index', & - ptr_patch=this%discomf_index_ref2m_u_patch, set_nourb=spval, default='inactive') - - this%discomf_index_ref2m_r_patch(begp:endp) = spval - call hist_addfld1d (fname='DISCOI_R', units='C', & - avgflag='A', long_name='Rural 2 m Discomfort Index', & - ptr_patch=this%discomf_index_ref2m_r_patch, set_spec=spval, default='inactive') - - this%discomf_index_ref2mS_patch(begp:endp) = spval - call hist_addfld1d (fname='DISCOIS', units='C', & - avgflag='A', long_name='2 m Stull Discomfort Index', & - ptr_patch=this%discomf_index_ref2mS_patch, default='inactive') - - this%discomf_index_ref2mS_u_patch(begp:endp) = spval - call hist_addfld1d (fname='DISCOIS_U', units='C', & - avgflag='A', long_name='Urban 2 m Stull Discomfort Index', & - ptr_patch=this%discomf_index_ref2mS_u_patch, set_nourb=spval, default='inactive') - - this%discomf_index_ref2mS_r_patch(begp:endp) = spval - call hist_addfld1d (fname='DISCOIS_R', units='C', & - avgflag='A', long_name='Rural 2 m Stull Discomfort Index', & - ptr_patch=this%discomf_index_ref2mS_r_patch, set_spec=spval, default='inactive') - this%nws_hi_ref2m_patch(begp:endp) = spval call hist_addfld1d (fname='HIA', units='C', & avgflag='A', long_name='2 m NWS Heat Index', & @@ -400,65 +330,159 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='Rural 2 m NWS Heat Index', & ptr_patch=this%nws_hi_ref2m_r_patch, set_spec=spval) - this%thip_ref2m_patch(begp:endp) = spval - call hist_addfld1d (fname='THIP', units='C', & - avgflag='A', long_name='2 m Temp Hum Index Physiology', & - ptr_patch=this%thip_ref2m_patch, default='inactive') - - this%thip_ref2m_u_patch(begp:endp) = spval - call hist_addfld1d (fname='THIP_U', units='C', & - avgflag='A', long_name='Urban 2 m Temp Hum Index Physiology', & - ptr_patch=this%thip_ref2m_u_patch, set_nourb=spval, default='inactive') - - this%thip_ref2m_r_patch(begp:endp) = spval - call hist_addfld1d (fname='THIP_R', units='C', & - avgflag='A', long_name='Rural 2 m Temp Hum Index Physiology', & - ptr_patch=this%thip_ref2m_r_patch, set_spec=spval, default='inactive') - - this%thic_ref2m_patch(begp:endp) = spval - call hist_addfld1d (fname='THIC', units='C', & - avgflag='A', long_name='2 m Temp Hum Index Comfort', & - ptr_patch=this%thic_ref2m_patch, default='inactive') - - this%thic_ref2m_u_patch(begp:endp) = spval - call hist_addfld1d (fname='THIC_U', units='C', & - avgflag='A', long_name='Urban 2 m Temp Hum Index Comfort', & - ptr_patch=this%thic_ref2m_u_patch, set_nourb=spval, default='inactive') - - this%thic_ref2m_r_patch(begp:endp) = spval - call hist_addfld1d (fname='THIC_R', units='C', & - avgflag='A', long_name='Rural 2 m Temp Hum Index Comfort', & - ptr_patch=this%thic_ref2m_r_patch, set_spec=spval, default='inactive') - - this%swmp65_ref2m_patch(begp:endp) = spval - call hist_addfld1d (fname='SWMP65', units='C', & - avgflag='A', long_name='2 m Swamp Cooler Temp 65% Eff', & - ptr_patch=this%swmp65_ref2m_patch, default='inactive') - - this%swmp65_ref2m_u_patch(begp:endp) = spval - call hist_addfld1d (fname='SWMP65_U', units='C', & - avgflag='A', long_name='Urban 2 m Swamp Cooler Temp 65% Eff', & - ptr_patch=this%swmp65_ref2m_u_patch, set_nourb=spval, default='inactive') - - this%swmp65_ref2m_r_patch(begp:endp) = spval - call hist_addfld1d (fname='SWMP65_R', units='C', & - avgflag='A', long_name='Rural 2 m Swamp Cooler Temp 65% Eff', & - ptr_patch=this%swmp65_ref2m_r_patch, set_spec=spval, default='inactive') - - this%swmp80_ref2m_patch(begp:endp) = spval - call hist_addfld1d (fname='SWMP80', units='C', & - avgflag='A', long_name='2 m Swamp Cooler Temp 80% Eff', & - ptr_patch=this%swmp80_ref2m_patch, default='inactive') - - this%swmp80_ref2m_u_patch(begp:endp) = spval - call hist_addfld1d (fname='SWMP80_U', units='C', & - avgflag='A', long_name='Urban 2 m Swamp Cooler Temp 80% Eff', & - ptr_patch=this%swmp80_ref2m_u_patch, set_nourb=spval, default='inactive') - - this%swmp80_ref2m_r_patch(begp:endp) = spval - call hist_addfld1d (fname='SWMP80_R', units='C', & - avgflag='A', long_name='Rural 2 m Swamp Cooler Temp 80% Eff', & - ptr_patch=this%swmp80_ref2m_r_patch, set_spec=spval, default='inactive') + if ( all_human_stress_indices )then + + this%appar_temp_ref2m_patch(begp:endp) = spval + call hist_addfld1d (fname='APPAR_TEMP', units='C', & + avgflag='A', long_name='2 m apparent temperature', & + ptr_patch=this%appar_temp_ref2m_patch) + + this%appar_temp_ref2m_u_patch(begp:endp) = spval + call hist_addfld1d (fname='APPAR_TEMP_U', units='C', & + avgflag='A', long_name='Urban 2 m apparent temperature', & + ptr_patch=this%appar_temp_ref2m_u_patch, set_nourb=spval) + + this%appar_temp_ref2m_r_patch(begp:endp) = spval + call hist_addfld1d (fname='APPAR_TEMP_R', units='C', & + avgflag='A', long_name='Rural 2 m apparent temperature', & + ptr_patch=this%appar_temp_ref2m_r_patch, set_spec=spval) + + this%wb_ref2m_patch(begp:endp) = spval + call hist_addfld1d (fname='WBA', units='C', & + avgflag='A', long_name='2 m Wet Bulb', & + ptr_patch=this%wb_ref2m_patch) + + this%wb_ref2m_u_patch(begp:endp) = spval + call hist_addfld1d (fname='WBA_U', units='C', & + avgflag='A', long_name='Urban 2 m Wet Bulb', & + ptr_patch=this%wb_ref2m_u_patch, set_nourb=spval) + + this%wb_ref2m_r_patch(begp:endp) = spval + call hist_addfld1d (fname='WBA_R', units='C', & + avgflag='A', long_name='Rural 2 m Wet Bulb', & + ptr_patch=this%wb_ref2m_r_patch, set_spec=spval) + + this%teq_ref2m_patch(begp:endp) = spval + call hist_addfld1d (fname='TEQ', units='K', & + avgflag='A', long_name='2 m Equiv Temp', & + ptr_patch=this%teq_ref2m_patch) + + this%teq_ref2m_u_patch(begp:endp) = spval + call hist_addfld1d (fname='TEQ_U', units='K', & + avgflag='A', long_name='Urban 2 m Equiv Temp', & + ptr_patch=this%teq_ref2m_u_patch, set_nourb=spval) + + this%teq_ref2m_r_patch(begp:endp) = spval + call hist_addfld1d (fname='TEQ_R', units='K', & + avgflag='A', long_name='Rural 2 m Equiv Temp', & + ptr_patch=this%teq_ref2m_r_patch, set_spec=spval) + + this%ept_ref2m_patch(begp:endp) = spval + call hist_addfld1d (fname='EPT', units='K', & + avgflag='A', long_name='2 m Equiv Pot Temp', & + ptr_patch=this%ept_ref2m_patch) + + this%ept_ref2m_u_patch(begp:endp) = spval + call hist_addfld1d (fname='EPT_U', units='K', & + avgflag='A', long_name='Urban 2 m Equiv Pot Temp', & + ptr_patch=this%ept_ref2m_u_patch, set_nourb=spval) + + this%ept_ref2m_r_patch(begp:endp) = spval + call hist_addfld1d (fname='EPT_R', units='K', & + avgflag='A', long_name='Rural 2 m Equiv Pot Temp', & + ptr_patch=this%ept_ref2m_r_patch, set_spec=spval) + + this%discomf_index_ref2m_patch(begp:endp) = spval + call hist_addfld1d (fname='DISCOI', units='C', & + avgflag='A', long_name='2 m Discomfort Index', & + ptr_patch=this%discomf_index_ref2m_patch) + + this%discomf_index_ref2m_u_patch(begp:endp) = spval + call hist_addfld1d (fname='DISCOI_U', units='C', & + avgflag='A', long_name='Urban 2 m Discomfort Index', & + ptr_patch=this%discomf_index_ref2m_u_patch, set_nourb=spval) + + this%discomf_index_ref2m_r_patch(begp:endp) = spval + call hist_addfld1d (fname='DISCOI_R', units='C', & + avgflag='A', long_name='Rural 2 m Discomfort Index', & + ptr_patch=this%discomf_index_ref2m_r_patch, set_spec=spval) + + this%thip_ref2m_patch(begp:endp) = spval + call hist_addfld1d (fname='THIP', units='C', & + avgflag='A', long_name='2 m Temp Hum Index Physiology', & + ptr_patch=this%thip_ref2m_patch) + + this%thip_ref2m_u_patch(begp:endp) = spval + call hist_addfld1d (fname='THIP_U', units='C', & + avgflag='A', long_name='Urban 2 m Temp Hum Index Physiology', & + ptr_patch=this%thip_ref2m_u_patch, set_nourb=spval) + + this%thip_ref2m_r_patch(begp:endp) = spval + call hist_addfld1d (fname='THIP_R', units='C', & + avgflag='A', long_name='Rural 2 m Temp Hum Index Physiology', & + ptr_patch=this%thip_ref2m_r_patch, set_spec=spval) + + this%thic_ref2m_patch(begp:endp) = spval + call hist_addfld1d (fname='THIC', units='C', & + avgflag='A', long_name='2 m Temp Hum Index Comfort', & + ptr_patch=this%thic_ref2m_patch) + + this%thic_ref2m_u_patch(begp:endp) = spval + call hist_addfld1d (fname='THIC_U', units='C', & + avgflag='A', long_name='Urban 2 m Temp Hum Index Comfort', & + ptr_patch=this%thic_ref2m_u_patch, set_nourb=spval) + + this%thic_ref2m_r_patch(begp:endp) = spval + call hist_addfld1d (fname='THIC_R', units='C', & + avgflag='A', long_name='Rural 2 m Temp Hum Index Comfort', & + ptr_patch=this%thic_ref2m_r_patch, set_spec=spval) + + this%swmp65_ref2m_patch(begp:endp) = spval + call hist_addfld1d (fname='SWMP65', units='C', & + avgflag='A', long_name='2 m Swamp Cooler Temp 65% Eff', & + ptr_patch=this%swmp65_ref2m_patch) + + this%swmp65_ref2m_u_patch(begp:endp) = spval + call hist_addfld1d (fname='SWMP65_U', units='C', & + avgflag='A', long_name='Urban 2 m Swamp Cooler Temp 65% Eff', & + ptr_patch=this%swmp65_ref2m_u_patch, set_nourb=spval) + + this%swmp65_ref2m_r_patch(begp:endp) = spval + call hist_addfld1d (fname='SWMP65_R', units='C', & + avgflag='A', long_name='Rural 2 m Swamp Cooler Temp 65% Eff', & + ptr_patch=this%swmp65_ref2m_r_patch, set_spec=spval) + + this%swmp80_ref2m_patch(begp:endp) = spval + call hist_addfld1d (fname='SWMP80', units='C', & + avgflag='A', long_name='2 m Swamp Cooler Temp 80% Eff', & + ptr_patch=this%swmp80_ref2m_patch) + + this%swmp80_ref2m_u_patch(begp:endp) = spval + call hist_addfld1d (fname='SWMP80_U', units='C', & + avgflag='A', long_name='Urban 2 m Swamp Cooler Temp 80% Eff', & + ptr_patch=this%swmp80_ref2m_u_patch, set_nourb=spval) + + this%swmp80_ref2m_r_patch(begp:endp) = spval + call hist_addfld1d (fname='SWMP80_R', units='C', & + avgflag='A', long_name='Rural 2 m Swamp Cooler Temp 80% Eff', & + ptr_patch=this%swmp80_ref2m_r_patch, set_spec=spval) + + this%discomf_index_ref2mS_patch(begp:endp) = spval + call hist_addfld1d (fname='DISCOIS', units='C', & + avgflag='A', long_name='2 m Stull Discomfort Index', & + ptr_patch=this%discomf_index_ref2mS_patch) + + this%discomf_index_ref2mS_u_patch(begp:endp) = spval + call hist_addfld1d (fname='DISCOIS_U', units='C', & + avgflag='A', long_name='Urban 2 m Stull Discomfort Index', & + ptr_patch=this%discomf_index_ref2mS_u_patch, set_nourb=spval) + + this%discomf_index_ref2mS_r_patch(begp:endp) = spval + call hist_addfld1d (fname='DISCOIS_R', units='C', & + avgflag='A', long_name='Rural 2 m Stull Discomfort Index', & + ptr_patch=this%discomf_index_ref2mS_r_patch, set_spec=spval) + + end if end subroutine InitHistory diff --git a/src/biogeophys/IrrigationMod.F90 b/src/biogeophys/IrrigationMod.F90 index cccb8dbc6b..ef12c1d063 100644 --- a/src/biogeophys/IrrigationMod.F90 +++ b/src/biogeophys/IrrigationMod.F90 @@ -51,7 +51,7 @@ module IrrigationMod use shr_log_mod , only : errMsg => shr_log_errMsg use abortutils , only : endrun use clm_varctl , only : iulog - use clm_varcon , only : isecspday, degpsec, denh2o, spval, namec + use clm_varcon , only : isecspday, denh2o, spval, namec use clm_varpar , only : nlevsoi, nlevgrnd use clm_time_manager , only : get_step_size use SoilWaterRetentionCurveMod, only : soil_water_retention_curve_type @@ -799,7 +799,7 @@ end subroutine ApplyIrrigation !----------------------------------------------------------------------- subroutine CalcIrrigationNeeded(this, bounds, num_exposedvegp, filter_exposedvegp, & - time_prev, elai, t_soisno, eff_porosity, h2osoi_liq, volr, rof_prognostic) + elai, t_soisno, eff_porosity, h2osoi_liq, volr, rof_prognostic) ! ! !DESCRIPTION: ! Calculate whether and how much irrigation is needed for each column. However, this @@ -812,9 +812,6 @@ subroutine CalcIrrigationNeeded(this, bounds, num_exposedvegp, filter_exposedveg class(irrigation_type) , intent(inout) :: this type(bounds_type) , intent(in) :: bounds - ! time of day (in seconds since 0Z) at start of timestep - integer, intent(in) :: time_prev - ! number of points in filter_exposedvegp integer, intent(in) :: num_exposedvegp @@ -911,7 +908,7 @@ subroutine CalcIrrigationNeeded(this, bounds, num_exposedvegp, filter_exposedveg check_for_irrig_patch(p) = this%PointNeedsCheckForIrrig( & pft_type=patch%itype(p), elai=elai(p), & - time_prev=time_prev, londeg=grc%londeg(g)) + londeg=grc%londeg(g)) if (check_for_irrig_patch(p)) then c = patch%column(p) check_for_irrig_col(c) = .true. @@ -1033,27 +1030,24 @@ subroutine CalcIrrigationNeeded(this, bounds, num_exposedvegp, filter_exposedveg end subroutine CalcIrrigationNeeded !----------------------------------------------------------------------- - function PointNeedsCheckForIrrig(this, pft_type, elai, time_prev, londeg) & + function PointNeedsCheckForIrrig(this, pft_type, elai, londeg) & result(check_for_irrig) ! ! !DESCRIPTION: ! Determine whether a given patch needs to be checked for irrigation now. ! ! !USES: - use pftconMod, only : pftcon + use clm_time_manager, only : get_local_time + use pftconMod , only : pftcon ! ! !ARGUMENTS: logical :: check_for_irrig ! function result class(irrigation_type), intent(in) :: this integer , intent(in) :: pft_type ! type of pft in this patch real(r8), intent(in) :: elai ! one-sided leaf area index with burying by snow - integer , intent(in) :: time_prev ! time of day (in seconds since 0Z) at start of timestep real(r8), intent(in) :: londeg ! longitude (degrees) ! ! !LOCAL VARIABLES: - ! local time at start of time step (seconds after solar midnight) - integer :: local_time - ! number of seconds since the prescribed irrigation start time integer :: seconds_since_irrig_start_time @@ -1063,8 +1057,7 @@ function PointNeedsCheckForIrrig(this, pft_type, elai, time_prev, londeg) & if (pftcon%irrigated(pft_type) == 1._r8 .and. & elai > this%params%irrig_min_lai) then ! see if it's the right time of day to start irrigating: - local_time = modulo(time_prev + nint(londeg/degpsec), isecspday) - seconds_since_irrig_start_time = modulo(local_time - this%params%irrig_start_time, isecspday) + seconds_since_irrig_start_time = get_local_time( londeg, starttime=this%params%irrig_start_time, offset=-this%dtime ) if (seconds_since_irrig_start_time < this%dtime) then check_for_irrig = .true. else diff --git a/src/biogeophys/LakeFluxesMod.F90 b/src/biogeophys/LakeFluxesMod.F90 index 3d447ec111..2c54cf3b69 100644 --- a/src/biogeophys/LakeFluxesMod.F90 +++ b/src/biogeophys/LakeFluxesMod.F90 @@ -55,7 +55,8 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep, use LakeCon , only : minz0lake, cur0, cus, curm, fcrit use QSatMod , only : QSat use FrictionVelocityMod , only : FrictionVelocity, MoninObukIni, frictionvel_parms_inst - use HumanIndexMod , only : calc_human_stress_indices, Wet_Bulb, Wet_BulbS, HeatIndex, AppTemp, & + use HumanIndexMod , only : all_human_stress_indices, fast_human_stress_indices, & + Wet_Bulb, Wet_BulbS, HeatIndex, AppTemp, & swbgt, hmdex, dis_coi, dis_coiS, THIndex, & SwampCoolEff, KtoC, VaporPres ! @@ -608,20 +609,22 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep, ! Human Heat Stress - if ( calc_human_stress_indices )then + if ( all_human_stress_indices .or. fast_human_stress_indices )then call KtoC(t_ref2m(p), tc_ref2m(p)) call VaporPres(rh_ref2m(p), e_ref2m, vap_ref2m(p)) - call Wet_Bulb(t_ref2m(p), vap_ref2m(p), forc_pbot(c), rh_ref2m(p), & - q_ref2m(p), teq_ref2m(p), ept_ref2m(p), wb_ref2m(p)) call Wet_BulbS(tc_ref2m(p), rh_ref2m(p), wbt_ref2m(p)) call HeatIndex(tc_ref2m(p), rh_ref2m(p), nws_hi_ref2m(p)) call AppTemp(tc_ref2m(p), vap_ref2m(p), u10_clm(p), appar_temp_ref2m(p)) call swbgt(tc_ref2m(p), vap_ref2m(p), swbgt_ref2m(p)) call hmdex(tc_ref2m(p), vap_ref2m(p), humidex_ref2m(p)) - call dis_coi(tc_ref2m(p), wb_ref2m(p), discomf_index_ref2m(p)) call dis_coiS(tc_ref2m(p), rh_ref2m(p), wbt_ref2m(p), discomf_index_ref2mS(p)) - call THIndex(tc_ref2m(p), wb_ref2m(p), thic_ref2m(p), thip_ref2m(p)) - call SwampCoolEff(tc_ref2m(p), wb_ref2m(p), swmp80_ref2m(p), swmp65_ref2m(p)) + if ( all_human_stress_indices ) then + call Wet_Bulb(t_ref2m(p), vap_ref2m(p), forc_pbot(c), rh_ref2m(p), & + q_ref2m(p), teq_ref2m(p), ept_ref2m(p), wb_ref2m(p)) + call dis_coi(tc_ref2m(p), wb_ref2m(p), discomf_index_ref2m(p)) + call THIndex(tc_ref2m(p), wb_ref2m(p), thic_ref2m(p), thip_ref2m(p)) + call SwampCoolEff(tc_ref2m(p), wb_ref2m(p), swmp80_ref2m(p), swmp65_ref2m(p)) + end if end if diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90 index e735d74ce0..3c6cef748d 100644 --- a/src/biogeophys/PhotosynthesisMod.F90 +++ b/src/biogeophys/PhotosynthesisMod.F90 @@ -955,9 +955,9 @@ subroutine Photosynthesis ( bounds, fn, filterp, & ! a multi-layer canopy ! ! !USES: - use clm_varcon , only : rgas, tfrz, spval, degpsec, isecspday + use clm_varcon , only : rgas, tfrz, spval use GridcellType , only : grc - use clm_time_manager , only : get_curr_date, get_step_size + use clm_time_manager , only : get_step_size, is_near_local_noon use clm_varctl , only : cnallocate_carbon_only use clm_varctl , only : lnc_opt, reduce_dayl_factor, vcmax_opt use pftconMod , only : nbrdlf_dcd_tmp_shrub, npcropmin @@ -1101,11 +1101,8 @@ subroutine Photosynthesis ( bounds, fn, filterp, & real(r8) :: total_lai integer :: nptreemax - integer :: local_secp1 ! seconds into current date in local time real(r8) :: dtime ! land model time step (sec) - integer :: year,month,day,secs ! calendar info for current time step integer :: g ! index - integer, parameter :: noonsec = isecspday / 2 ! seconds at local noon !------------------------------------------------------------------------------ ! Temperature and soil water response functions @@ -1220,7 +1217,6 @@ subroutine Photosynthesis ( bounds, fn, filterp, & ! Determine seconds of current time step dtime = get_step_size() - call get_curr_date (year, month, day, secs) ! vcmax25 parameters, from CN @@ -1654,12 +1650,8 @@ subroutine Photosynthesis ( bounds, fn, filterp, & if (an(p,iv) < 0._r8) gs_mol(p,iv) = bbb(p) - ! Get local noon sunlit and shaded stomatal conductance - local_secp1 = secs + nint((grc%londeg(g)/degpsec)/dtime)*dtime - local_secp1 = mod(local_secp1,isecspday) - ! Use time period 1 hour before and 1 hour after local noon inclusive (11AM-1PM) - if (local_secp1 >= (isecspday/2 - 3600) .and. local_secp1 <= (isecspday/2 + 3600)) then + if ( is_near_local_noon( grc%londeg(g), deltasec=3600 ) )then if (phase == 'sun') then gs_mol_sun_ln(p,iv) = gs_mol(p,iv) else if (phase == 'sha') then @@ -2439,9 +2431,9 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, & ! method ! ! !USES: - use clm_varcon , only : rgas, tfrz, rpi, spval, degpsec, isecspday + use clm_varcon , only : rgas, tfrz, rpi, spval use GridcellType , only : grc - use clm_time_manager , only : get_curr_date, get_step_size + use clm_time_manager , only : get_step_size, is_near_local_noon use clm_varctl , only : cnallocate_carbon_only use clm_varctl , only : lnc_opt, reduce_dayl_factor, vcmax_opt use clm_varpar , only : nlevsoi @@ -2631,9 +2623,7 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, & real(r8) :: sum_nscaler real(r8) :: total_lai integer :: nptreemax - integer :: local_secp1 ! seconds into current date in local time real(r8) :: dtime ! land model time step (sec) - integer :: year,month,day,secs ! calendar info for current time step integer :: j,g ! index real(r8) :: rs_resis ! combined soil-root resistance [s] real(r8) :: r_soil ! root spacing [m] @@ -2653,7 +2643,6 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, & real(r8), parameter :: croot_lateral_length = 0.25_r8 ! specified lateral coarse root length [m] real(r8), parameter :: c_to_b = 2.0_r8 !(g biomass /g C) !Note that root density is for dry biomass not carbon. CLM provides root biomass as carbon. The conversion is 0.5 g C / g biomass - integer, parameter :: noonsec = isecspday / 2 ! seconds at local noon !------------------------------------------------------------------------------ @@ -2791,7 +2780,6 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, & ! Determine seconds off current time step dtime = get_step_size() - call get_curr_date (year, month, day, secs) ! vcmax25 parameters, from CN @@ -3356,11 +3344,8 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, & if (an_sun(p,iv) < 0._r8) gs_mol_sun(p,iv) = max( bsun(p)*gsminsun, 1._r8 ) if (an_sha(p,iv) < 0._r8) gs_mol_sha(p,iv) = max( bsha(p)*gsminsha, 1._r8 ) - ! Get local noon sunlit and shaded stomatal conductance - local_secp1 = secs + nint((grc%londeg(g)/degpsec)/dtime)*dtime - local_secp1 = mod(local_secp1,isecspday) ! Use time period 1 hour before and 1 hour after local noon inclusive (11AM-1PM) - if (local_secp1 >= (isecspday/2 - 3600) .and. local_secp1 <= (isecspday/2 + 3600)) then + if ( is_near_local_noon( grc%londeg(g), deltasec=3600 ) )then gs_mol_sun_ln(p,iv) = gs_mol_sun(p,iv) gs_mol_sha_ln(p,iv) = gs_mol_sha(p,iv) else diff --git a/src/biogeophys/SolarAbsorbedType.F90 b/src/biogeophys/SolarAbsorbedType.F90 index cff8cf9956..2a9c9a7543 100644 --- a/src/biogeophys/SolarAbsorbedType.F90 +++ b/src/biogeophys/SolarAbsorbedType.F90 @@ -17,8 +17,9 @@ module SolarAbsorbedType type, public :: solarabs_type ! Solar reflected - real(r8), pointer :: fsr_patch (:) ! patch solar radiation reflected (W/m**2) - + real(r8), pointer :: fsr_patch (:) ! patch solar radiation reflected (W/m**2) + real(r8), pointer :: fsrSF_patch (:) ! diagnostic snow-free patch solar radiation reflected (W/m**2) + real(r8), pointer :: ssre_fsr_patch (:) ! snow radiative effect on patch solar radiation reflected (W/m**2) ! Solar Absorbed real(r8), pointer :: fsa_patch (:) ! patch solar radiation absorbed (total) (W/m**2) real(r8), pointer :: fsa_u_patch (:) ! patch urban solar radiation absorbed (total) (W/m**2) @@ -53,12 +54,18 @@ module SolarAbsorbedType ! Currently needed by lake code ! TODO (MV 8/20/2014) should be moved in the future real(r8), pointer :: fsds_nir_d_patch (:) ! patch incident direct beam nir solar radiation (W/m**2) - real(r8), pointer :: fsds_nir_i_patch (:) ! patch incident diffuse nir solar radiation (W/m**2) + real(r8), pointer :: fsds_nir_i_patch (:) ! patch incident diffuse nir solar radiation (W/m**2) real(r8), pointer :: fsds_nir_d_ln_patch (:) ! patch incident direct beam nir solar radiation at local noon (W/m**2) - real(r8), pointer :: fsr_nir_d_patch (:) ! patch reflected direct beam nir solar radiation (W/m**2) - real(r8), pointer :: fsr_nir_i_patch (:) ! patch reflected diffuse nir solar radiation (W/m**2) + real(r8), pointer :: fsr_nir_d_patch (:) ! patch reflected direct beam nir solar radiation (W/m**2) + real(r8), pointer :: fsr_nir_i_patch (:) ! patch reflected diffuse nir solar radiation (W/m**2) real(r8), pointer :: fsr_nir_d_ln_patch (:) ! patch reflected direct beam nir solar radiation at local noon (W/m**2) - + ! optional diagnostic fluxes: + real(r8), pointer :: fsrSF_nir_d_patch (:) ! snow-free patch reflected direct beam nir solar radiation (W/m**2) + real(r8), pointer :: fsrSF_nir_i_patch (:) ! snow-free patch reflected diffuse nir solar radiation (W/m**2) + real(r8), pointer :: fsrSF_nir_d_ln_patch (:) ! snow-free patch reflected direct beam nir solar radiation at local noon (W/m**2) + real(r8), pointer :: ssre_fsr_nir_d_patch (:) ! snow-free patch reflected direct beam nir solar radiation (W/m**2) + real(r8), pointer :: ssre_fsr_nir_i_patch (:) ! snow-free patch reflected diffuse nir solar radiation (W/m**2) + real(r8), pointer :: ssre_fsr_nir_d_ln_patch(:) ! snow-free patch reflected direct beam nir solar radiation at local noon (W/m**2) contains procedure, public :: Init @@ -140,6 +147,14 @@ subroutine InitAllocate(this, bounds) allocate(this%fsr_nir_d_patch (begp:endp)) ; this%fsr_nir_d_patch (:) = nan allocate(this%fsr_nir_i_patch (begp:endp)) ; this%fsr_nir_i_patch (:) = nan allocate(this%fsr_nir_d_ln_patch (begp:endp)) ; this%fsr_nir_d_ln_patch (:) = nan + allocate(this%fsrSF_patch (begp:endp)) ; this%fsrSF_patch (:) = nan + allocate(this%fsrSF_nir_d_patch (begp:endp)) ; this%fsrSF_nir_d_patch (:) = nan + allocate(this%fsrSF_nir_i_patch (begp:endp)) ; this%fsrSF_nir_i_patch (:) = nan + allocate(this%fsrSF_nir_d_ln_patch (begp:endp)) ; this%fsrSF_nir_d_ln_patch (:) = nan + allocate(this%ssre_fsr_patch (begp:endp)) ; this%ssre_fsr_patch (:) = nan + allocate(this%ssre_fsr_nir_d_patch (begp:endp)) ; this%ssre_fsr_nir_d_patch (:) = nan + allocate(this%ssre_fsr_nir_i_patch (begp:endp)) ; this%ssre_fsr_nir_i_patch (:) = nan + allocate(this%ssre_fsr_nir_d_ln_patch(begp:endp)) ; this%ssre_fsr_nir_d_ln_patch(:) = nan allocate(this%fsds_nir_d_patch (begp:endp)) ; this%fsds_nir_d_patch (:) = nan allocate(this%fsds_nir_i_patch (begp:endp)) ; this%fsds_nir_i_patch (:) = nan allocate(this%fsds_nir_d_ln_patch (begp:endp)) ; this%fsds_nir_d_ln_patch (:) = nan @@ -153,7 +168,7 @@ subroutine InitHistory(this, bounds) ! ! !USES: use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) - use clm_varctl , only : use_snicar_frc + use clm_varctl , only : use_snicar_frc , use_SSRE use clm_varpar , only : nlevsno use histFileMod , only : hist_addfld1d, hist_addfld2d use histFileMod , only : no_snow_normal @@ -201,7 +216,6 @@ subroutine InitHistory(this, bounds) call hist_addfld1d (fname='SWup', units='W/m^2', & avgflag='A', long_name='upwelling shortwave radiation', & ptr_patch=this%fsr_patch, c2l_scale_type='urbanf', default='inactive') - call hist_addfld1d (fname='FSR_ICE', units='W/m^2', & avgflag='A', long_name='reflected solar radiation (ice landunits only)', & ptr_patch=this%fsr_patch, c2l_scale_type='urbanf', l2g_scale_type='ice', & @@ -263,7 +277,47 @@ subroutine InitHistory(this, bounds) call hist_addfld1d (fname='FSRNDLN', units='W/m^2', & avgflag='A', long_name='direct nir reflected solar radiation at local noon', & ptr_patch=this%fsr_nir_d_ln_patch, c2l_scale_type='urbanf') - + ! diagnostic fluxes for ESM-SnowMIP + if (use_SSRE) then + this%fsrSF_patch(begp:endp) = spval + call hist_addfld1d (fname='FSRSF', units='W/m^2', & + avgflag='A', long_name='reflected solar radiation', & + ptr_patch=this%fsrSF_patch, c2l_scale_type='urbanf') + + this%ssre_fsr_patch(begp:endp) = spval + call hist_addfld1d (fname='SSRE_FSR', units='W/m^2', & + avgflag='A', long_name='surface snow effect on reflected solar radiation', & + ptr_patch=this%ssre_fsr_patch, c2l_scale_type='urbanf') + this%fsrSF_nir_d_patch(begp:endp) = spval + call hist_addfld1d (fname='FSRSFND', units='W/m^2', & + avgflag='A', long_name='direct nir reflected solar radiation', & + ptr_patch=this%fsrSF_nir_d_patch, c2l_scale_type='urbanf') + + this%fsrSF_nir_i_patch(begp:endp) = spval + call hist_addfld1d (fname='FSRSFNI', units='W/m^2', & + avgflag='A', long_name='diffuse nir reflected solar radiation', & + ptr_patch=this%fsrSF_nir_i_patch, c2l_scale_type='urbanf') + + this%fsrSF_nir_d_ln_patch(begp:endp) = spval + call hist_addfld1d (fname='FSRSFNDLN', units='W/m^2', & + avgflag='A', long_name='direct nir reflected solar radiation at local noon', & + ptr_patch=this%fsrSF_nir_d_ln_patch, c2l_scale_type='urbanf') + + this%ssre_fsr_nir_d_patch(begp:endp) = spval + call hist_addfld1d (fname='SSRE_FSRND', units='W/m^2', & + avgflag='A', long_name='surface snow effect on direct nir reflected solar radiation', & + ptr_patch=this%ssre_fsr_nir_d_patch, c2l_scale_type='urbanf') + + this%ssre_fsr_nir_i_patch(begp:endp) = spval + call hist_addfld1d (fname='SSRE_FSRNI', units='W/m^2', & + avgflag='A', long_name='surface snow effect on diffuse nir reflected solar radiation', & + ptr_patch=this%ssre_fsr_nir_i_patch, c2l_scale_type='urbanf') + + this%ssre_fsr_nir_d_ln_patch(begp:endp) = spval + call hist_addfld1d (fname='SSRE_FSRNDLN', units='W/m^2', & + avgflag='A', long_name='surface snow effect on direct nir reflected solar radiation at local noon', & + ptr_patch=this%ssre_fsr_nir_d_ln_patch, c2l_scale_type='urbanf') + end if this%sub_surf_abs_SW_patch(begp:endp) = spval call hist_addfld1d (fname='SNOINTABS', units='-', & avgflag='A', long_name='Fraction of incoming solar absorbed by lower snow layers', & diff --git a/src/biogeophys/SurfaceAlbedoMod.F90 b/src/biogeophys/SurfaceAlbedoMod.F90 index 9ad3229543..7d1763041f 100644 --- a/src/biogeophys/SurfaceAlbedoMod.F90 +++ b/src/biogeophys/SurfaceAlbedoMod.F90 @@ -13,7 +13,7 @@ module SurfaceAlbedoMod use landunit_varcon , only : istsoil, istcrop, istdlak use clm_varcon , only : grlnd, namep use clm_varpar , only : numrad, nlevcan, nlevsno, nlevcan - use clm_varctl , only : fsurdat, iulog, use_snicar_frc + use clm_varctl , only : fsurdat, iulog, use_snicar_frc, use_SSRE use pftconMod , only : pftcon use SnowSnicarMod , only : sno_nbr_aer, SNICAR_RT, DO_SNO_AER, DO_SNO_OC use AerosolMod , only : aerosol_type @@ -330,8 +330,10 @@ subroutine SurfaceAlbedo(bounds,nc, & albgri_dst => surfalb_inst%albgri_dst_col , & ! Output: [real(r8) (:,:) ] ground albedo without dust (diffuse) albsnd_hst => surfalb_inst%albsnd_hst_col , & ! Output: [real(r8) (:,:) ] snow albedo, direct, for history files (col,bnd) [frc] albsni_hst => surfalb_inst%albsni_hst_col , & ! Output: [real(r8) (:,:) ] snow ground albedo, diffuse, for history files (col,bnd) [frc] - albd => surfalb_inst%albd_patch , & ! Output: [real(r8) (:,:) ] surface albedo (direct) - albi => surfalb_inst%albi_patch , & ! Output: [real(r8) (:,:) ] surface albedo (diffuse) + albd => surfalb_inst%albd_patch , & ! Output: [real(r8) (:,:) ] surface albedo (direct) + albi => surfalb_inst%albi_patch , & ! Output: [real(r8) (:,:) ] surface albedo (diffuse) + albdSF => surfalb_inst%albdSF_patch , & ! Output: [real(r8) (:,:) ] diagnostic snow-free surface albedo (direct) + albiSF => surfalb_inst%albiSF_patch , & ! Output: [real(r8) (:,:) ] diagnostic snow-free surface albedo (diffuse) fabd => surfalb_inst%fabd_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by canopy per unit direct flux fabd_sun => surfalb_inst%fabd_sun_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by sunlit canopy per unit direct flux fabd_sha => surfalb_inst%fabd_sha_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by shaded canopy per unit direct flux @@ -395,6 +397,10 @@ subroutine SurfaceAlbedo(bounds,nc, & p = filter_nourbanp(fp) albd(p,ib) = 1._r8 albi(p,ib) = 1._r8 + if (use_SSRE) then + albdSF(p,ib) = 1._r8 + albiSF(p,ib) = 1._r8 + end if fabd(p,ib) = 0._r8 fabd_sun(p,ib) = 0._r8 fabd_sha(p,ib) = 0._r8 @@ -932,7 +938,19 @@ subroutine SurfaceAlbedo(bounds,nc, & rho(bounds%begp:bounds%endp, :), & tau(bounds%begp:bounds%endp, :), & canopystate_inst, temperature_inst, waterstate_inst, surfalb_inst) - + ! Run TwoStream again just to calculate the Snow Free (SF) albedo's + if (use_SSRE) then + if ( nlevcan > 1 )then + call endrun( 'ERROR: use_ssre option was NOT developed with allowance for multi-layer canopy: '// & + 'nlevcan can ONLY be 1 in when use_ssre is on') + end if + call TwoStream (bounds, filter_vegsol, num_vegsol, & + coszen_patch(bounds%begp:bounds%endp), & + rho(bounds%begp:bounds%endp, :), & + tau(bounds%begp:bounds%endp, :), & + canopystate_inst, temperature_inst, waterstate_inst, surfalb_inst, & + SFonly=.true.) + end if endif ! Determine values for non-vegetated patches where coszen > 0 @@ -952,6 +970,10 @@ subroutine SurfaceAlbedo(bounds,nc, & ftii(p,ib) = 1._r8 albd(p,ib) = albgrd(c,ib) albi(p,ib) = albgri(c,ib) + if (use_SSRE) then + albdSF(p,ib) = albsod(c,ib) + albiSF(p,ib) = albsoi(c,ib) + end if end do end do @@ -1094,7 +1116,8 @@ end subroutine SoilAlbedo subroutine TwoStream (bounds, & filter_vegsol, num_vegsol, & coszen, rho, tau, & - canopystate_inst, temperature_inst, waterstate_inst, surfalb_inst) + canopystate_inst, temperature_inst, waterstate_inst, surfalb_inst, & + SFonly) ! ! !DESCRIPTION: ! Two-stream fluxes for canopy radiative transfer @@ -1123,6 +1146,7 @@ subroutine TwoStream (bounds, & type(temperature_type) , intent(in) :: temperature_inst type(waterstate_type) , intent(in) :: waterstate_inst type(surfalb_type) , intent(inout) :: surfalb_inst + logical, optional , intent(in) :: SFonly ! If should just calculate the Snow Free albedos ! ! !LOCAL VARIABLES: integer :: fp,p,c,iv ! array indices @@ -1158,6 +1182,7 @@ subroutine TwoStream (bounds, & real(r8) :: laisum ! cumulative lai+sai for canopy layer (at middle of layer) real(r8) :: extkb ! direct beam extinction coefficient real(r8) :: extkn ! nitrogen allocation coefficient + logical :: lSFonly ! Local version of SFonly (Snow Free) flag !----------------------------------------------------------------------- ! Enforce expected array sizes @@ -1165,6 +1190,12 @@ subroutine TwoStream (bounds, & SHR_ASSERT_ALL((ubound(rho) == (/bounds%endp, numrad/)), errMsg(sourcefile, __LINE__)) SHR_ASSERT_ALL((ubound(tau) == (/bounds%endp, numrad/)), errMsg(sourcefile, __LINE__)) + if ( present(SFonly) )then + lSFonly = SFonly + else + lSFonly = .false. + end if + associate(& xl => pftcon%xl , & ! Input: ecophys const - leaf/stem orientation index @@ -1181,6 +1212,8 @@ subroutine TwoStream (bounds, & nrad => surfalb_inst%nrad_patch , & ! Input: [integer (:) ] number of canopy layers, above snow for radiative transfer albgrd => surfalb_inst%albgrd_col , & ! Input: [real(r8) (:,:) ] ground albedo (direct) (column-level) albgri => surfalb_inst%albgri_col , & ! Input: [real(r8) (:,:) ] ground albedo (diffuse)(column-level) + + ! For non-Snow Free fsun_z => surfalb_inst%fsun_z_patch , & ! Output: [real(r8) (:,:) ] sunlit fraction of canopy layer vcmaxcintsun => surfalb_inst%vcmaxcintsun_patch , & ! Output: [real(r8) (:) ] leaf to canopy scaling coefficient, sunlit leaf vcmax vcmaxcintsha => surfalb_inst%vcmaxcintsha_patch , & ! Output: [real(r8) (:) ] leaf to canopy scaling coefficient, shaded leaf vcmax @@ -1198,7 +1231,13 @@ subroutine TwoStream (bounds, & fabi_sha => surfalb_inst%fabi_sha_patch , & ! Output: [real(r8) (:,:) ] flux absorbed by shaded canopy per unit diffuse flux ftdd => surfalb_inst%ftdd_patch , & ! Output: [real(r8) (:,:) ] down direct flux below canopy per unit direct flx ftid => surfalb_inst%ftid_patch , & ! Output: [real(r8) (:,:) ] down diffuse flux below canopy per unit direct flx - ftii => surfalb_inst%ftii_patch & ! Output: [real(r8) (:,:) ] down diffuse flux below canopy per unit diffuse flx + ftii => surfalb_inst%ftii_patch , & ! Output: [real(r8) (:,:) ] down diffuse flux below canopy per unit diffuse flx + + ! Needed for SF Snow free case + albsod => surfalb_inst%albsod_col , & ! Input: [real(r8) (:,:) ] soil albedo (direct) + albsoi => surfalb_inst%albsoi_col , & ! Input: [real(r8) (:,:) ] soil albedo (diffuse) + albdSF => surfalb_inst%albdSF_patch , & ! Output: [real(r8) (:,:) ] Snow Free surface albedo (direct) + albiSF => surfalb_inst%albiSF_patch & ! Output: [real(r8) (:,:) ] Snow Free surface albedo (diffuse) ) ! Calculate two-stream parameters that are independent of waveband: @@ -1277,23 +1316,24 @@ subroutine TwoStream (bounds, & betail = 0.5_r8 * ((rho(p,ib)+tau(p,ib)) + (rho(p,ib)-tau(p,ib)) & * ((1._r8+chil(p))/2._r8)**2) / omegal - ! Adjust omega, betad, and betai for intercepted snow - - if (snowveg_onrad) then - tmp0 = (1._r8-fcansno(p))*omegal + fcansno(p)*omegas(ib) - tmp1 = ( (1._r8-fcansno(p))*omegal*betadl + fcansno(p)*omegas(ib)*betads ) / tmp0 - tmp2 = ( (1._r8-fcansno(p))*omegal*betail + fcansno(p)*omegas(ib)*betais ) / tmp0 + if ( lSFonly .or. ( (.not. snowveg_onrad) .and. (t_veg(p) > tfrz) ) ) then + ! Keep omega, betad, and betai as they are (for Snow free case or + ! when there is no snow + tmp0 = omegal + tmp1 = betadl + tmp2 = betail else - if (t_veg(p) > tfrz) then !no snow - tmp0 = omegal - tmp1 = betadl - tmp2 = betail + ! Adjust omega, betad, and betai for intercepted snow + if (snowveg_onrad) then + tmp0 = (1._r8-fcansno(p))*omegal + fcansno(p)*omegas(ib) + tmp1 = ( (1._r8-fcansno(p))*omegal*betadl + fcansno(p)*omegas(ib)*betads ) / tmp0 + tmp2 = ( (1._r8-fcansno(p))*omegal*betail + fcansno(p)*omegas(ib)*betais ) / tmp0 else tmp0 = (1._r8-fwet(p))*omegal + fwet(p)*omegas(ib) tmp1 = ( (1._r8-fwet(p))*omegal*betadl + fwet(p)*omegas(ib)*betads ) / tmp0 tmp2 = ( (1._r8-fwet(p))*omegal*betail + fwet(p)*omegas(ib)*betais ) / tmp0 end if - end if + end if ! end Snow free omega(p,ib) = tmp0 betad = tmp1 @@ -1323,10 +1363,17 @@ subroutine TwoStream (bounds, & s2 = exp(-t1) ! Direct beam - - u1 = b - c1/albgrd(c,ib) - u2 = b - c1*albgrd(c,ib) - u3 = f + c1*albgrd(c,ib) + if ( .not. lSFonly )then + u1 = b - c1/albgrd(c,ib) + u2 = b - c1*albgrd(c,ib) + u3 = f + c1*albgrd(c,ib) + else + ! Snow Free (SF) only + ! albsod instead of albgrd here: + u1 = b - c1/albsod(c,ib) + u2 = b - c1*albsod(c,ib) + u3 = f + c1*albsod(c,ib) + end if tmp2 = u1 - avmu(p)*h tmp3 = u1 + avmu(p)*h d1 = p1*tmp2/s1 - p2*tmp3*s1 @@ -1343,11 +1390,15 @@ subroutine TwoStream (bounds, & tmp9 = ( u3 - tmp8*(u2-tmp0) ) * s2 h5 = - ( tmp8*tmp4/s1 + tmp9 ) / d2 h6 = ( tmp8*tmp5*s1 + tmp9 ) / d2 - - albd(p,ib) = h1/sigma + h2 + h3 - ftid(p,ib) = h4*s2/sigma + h5*s1 + h6/s1 - ftdd(p,ib) = s2 - fabd(p,ib) = 1._r8 - albd(p,ib) - (1._r8-albgrd(c,ib))*ftdd(p,ib) - (1._r8-albgri(c,ib))*ftid(p,ib) + if ( .not. lSFonly )then + albd(p,ib) = h1/sigma + h2 + h3 + ftid(p,ib) = h4*s2/sigma + h5*s1 + h6/s1 + ftdd(p,ib) = s2 + fabd(p,ib) = 1._r8 - albd(p,ib) - (1._r8-albgrd(c,ib))*ftdd(p,ib) - (1._r8-albgri(c,ib))*ftid(p,ib) + else + albdSF(p,ib) = h1/sigma + h2 + h3 + end if + a1 = h1 / sigma * (1._r8 - s2*s2) / (2._r8 * twostext(p)) & + h2 * (1._r8 - s2*s1) / (twostext(p) + h) & @@ -1356,14 +1407,21 @@ subroutine TwoStream (bounds, & a2 = h4 / sigma * (1._r8 - s2*s2) / (2._r8 * twostext(p)) & + h5 * (1._r8 - s2*s1) / (twostext(p) + h) & + h6 * (1._r8 - s2/s1) / (twostext(p) - h) - - fabd_sun(p,ib) = (1._r8 - omega(p,ib)) * ( 1._r8 - s2 + 1._r8 / avmu(p) * (a1 + a2) ) - fabd_sha(p,ib) = fabd(p,ib) - fabd_sun(p,ib) + if ( .not. lSFonly )then + fabd_sun(p,ib) = (1._r8 - omega(p,ib)) * ( 1._r8 - s2 + 1._r8 / avmu(p) * (a1 + a2) ) + fabd_sha(p,ib) = fabd(p,ib) - fabd_sun(p,ib) + end if ! Diffuse - - u1 = b - c1/albgri(c,ib) - u2 = b - c1*albgri(c,ib) + if ( .not. lSFonly )then + u1 = b - c1/albgri(c,ib) + u2 = b - c1*albgri(c,ib) + else + ! Snow Free (SF) only + ! albsoi instead of albgri here: + u1 = b - c1/albsoi(c,ib) + u2 = b - c1*albsoi(c,ib) + end if tmp2 = u1 - avmu(p)*h tmp3 = u1 + avmu(p)*h d1 = p1*tmp2/s1 - p2*tmp3*s1 @@ -1375,222 +1433,195 @@ subroutine TwoStream (bounds, & h9 = tmp4 / (d2*s1) h10 = (-tmp5*s1) / d2 - albi(p,ib) = h7 + h8 - ftii(p,ib) = h9*s1 + h10/s1 - fabi(p,ib) = 1._r8 - albi(p,ib) - (1._r8-albgri(c,ib))*ftii(p,ib) - - a1 = h7 * (1._r8 - s2*s1) / (twostext(p) + h) + h8 * (1._r8 - s2/s1) / (twostext(p) - h) - a2 = h9 * (1._r8 - s2*s1) / (twostext(p) + h) + h10 * (1._r8 - s2/s1) / (twostext(p) - h) - - fabi_sun(p,ib) = (1._r8 - omega(p,ib)) / avmu(p) * (a1 + a2) - fabi_sha(p,ib) = fabi(p,ib) - fabi_sun(p,ib) - - ! Repeat two-stream calculations for each canopy layer to calculate derivatives. - ! tlai_z and tsai_z are the leaf+stem area increment for a layer. Derivatives are - ! calculated at the center of the layer. Derivatives are needed only for the - ! visible waveband to calculate absorbed PAR (per unit lai+sai) for each canopy layer. - ! Derivatives are calculated first per unit lai+sai and then normalized for sunlit - ! or shaded fraction of canopy layer. - - ! Sun/shade big leaf code uses only one layer, with canopy integrated values from above - ! and also canopy-integrated scaling coefficients - - if (ib == 1) then - if (nlevcan == 1) then - - ! sunlit fraction of canopy - fsun_z(p,1) = (1._r8 - s2) / t1 - - ! absorbed PAR (per unit sun/shade lai+sai) - laisum = elai(p)+esai(p) - fabd_sun_z(p,1) = fabd_sun(p,ib) / (fsun_z(p,1)*laisum) - fabi_sun_z(p,1) = fabi_sun(p,ib) / (fsun_z(p,1)*laisum) - fabd_sha_z(p,1) = fabd_sha(p,ib) / ((1._r8 - fsun_z(p,1))*laisum) - fabi_sha_z(p,1) = fabi_sha(p,ib) / ((1._r8 - fsun_z(p,1))*laisum) - - ! leaf to canopy scaling coefficients - extkn = 0.30_r8 - extkb = twostext(p) - vcmaxcintsun(p) = (1._r8 - exp(-(extkn+extkb)*elai(p))) / (extkn + extkb) - vcmaxcintsha(p) = (1._r8 - exp(-extkn*elai(p))) / extkn - vcmaxcintsun(p) - if (elai(p) > 0._r8) then - vcmaxcintsun(p) = vcmaxcintsun(p) / (fsun_z(p,1)*elai(p)) - vcmaxcintsha(p) = vcmaxcintsha(p) / ((1._r8 - fsun_z(p,1))*elai(p)) - else - vcmaxcintsun(p) = 0._r8 - vcmaxcintsha(p) = 0._r8 - end if - - else if (nlevcan > 1) then - do iv = 1, nrad(p) - - ! Cumulative lai+sai at center of layer - - if (iv == 1) then - laisum = 0.5_r8 * (tlai_z(p,iv)+tsai_z(p,iv)) - else - laisum = laisum + 0.5_r8 * ((tlai_z(p,iv-1)+tsai_z(p,iv-1))+(tlai_z(p,iv)+tsai_z(p,iv))) - end if - - ! Coefficients s1 and s2 depend on cumulative lai+sai. s2 is the sunlit fraction - - t1 = min(h*laisum, 40._r8) - s1 = exp(-t1) - t1 = min(twostext(p)*laisum, 40._r8) - s2 = exp(-t1) - fsun_z(p,iv) = s2 - - ! =============== - ! Direct beam - ! =============== - - ! Coefficients h1-h6 and a1,a2 depend of cumulative lai+sai - - u1 = b - c1/albgrd(c,ib) - u2 = b - c1*albgrd(c,ib) - u3 = f + c1*albgrd(c,ib) - tmp2 = u1 - avmu(p)*h - tmp3 = u1 + avmu(p)*h - d1 = p1*tmp2/s1 - p2*tmp3*s1 - tmp4 = u2 + avmu(p)*h - tmp5 = u2 - avmu(p)*h - d2 = tmp4/s1 - tmp5*s1 - h1 = -d*p4 - c1*f - tmp6 = d - h1*p3/sigma - tmp7 = ( d - c1 - h1/sigma*(u1+tmp0) ) * s2 - h2 = ( tmp6*tmp2/s1 - p2*tmp7 ) / d1 - h3 = - ( tmp6*tmp3*s1 - p1*tmp7 ) / d1 - h4 = -f*p3 - c1*d - tmp8 = h4/sigma - tmp9 = ( u3 - tmp8*(u2-tmp0) ) * s2 - h5 = - ( tmp8*tmp4/s1 + tmp9 ) / d2 - h6 = ( tmp8*tmp5*s1 + tmp9 ) / d2 - - a1 = h1 / sigma * (1._r8 - s2*s2) / (2._r8 * twostext(p)) & - + h2 * (1._r8 - s2*s1) / (twostext(p) + h) & - + h3 * (1._r8 - s2/s1) / (twostext(p) - h) - - a2 = h4 / sigma * (1._r8 - s2*s2) / (2._r8 * twostext(p)) & - + h5 * (1._r8 - s2*s1) / (twostext(p) + h) & - + h6 * (1._r8 - s2/s1) / (twostext(p) - h) - - ! Derivatives for h2, h3, h5, h6 and a1, a2 - - v = d1 - dv = h * p1 * tmp2 / s1 + h * p2 * tmp3 * s1 - - u = tmp6 * tmp2 / s1 - p2 * tmp7 - du = h * tmp6 * tmp2 / s1 + twostext(p) * p2 * tmp7 - dh2 = (v * du - u * dv) / (v * v) - - u = -tmp6 * tmp3 * s1 + p1 * tmp7 - du = h * tmp6 * tmp3 * s1 - twostext(p) * p1 * tmp7 - dh3 = (v * du - u * dv) / (v * v) - - v = d2 - dv = h * tmp4 / s1 + h * tmp5 * s1 - - u = -h4/sigma * tmp4 / s1 - tmp9 - du = -h * h4/sigma * tmp4 / s1 + twostext(p) * tmp9 - dh5 = (v * du - u * dv) / (v * v) - - u = h4/sigma * tmp5 * s1 + tmp9 - du = -h * h4/sigma * tmp5 * s1 - twostext(p) * tmp9 - dh6 = (v * du - u * dv) / (v * v) - - da1 = h1/sigma * s2*s2 + h2 * s2*s1 + h3 * s2/s1 & - + (1._r8 - s2*s1) / (twostext(p) + h) * dh2 & - + (1._r8 - s2/s1) / (twostext(p) - h) * dh3 - da2 = h4/sigma * s2*s2 + h5 * s2*s1 + h6 * s2/s1 & - + (1._r8 - s2*s1) / (twostext(p) + h) * dh5 & - + (1._r8 - s2/s1) / (twostext(p) - h) * dh6 - - ! Flux derivatives - - d_ftid = -twostext(p)*h4/sigma*s2 - h*h5*s1 + h*h6/s1 + dh5*s1 + dh6/s1 - d_fabd = -(dh2+dh3) + (1._r8-albgrd(c,ib))*twostext(p)*s2 - (1._r8-albgri(c,ib))*d_ftid - d_fabd_sun = (1._r8 - omega(p,ib)) * (twostext(p)*s2 + 1._r8 / avmu(p) * (da1 + da2)) - d_fabd_sha = d_fabd - d_fabd_sun - - fabd_sun_z(p,iv) = max(d_fabd_sun, 0._r8) - fabd_sha_z(p,iv) = max(d_fabd_sha, 0._r8) - - ! Flux derivatives are APARsun and APARsha per unit (LAI+SAI). Need - ! to normalize derivatives by sunlit or shaded fraction to get - ! APARsun per unit (LAI+SAI)sun and APARsha per unit (LAI+SAI)sha - - fabd_sun_z(p,iv) = fabd_sun_z(p,iv) / fsun_z(p,iv) - fabd_sha_z(p,iv) = fabd_sha_z(p,iv) / (1._r8 - fsun_z(p,iv)) - - ! =============== - ! Diffuse - ! =============== - - ! Coefficients h7-h10 and a1,a2 depend of cumulative lai+sai - - u1 = b - c1/albgri(c,ib) - u2 = b - c1*albgri(c,ib) - tmp2 = u1 - avmu(p)*h - tmp3 = u1 + avmu(p)*h - d1 = p1*tmp2/s1 - p2*tmp3*s1 - tmp4 = u2 + avmu(p)*h - tmp5 = u2 - avmu(p)*h - d2 = tmp4/s1 - tmp5*s1 - h7 = (c1*tmp2) / (d1*s1) - h8 = (-c1*tmp3*s1) / d1 - h9 = tmp4 / (d2*s1) - h10 = (-tmp5*s1) / d2 - - a1 = h7 * (1._r8 - s2*s1) / (twostext(p) + h) + h8 * (1._r8 - s2/s1) / (twostext(p) - h) - a2 = h9 * (1._r8 - s2*s1) / (twostext(p) + h) + h10 * (1._r8 - s2/s1) / (twostext(p) - h) - - ! Derivatives for h7, h8, h9, h10 and a1, a2 - - v = d1 - dv = h * p1 * tmp2 / s1 + h * p2 * tmp3 * s1 - - u = c1 * tmp2 / s1 - du = h * c1 * tmp2 / s1 - dh7 = (v * du - u * dv) / (v * v) - - u = -c1 * tmp3 * s1 - du = h * c1 * tmp3 * s1 - dh8 = (v * du - u * dv) / (v * v) - - v = d2 - dv = h * tmp4 / s1 + h * tmp5 * s1 - - u = tmp4 / s1 - du = h * tmp4 / s1 - dh9 = (v * du - u * dv) / (v * v) - - u = -tmp5 * s1 - du = h * tmp5 * s1 - dh10 = (v * du - u * dv) / (v * v) - - da1 = h7*s2*s1 + h8*s2/s1 + (1._r8-s2*s1)/(twostext(p)+h)*dh7 + (1._r8-s2/s1)/(twostext(p)-h)*dh8 - da2 = h9*s2*s1 + h10*s2/s1 + (1._r8-s2*s1)/(twostext(p)+h)*dh9 + (1._r8-s2/s1)/(twostext(p)-h)*dh10 - - ! Flux derivatives - - d_ftii = -h * h9 * s1 + h * h10 / s1 + dh9 * s1 + dh10 / s1 - d_fabi = -(dh7+dh8) - (1._r8-albgri(c,ib))*d_ftii - d_fabi_sun = (1._r8 - omega(p,ib)) / avmu(p) * (da1 + da2) - d_fabi_sha = d_fabi - d_fabi_sun - - fabi_sun_z(p,iv) = max(d_fabi_sun, 0._r8) - fabi_sha_z(p,iv) = max(d_fabi_sha, 0._r8) + + ! Final Snow Free albedo + if ( lSFonly )then + albiSF(p,ib) = h7 + h8 + else + ! For non snow Free case, adjustments continue + albi(p,ib) = h7 + h8 + ftii(p,ib) = h9*s1 + h10/s1 + fabi(p,ib) = 1._r8 - albi(p,ib) - (1._r8-albgri(c,ib))*ftii(p,ib) - ! Flux derivatives are APARsun and APARsha per unit (LAI+SAI). Need - ! to normalize derivatives by sunlit or shaded fraction to get - ! APARsun per unit (LAI+SAI)sun and APARsha per unit (LAI+SAI)sha + a1 = h7 * (1._r8 - s2*s1) / (twostext(p) + h) + h8 * (1._r8 - s2/s1) / (twostext(p) - h) + a2 = h9 * (1._r8 - s2*s1) / (twostext(p) + h) + h10 * (1._r8 - s2/s1) / (twostext(p) - h) - fabi_sun_z(p,iv) = fabi_sun_z(p,iv) / fsun_z(p,iv) - fabi_sha_z(p,iv) = fabi_sha_z(p,iv) / (1._r8 - fsun_z(p,iv)) + fabi_sun(p,ib) = (1._r8 - omega(p,ib)) / avmu(p) * (a1 + a2) + fabi_sha(p,ib) = fabi(p,ib) - fabi_sun(p,ib) + + ! Repeat two-stream calculations for each canopy layer to calculate derivatives. + ! tlai_z and tsai_z are the leaf+stem area increment for a layer. Derivatives are + ! calculated at the center of the layer. Derivatives are needed only for the + ! visible waveband to calculate absorbed PAR (per unit lai+sai) for each canopy layer. + ! Derivatives are calculated first per unit lai+sai and then normalized for sunlit + ! or shaded fraction of canopy layer. + + ! Sun/shade big leaf code uses only one layer, with canopy integrated values from above + ! and also canopy-integrated scaling coefficients + + if (ib == 1) then + if (nlevcan == 1) then + + ! sunlit fraction of canopy + fsun_z(p,1) = (1._r8 - s2) / t1 + + ! absorbed PAR (per unit sun/shade lai+sai) + laisum = elai(p)+esai(p) + fabd_sun_z(p,1) = fabd_sun(p,ib) / (fsun_z(p,1)*laisum) + fabi_sun_z(p,1) = fabi_sun(p,ib) / (fsun_z(p,1)*laisum) + fabd_sha_z(p,1) = fabd_sha(p,ib) / ((1._r8 - fsun_z(p,1))*laisum) + fabi_sha_z(p,1) = fabi_sha(p,ib) / ((1._r8 - fsun_z(p,1))*laisum) + + ! leaf to canopy scaling coefficients + extkn = 0.30_r8 + extkb = twostext(p) + vcmaxcintsun(p) = (1._r8 - exp(-(extkn+extkb)*elai(p))) / (extkn + extkb) + vcmaxcintsha(p) = (1._r8 - exp(-extkn*elai(p))) / extkn - vcmaxcintsun(p) + if (elai(p) > 0._r8) then + vcmaxcintsun(p) = vcmaxcintsun(p) / (fsun_z(p,1)*elai(p)) + vcmaxcintsha(p) = vcmaxcintsha(p) / ((1._r8 - fsun_z(p,1))*elai(p)) + else + vcmaxcintsun(p) = 0._r8 + vcmaxcintsha(p) = 0._r8 + end if + + else if (nlevcan > 1)then + do iv = 1, nrad(p) + + ! Cumulative lai+sai at center of layer + + if (iv == 1) then + laisum = 0.5_r8 * (tlai_z(p,iv)+tsai_z(p,iv)) + else + laisum = laisum + 0.5_r8 * ((tlai_z(p,iv-1)+tsai_z(p,iv-1))+(tlai_z(p,iv)+tsai_z(p,iv))) + end if + + ! Coefficients s1 and s2 depend on cumulative lai+sai. s2 is the sunlit fraction + + t1 = min(h*laisum, 40._r8) + s1 = exp(-t1) + t1 = min(twostext(p)*laisum, 40._r8) + s2 = exp(-t1) + fsun_z(p,iv) = s2 + + ! =============== + ! Direct beam + ! =============== + + ! Coefficients h1-h6 and a1,a2 depend of cumulative lai+sai + + u1 = b - c1/albgrd(c,ib) + u2 = b - c1*albgrd(c,ib) + u3 = f + c1*albgrd(c,ib) + + ! Derivatives for h2, h3, h5, h6 and a1, a2 + + v = d1 + dv = h * p1 * tmp2 / s1 + h * p2 * tmp3 * s1 + + u = tmp6 * tmp2 / s1 - p2 * tmp7 + du = h * tmp6 * tmp2 / s1 + twostext(p) * p2 * tmp7 + dh2 = (v * du - u * dv) / (v * v) + + u = -tmp6 * tmp3 * s1 + p1 * tmp7 + du = h * tmp6 * tmp3 * s1 - twostext(p) * p1 * tmp7 + dh3 = (v * du - u * dv) / (v * v) + + v = d2 + dv = h * tmp4 / s1 + h * tmp5 * s1 + + u = -h4/sigma * tmp4 / s1 - tmp9 + du = -h * h4/sigma * tmp4 / s1 + twostext(p) * tmp9 + dh5 = (v * du - u * dv) / (v * v) + + u = h4/sigma * tmp5 * s1 + tmp9 + du = -h * h4/sigma * tmp5 * s1 - twostext(p) * tmp9 + dh6 = (v * du - u * dv) / (v * v) + + da1 = h1/sigma * s2*s2 + h2 * s2*s1 + h3 * s2/s1 & + + (1._r8 - s2*s1) / (twostext(p) + h) * dh2 & + + (1._r8 - s2/s1) / (twostext(p) - h) * dh3 + da2 = h4/sigma * s2*s2 + h5 * s2*s1 + h6 * s2/s1 & + + (1._r8 - s2*s1) / (twostext(p) + h) * dh5 & + + (1._r8 - s2/s1) / (twostext(p) - h) * dh6 + + ! Flux derivatives + + d_ftid = -twostext(p)*h4/sigma*s2 - h*h5*s1 + h*h6/s1 + dh5*s1 + dh6/s1 + d_fabd = -(dh2+dh3) + (1._r8-albgrd(c,ib))*twostext(p)*s2 - (1._r8-albgri(c,ib))*d_ftid + d_fabd_sun = (1._r8 - omega(p,ib)) * (twostext(p)*s2 + 1._r8 / avmu(p) * (da1 + da2)) + d_fabd_sha = d_fabd - d_fabd_sun + + fabd_sun_z(p,iv) = max(d_fabd_sun, 0._r8) + fabd_sha_z(p,iv) = max(d_fabd_sha, 0._r8) + + ! Flux derivatives are APARsun and APARsha per unit (LAI+SAI). Need + ! to normalize derivatives by sunlit or shaded fraction to get + ! APARsun per unit (LAI+SAI)sun and APARsha per unit (LAI+SAI)sha + + fabd_sun_z(p,iv) = fabd_sun_z(p,iv) / fsun_z(p,iv) + fabd_sha_z(p,iv) = fabd_sha_z(p,iv) / (1._r8 - fsun_z(p,iv)) + + ! =============== + ! Diffuse + ! =============== + + ! Coefficients h7-h10 and a1,a2 depend of cumulative lai+sai + + u1 = b - c1/albgri(c,ib) + u2 = b - c1*albgri(c,ib) - end do ! end of canopy layer loop - end if - end if + a1 = h7 * (1._r8 - s2*s1) / (twostext(p) + h) + h8 * (1._r8 - s2/s1) / (twostext(p) - h) + a2 = h9 * (1._r8 - s2*s1) / (twostext(p) + h) + h10 * (1._r8 - s2/s1) / (twostext(p) - h) + + ! Derivatives for h7, h8, h9, h10 and a1, a2 + + v = d1 + dv = h * p1 * tmp2 / s1 + h * p2 * tmp3 * s1 + + u = c1 * tmp2 / s1 + du = h * c1 * tmp2 / s1 + dh7 = (v * du - u * dv) / (v * v) + + u = -c1 * tmp3 * s1 + du = h * c1 * tmp3 * s1 + dh8 = (v * du - u * dv) / (v * v) + + v = d2 + dv = h * tmp4 / s1 + h * tmp5 * s1 + + u = tmp4 / s1 + du = h * tmp4 / s1 + dh9 = (v * du - u * dv) / (v * v) + + u = -tmp5 * s1 + du = h * tmp5 * s1 + dh10 = (v * du - u * dv) / (v * v) + + da1 = h7*s2*s1 + h8*s2/s1 + (1._r8-s2*s1)/(twostext(p)+h)*dh7 + (1._r8-s2/s1)/(twostext(p)-h)*dh8 + da2 = h9*s2*s1 + h10*s2/s1 + (1._r8-s2*s1)/(twostext(p)+h)*dh9 + (1._r8-s2/s1)/(twostext(p)-h)*dh10 + + ! Flux derivatives + + d_ftii = -h * h9 * s1 + h * h10 / s1 + dh9 * s1 + dh10 / s1 + d_fabi = -(dh7+dh8) - (1._r8-albgri(c,ib))*d_ftii + d_fabi_sun = (1._r8 - omega(p,ib)) / avmu(p) * (da1 + da2) + d_fabi_sha = d_fabi - d_fabi_sun + + fabi_sun_z(p,iv) = max(d_fabi_sun, 0._r8) + fabi_sha_z(p,iv) = max(d_fabi_sha, 0._r8) + + ! Flux derivatives are APARsun and APARsha per unit (LAI+SAI). Need + ! to normalize derivatives by sunlit or shaded fraction to get + ! APARsun per unit (LAI+SAI)sun and APARsha per unit (LAI+SAI)sha + + fabi_sun_z(p,iv) = fabi_sun_z(p,iv) / fsun_z(p,iv) + fabi_sha_z(p,iv) = fabi_sha_z(p,iv) / (1._r8 - fsun_z(p,iv)) + + end do ! end of iv loop + end if ! nlevcan + end if ! first band + end if ! NOT lSFonly end do ! end of pft loop end do ! end of radiation band loop diff --git a/src/biogeophys/SurfaceAlbedoType.F90 b/src/biogeophys/SurfaceAlbedoType.F90 index 1540d9f991..0319c50b47 100644 --- a/src/biogeophys/SurfaceAlbedoType.F90 +++ b/src/biogeophys/SurfaceAlbedoType.F90 @@ -8,6 +8,7 @@ module SurfaceAlbedoType use decompMod , only : bounds_type use clm_varpar , only : numrad, nlevcan, nlevsno use abortutils , only : endrun + use clm_varctl , only : use_SSRE ! ! !PUBLIC TYPES: implicit none @@ -19,6 +20,8 @@ module SurfaceAlbedoType real(r8), pointer :: coszen_col (:) ! col cosine of solar zenith angle real(r8), pointer :: albd_patch (:,:) ! patch surface albedo (direct) (numrad) real(r8), pointer :: albi_patch (:,:) ! patch surface albedo (diffuse) (numrad) + real(r8), pointer :: albdSF_patch (:,:) ! patch snow-free surface albedo (direct) (numrad) + real(r8), pointer :: albiSF_patch (:,:) ! patch snow-free surface albedo (diffuse) (numrad) real(r8), pointer :: albgrd_pur_col (:,:) ! col pure snow ground direct albedo (numrad) real(r8), pointer :: albgri_pur_col (:,:) ! col pure snow ground diffuse albedo (numrad) real(r8), pointer :: albgrd_bc_col (:,:) ! col ground direct albedo without BC (numrad) @@ -96,6 +99,7 @@ subroutine InitAllocate(this, bounds) ! !USES: use shr_infnan_mod, only: nan => shr_infnan_nan, assignment(=) use clm_varcon , only: spval, ispval + use clm_varctl , only: use_SSRE ! ! !ARGUMENTS: class(surfalb_type) :: this @@ -126,7 +130,8 @@ subroutine InitAllocate(this, bounds) allocate(this%albgri_dst_col (begc:endc,numrad)) ; this%albgri_dst_col (:,:) = nan allocate(this%albd_patch (begp:endp,numrad)) ; this%albd_patch (:,:) = nan allocate(this%albi_patch (begp:endp,numrad)) ; this%albi_patch (:,:) = nan - + allocate(this%albdSF_patch (begp:endp,numrad)) ; this%albdSF_patch (:,:) = nan + allocate(this%albiSF_patch (begp:endp,numrad)) ; this%albiSF_patch (:,:) = nan allocate(this%ftdd_patch (begp:endp,numrad)) ; this%ftdd_patch (:,:) = nan allocate(this%ftid_patch (begp:endp,numrad)) ; this%ftid_patch (:,:) = nan allocate(this%ftii_patch (begp:endp,numrad)) ; this%ftii_patch (:,:) = nan @@ -161,6 +166,7 @@ subroutine InitHistory(this, bounds) ! History fields initialization ! ! !USES: + use shr_kind_mod , only: cs => shr_kind_CS use shr_infnan_mod, only: nan => shr_infnan_nan, assignment(=) use clm_varcon , only: spval use histFileMod , only: hist_addfld1d, hist_addfld2d @@ -172,6 +178,7 @@ subroutine InitHistory(this, bounds) ! !LOCAL VARIABLES: integer :: begp, endp integer :: begc, endc + character(len=cs) :: defaultoutput !--------------------------------------------------------------------- begp = bounds%begp; endp = bounds%endp @@ -192,15 +199,28 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='ground albedo (indirect)', & ptr_col=this%albgri_col, default='inactive') + if (use_SSRE) then + this%albdSF_patch(begp:endp,:) = spval + call hist_addfld2d (fname='ALBDSF', units='proportion', type2d='numrad', & + avgflag='A', long_name='diagnostic snow-free surface albedo (direct)', & + ptr_patch=this%albdSF_patch, default='active', c2l_scale_type='urbanf') + this%albiSF_patch(begp:endp,:) = spval + call hist_addfld2d (fname='ALBISF', units='proportion', type2d='numrad', & + avgflag='A', long_name='diagnostic snow-free surface albedo (indirect)', & + ptr_patch=this%albiSF_patch, default='active', c2l_scale_type='urbanf') + defaultoutput = "active" + else + defaultoutput = "inactive" + end if this%albd_patch(begp:endp,:) = spval call hist_addfld2d (fname='ALBD', units='proportion', type2d='numrad', & avgflag='A', long_name='surface albedo (direct)', & - ptr_patch=this%albd_patch, default='inactive', c2l_scale_type='urbanf') + ptr_patch=this%albd_patch, default=defaultoutput, c2l_scale_type='urbanf') this%albi_patch(begp:endp,:) = spval call hist_addfld2d (fname='ALBI', units='proportion', type2d='numrad', & avgflag='A', long_name='surface albedo (indirect)', & - ptr_patch=this%albi_patch, default='inactive', c2l_scale_type='urbanf') + ptr_patch=this%albi_patch, default=defaultoutput, c2l_scale_type='urbanf') end subroutine InitHistory @@ -229,7 +249,10 @@ subroutine InitCold(this, bounds) this%albsni_hst_col (begc:endc, :) = 0.6_r8 this%albd_patch (begp:endp, :) = 0.2_r8 this%albi_patch (begp:endp, :) = 0.2_r8 - + if (use_SSRE) then + this%albdSF_patch (begp:endp, :) = 0.2_r8 + this%albiSF_patch (begp:endp, :) = 0.2_r8 + end if this%albgrd_pur_col (begc:endc, :) = 0.2_r8 this%albgri_pur_col (begc:endc, :) = 0.2_r8 this%albgrd_bc_col (begc:endc, :) = 0.2_r8 @@ -301,7 +324,17 @@ subroutine Restart(this, bounds, ncid, flag, & dim1name='pft', dim2name='numrad', switchdim=.true., & long_name='surface albedo (diffuse) (0 to 1)', units='', & interpinic_flag='interp', readvar=readvar, data=this%albi_patch) - + if (use_SSRE) then + call restartvar(ncid=ncid, flag=flag, varname='albdSF', xtype=ncd_double, & + dim1name='pft', dim2name='numrad', switchdim=.true., & + long_name='diagnostic snow-free surface albedo (direct) (0 to 1)', units='', & + interpinic_flag='interp', readvar=readvar, data=this%albdSF_patch) + + call restartvar(ncid=ncid, flag=flag, varname='albiSF', xtype=ncd_double, & + dim1name='pft', dim2name='numrad', switchdim=.true., & + long_name='diagnostic snow-free surface albedo (diffuse) (0 to 1)', units='', & + interpinic_flag='interp', readvar=readvar, data=this%albiSF_patch) + end if call restartvar(ncid=ncid, flag=flag, varname='albgrd', xtype=ncd_double, & dim1name='column', dim2name='numrad', switchdim=.true., & long_name='ground albedo (direct) (0 to 1)', units='', & diff --git a/src/biogeophys/SurfaceRadiationMod.F90 b/src/biogeophys/SurfaceRadiationMod.F90 index d0b56f6a27..8af0aaa633 100644 --- a/src/biogeophys/SurfaceRadiationMod.F90 +++ b/src/biogeophys/SurfaceRadiationMod.F90 @@ -54,7 +54,13 @@ module SurfaceRadiationMod real(r8), pointer, private :: fsr_vis_d_patch (:) ! patch reflected direct beam vis solar radiation (W/m**2) real(r8), pointer, private :: fsr_vis_i_patch (:) ! patch reflected diffuse vis solar radiation (W/m**2) real(r8), pointer, private :: fsr_vis_d_ln_patch (:) ! patch reflected direct beam vis solar radiation at local noon (W/m**2) - + ! diagnostic fluxes: + real(r8), pointer, private :: fsrSF_vis_d_patch (:) ! snow-free patch reflected direct beam vis solar radiation (W/m**2) + real(r8), pointer, private :: fsrSF_vis_i_patch (:) ! snow-free patch reflected diffuse vis solar radiation (W/m**2) + real(r8), pointer, private :: fsrSF_vis_d_ln_patch (:) ! snow-free patch reflected direct beam vis solar radiation at local noon (W/m**2) + real(r8), pointer, private :: ssre_fsr_vis_d_patch (:) ! snow radiative effect + real(r8), pointer, private :: ssre_fsr_vis_i_patch (:) ! snow radiative effect + real(r8), pointer, private :: ssre_fsr_vis_d_ln_patch(:)! snow radiative effect real(r8), pointer, private :: fsds_sno_vd_patch (:) ! patch incident visible, direct radiation on snow (for history files) [W/m2] real(r8), pointer, private :: fsds_sno_nd_patch (:) ! patch incident near-IR, direct radiation on snow (for history files) [W/m2] real(r8), pointer, private :: fsds_sno_vi_patch (:) ! patch incident visible, diffuse radiation on snow (for history files) [W/m2] @@ -122,6 +128,12 @@ subroutine InitAllocate(this, bounds) allocate(this%fsr_vis_d_patch (begp:endp)) ; this%fsr_vis_d_patch (:) = nan allocate(this%fsr_vis_d_ln_patch (begp:endp)) ; this%fsr_vis_d_ln_patch (:) = nan allocate(this%fsr_vis_i_patch (begp:endp)) ; this%fsr_vis_i_patch (:) = nan + allocate(this%fsrSF_vis_d_patch (begp:endp)) ; this%fsrSF_vis_d_patch (:) = nan + allocate(this%fsrSF_vis_d_ln_patch (begp:endp)) ; this%fsrSF_vis_d_ln_patch (:) = nan + allocate(this%fsrSF_vis_i_patch (begp:endp)) ; this%fsrSF_vis_i_patch (:) = nan + allocate(this%ssre_fsr_vis_d_patch (begp:endp)) ; this%ssre_fsr_vis_d_patch (:) = nan + allocate(this%ssre_fsr_vis_d_ln_patch(begp:endp)) ; this%ssre_fsr_vis_d_ln_patch(:) = nan + allocate(this%ssre_fsr_vis_i_patch (begp:endp)) ; this%ssre_fsr_vis_i_patch (:) = nan allocate(this%fsr_sno_vd_patch (begp:endp)) ; this%fsr_sno_vd_patch (:) = nan allocate(this%fsr_sno_nd_patch (begp:endp)) ; this%fsr_sno_nd_patch (:) = nan allocate(this%fsr_sno_vi_patch (begp:endp)) ; this%fsr_sno_vi_patch (:) = nan @@ -147,6 +159,7 @@ subroutine InitHistory(this, bounds) use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) use clm_varcon , only : spval use histFileMod , only : hist_addfld1d, hist_addfld2d + use clm_varctl , only : use_SSRE ! ! !ARGUMENTS: class(surfrad_type) :: this @@ -217,12 +230,30 @@ subroutine InitHistory(this, bounds) call hist_addfld1d (fname='FSRVD', units='W/m^2', & avgflag='A', long_name='direct vis reflected solar radiation', & ptr_patch=this%fsr_vis_d_patch, c2l_scale_type='urbanf') - this%fsr_vis_i_patch(begp:endp) = spval call hist_addfld1d (fname='FSRVI', units='W/m^2', & avgflag='A', long_name='diffuse vis reflected solar radiation', & ptr_patch=this%fsr_vis_i_patch, c2l_scale_type='urbanf') - + ! diagnostic fluxes + if (use_SSRE) then + this%fsrSF_vis_d_patch(begp:endp) = spval + call hist_addfld1d (fname='FSRSFVD', units='W/m^2', & + avgflag='A', long_name='direct vis reflected solar radiation', & + ptr_patch=this%fsrSF_vis_d_patch, c2l_scale_type='urbanf') + this%fsrSF_vis_i_patch(begp:endp) = spval + call hist_addfld1d (fname='FSRSFVI', units='W/m^2', & + avgflag='A', long_name='diffuse vis reflected solar radiation', & + ptr_patch=this%fsrSF_vis_i_patch, c2l_scale_type='urbanf') + + this%ssre_fsr_vis_d_patch(begp:endp) = spval + call hist_addfld1d (fname='SSRE_FSRVD', units='W/m^2', & + avgflag='A', long_name='surface snow radiatve effect on direct vis reflected solar radiation', & + ptr_patch=this%ssre_fsr_vis_d_patch, c2l_scale_type='urbanf') + this%ssre_fsr_vis_i_patch(begp:endp) = spval + call hist_addfld1d (fname='SSRE_FSRVI', units='W/m^2', & + avgflag='A', long_name='surface snow radiatve effect on diffuse vis reflected solar radiation', & + ptr_patch=this%ssre_fsr_vis_i_patch, c2l_scale_type='urbanf') + end if this%fsds_vis_d_ln_patch(begp:endp) = spval call hist_addfld1d (fname='FSDSVDLN', units='W/m^2', & avgflag='A', long_name='direct vis incident solar radiation at local noon', & @@ -242,7 +273,17 @@ subroutine InitHistory(this, bounds) call hist_addfld1d (fname='FSRVDLN', units='W/m^2', & avgflag='A', long_name='direct vis reflected solar radiation at local noon', & ptr_patch=this%fsr_vis_d_ln_patch, c2l_scale_type='urbanf') - + ! diagnostic flux + if (use_SSRE) then + this%fsrSF_vis_d_ln_patch(begp:endp) = spval + call hist_addfld1d (fname='FSRSFVDLN', units='W/m^2', & + avgflag='A', long_name='direct vis reflected solar radiation at local noon', & + ptr_patch=this%fsrSF_vis_d_ln_patch, c2l_scale_type='urbanf') + this%ssre_fsr_vis_d_ln_patch(begp:endp) = spval + call hist_addfld1d (fname='SSRE_FSRVDLN', units='W/m^2', & + avgflag='A', long_name='surface snow radiatve effect on direct vis reflected solar radiation at local noon', & + ptr_patch=this%ssre_fsr_vis_d_ln_patch, c2l_scale_type='urbanf') + end if this%fsds_sno_vd_patch(begp:endp) = spval call hist_addfld1d (fname='SNOFSDSVD', units='W/m^2', & avgflag='A', long_name='direct vis incident solar radiation on snow', & @@ -435,10 +476,10 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & ! ! !USES: use clm_varpar , only : numrad, nlevsno - use clm_varcon , only : spval, degpsec, isecspday + use clm_varcon , only : spval use landunit_varcon , only : istsoil, istcrop - use clm_varctl , only : subgridflag, use_snicar_frc, iulog - use clm_time_manager , only : get_curr_date, get_step_size + use clm_varctl , only : subgridflag, use_snicar_frc, iulog, use_SSRE + use clm_time_manager , only : get_step_size, is_near_local_noon use SnowSnicarMod , only : DO_SNO_OC use abortutils , only : endrun ! @@ -471,13 +512,13 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & integer :: i ! layer index [idx] real(r8) :: rnir ! reflected solar radiation [nir] (W/m**2) real(r8) :: rvis ! reflected solar radiation [vis] (W/m**2) + real(r8) :: rnirSF ! snow-free reflected solar radiation [nir] (W/m**2) + real(r8) :: rvisSF ! snow-free reflected solar radiation [vis] (W/m**2) real(r8) :: trd(bounds%begp:bounds%endp,numrad) ! transmitted solar radiation: direct (W/m**2) real(r8) :: tri(bounds%begp:bounds%endp,numrad) ! transmitted solar radiation: diffuse (W/m**2) real(r8) :: cad(bounds%begp:bounds%endp,numrad) ! direct beam absorbed by canopy (W/m**2) real(r8) :: cai(bounds%begp:bounds%endp,numrad) ! diffuse radiation absorbed by canopy (W/m**2) - integer :: local_secp1 ! seconds into current date in local time real(r8) :: dtime ! land model time step (sec) - integer :: year,month,day,secs ! calendar info for current time step real(r8) :: sabg_snl_sum ! temporary, absorbed energy in all active snow layers [W/m2] real(r8) :: absrad_pur ! temp: absorbed solar radiation by pure snow [W/m2] real(r8) :: absrad_bc ! temp: absorbed solar radiation without BC [W/m2] @@ -489,8 +530,6 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & real(r8) :: sabg_dst(bounds%begp:bounds%endp) ! solar radiation absorbed by ground without dust [W/m2] real(r8) :: parveg(bounds%begp:bounds%endp) ! absorbed par by vegetation (W/m**2) ! - integer, parameter :: noonsec = isecspday / 2 ! seconds at local noon - ! !------------------------------------------------------------------------------ associate( & @@ -520,6 +559,8 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & albsoi => surfalb_inst%albsoi_col , & ! Input: [real(r8) (:,:) ] diffuse soil albedo (col,bnd) [frc] albd => surfalb_inst%albd_patch , & ! Input: [real(r8) (:,:) ] surface albedo (direct) albi => surfalb_inst%albi_patch , & ! Input: [real(r8) (:,:) ] surface albedo (diffuse) + albdSF => surfalb_inst%albdSF_patch , & ! Input: [real(r8) (:,:) ] snow-free surface albedo (direct) + albiSF => surfalb_inst%albiSF_patch , & ! Input: [real(r8) (:,:) ] snow-free surface albedo (diffuse) fabd => surfalb_inst%fabd_patch , & ! Input: [real(r8) (:,:) ] flux absorbed by canopy per unit direct flux fabd_sun => surfalb_inst%fabd_sun_patch , & ! Input: [real(r8) (:,:) ] flux absorbed by sunlit canopy per unit direct flux fabd_sha => surfalb_inst%fabd_sha_patch , & ! Input: [real(r8) (:,:) ] flux absorbed by shaded canopy per unit direct flux @@ -543,6 +584,8 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & fsun => canopystate_inst%fsun_patch , & ! Output: [real(r8) (:) ] sunlit fraction of canopy fsa => solarabs_inst%fsa_patch , & ! Output: [real(r8) (:) ] solar radiation absorbed (total) (W/m**2) fsr => solarabs_inst%fsr_patch , & ! Output: [real(r8) (:) ] solar radiation reflected (W/m**2) + fsrSF => solarabs_inst%fsrSF_patch , & ! Output: [real(r8) (:) ] diagnostic snow-free solar radiation reflected (W/m**2) + ssre_fsr => solarabs_inst%ssre_fsr_patch , & ! Output: [real(r8) (:) ] diagnostic snow-free solar radiation reflected (W/m**2) sabv => solarabs_inst%sabv_patch , & ! Output: [real(r8) (:) ] solar radiation absorbed by vegetation (W/m**2) sabg => solarabs_inst%sabg_patch , & ! Output: [real(r8) (:) ] solar radiation absorbed by ground (W/m**2) sabg_pen => solarabs_inst%sabg_pen_patch , & ! Output: [real(r8) (:) ] solar (rural) radiation penetrating top soisno layer (W/m**2) @@ -555,14 +598,25 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & fsds_nir_d => solarabs_inst%fsds_nir_d_patch , & ! Output: [real(r8) (:) ] incident direct beam nir solar radiation (W/m**2) fsds_nir_d_ln => solarabs_inst%fsds_nir_d_ln_patch , & ! Output: [real(r8) (:) ] incident direct beam nir solar rad at local noon (W/m**2) fsds_nir_i => solarabs_inst%fsds_nir_i_patch , & ! Output: [real(r8) (:) ] incident diffuse nir solar radiation (W/m**2) + fsrSF_nir_d => solarabs_inst%fsrSF_nir_d_patch , & ! Output: [real(r8) (:) ] snow-free reflected direct beam nir solar radiation (W/m**2) + fsrSF_nir_i => solarabs_inst%fsrSF_nir_i_patch , & ! Output: [real(r8) (:) ] snow-free reflected diffuse nir solar radiation (W/m**2) + fsrSF_nir_d_ln => solarabs_inst%fsrSF_nir_d_ln_patch, & ! Output: [real(r8) (:) ] snow-free reflected direct beam nir solar rad at local noon (W/m**2) + ssre_fsr_nir_d => solarabs_inst%ssre_fsr_nir_d_patch, & ! Output: [real(r8) (:) ] snow-free reflected direct beam nir solar radiation (W/m**2) + ssre_fsr_nir_i => solarabs_inst%ssre_fsr_nir_i_patch, & ! Output: [real(r8) (:) ] snow-free reflected diffuse nir solar radiation (W/m**2) + ssre_fsr_nir_d_ln=> solarabs_inst%ssre_fsr_nir_d_ln_patch,&!Output: [real(r8) (:) ] snow-free reflected direct beam nir solar rad at local noon (W/m**2) fsa_r => solarabs_inst%fsa_r_patch , & ! Output: [real(r8) (:) ] rural solar radiation absorbed (total) (W/m**2) - sub_surf_abs_SW => solarabs_inst%sub_surf_abs_SW_patch , & ! Output: [real(r8) (:) ] fraction of solar radiation absorbed below first snow layer (W/M**2) + sub_surf_abs_SW => solarabs_inst%sub_surf_abs_SW_patch,& ! Output: [real(r8) (:) ] fraction of solar radiation absorbed below first snow layer (W/M**2) parveg_ln => surfrad_inst%parveg_ln_patch , & ! Output: [real(r8) (:) ] absorbed par by vegetation at local noon (W/m**2) fsr_vis_d => surfrad_inst%fsr_vis_d_patch , & ! Output: [real(r8) (:) ] reflected direct beam vis solar radiation (W/m**2) fsr_vis_i => surfrad_inst%fsr_vis_i_patch , & ! Output: [real(r8) (:) ] reflected diffuse vis solar radiation (W/m**2) + fsrSF_vis_d => surfrad_inst%fsrSF_vis_d_patch , & ! Output: [real(r8) (:) ] snow-free reflected direct beam vis solar radiation (W/m**2) + fsrSF_vis_i => surfrad_inst%fsrSF_vis_i_patch , & ! Output: [real(r8) (:) ] snow-free reflected diffuse vis solar radiation (W/m**2) + ssre_fsr_vis_d => surfrad_inst%ssre_fsr_vis_d_patch , & ! Output: [real(r8) (:) ] snow-free reflected direct beam vis solar radiation (W/m**2) + ssre_fsr_vis_i => surfrad_inst%ssre_fsr_vis_i_patch , & ! Output: [real(r8) (:) ] snow-free reflected diffuse vis solar radiation (W/m**2) fsds_vis_i_ln => surfrad_inst%fsds_vis_i_ln_patch , & ! Output: [real(r8) (:) ] incident diffuse beam vis solar rad at local noon (W/m**2) fsr_vis_d_ln => surfrad_inst%fsr_vis_d_ln_patch , & ! Output: [real(r8) (:) ] reflected direct beam vis solar rad at local noon (W/m**2) + fsrSF_vis_d_ln => surfrad_inst%fsrSF_vis_d_ln_patch , & ! Output: [real(r8) (:) ] snow-free reflected direct beam vis solar rad at local noon (W/m**2) fsds_vis_d => surfrad_inst%fsds_vis_d_patch , & ! Output: [real(r8) (:) ] incident direct beam vis solar radiation (W/m**2) fsds_vis_i => surfrad_inst%fsds_vis_i_patch , & ! Output: [real(r8) (:) ] incident diffuse vis solar radiation (W/m**2) fsds_vis_d_ln => surfrad_inst%fsds_vis_d_ln_patch , & ! Output: [real(r8) (:) ] incident direct beam vis solar rad at local noon (W/m**2) @@ -589,7 +643,6 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & ! Determine seconds off current time step dtime = get_step_size() - call get_curr_date (year, month, day, secs) ! Initialize fluxes @@ -843,7 +896,12 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & rvis = albd(p,1)*forc_solad(g,1) + albi(p,1)*forc_solai(g,1) rnir = albd(p,2)*forc_solad(g,2) + albi(p,2)*forc_solai(g,2) fsr(p) = rvis + rnir - + if (use_SSRE) then + rvisSF = albdSF(p,1)*forc_solad(g,1) + albiSF(p,1)*forc_solai(g,1) + rnirSF = albdSF(p,2)*forc_solad(g,2) + albiSF(p,2)*forc_solai(g,2) + fsrSF(p) = rvisSF + rnirSF + ssre_fsr(p) = fsr(p)-fsrSF(p) + end if fsds_vis_d(p) = forc_solad(g,1) fsds_nir_d(p) = forc_solad(g,2) fsds_vis_i(p) = forc_solai(g,1) @@ -852,10 +910,18 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & fsr_nir_d(p) = albd(p,2)*forc_solad(g,2) fsr_vis_i(p) = albi(p,1)*forc_solai(g,1) fsr_nir_i(p) = albi(p,2)*forc_solai(g,2) - - local_secp1 = secs + nint((grc%londeg(g)/degpsec)/dtime)*dtime - local_secp1 = mod(local_secp1,isecspday) - if (local_secp1 == isecspday/2) then + if (use_SSRE) then + fsrSF_vis_d(p) = albdSF(p,1)*forc_solad(g,1) + fsrSF_nir_d(p) = albdSF(p,2)*forc_solad(g,2) + fsrSF_vis_i(p) = albiSF(p,1)*forc_solai(g,1) + fsrSF_nir_i(p) = albiSF(p,2)*forc_solai(g,2) + + ssre_fsr_vis_d(p) = fsrSF_vis_d(p)-fsr_vis_d(p) + ssre_fsr_nir_d(p) = fsrSF_nir_d(p)-fsr_nir_d(p) + ssre_fsr_vis_i(p) = fsrSF_vis_i(p)-fsr_vis_i(p) + ssre_fsr_nir_i(p) = fsrSF_nir_i(p)-fsr_nir_i(p) + end if + if ( is_near_local_noon( grc%londeg(g), deltasec=nint(dtime)/2 ) )then fsds_vis_d_ln(p) = forc_solad(g,1) fsds_nir_d_ln(p) = forc_solad(g,2) fsr_vis_d_ln(p) = albd(p,1)*forc_solad(g,1) @@ -870,7 +936,15 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & fsds_vis_i_ln(p) = spval parveg_ln(p) = spval end if - + if (use_SSRE) then + if ( is_near_local_noon( grc%londeg(g), deltasec=nint(dtime)/2 ) )then + fsrSF_vis_d_ln(p) = albdSF(p,1)*forc_solad(g,1) + fsrSF_nir_d_ln(p) = albdSF(p,2)*forc_solad(g,2) + else + fsrSF_vis_d_ln(p) = spval + fsrSF_nir_d_ln(p) = spval + end if + end if ! diagnostic variables (downwelling and absorbed radiation partitioning) for history files ! (OPTIONAL) c = patch%column(p) @@ -897,16 +971,14 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & endif end do + ! TODO: urban snow-free albedos: do fp = 1,num_urbanp p = filter_urbanp(fp) g = patch%gridcell(p) - local_secp1 = secs + nint((grc%londeg(g)/degpsec)/dtime)*dtime - local_secp1 = mod(local_secp1,isecspday) - - if(elai(p)==0.0_r8.and.fabd(p,1)>0._r8)then - if ( DEBUG ) write(iulog,*) 'absorption without LAI',elai(p),tlai(p),fabd(p,1),p - endif + if(elai(p)==0.0_r8.and.fabd(p,1)>0._r8)then + if ( DEBUG ) write(iulog,*) 'absorption without LAI',elai(p),tlai(p),fabd(p,1),p + endif ! Solar incident fsds_vis_d(p) = forc_solad(g,1) @@ -915,7 +987,7 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & fsds_nir_i(p) = forc_solai(g,2) ! Determine local noon incident solar - if (local_secp1 == noonsec) then + if ( is_near_local_noon( grc%londeg(g), deltasec=nint(dtime)/2 ) )then fsds_vis_d_ln(p) = forc_solad(g,1) fsds_nir_d_ln(p) = forc_solad(g,2) fsds_vis_i_ln(p) = forc_solai(g,1) @@ -936,7 +1008,7 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, & fsr_nir_i(p) = albi(p,2) * forc_solai(g,2) ! Determine local noon reflected solar - if (local_secp1 == noonsec) then + if ( is_near_local_noon( grc%londeg(g), deltasec=nint(dtime)/2 ) )then fsr_vis_d_ln(p) = fsr_vis_d(p) fsr_nir_d_ln(p) = fsr_nir_d(p) else diff --git a/src/biogeophys/UrbanFluxesMod.F90 b/src/biogeophys/UrbanFluxesMod.F90 index 0a9702f295..29f6b71aba 100644 --- a/src/biogeophys/UrbanFluxesMod.F90 +++ b/src/biogeophys/UrbanFluxesMod.F90 @@ -10,7 +10,7 @@ module UrbanFluxesMod use shr_log_mod , only : errMsg => shr_log_errMsg use decompMod , only : bounds_type use clm_varpar , only : numrad - use clm_varcon , only : isecspday, degpsec, namel + use clm_varcon , only : namel use clm_varctl , only : iulog use abortutils , only : endrun use UrbanParamsType , only : urbanparams_type @@ -68,8 +68,9 @@ subroutine UrbanFluxes (bounds, num_nourbanl, filter_nourbanl, use FrictionVelocityMod , only : FrictionVelocity, MoninObukIni, frictionvel_parms_inst use QSatMod , only : QSat use clm_varpar , only : maxpatch_urb, nlevurb, nlevgrnd - use clm_time_manager , only : get_curr_date, get_step_size, get_nstep - use HumanIndexMod , only : calc_human_stress_indices, Wet_Bulb, Wet_BulbS, HeatIndex, AppTemp, & + use clm_time_manager , only : get_step_size, get_nstep, is_near_local_noon + use HumanIndexMod , only : all_human_stress_indices, fast_human_stress_indices, & + Wet_Bulb, Wet_BulbS, HeatIndex, AppTemp, & swbgt, hmdex, dis_coi, dis_coiS, THIndex, & SwampCoolEff, KtoC, VaporPres ! @@ -174,9 +175,7 @@ subroutine UrbanFluxes (bounds, num_nourbanl, filter_nourbanl, real(r8) :: qflx_err(bounds%begl:bounds%endl) ! water vapor flux error (kg/m**2/s) real(r8) :: fwet_roof ! fraction of roof surface that is wet (-) real(r8) :: fwet_road_imperv ! fraction of impervious road surface that is wet (-) - integer :: local_secp1(bounds%begl:bounds%endl) ! seconds into current date in local time (sec) real(r8) :: dtime ! land model time step (sec) - integer :: year,month,day,secs ! calendar info for current time step logical :: found ! flag in search loop integer :: indexl ! index of first found in search loop integer :: nstep ! time step number @@ -318,7 +317,6 @@ subroutine UrbanFluxes (bounds, num_nourbanl, filter_nourbanl, ! Get current date dtime = get_step_size() - call get_curr_date (year, month, day, secs) ! Compute canyontop wind using Masson (2000) @@ -326,9 +324,6 @@ subroutine UrbanFluxes (bounds, num_nourbanl, filter_nourbanl, l = filter_urbanl(fl) g = lun%gridcell(l) - local_secp1(l) = secs + nint((grc%londeg(g)/degpsec)/dtime)*dtime - local_secp1(l) = mod(local_secp1(l),isecspday) - ! Error checks if (ht_roof(l) - z_d_town(l) <= z_0_town(l)) then @@ -868,36 +863,39 @@ subroutine UrbanFluxes (bounds, num_nourbanl, filter_nourbanl, rh_ref2m_u(p) = rh_ref2m(p) ! Human Heat Stress - if ( calc_human_stress_indices )then - + if ( all_human_stress_indices .or. fast_human_stress_indices )then call KtoC(t_ref2m(p), tc_ref2m(p)) call VaporPres(rh_ref2m(p), e_ref2m, vap_ref2m(p)) - call Wet_Bulb(t_ref2m(p), vap_ref2m(p), forc_pbot(g), rh_ref2m(p), q_ref2m(p), & - teq_ref2m(p), ept_ref2m(p), wb_ref2m(p)) call Wet_BulbS(tc_ref2m(p), rh_ref2m(p), wbt_ref2m(p)) call HeatIndex(tc_ref2m(p), rh_ref2m(p), nws_hi_ref2m(p)) call AppTemp(tc_ref2m(p), vap_ref2m(p), u10_clm(p), appar_temp_ref2m(p)) call swbgt(tc_ref2m(p), vap_ref2m(p), swbgt_ref2m(p)) call hmdex(tc_ref2m(p), vap_ref2m(p), humidex_ref2m(p)) - call dis_coi(tc_ref2m(p), wb_ref2m(p), discomf_index_ref2m(p)) call dis_coiS(tc_ref2m(p), rh_ref2m(p), wbt_ref2m(p), discomf_index_ref2mS(p)) - call THIndex(tc_ref2m(p), wb_ref2m(p), thic_ref2m(p), thip_ref2m(p)) - call SwampCoolEff(tc_ref2m(p), wb_ref2m(p), swmp80_ref2m(p), swmp65_ref2m(p)) + if ( all_human_stress_indices ) then + call Wet_Bulb(t_ref2m(p), vap_ref2m(p), forc_pbot(g), rh_ref2m(p), q_ref2m(p), & + teq_ref2m(p), ept_ref2m(p), wb_ref2m(p)) + call dis_coi(tc_ref2m(p), wb_ref2m(p), discomf_index_ref2m(p)) + call THIndex(tc_ref2m(p), wb_ref2m(p), thic_ref2m(p), thip_ref2m(p)) + call SwampCoolEff(tc_ref2m(p), wb_ref2m(p), swmp80_ref2m(p), swmp65_ref2m(p)) + end if - teq_ref2m_u(p) = teq_ref2m(p) - ept_ref2m_u(p) = ept_ref2m(p) - wb_ref2m_u(p) = wb_ref2m(p) wbt_ref2m_u(p) = wbt_ref2m(p) nws_hi_ref2m_u(p) = nws_hi_ref2m(p) appar_temp_ref2m_u(p) = appar_temp_ref2m(p) swbgt_ref2m_u(p) = swbgt_ref2m(p) humidex_ref2m_u(p) = humidex_ref2m(p) - discomf_index_ref2m_u(p) = discomf_index_ref2m(p) discomf_index_ref2mS_u(p) = discomf_index_ref2mS(p) - thic_ref2m_u(p) = thic_ref2m(p) - thip_ref2m_u(p) = thip_ref2m(p) - swmp80_ref2m_u(p) = swmp80_ref2m(p) - swmp65_ref2m_u(p) = swmp65_ref2m(p) + if ( all_human_stress_indices ) then + teq_ref2m_u(p) = teq_ref2m(p) + ept_ref2m_u(p) = ept_ref2m(p) + wb_ref2m_u(p) = wb_ref2m(p) + discomf_index_ref2m_u(p) = discomf_index_ref2m(p) + thic_ref2m_u(p) = thic_ref2m(p) + thip_ref2m_u(p) = thip_ref2m(p) + swmp80_ref2m_u(p) = swmp80_ref2m(p) + swmp65_ref2m_u(p) = swmp65_ref2m(p) + end if end if ! Variables needed by history tape diff --git a/src/biogeophys/UrbanRadiationMod.F90 b/src/biogeophys/UrbanRadiationMod.F90 index 102d2b6803..4e48693f2d 100644 --- a/src/biogeophys/UrbanRadiationMod.F90 +++ b/src/biogeophys/UrbanRadiationMod.F90 @@ -12,7 +12,7 @@ module UrbanRadiationMod use shr_log_mod , only : errMsg => shr_log_errMsg use decompMod , only : bounds_type use clm_varpar , only : numrad - use clm_varcon , only : isecspday, degpsec, namel + use clm_varcon , only : namel use clm_varctl , only : iulog use abortutils , only : endrun use UrbanParamsType , only : urbanparams_type @@ -61,7 +61,7 @@ subroutine UrbanRadiation (bounds , & use clm_varcon , only : spval, sb, tfrz use column_varcon , only : icol_road_perv, icol_road_imperv use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall - use clm_time_manager , only : get_curr_date, get_step_size + use clm_time_manager , only : get_step_size ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds @@ -83,10 +83,7 @@ subroutine UrbanRadiation (bounds , & ! ! !LOCAL VARIABLES: integer :: fp,fl,p,c,l,g ! indices - integer :: local_secp1 ! seconds into current date in local time real(r8) :: dtime ! land model time step (sec) - integer :: year,month,day ! temporaries (not used) - integer :: secs ! seconds into current date real(r8), parameter :: mpe = 1.e-06_r8 ! prevents overflow for division by zero real(r8), parameter :: snoem = 0.97_r8 ! snow emissivity (should use value from Biogeophysics1) @@ -250,7 +247,6 @@ subroutine UrbanRadiation (bounds , & end if dtime = get_step_size() - call get_curr_date (year, month, day, secs) ! Determine variables needed for history output and communication with atm ! Loop over urban patches in clump diff --git a/src/biogeophys/WaterfluxType.F90 b/src/biogeophys/WaterfluxType.F90 index 96e8c92a09..905fe9471f 100644 --- a/src/biogeophys/WaterfluxType.F90 +++ b/src/biogeophys/WaterfluxType.F90 @@ -599,7 +599,7 @@ subroutine InitAccVars (this, bounds) if (use_fun) then call extract_accum_field ('AnnET', rbufslp, nstep) - this%qflx_evap_tot_col(begc:endc) = rbufslp(begc:endc) + this%AnnET(begc:endc) = rbufslp(begc:endc) end if deallocate(rbufslp) @@ -681,8 +681,6 @@ subroutine InitCold(this, bounds) ! the sake of columns outside this filter this%qflx_ice_runoff_xs_col(bounds%begc:bounds%endc) = 0._r8 - this%AnnEt(bounds%begc:bounds%endc) = 0._r8 - ! needed for CNNLeaching do c = bounds%begc, bounds%endc l = col%landunit(c) @@ -730,16 +728,6 @@ subroutine Restart(this, bounds, ncid, flag) this%qflx_snow_drain_col(bounds%begc:bounds%endc) = 0._r8 endif - - call restartvar(ncid=ncid, flag=flag, varname='AnnET', xtype=ncd_double, & - dim1name='column', & - long_name='Annual ET ', units='mm/s', & - interpinic_flag='interp', readvar=readvar, data=this%AnnET) - if (flag == 'read' .and. .not. readvar) then - ! initial run, not restart: initialize qflx_snow_drain to zero - this%AnnET(bounds%begc:bounds%endc) = 0._r8 - endif - call this%qflx_liq_dynbal_dribbler%Restart(bounds, ncid, flag) call this%qflx_ice_dynbal_dribbler%Restart(bounds, ncid, flag) diff --git a/src/biogeophys/test/Balance_test/CMakeLists.txt b/src/biogeophys/test/Balance_test/CMakeLists.txt new file mode 100644 index 0000000000..2d0cd75c00 --- /dev/null +++ b/src/biogeophys/test/Balance_test/CMakeLists.txt @@ -0,0 +1,4 @@ +create_pFUnit_test(balance test_balance_exe + "test_Balance.pf" "") + +target_link_libraries(test_balance_exe clm csm_share esmf_wrf_timemgr) diff --git a/src/biogeophys/test/Balance_test/test_Balance.pf b/src/biogeophys/test/Balance_test/test_Balance.pf new file mode 100644 index 0000000000..3d07385ffb --- /dev/null +++ b/src/biogeophys/test/Balance_test/test_Balance.pf @@ -0,0 +1,107 @@ +module test_balance + + ! Some tests of the balance check system + + use pfunit_mod + + use shr_kind_mod , only : r8 => shr_kind_r8 + use unittestTimeManagerMod, only : unittest_timemgr_setup, unittest_timemgr_teardown + use unittestSubgridMod + use ncdio_pio ! use the fake version of this module + use BalanceCheckMod + use unittestUtils , only : endrun_msg + + implicit none + save + + @TestCase + type, extends(TestCase) :: TestBalance + contains + procedure :: setUp + procedure :: tearDown + end type TestBalance + +contains + + subroutine setUp(this) + class(TestBalance), intent(inout) :: this + + end subroutine setUp + + subroutine tearDown(this) + class(TestBalance), intent(inout) :: this + + call unittest_timemgr_teardown() + call BalanceCheckClean() + + end subroutine tearDown + + @Test + subroutine test_balance_init( this ) + class(TestBalance), intent(inout) :: this + integer :: dtime + integer :: nskip + + dtime = 1800 + call unittest_timemgr_setup(dtime=dtime) + call BalanceCheckInit() + nskip = GetBalanceCheckSkipSteps() + @assertEqual( 2, nskip, message="Ensure standard balance check is 2 time-steps" ) + end subroutine test_balance_init + + @Test + subroutine test_balance_longstep( this ) + class(TestBalance), intent(inout) :: this + integer :: dtime + integer :: nskip + + dtime = 7200 + call unittest_timemgr_setup(dtime=dtime) + call BalanceCheckInit() + nskip = GetBalanceCheckSkipSteps() + @assertEqual( 2, nskip, message="Ensure even with a long time-step skip is 2 time-steps" ) + end subroutine test_balance_longstep + + @Test + subroutine test_balance_300sec( this ) + class(TestBalance), intent(inout) :: this + integer :: dtime + integer :: nskip + + dtime = 300 + call unittest_timemgr_setup(dtime=dtime) + call BalanceCheckInit() + nskip = GetBalanceCheckSkipSteps() + @assertEqual( 12, nskip, message="Check skip length for 300 sec time-step" ) + end subroutine test_balance_300sec + + @Test + subroutine test_balance_Fail( this ) + class(TestBalance), intent(inout) :: this + integer :: dtime + integer :: nskip + character(len=256) :: expected_msg + + dtime = 1800 + call unittest_timemgr_setup(dtime=dtime) + nskip = GetBalanceCheckSkipSteps() + expected_msg = endrun_msg( & + "ERROR:: GetBalanceCheckSkipSteps called before BalanceCheckInit" ) + @assertExceptionRaised(expected_msg) + call BalanceCheckInit() + end subroutine test_balance_Fail + + @Test + subroutine test_balance_shortstep( this ) + class(TestBalance), intent(inout) :: this + integer :: dtime + integer :: nskip + + dtime = 36 + call unittest_timemgr_setup(dtime=dtime) + call BalanceCheckInit() + nskip = GetBalanceCheckSkipSteps() + @assertEqual( 100, nskip, message="Ensure with a short step correct number of skip steps is done" ) + end subroutine test_balance_shortstep + +end module test_balance diff --git a/src/biogeophys/test/CMakeLists.txt b/src/biogeophys/test/CMakeLists.txt index 75b77d5391..16007e5e70 100644 --- a/src/biogeophys/test/CMakeLists.txt +++ b/src/biogeophys/test/CMakeLists.txt @@ -2,4 +2,5 @@ add_subdirectory(Daylength_test) add_subdirectory(Irrigation_test) add_subdirectory(HumanStress_test) add_subdirectory(SnowHydrology_test) +add_subdirectory(Balance_test) add_subdirectory(TotalWaterAndHeat_test) diff --git a/src/biogeophys/test/Irrigation_test/test_irrigation.pf b/src/biogeophys/test/Irrigation_test/test_irrigation.pf index 3c599aac02..742edf19d8 100644 --- a/src/biogeophys/test/Irrigation_test/test_irrigation.pf +++ b/src/biogeophys/test/Irrigation_test/test_irrigation.pf @@ -4,6 +4,9 @@ module test_irrigation use pfunit_mod use unittestSubgridMod + use unittestTimeManagerMod, only : unittest_timemgr_setup, unittest_timemgr_teardown + use unittestTimeManagerMod, only : unittest_timemgr_set_curr_date + use clm_time_manager, only: advance_timestep use IrrigationMod, only : irrigation_type, irrigation_params_type use shr_kind_mod, only : r8 => shr_kind_r8 use clm_varpar, only : nlevsoi, nlevgrnd @@ -20,6 +23,7 @@ module test_irrigation real(r8), parameter :: tol = 1.e-13_r8 integer , parameter :: dtime = 1800 ! model time step, seconds + integer , parameter :: irrig_start = 21600 @TestCase type, extends(TestCase) :: TestIrrigation @@ -39,8 +43,6 @@ module test_irrigation real(r8), allocatable :: relsat_target(:,:) real(r8), allocatable :: volr(:) - ! Previous model time - integer :: time_prev contains procedure :: setUp procedure :: tearDown @@ -68,6 +70,11 @@ contains subroutine setUp(this) class(TestIrrigation), intent(inout) :: this + ! Setup time manager + call unittest_timemgr_setup(dtime=dtime) + ! Set time to one time-step ahead of start time, as irrigation uses the previous time step + call unittest_timemgr_set_curr_date(yr=5, mon=1, day=1, tod=irrig_start+dtime) + ! Need to set nlevgrnd before doing the subgrid setup (because it is needed when ! allocating the col object). So we must do this before setupEnvironment, because ! that assumes that the subgrid setup has already been done. @@ -83,6 +90,7 @@ contains subroutine tearDown(this) class(TestIrrigation), intent(inout) :: this + call unittest_timemgr_teardown() call this%irrigation%Clean() call this%teardownEnvironment() call unittest_subgrid_teardown() @@ -170,7 +178,7 @@ contains ! Set parameters this%irrigation_params = irrigation_params_type( & irrig_min_lai = 0.0_r8, & - irrig_start_time = 21600, & + irrig_start_time = irrig_start, & irrig_length = 14400, & irrig_target_smp = -3400._r8, & irrig_depth = irrig_depth, & @@ -196,9 +204,8 @@ contains end do end do - ! Set time_prev to the irrig_start_time minus 1 hour (since we're using a longitude - ! about 1 hour east of 0Z) - this%time_prev = this%irrigation_params%irrig_start_time - 3600 + ! Set time to one time-step ahead of start time, as irrigation uses the previous time step + call unittest_timemgr_set_curr_date(yr=5, mon=1, day=1, tod=irrig_start+dtime) call this%irrigation%InitForTesting(bounds, this%irrigation_params, dtime, & this%relsat_wilting_point, this%relsat_target) @@ -238,8 +245,8 @@ contains col%nbedrock(c) = nlevsoi end do - ! slightly greater than 1 hour offset - grc%londeg(:) = 15.1_r8 + ! Use longitude along Greenich so don't have to calculate offsets for longitudes (that's calculated in clm_time_manager) + grc%londeg(:) = 0.0_r8 grc%area(:) = 10.0_r8 @@ -325,7 +332,6 @@ contains bounds=bounds, & num_exposedvegp = this%numf, & filter_exposedvegp = this%filter, & - time_prev = this%time_prev, & elai = this%elai, & t_soisno = this%t_soisno, & eff_porosity = this%eff_porosity, & @@ -480,8 +486,8 @@ contains ! Setup call this%setupSinglePatch() call this%setupIrrigation() - ! Set previous time to be one time step before the time when we would start irrigating - this%time_prev = this%time_prev - dtime + ! Set current time to start time, as irrigation uses the time step befor ethat + call unittest_timemgr_set_curr_date(yr=5, mon=1, day=1, tod=irrig_start) ! Call irrigation routines call this%calculateAndApplyIrrigation() @@ -555,27 +561,6 @@ contains @assertEqual(expected, this%irrigation%qflx_irrig_patch(bounds%begp), tolerance=tol) end subroutine limited_irrigation_for_limiting_volr - @Test - subroutine irrigation_should_happen_for_big_longitude(this) - use GridcellType, only : grc - class(TestIrrigation), intent(inout) :: this - - ! Setup - call this%setupSinglePatch() - call this%setupIrrigation() - ! Use a big longitude and a time_prev that should lead to irrigation at that longitude - ! The main point of this is to test the modulo in the local_time calculation - grc%londeg(:) = 359.9_r8 - this%time_prev = this%irrigation_params%irrig_start_time + dtime - - ! Call irrigation routines - call this%calculateAndApplyIrrigation() - - ! Check result - @assertTrue(this%irrigation%qflx_irrig_patch(bounds%begp) > 0._r8) - - end subroutine irrigation_should_happen_for_big_longitude - @Test subroutine irrigation_continues_at_same_rate_for_multiple_time_steps(this) class(TestIrrigation), intent(inout) :: this @@ -592,7 +577,7 @@ contains call this%calculateAndApplyIrrigation() call this%computeDeficits(deficits) expected = sum(deficits(bounds%begp,1:nlevsoi)) / this%irrigation_params%irrig_length - this%time_prev = this%time_prev + dtime + call advance_timestep() this%h2osoi_liq = 100._r8 call this%calculateAndApplyIrrigation() @@ -618,7 +603,7 @@ contains ! steps do time = 1, expected_num_time_steps call this%calculateAndApplyIrrigation() - this%time_prev = this%time_prev + dtime + call advance_timestep() end do @assertTrue(this%irrigation%qflx_irrig_patch(bounds%begp) > 0._r8) @@ -637,26 +622,24 @@ contains class(TestIrrigation), intent(inout) :: this real(r8), allocatable :: deficits(:,:) real(r8) :: expected - integer :: time_prev_orig integer :: time integer :: expected_num_time_steps ! Setup call this%setupSinglePatch() call this%setupIrrigation() - time_prev_orig = this%time_prev ! Call irrigation routines for long enough that irrigation should go to 0 expected_num_time_steps = this%irrigation_params%irrig_length / dtime do time = 1, expected_num_time_steps + 1 call this%calculateAndApplyIrrigation() - this%time_prev = this%time_prev + dtime + call advance_timestep() end do ! The following assertion is mainly here to make sure the test is working as intended @assertEqual(0._r8, this%irrigation%qflx_irrig_patch(bounds%begp)) ! Now reset time, change soil moisture, and make sure that irrigation happens as expected - this%time_prev = time_prev_orig + call unittest_timemgr_set_curr_date(yr=5, mon=1, day=1, tod=irrig_start+dtime) this%h2osoi_liq(:,:) = 100._r8 call this%calculateAndApplyIrrigation() call this%computeDeficits(deficits) diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90 index b1a4a9beb8..c66112c0bc 100644 --- a/src/main/clm_driver.F90 +++ b/src/main/clm_driver.F90 @@ -131,10 +131,6 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro integer :: mon ! month (1, ..., 12) integer :: day ! day of month (1, ..., 31) integer :: sec ! seconds of the day - integer :: yr_prev ! year (0, ...) at start of timestep - integer :: mon_prev ! month (1, ..., 12) at start of timestep - integer :: day_prev ! day of month (1, ..., 31) at start of timestep - integer :: sec_prev ! seconds of the day at start of timestep character(len=256) :: filer ! restart file name integer :: ier ! error code logical :: need_glacier_initialization ! true if we need to initialize glacier areas in this time step @@ -393,9 +389,6 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro ! snow accumulation exceeds 10 mm. ! ============================================================================ - ! Get time as of beginning of time step - call get_prev_date(yr_prev, mon_prev, day_prev, sec_prev) - !$OMP PARALLEL DO PRIVATE (nc,l,c, bounds_clump, downreg_patch, leafn_patch, agnpp_patch, bgnpp_patch, annsum_npp_patch, rr_patch, froot_carbon, croot_carbon) do nc = 1,nclumps call get_clump_bounds(nc, bounds_clump) @@ -596,7 +589,6 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro bounds = bounds_clump, & num_exposedvegp = filter(nc)%num_exposedvegp, & filter_exposedvegp = filter(nc)%exposedvegp, & - time_prev = sec_prev, & elai = canopystate_inst%elai_patch(bounds_clump%begp:bounds_clump%endp), & t_soisno = temperature_inst%t_soisno_col(bounds_clump%begc:bounds_clump%endc , 1:nlevgrnd), & eff_porosity = soilstate_inst%eff_porosity_col(bounds_clump%begc:bounds_clump%endc, 1:nlevgrnd), & @@ -929,6 +921,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro call BalanceCheck(bounds_clump, & atm2lnd_inst, solarabs_inst, waterflux_inst, & waterstate_inst, irrigation_inst, glacier_smb_inst, & + surfalb_inst, & energyflux_inst, canopystate_inst) call t_stopf('balchk') diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index 42f71ec526..692ba83ab8 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -289,6 +289,7 @@ subroutine initialize2( ) use NutrientCompetitionFactoryMod, only : create_nutrient_competition_method use controlMod , only : NLFilename use clm_instMod , only : clm_fates + use BalanceCheckMod , only : BalanceCheckInit ! ! !ARGUMENTS ! @@ -372,6 +373,9 @@ subroutine initialize2( ) call t_stopf('init_orbd') call InitDaylength(bounds_proc, declin=declin, declinm1=declinm1, obliquity=obliqr) + + ! Initialize Balance checking (after time-manager) + call BalanceCheckInit() ! History file variables diff --git a/src/main/clm_instMod.F90 b/src/main/clm_instMod.F90 index 1b1b44adba..c798e330ef 100644 --- a/src/main/clm_instMod.F90 +++ b/src/main/clm_instMod.F90 @@ -186,6 +186,7 @@ subroutine clm_instInit(bounds) use accumulMod , only : print_accum_fields use SoilWaterRetentionCurveFactoryMod , only : create_soil_water_retention_curve use decompMod , only : get_proc_bounds + use BalanceCheckMod , only : GetBalanceCheckSkipSteps ! ! !ARGUMENTS type(bounds_type), intent(in) :: bounds ! processor bounds @@ -391,7 +392,7 @@ subroutine clm_instInit(bounds) end if ! end of if use_cn ! Note - always call Init for bgc_vegetation_inst: some pieces need to be initialized always - call bgc_vegetation_inst%Init(bounds, nlfilename) + call bgc_vegetation_inst%Init(bounds, nlfilename, GetBalanceCheckSkipSteps() ) if (use_cn .or. use_fates) then call crop_inst%Init(bounds) diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index 938155c5dd..f7e5a793ca 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -308,6 +308,10 @@ module clm_varctl ! FATES !---------------------------------------------------------- character(len=fname_len), public :: fates_paramfile = ' ' + !---------------------------------------------------------- + ! SSRE diagnostic + !---------------------------------------------------------- + logical, public :: use_SSRE = .false. ! flag for SSRE diagnostic !---------------------------------------------------------- ! Migration of CPP variables diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index ae9c2fcafe..a78f56785e 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -250,6 +250,9 @@ subroutine control_init( ) ! max number of plant functional types in naturally vegetated landunit namelist /clm_inparm/ maxpatch_pft + ! flag for SSRE diagnostic + namelist /clm_inparm/ use_SSRE + namelist /clm_inparm/ & use_lch4, use_nitrif_denitrif, use_vertsoilc, use_extralakelayers, & use_vichydro, use_century_decomp, use_cn, use_cndv, use_crop, use_fertilizer, use_ozone, & @@ -587,6 +590,7 @@ subroutine control_spmd() call mpi_bcast (use_vancouver, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (use_mexicocity, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (use_noio, 1, MPI_LOGICAL, 0, mpicom, ier) + call mpi_bcast (use_SSRE, 1, MPI_LOGICAL, 0, mpicom, ier) ! initial file variables call mpi_bcast (nrevsn, len(nrevsn), MPI_CHARACTER, 0, mpicom, ier) @@ -799,7 +803,7 @@ subroutine control_print () write(iulog,*) ' use_vancouver = ', use_vancouver write(iulog,*) ' use_mexicocity = ', use_mexicocity write(iulog,*) ' use_noio = ', use_noio - + write(iulog,*) ' use_SSRE = ', use_SSRE write(iulog,*) 'input data files:' write(iulog,*) ' PFT physiology and parameters file = ',trim(paramfile) if (fsurdat == ' ') then diff --git a/src/main/surfrdMod.F90 b/src/main/surfrdMod.F90 index ef593b1a6e..a7808f1e6e 100644 --- a/src/main/surfrdMod.F90 +++ b/src/main/surfrdMod.F90 @@ -188,11 +188,12 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename) call domain_init(ldomain, isgrid2d=isgrid2d, ni=ni, nj=nj, nbeg=begg, nend=endg) ! Determine type of file - old style grid file or new style domain file - call check_var(ncid=ncid, varname='LONGXY', vardesc=vardesc, readvar=readvar) - if (readvar) istype_domain = .false. - call check_var(ncid=ncid, varname='xc', vardesc=vardesc, readvar=readvar) - if (readvar) istype_domain = .true. + if (readvar)then + istype_domain = .true. + else + istype_domain = .false. + end if ! Read in area, lon, lat @@ -211,33 +212,15 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename) dim1name=grlnd, readvar=readvar) if (.not. readvar) call endrun( msg=' ERROR: yc NOT on file'//errMsg(sourcefile, __LINE__)) else - call ncd_io(ncid=ncid, varname= 'AREA', flag='read', data=ldomain%area, & - dim1name=grlnd, readvar=readvar) - if (.not. readvar) call endrun( msg=' ERROR: AREA NOT on file'//errMsg(sourcefile, __LINE__)) - - call ncd_io(ncid=ncid, varname= 'LONGXY', flag='read', data=ldomain%lonc, & - dim1name=grlnd, readvar=readvar) - if (.not. readvar) call endrun( msg=' ERROR: LONGXY NOT on file'//errMsg(sourcefile, __LINE__)) - - call ncd_io(ncid=ncid, varname= 'LATIXY', flag='read', data=ldomain%latc, & - dim1name=grlnd, readvar=readvar) - if (.not. readvar) call endrun( msg=' ERROR: LATIXY NOT on file'//errMsg(sourcefile, __LINE__)) + call endrun( msg=" ERROR: can no longer read non domain files" ) end if if (isgrid2d) then allocate(rdata2d(ni,nj), lon1d(ni), lat1d(nj)) - if (istype_domain) then - vname = 'xc' - else - vname = 'LONGXY' - end if + if (istype_domain) vname = 'xc' call ncd_io(ncid=ncid, varname=trim(vname), data=rdata2d, flag='read', readvar=readvar) lon1d(:) = rdata2d(:,1) - if (istype_domain) then - vname = 'yc' - else - vname = 'LATIXY' - end if + if (istype_domain) vname = 'yc' call ncd_io(ncid=ncid, varname=trim(vname), data=rdata2d, flag='read', readvar=readvar) lat1d(:) = rdata2d(1,:) deallocate(rdata2d) @@ -255,25 +238,21 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename) ! where (ldomain%latc < -90.0_r8) ldomain%latc = -90.0_r8 ! where (ldomain%latc > 90.0_r8) ldomain%latc = 90.0_r8 endif + if ( any(ldomain%lonc < 0.0_r8) )then + call endrun( msg=' ERROR: lonc is negative and currently can NOT be (see https://github.com/ESCOMP/ctsm/issues/507)' & + //errMsg(sourcefile, __LINE__)) + endif - call ncd_io(ncid=ncid, varname='LANDMASK', flag='read', data=ldomain%mask, & + call ncd_io(ncid=ncid, varname='mask', flag='read', data=ldomain%mask, & dim1name=grlnd, readvar=readvar) if (.not. readvar) then - call ncd_io(ncid=ncid, varname='mask', flag='read', data=ldomain%mask, & - dim1name=grlnd, readvar=readvar) - if (.not. readvar) then - call endrun( msg=' ERROR: LANDMASK NOT on fracdata file'//errMsg(sourcefile, __LINE__)) - end if + call endrun( msg=' ERROR: LANDMASK NOT on fracdata file'//errMsg(sourcefile, __LINE__)) end if - call ncd_io(ncid=ncid, varname='LANDFRAC', flag='read', data=ldomain%frac, & + call ncd_io(ncid=ncid, varname='frac', flag='read', data=ldomain%frac, & dim1name=grlnd, readvar=readvar) if (.not. readvar) then - call ncd_io(ncid=ncid, varname='frac', flag='read', data=ldomain%frac, & - dim1name=grlnd, readvar=readvar) - if (.not. readvar) then - call endrun( msg=' ERROR: LANDFRAC NOT on fracdata file'//errMsg(sourcefile, __LINE__)) - end if + call endrun( msg=' ERROR: LANDFRAC NOT on fracdata file'//errMsg(sourcefile, __LINE__)) end if call ncd_pio_closefile(ncid) diff --git a/src/soilbiogeochem/SoilBiogeochemNitrogenFluxType.F90 b/src/soilbiogeochem/SoilBiogeochemNitrogenFluxType.F90 index 49371a8c9d..d60b66595d 100644 --- a/src/soilbiogeochem/SoilBiogeochemNitrogenFluxType.F90 +++ b/src/soilbiogeochem/SoilBiogeochemNitrogenFluxType.F90 @@ -326,10 +326,12 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='symbiotic/asymbiotic N fixation to soil mineral N', & ptr_col=this%nfix_to_sminn_col, default=default) - this%ffix_to_sminn_col(begc:endc) = spval - call hist_addfld1d (fname='FFIX_TO_SMINN', units='gN/m^2/s', & - avgflag='A', long_name='free living N fixation to soil mineral N', & - ptr_col=this%ffix_to_sminn_col, default=default) + if ( use_fun )then + this%ffix_to_sminn_col(begc:endc) = spval + call hist_addfld1d (fname='FFIX_TO_SMINN', units='gN/m^2/s', & + avgflag='A', long_name='free living N fixation to soil mineral N', & + ptr_col=this%ffix_to_sminn_col, default='active') + end if do l = 1, ndecomp_cascade_transitions ! vertically integrated fluxes diff --git a/src/unit_test_stubs/main/GetGlobalValuesMod_stub.F90 b/src/unit_test_stubs/main/GetGlobalValuesMod_stub.F90 index 1e61fc17be..d5b7258acb 100644 --- a/src/unit_test_stubs/main/GetGlobalValuesMod_stub.F90 +++ b/src/unit_test_stubs/main/GetGlobalValuesMod_stub.F90 @@ -6,6 +6,7 @@ module GetGlobalValuesMod implicit none public :: GetGlobalWrite + public :: GetGlobalIndex contains @@ -16,4 +17,21 @@ subroutine GetGlobalWrite(decomp_index, clmlevel) ! do nothing end subroutine GetGlobalWrite + !----------------------------------------------------------------------- + integer function GetGlobalIndex(decomp_index, clmlevel) + + !---------------------------------------------------------------- + ! Description + ! Determine global index space value for target point at given clmlevel + ! + ! Uses: + ! + ! Arguments + integer , intent(in) :: decomp_index + character(len=*) , intent(in) :: clmlevel + + ! De essentially nothing, just set the index to a negative value to signal it's not real + GetGlobalIndex = -1 + end function GetGlobalIndex + end module GetGlobalValuesMod diff --git a/src/unit_test_stubs/main/ncdio_pio_fake.F90.in b/src/unit_test_stubs/main/ncdio_pio_fake.F90.in index 0d88830b88..47f7fabec1 100644 --- a/src/unit_test_stubs/main/ncdio_pio_fake.F90.in +++ b/src/unit_test_stubs/main/ncdio_pio_fake.F90.in @@ -10,6 +10,7 @@ module ncdio_pio use shr_assert_mod , only : shr_assert use ncdio_var, only : ncdio_var_type use ncdio_dim, only : ncdio_dim_type + use abortutils,only : endrun ! !PUBLIC TYPES: implicit none @@ -55,6 +56,7 @@ module ncdio_pio public :: ncd_inqdid ! stub: inquire dimension id public :: ncd_inqvdlen ! stub: inquire size of a dimension public :: ncd_inqdlen ! stub: inquire size of a dimension + public :: ncd_defvar ! define variables public :: ncd_inqfdims ! stub: inquire file dimensions public :: ncd_getatt ! stub: get attribute public :: ncd_putatt ! stub: put attribute @@ -70,6 +72,11 @@ module ncdio_pio ! ! !PRIVATE MEMBER FUNCTIONS: + + interface ncd_defvar + module procedure ncd_defvar_bynf + module procedure ncd_defvar_bygrid + end interface private :: ncd_get_variable_index ! return the index of a given variable @@ -326,6 +333,77 @@ contains ! Stubs for the actual ncdio_pio functionality (do nothing) ! ====================================================================== + !----------------------------------------------------------------------- + subroutine ncd_defvar_bynf(ncid, varname, xtype, ndims, dimid, varid, & + long_name, units, cell_method, missing_value, fill_value, & + imissing_value, ifill_value, comment, flag_meanings, & + flag_values, nvalid_range ) + ! + ! !DESCRIPTION: + ! Define a netcdf variable + ! + ! !ARGUMENTS: + class(file_desc_t) , intent(inout) :: ncid ! netcdf file id + character(len=*) , intent(in) :: varname ! variable name + integer , intent(in) :: xtype ! external type + integer , intent(in) :: ndims ! number of dims + integer , intent(inout) :: varid ! returned var id + integer , intent(in), optional :: dimid(:) ! dimids + character(len=*) , intent(in), optional :: long_name ! attribute + character(len=*) , intent(in), optional :: units ! attribute + character(len=*) , intent(in), optional :: cell_method ! attribute + character(len=*) , intent(in), optional :: comment ! attribute + character(len=*) , intent(in), optional :: flag_meanings(:) ! attribute + real(r8) , intent(in), optional :: missing_value ! attribute for real + real(r8) , intent(in), optional :: fill_value ! attribute for real + integer , intent(in), optional :: imissing_value ! attribute for int + integer , intent(in), optional :: ifill_value ! attribute for int + integer , intent(in), optional :: flag_values(:) ! attribute for int + integer , intent(in), optional :: nvalid_range(2) ! attribute for int + + ! Do nothing + call endrun( 'ncd_defvar_bynf should not be actually called for unit tests' ) + + end subroutine ncd_defvar_bynf + + !----------------------------------------------------------------------- + subroutine ncd_defvar_bygrid(ncid, varname, xtype, & + dim1name, dim2name, dim3name, dim4name, dim5name, & + long_name, units, cell_method, missing_value, fill_value, & + imissing_value, ifill_value, switchdim, comment, & + flag_meanings, flag_values, nvalid_range ) + ! + ! !DESCRIPTION: + ! Define a netcdf variable + ! + ! !ARGUMENTS: + class(file_desc_t) , intent(inout) :: ncid ! netcdf file id + character(len=*) , intent(in) :: varname ! variable name + integer , intent(in) :: xtype ! external type + character(len=*) , intent(in), optional :: dim1name ! dimension name + character(len=*) , intent(in), optional :: dim2name ! dimension name + character(len=*) , intent(in), optional :: dim3name ! dimension name + character(len=*) , intent(in), optional :: dim4name ! dimension name + character(len=*) , intent(in), optional :: dim5name ! dimension name + character(len=*) , intent(in), optional :: long_name ! attribute + character(len=*) , intent(in), optional :: units ! attribute + character(len=*) , intent(in), optional :: cell_method ! attribute + character(len=*) , intent(in), optional :: comment ! attribute + character(len=*) , intent(in), optional :: flag_meanings(:) ! attribute + real(r8) , intent(in), optional :: missing_value ! attribute for real + real(r8) , intent(in), optional :: fill_value ! attribute for real + integer , intent(in), optional :: imissing_value ! attribute for int + integer , intent(in), optional :: ifill_value ! attribute for int + logical , intent(in), optional :: switchdim ! true=> permute dim1 and dim2 for output + integer , intent(in), optional :: flag_values(:) ! attribute for int + integer , intent(in), optional :: nvalid_range(2) ! attribute for int + + ! Do nothing + call endrun( 'ncd_defvar_bygrid should not be actually called for unit tests' ) + + end subroutine ncd_defvar_bygrid + + !----------------------------------------------------------------------- subroutine ncd_pio_openfile(file, fname, mode) ! diff --git a/src/utils/clm_time_manager.F90 b/src/utils/clm_time_manager.F90 index ec62c9f071..889258fd06 100644 --- a/src/utils/clm_time_manager.F90 +++ b/src/utils/clm_time_manager.F90 @@ -1,5 +1,6 @@ module clm_time_manager +#include "shr_assert.h" use shr_kind_mod, only: r8 => shr_kind_r8 use shr_sys_mod , only: shr_sys_abort use spmdMod , only: masterproc @@ -42,6 +43,8 @@ module clm_time_manager get_curr_yearfrac, &! return the fractional position in the current year, as of the end of the current timestep get_prev_yearfrac, &! return the fractional position in the current year, as of the beginning of the current timestep get_rest_date, &! return the date from the restart file + get_local_timestep_time, &! return the local time for the input longitude to the nearest time-step + get_local_time, &! return the local time for the input longitude set_nextsw_cday, &! set the next radiation calendar day is_first_step, &! return true on first step of initial run is_first_restart_step, &! return true on first step of restart or branch run @@ -53,6 +56,7 @@ module clm_time_manager is_end_curr_year, &! return true on last timestep in current year is_last_step, &! return true on last timestep is_perpetual, &! return true if perpetual calendar is in use + is_near_local_noon, &! return true if near local noon is_restart, &! return true if this is a restart run update_rad_dtime, &! track radiation interval via nstep update_DA_nstep, &! update the Data Assimulation time step @@ -821,7 +825,7 @@ subroutine get_clock( clock ) type(ESMF_Time) :: start_date, stop_date, ref_date integer :: rc - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet( tm_clock, timeStep=step_size, startTime=start_date, & stoptime=stop_date, reftime=ref_date, rc=rc ) @@ -842,7 +846,7 @@ function get_curr_ESMF_Time( ) character(len=*), parameter :: sub = 'clm::get_curr_ESMF_Time' integer :: rc - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet( tm_clock, currTime=get_curr_ESMF_Time, rc=rc ) call chkrc(rc, sub//': error return from ESMF_ClockGet') @@ -859,7 +863,7 @@ integer function get_step_size() type(ESMF_TimeInterval) :: step_size ! timestep size integer :: rc - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet(tm_clock, timeStep=step_size, rc=rc) call chkrc(rc, sub//': error return from ESMF_ClockGet') @@ -917,7 +921,7 @@ integer function get_rad_step_size() character(len=*), parameter :: sub = 'clm::get_rad_step_size' - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return if (nstep_rad_prev == uninit_int ) then get_rad_step_size=get_step_size() @@ -937,7 +941,7 @@ integer function get_nstep() integer :: rc integer(ESMF_KIND_I8) :: step_no - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet(tm_clock, advanceCount=step_no, rc=rc) call chkrc(rc, sub//': error return from ESMF_ClockGet') @@ -982,7 +986,7 @@ subroutine get_curr_date(yr, mon, day, tod, offset) type(ESMF_TimeInterval) :: off !----------------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet( tm_clock, currTime=date, rc=rc ) call chkrc(rc, sub//': error return from ESMF_ClockGet') @@ -1028,7 +1032,7 @@ subroutine get_perp_date(yr, mon, day, tod, offset) type(ESMF_TimeInterval) :: DelTime !----------------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet( tm_clock, currTime=date, rc=rc ) ! Get time of day add it to perpetual date @@ -1072,7 +1076,7 @@ subroutine get_prev_date(yr, mon, day, tod) type(ESMF_Time) :: date !----------------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet(tm_clock, prevTime=date, rc=rc ) call chkrc(rc, sub//': error return from ESMF_ClockGet') @@ -1101,7 +1105,7 @@ subroutine get_start_date(yr, mon, day, tod) type(ESMF_Time) :: date !----------------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet(tm_clock, startTime=date, rc=rc) call chkrc(rc, sub//': error return from ESMF_ClockGet') @@ -1127,7 +1131,7 @@ integer function get_driver_start_ymd( tod ) character(len=*), parameter :: sub = 'clm::get_driver_start_ymd' !----------------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return if ( start_ymd == uninit_int )then call shr_sys_abort( sub//': error driver start date is NOT set yet' ) @@ -1164,7 +1168,7 @@ subroutine get_ref_date(yr, mon, day, tod) type(ESMF_Time) :: date !----------------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet(tm_clock, refTime=date, rc=rc) call chkrc(rc, sub//': error return from ESMF_ClockGet') @@ -1193,7 +1197,7 @@ subroutine get_curr_time(days, seconds) type(ESMF_TimeInterval) :: diff !----------------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet( tm_clock, currTime=cdate, rc=rc ) call chkrc(rc, sub//': error return from ESMF_ClockGet') @@ -1227,7 +1231,7 @@ subroutine get_prev_time(days, seconds) type(ESMF_TimeInterval) :: diff !----------------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet(tm_clock, prevTime=date, rc=rc ) call chkrc(rc, sub//': error return from ESMF_ClockGet for prevTime') @@ -1261,7 +1265,7 @@ function get_curr_calday(offset) integer :: year, month, day, tod !----------------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet( tm_clock, currTime=date, rc=rc ) call chkrc(rc, sub//': error return from ESMF_ClockGet') @@ -1383,7 +1387,7 @@ integer function get_days_per_year( offset ) integer :: rc ! ESMF return code !--------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return if ( present(offset) )then call get_curr_date(yr, mon, day, tod, offset ) @@ -1416,7 +1420,7 @@ function get_curr_yearfrac( offset ) real(r8) :: cday ! current calendar day (1.0 = 0Z on Jan 1) real(r8) :: days_per_year ! days per year - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return cday = get_curr_calday(offset=offset) days_per_year = get_days_per_year() @@ -1439,7 +1443,7 @@ function get_prev_yearfrac() character(len=*), parameter :: sub = 'clm::get_curr_yearfrac' - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return get_prev_yearfrac = get_curr_yearfrac(offset = -dtime) @@ -1484,6 +1488,111 @@ end subroutine get_rest_date !========================================================================================= + integer function get_local_timestep_time( londeg, offset ) + + !--------------------------------------------------------------------------------- + ! Get the local time for this longitude that is evenly divisible by the time-step + ! + ! uses + use clm_varcon, only: degpsec, isecspday + ! Arguments + real(r8) , intent(in) :: londeg ! Longitude in degrees + integer, optional, intent(in) :: offset ! Offset from current time in seconds (either sign) + + ! Local variables + integer :: yr, mon, day ! year, month, day, unused + integer :: secs ! seconds into the day + real(r8) :: lon ! positive longitude + integer :: offset_sec ! offset seconds (either 0 for current time or -dtime for previous time) + !--------------------------------------------------------------------------------- + if ( present(offset) ) then + offset_sec = offset + else + offset_sec = 0 + end if + SHR_ASSERT( londeg >= -180.0_r8, "londeg must be greater than -180" ) + SHR_ASSERT( londeg <= 360.0_r8, "londeg must be less than 360" ) + call get_curr_date(yr, mon, day, secs, offset=offset_sec ) + lon = londeg + if ( lon < 0.0_r8 ) lon = lon + 360.0_r8 + get_local_timestep_time = secs + nint((lon/degpsec)/real(dtime,r8))*dtime + get_local_timestep_time = mod(get_local_timestep_time,isecspday) + end function get_local_timestep_time + + !========================================================================================= + + integer function get_local_time( londeg, starttime, offset ) + + !--------------------------------------------------------------------------------- + ! Get the local time for this longitude + ! + ! uses + use clm_varcon, only: degpsec, isecspday + ! Arguments + real(r8) , intent(in) :: londeg ! Longitude in degrees + integer, optional, intent(in) :: starttime ! Start time (sec) + integer, optional, intent(in) :: offset ! Offset from current time in seconds (either sign) + + ! Local variables + integer :: yr, mon, day ! year, month, day, unused + integer :: secs ! seconds into the day + integer :: start ! start seconds + integer :: offset_sec ! offset seconds (either 0 for current time or -dtime for previous time) + real(r8) :: lon ! positive longitude + !--------------------------------------------------------------------------------- + if ( present(starttime) ) then + start = starttime + else + start = 0 + end if + if ( present(offset) ) then + offset_sec = offset + else + offset_sec = 0 + end if + SHR_ASSERT( start >= 0, "starttime must be greater than or equal to zero" ) + SHR_ASSERT( start <= isecspday, "starttime must be less than or equal to number of seconds in a day" ) + SHR_ASSERT( londeg >= -180.0_r8, "londeg must be greater than -180" ) + SHR_ASSERT( londeg <= 360.0_r8, "londeg must be less than 360" ) + SHR_ASSERT( (offset_sec == 0) .or. (offset_sec == -dtime), "offset must be zero or negative time-step" ) + call get_curr_date(yr, mon, day, secs, offset=offset_sec ) + lon = londeg + if ( lon < 0.0_r8 ) lon = lon + 360.0_r8 + get_local_time = modulo(secs + nint(londeg/degpsec), isecspday) + get_local_time = modulo(get_local_time - start,isecspday) + end function get_local_time + + !========================================================================================= + + logical function is_near_local_noon( londeg, deltasec ) + + !--------------------------------------------------------------------------------- + ! Is this longitude near it's local noon? + ! + ! uses + use clm_varcon, only: degpsec, isecspday + ! Arguments + real(r8), intent(in) :: londeg ! Longitude in degrees + integer , intent(in) :: deltasec ! Number of seconds before or after local noon + + ! Local variables + integer :: local_secs ! Local time in seconds + integer, parameter :: noonsec = isecspday / 2 ! seconds at local noon + !--------------------------------------------------------------------------------- + SHR_ASSERT( deltasec < noonsec, "deltasec must be less than 12 hours" ) + local_secs = get_local_timestep_time( londeg ) + + if ( local_secs >= (noonsec - deltasec) .and. local_secs <= (noonsec + deltasec)) then + is_near_local_noon = .true. + else + is_near_local_noon = .false. + end if + + !--------------------------------------------------------------------------------- + end function is_near_local_noon + + !========================================================================================= + subroutine set_nextsw_cday( nextsw_cday_in ) ! Set the next radiation calendar day, so that radiation step can be calculated @@ -1515,7 +1624,7 @@ function is_beg_curr_day() character(len=*), parameter :: sub = 'clm::is_beg_curr_day' - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call get_curr_date(yr, mon, day, tod) is_beg_curr_day = ( tod == dtime ) @@ -1542,7 +1651,7 @@ function is_end_curr_day() character(len=*), parameter :: sub = 'clm::is_end_curr_day' !--------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call get_curr_date(yr, mon, day, tod) is_end_curr_day = (tod == 0) @@ -1566,7 +1675,7 @@ logical function is_end_curr_month() character(len=*), parameter :: sub = 'clm::is_end_curr_month' !--------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call get_curr_date(yr, mon, day, tod) is_end_curr_month = (day == 1 .and. tod == 0) @@ -1589,7 +1698,7 @@ logical function is_beg_curr_year() character(len=*), parameter :: subname = 'is_beg_curr_year' !----------------------------------------------------------------------- - call check_timemgr_initialized(subname) + if ( .not. check_timemgr_initialized(subname) ) return call get_curr_date(yr, mon, day, tod) is_beg_curr_year = (mon == 1 .and. day == 1 .and. tod == dtime) @@ -1612,7 +1721,7 @@ logical function is_end_curr_year() character(len=*), parameter :: subname = 'is_end_curr_year' !----------------------------------------------------------------------- - call check_timemgr_initialized(subname) + if ( .not. check_timemgr_initialized(subname) ) return call get_curr_date(yr, mon, day, tod) is_end_curr_year = (mon == 1 .and. day == 1 .and. tod == 0) @@ -1634,7 +1743,7 @@ logical function is_first_step() integer(ESMF_KIND_I8) :: step_no !--------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet( tm_clock, advanceCount=step_no, rc=rc ) call chkrc(rc, sub//': error return from ESMF_ClockGet') @@ -1649,7 +1758,7 @@ logical function is_first_restart_step() ! Return true on first step of restart or branch run only. character(len=*), parameter :: sub = 'clm::is_first_restart_step' - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return is_first_restart_step = tm_first_restart_step @@ -1663,7 +1772,7 @@ logical function is_first_step_of_this_run_segment() ! the first step of a startup, restart or branch run. character(len=*), parameter :: sub = 'clm::is_first_step_of_this_run_segment' - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return is_first_step_of_this_run_segment = (is_first_step() .or. is_first_restart_step()) @@ -1684,7 +1793,7 @@ logical function is_last_step() integer :: rc !--------------------------------------------------------------------------------- - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return call ESMF_ClockGet( tm_clock, stopTime=stop_date, & currTime=curr_date, TimeStep=time_step, rc=rc ) @@ -1704,7 +1813,7 @@ logical function is_perpetual() ! Return true on last timestep. character(len=*), parameter :: sub = 'clm::is_perpetual' - call check_timemgr_initialized(sub) + if ( .not. check_timemgr_initialized(sub) ) return is_perpetual = tm_perp_calendar @@ -1806,7 +1915,7 @@ end subroutine timemgr_spmdbcast !========================================================================================= - subroutine check_timemgr_initialized(caller) + logical function check_timemgr_initialized(caller) ! ! !DESCRIPTION: ! Checks if the time manager has been initialized. If not, aborts with an error @@ -1824,9 +1933,12 @@ subroutine check_timemgr_initialized(caller) if (.not. timemgr_set) then call shr_sys_abort(trim(caller)//":: Time manager has not been initialized") + check_timemgr_initialized = .false. + else + check_timemgr_initialized = .true. end if - end subroutine check_timemgr_initialized + end function check_timemgr_initialized !----------------------------------------------------------------------- subroutine timemgr_reset() diff --git a/src/utils/test/clm_time_manager_test/test_clm_time_manager.pf b/src/utils/test/clm_time_manager_test/test_clm_time_manager.pf index b27a16f912..8a13a9c17c 100644 --- a/src/utils/test/clm_time_manager_test/test_clm_time_manager.pf +++ b/src/utils/test/clm_time_manager_test/test_clm_time_manager.pf @@ -9,6 +9,7 @@ module test_clm_time_manager use unittestTimeManagerMod, only : & set_date => unittest_timemgr_set_curr_date, & set_nstep => unittest_timemgr_set_nstep + use unittestUtils , only : endrun_msg implicit none save @@ -220,4 +221,156 @@ contains end subroutine check_DA_nstep + @Test + subroutine check_local_time(this) + class(TestTimeManager), intent(inout) :: this + + integer :: secs + real(r8) :: londeg + integer :: expected + + ! Check for local noon at Greenich + londeg = 0.0_r8 + secs = 3600*12 + call set_date(yr=2, mon=3, day=1, tod=secs) + expected = secs + @assertEqual( expected, get_local_time( londeg ) ) + londeg = 360.0_r8 + @assertEqual( expected, get_local_time( londeg ) ) + + ! Check for local noon on otherside of the world + secs = 0 + londeg = 180.0_r8 + call set_date(yr=2, mon=3, day=1, tod=secs ) + @assertEqual( expected, get_local_time( londeg ) ) + londeg = -180.0_r8 + @assertEqual( expected, get_local_time( londeg ) ) + end subroutine check_local_time + + @Test + subroutine compare_local_time(this) + class(TestTimeManager), intent(inout) :: this + + integer :: secs + real(r8) :: londeg + + londeg = 0.0_r8 + do while ( londeg <= 360.0_r8 ) + londeg = londeg + 0.1_r8 + ! Start at 0 Z + secs = 0 + call set_date(yr=1, mon=1, day=1, tod=secs) + do while ( .not. is_end_curr_year() ) + @assertEqual( get_local_time( londeg, starttime=0 ), get_local_timestep_time( londeg ) ) + call advance_timestep() + end do + end do + + end subroutine compare_local_time + + @Test + subroutine check_local_time_woffset(this) + class(TestTimeManager), intent(inout) :: this + + integer :: secs + real(r8) :: londeg + integer :: expected + + ! Check for local noon at Greenich for 1 time step ahead + londeg = 0.0_r8 + secs = 3600*12 + dtime + call set_date(yr=2018, mon=9, day=3, tod=secs) + expected = secs - dtime + @assertEqual( expected, get_local_time( londeg, offset=-dtime ) ) + londeg = 360.0_r8 + @assertEqual( expected, get_local_time( londeg, offset=-dtime ) ) + end subroutine check_local_time_woffset + + + @Test + subroutine check_near_local_noon(this) + class(TestTimeManager), intent(inout) :: this + + integer :: secs + real(r8) :: londeg + + ! Check for local noon at Greenich will be true from 11 to 1pm + londeg = 0.0_r8 + secs = 3600*11 + call set_date(yr=2, mon=3, day=1, tod=secs) + @assertTrue( is_near_local_noon( londeg, deltasec=3600 ) ) + secs = 3600*12 + call set_date(yr=2, mon=2, day=2, tod=secs) + @assertTrue( is_near_local_noon( londeg, deltasec=3600 ) ) + secs = 3600*13 + call set_date(yr=2, mon=1, day=10, tod=secs) + @assertTrue( is_near_local_noon( londeg, deltasec=3600 ) ) + + ! anything before or after will be False + secs =0 + call set_date(yr=2, mon=2, day=19, tod=secs) + @assertFalse( is_near_local_noon( londeg, deltasec=3600 ) ) + secs = 3600*11 - 1 + call set_date(yr=2, mon=4, day=23, tod=secs) + @assertFalse( is_near_local_noon( londeg, deltasec=3600 ) ) + secs = 3600*13 + 1 + call set_date(yr=2, mon=6, day=30, tod=secs) + @assertFalse( is_near_local_noon( londeg, deltasec=3600 ) ) + secs = 3600*24 - 1 + call set_date(yr=2, mon=12, day=31, tod=secs) + @assertFalse( is_near_local_noon( londeg, deltasec=3600 ) ) + end subroutine check_near_local_noon + + @Test + subroutine bad_deltasec(this) + class(TestTimeManager), intent(inout) :: this + + character(len=256) :: expected_msg + real(r8) :: londeg + integer :: secs + logical :: check + + londeg = 0.0_r8 + secs = get_local_time( londeg ) + secs = 3600*24 - 1 + call set_date(yr=2, mon=12, day=31, tod=secs) + check = is_near_local_noon( londeg, deltasec=3600*12+1 ) + expected_msg = endrun_msg( & + "ERROR: deltasec must be less than 12 hours" ) + @assertExceptionRaised(expected_msg) + + end subroutine bad_deltasec + + @Test + subroutine bad_neglontolocal_time(this) + class(TestTimeManager), intent(inout) :: this + + character(len=256) :: expected_msg + real(r8) :: londeg + integer :: secs + + londeg = -200.0_r8 + secs = get_local_time( londeg ) + expected_msg = endrun_msg( & + "ERROR: londeg must be greater than -180" ) + @assertExceptionRaised(expected_msg) + + end subroutine bad_neglontolocal_time + + @Test + subroutine bad_hilontolocal_time(this) + class(TestTimeManager), intent(inout) :: this + + character(len=256) :: expected_msg + real(r8) :: londeg + integer :: secs + + londeg = 400.0_r8 + secs = get_local_time( londeg ) + expected_msg = endrun_msg( & + "ERROR: londeg must be less than 360" ) + @assertExceptionRaised(expected_msg) + + end subroutine bad_hilontolocal_time + end module test_clm_time_manager diff --git a/tools/mkmapdata/mkmapdata.sh b/tools/mkmapdata/mkmapdata.sh index 995bde9e63..2aad9f24fe 100755 --- a/tools/mkmapdata/mkmapdata.sh +++ b/tools/mkmapdata/mkmapdata.sh @@ -77,6 +77,8 @@ usage() { echo " Displays this help message" echo "[-v|--verbose]" echo " Toggle verbose usage -- log more information on what is happening " + echo "[--fast]" + echo " Toggle fast maps only -- only create the maps that can be done quickly " echo "" echo " You can also set the following env variables:" echo " ESMFBIN_PATH - Path to ESMF binaries " @@ -134,6 +136,7 @@ verbose="no" list="no" outgrid="" gridfile="default" +fast="no" while [ $# -gt 0 ]; do case $1 in @@ -146,6 +149,9 @@ while [ $# -gt 0 ]; do -d|--debug) debug="YES" ;; + --fast) + fast="YES" + ;; -l|--list) debug="YES" list="YES" @@ -498,6 +504,9 @@ until ((nfile>${#INGRID[*]})); do # Skip if file already exists if [ -f "${OUTFILE[nfile]}" ]; then echo "Skipping creation of ${OUTFILE[nfile]} as already exists" + # Skip if large file and Fast mode is on + elif [ "$fast" = "YES" ] && [ "${SRC_LRGFIL[nfile]}" = "netcdf4" ]; then + echo "Skipping creation of ${OUTFILE[nfile]} as fast mode is on so skipping large files in NetCDF4 format" else cmd="$mpirun $ESMF_REGRID --ignore_unmapped -s ${INGRID[nfile]} " diff --git a/tools/mkmapdata/regridbatch.sh b/tools/mkmapdata/regridbatch.sh index af31db30de..7f266d4b66 100755 --- a/tools/mkmapdata/regridbatch.sh +++ b/tools/mkmapdata/regridbatch.sh @@ -29,9 +29,30 @@ if [ -z "$RES" ]; then echo "Run for all valid resolutions" resols=`../../bld/queryDefaultNamelist.pl -res list -silent` + if [ ! -z "$GRIDFILE" ]; then + echo "When GRIDFILE set RES also needs to be set for a single resolution" + exit 1 + fi else resols="$RES" fi +if [ -z "$GRIDFILE" ]; then + grid="" +else + if [[ ${#resols[@]} > 1 ]]; then + echo "When GRIDFILE is specificed only one resolution can also be given (# resolutions ${#resols[@]})" + echo "Resolutions input is: $resols" + exit 1 + fi + grid="-f $GRIDFILE" +fi + +if [ -z "$MKMAPDATA_OPTIONS" ]; then + echo "Run with standard options" + options=" " +else + options="$MKMAPDATA_OPTIONS" +fi echo "Create mapping files for this list of resolutions: $resols" #---------------------------------------------------------------------- @@ -39,7 +60,7 @@ echo "Create mapping files for this list of resolutions: $resols" for res in $resols; do echo "Create mapping files for: $res" #---------------------------------------------------------------------- - cmdargs="-r $res" + cmdargs="-r $res $grid $options" # For single-point and regional resolutions, tell mkmapdata that # output type is regional @@ -48,6 +69,10 @@ for res in $resols; do else res_type="global" fi + # Assume if you are providing a gridfile that the grid is regional + if [ $grid != "" ];then + res_type="regional" + fi cmdargs="$cmdargs -t $res_type" diff --git a/tools/mksurfdata_map/mksurfdata.pl b/tools/mksurfdata_map/mksurfdata.pl index 6f7c784388..7a2955421f 100755 --- a/tools/mksurfdata_map/mksurfdata.pl +++ b/tools/mksurfdata_map/mksurfdata.pl @@ -56,6 +56,7 @@ exedir=>undef, allownofile=>undef, crop=>1, + fast_maps=>0, hirespft=>undef, years=>"1850,2000", glc_nec=>10, @@ -115,6 +116,7 @@ sub usage { -dynpft "filename" Dynamic PFT/harvesting file to use (rather than create it on the fly) (must be consistent with first year) + -fast_maps Toggle fast mode which doesn't use the large mapping files -glc_nec "number" Number of glacier elevation classes to use (by default $opts{'glc_nec'}) -merge_gis If you want to use the glacier dataset that merges in the Greenland Ice Sheet data that CISM uses (typically @@ -330,7 +332,6 @@ sub write_namelist_file { map_fpeat = '$map->{'peat'}' map_fsoildepth = '$map->{'soildepth'}' map_fabm = '$map->{'abm'}' - map_ftopostats = '$map->{'topostats'}' map_fvic = '$map->{'vic'}' map_fch4 = '$map->{'ch4'}' mksrf_fsoitex = '$datfil->{'tex'}' @@ -346,7 +347,6 @@ sub write_namelist_file { mksrf_fpeat = '$datfil->{'peat'}' mksrf_fsoildepth = '$datfil->{'soildepth'}' mksrf_fabm = '$datfil->{'abm'}' - mksrf_ftopostats = '$datfil->{'topostats'}' mksrf_fvic = '$datfil->{'vic'}' mksrf_fch4 = '$datfil->{'ch4'}' outnc_double = $double @@ -355,6 +355,16 @@ sub write_namelist_file { mksrf_furban = '$datfil->{'urb'}' gitdescribe = '$gitdescribe' EOF + if ( ! $opts{'fast_maps'} ) { + print $fh <<"EOF"; + map_ftopostats = '$map->{'topostats'}' + mksrf_ftopostats = '$datfil->{'topostats'}' +EOF + } else { + print $fh <<"EOF"; + std_elev = 371.0d00 +EOF + } if ( defined($opts{'soil_override'}) ) { print $fh <<"EOF"; soil_clay = $opts{'soil_cly'} @@ -423,6 +433,7 @@ sub write_namelist_file { "hirespft" => \$opts{'hirespft'}, "l|dinlc=s" => \$opts{'csmdata'}, "d|debug" => \$opts{'debug'}, + "fast_maps" => \$opts{'fast_maps'}, "dynpft=s" => \$opts{'dynpft'}, "y|years=s" => \$opts{'years'}, "exedir=s" => \$opts{'exedir'}, @@ -603,9 +614,13 @@ sub write_namelist_file { } my $mopts = "$queryopts -namelist default_settings $usrnam"; my $mkopts = "-csmdata $CSMDATA -silent -justvalue -namelist clmexp $usrnam"; - foreach my $typ ( "lak", "veg", "voc", "tex", "col", "hrv", + my @typlist = ( "lak", "veg", "voc", "tex", "col", "hrv", "fmx", "lai", "urb", "org", "glc", "glcregion", "utp", "wet", - "gdp", "peat","soildepth","abm", "topostats" , "vic", "ch4") { + "gdp", "peat","soildepth","abm", "vic", "ch4"); + if ( ! $opts{'fast_maps'} ) { + push( @typlist, "topostats" ); + } + foreach my $typ ( @typlist ) { my $lmask = `$scrdir/../../bld/queryDefaultNamelist.pl $mopts -options type=$typ,mergeGIS=$merge_gis,hirespft=$hirespft -var lmask`; $lmask = trim($lmask); my $hgrid_cmd = "$scrdir/../../bld/queryDefaultNamelist.pl $mopts -options type=$typ,hirespft=$hirespft -var hgrid"; diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90 index fbbf57171a..5089fabf7d 100644 --- a/tools/mksurfdata_map/src/mksurfdat.F90 +++ b/tools/mksurfdata_map/src/mksurfdat.F90 @@ -139,6 +139,8 @@ program mksurfdat real(r8), allocatable :: p3(:) ! coefficient for qflx_surf_lag for finundated (s/mm) real(r8), allocatable :: zwt0(:) ! decay factor for finundated (m) + real(r8) :: std_elev = -999.99_r8 ! Standard deviation of elevation (m) to use for entire grid + integer, allocatable :: harvind1D(:) ! Indices of 1D harvest fields integer, allocatable :: harvind2D(:) ! Indices of 2D harvest fields @@ -214,6 +216,7 @@ program mksurfdat fsurdat, & fdyndat, & fsurlog, & + std_elev, & urban_skip_abort_on_invalid_data_check !----------------------------------------------------------------------- @@ -655,7 +658,7 @@ program mksurfdat ! Compute topography statistics [topo_stddev, slope] from [ftopostats] call mktopostats (ldomain, mapfname=map_ftopostats, datfname=mksrf_ftopostats, & - ndiag=ndiag, topo_stddev_o=topo_stddev, slope_o=slope) + ndiag=ndiag, topo_stddev_o=topo_stddev, slope_o=slope, std_elev=std_elev) ! Make VIC parameters [binfl, ws, dsmax, ds] from [fvic] call mkVICparams (ldomain, mapfname=map_fvic, datfname=mksrf_fvic, ndiag=ndiag, & diff --git a/tools/mksurfdata_map/src/mktopostatsMod.F90 b/tools/mksurfdata_map/src/mktopostatsMod.F90 index 2348a9f6f3..cc1f541354 100644 --- a/tools/mksurfdata_map/src/mktopostatsMod.F90 +++ b/tools/mksurfdata_map/src/mktopostatsMod.F90 @@ -36,7 +36,7 @@ module mktopostatsMod ! !IROUTINE: mktopostats ! ! !INTERFACE: -subroutine mktopostats(ldomain, mapfname, datfname, ndiag, topo_stddev_o, slope_o) +subroutine mktopostats(ldomain, mapfname, datfname, ndiag, topo_stddev_o, slope_o, std_elev) ! ! !DESCRIPTION: ! make various topography statistics @@ -55,6 +55,7 @@ subroutine mktopostats(ldomain, mapfname, datfname, ndiag, topo_stddev_o, slope_ character(len=*) , intent(in) :: mapfname ! input mapping file name character(len=*) , intent(in) :: datfname ! input data file name integer , intent(in) :: ndiag ! unit number for diag out + real(r8) , intent(in) :: std_elev ! standard deviation of elevation (m) to use when not using input file real(r8) , intent(out):: topo_stddev_o(:) ! output grid: standard deviation of elevation (m) real(r8) , intent(out):: slope_o(:) ! output grid: slope (degrees) ! @@ -72,6 +73,7 @@ subroutine mktopostats(ldomain, mapfname, datfname, ndiag, topo_stddev_o, slope_ real(r8), allocatable :: data_i(:) ! data on input grid integer :: ncid,varid ! input netCDF id's integer :: ier ! error status + logical :: bypass_reading ! If should bypass reading dataset and just use a global value real(r8), parameter :: min_valid_topo_stddev = 0._r8 @@ -82,68 +84,87 @@ subroutine mktopostats(ldomain, mapfname, datfname, ndiag, topo_stddev_o, slope_ !----------------------------------------------------------------------- write (6,*) 'Attempting to make Topography statistics.....' + if ( std_elev >= 0.0_r8 )then + bypass_reading = .true. + write (6,*) ' By pass the reading and just use global values' + else + bypass_reading = .false. + end if call shr_sys_flush(6) ! ----------------------------------------------------------------- ! Read domain and mapping information, check for consistency ! ----------------------------------------------------------------- - call domain_read(tdomain,datfname) + if ( .not. bypass_reading )then + call domain_read(tdomain,datfname) - call gridmap_mapread(tgridmap, mapfname ) - call gridmap_check( tgridmap, subname ) + call gridmap_mapread(tgridmap, mapfname ) + call gridmap_check( tgridmap, subname ) - call domain_checksame( tdomain, ldomain, tgridmap ) + call domain_checksame( tdomain, ldomain, tgridmap ) - ! ----------------------------------------------------------------- - ! Open input file, allocate memory for input data - ! ----------------------------------------------------------------- + ! ----------------------------------------------------------------- + ! Open input file, allocate memory for input data + ! ----------------------------------------------------------------- - write(6,*)'Open Topography file: ', trim(datfname) - call check_ret(nf_open(datfname, 0, ncid), subname) + write(6,*)'Open Topography file: ', trim(datfname) + call check_ret(nf_open(datfname, 0, ncid), subname) - allocate(data_i(tdomain%ns), stat=ier) - if (ier/=0) call abort() + allocate(data_i(tdomain%ns), stat=ier) + if (ier/=0) call abort() - ! ----------------------------------------------------------------- - ! Make topography standard deviation - ! ----------------------------------------------------------------- + ! ----------------------------------------------------------------- + ! Make topography standard deviation + ! ----------------------------------------------------------------- - call check_ret(nf_inq_varid (ncid, 'ELEVATION', varid), subname) - call check_ret(nf_get_var_double (ncid, varid, data_i), subname) - call gridmap_areastddev(tgridmap, data_i, topo_stddev_o, nodata=0._r8) + call check_ret(nf_inq_varid (ncid, 'ELEVATION', varid), subname) + call check_ret(nf_get_var_double (ncid, varid, data_i), subname) + call gridmap_areastddev(tgridmap, data_i, topo_stddev_o, nodata=0._r8) + + call output_diagnostics_continuous_outonly(topo_stddev_o, tgridmap, "Topo Std Dev", "m", ndiag) + else + write (6,*) ' Set std deviation of topography to ', std_elev + topo_stddev_o = std_elev + end if ! Check validity of output data if (min_bad(topo_stddev_o, min_valid_topo_stddev, 'topo_stddev')) then stop end if - call output_diagnostics_continuous_outonly(topo_stddev_o, tgridmap, "Topo Std Dev", "m", ndiag) ! ----------------------------------------------------------------- ! Regrid slope ! ----------------------------------------------------------------- - call check_ret(nf_inq_varid (ncid, 'SLOPE', varid), subname) - call check_ret(nf_get_var_double (ncid, varid, data_i), subname) - call gridmap_areaave(tgridmap, data_i, slope_o, nodata=0._r8) + if ( .not. bypass_reading )then + call check_ret(nf_inq_varid (ncid, 'SLOPE', varid), subname) + call check_ret(nf_get_var_double (ncid, varid, data_i), subname) + call gridmap_areaave(tgridmap, data_i, slope_o, nodata=0._r8) + call output_diagnostics_continuous(data_i, slope_o, tgridmap, "Slope", "degrees", ndiag) + else + write (6,*) ' Set slope of topography to ', 0.0_r8 + slope_o = 0.0_r8 + end if ! Check validity of output data if (min_bad(slope_o, min_valid_slope, 'slope') .or. & max_bad(slope_o, max_valid_slope, 'slope')) then stop end if - call output_diagnostics_continuous(data_i, slope_o, tgridmap, "Slope", "degrees", ndiag) ! ----------------------------------------------------------------- ! Close files and deallocate dynamic memory ! ----------------------------------------------------------------- - call check_ret(nf_close(ncid), subname) - call domain_clean(tdomain) - call gridmap_clean(tgridmap) - deallocate (data_i) + if ( .not. bypass_reading )then + call check_ret(nf_close(ncid), subname) + call domain_clean(tdomain) + call gridmap_clean(tgridmap) + deallocate (data_i) + end if write (6,*) 'Successfully made Topography statistics' write (6,*) diff --git a/tools/ncl_scripts/getco2_historical.ncl b/tools/ncl_scripts/getco2_historical.ncl index 13630de4f9..7284dc88b7 100644 --- a/tools/ncl_scripts/getco2_historical.ncl +++ b/tools/ncl_scripts/getco2_historical.ncl @@ -1,8 +1,8 @@ ; ; Take the greenhouse gas file used by CAM for historical (and future) representations of ; greenhouse gases, and convert it to a format that can be used by streams. -; So include domain data for a single point that covers the globe, as well -; as CO2 data over that single point. In the process we also discard the other +; So include domain data for a single point (or latitude bands) that covers the globe, as well +; as CO2 data over those latitude bands. In the process we also discard the other ; greenhouse gases, as the datm can only pass CO2. ; ; Erik Kluzek @@ -18,7 +18,7 @@ begin ; csmdata = getenv("CSMDATA"); clmroot = getenv("CLM_ROOT"); - rcp = getenv("RCP"); ; Get representative concentration pathway from env variable + hgrid = getenv("HGRID"); ; Get horizontal grid to use from env variable querynml = "bld/queryDefaultNamelist.pl -silent -justvalue "; if ( .not. ismissing(csmdata) )then querynml = querynml+" -csmdata "+csmdata; @@ -28,15 +28,15 @@ begin else querynml = clmroot+"/"+querynml; end if - histrcp = -999.9 - if ( ismissing(rcp) )then - rcp = histrcp + if ( ismissing(hgrid) )then + hgrid = "lat-bands" end if ; ; Get input Greenhouse gas file and open it ; filetype = "mkghg_bndtvghg"; - ghgfile = systemfunc( querynml+" -namelist clmexp -var "+filetype+" -options rcp="+rcp ); + print( querynml+" -namelist clmexp -var "+filetype+" -options hgrid="+hgrid ); + ghgfile = systemfunc( querynml+" -namelist clmexp -var "+filetype+" -options hgrid="+hgrid ); print( "Use "+filetype+" file: "+ghgfile ); if ( systemfunc("test -f "+ghgfile+"; echo $?" ) .ne. 0 )then print( "Input "+filetype+" file does not exist or not found: "+ghgfile ); @@ -51,35 +51,41 @@ begin ldate = systemfunc( "date" ); sim_yr0 = ncg->date(0) / 10000; - nyrs = dimsizes( ncg->date ); - sim_yr2 = ncg->date(nyrs-1) / 10000; + ntime = dimsizes( ncg->date ); + sim_yr2 = ncg->date(ntime-1) / 10000; - sim_yr_rng = sim_yr0 + "-" + sim_yr2; + sim_yr_rng = "_simyr_"+sim_yr0 + "-" + sim_yr2; - if ( rcp .eq. histrcp )then - outco2filename = "fco2_datm_"+sim_yr_rng+"_c"+sdate+".nc"; - else - outco2filename = "fco2_datm_rcp"+rcp+"_"+sim_yr_rng+"_c"+sdate+".nc"; - end if + cmip_vers = "_CMIP6_"; + outco2filename = "fco2_datm_"+hgrid+sim_yr_rng+cmip_vers+"c"+sdate+".nc"; system( "/bin/rm -f "+outco2filename ); print( "output file: "+outco2filename ); nco = addfile( outco2filename, "c" ); ; ; Define dimensions ; - nlat = 1; + if ( hgrid .eq. "lat-bands" )then + nlat = dimsizes(ncg->lat); + else + if ( hgrid .eq. "global" )then + nlat = 1 + else + print( "hgrid type can only be global or lat-bands: "+hgrid ) + exit + end if + end if nlon = 1; nv = 4; - dimnames = (/ "time", "lat", "lon", "nv" /); - dsizes = (/ nyrs, nlat, nlon, nv /); - is_unlim = (/ True, False, False, False /); + dimnames = (/ "time", "lat", "lon", "nv", "bounds" /); + dsizes = (/ ntime, nlat, nlon, nv, 2 /); + is_unlim = (/ True, False, False, False, False /); filedimdef( nco, dimnames, dsizes, is_unlim ); ; ; Define variables ; vars = (/ "lonc", "latc", "lonv", "latv", "mask", "frac", "area", "CO2" /); units= (/ "degrees_east", "degrees_north", "degree_east", "degrees_north", "unitless", "unitless", "radians^2", "ppmv" /); - lname= (/ "Longitude of grid cell center", "Latitude of grid cell center", "Longitudesof grid cell vertices", "Latitudes of grid cell vertices", "Mask of active cells: 1=active", "Fraction of grid cell that is active", "Area of grid cell", "CO2 concentration" /); + lname= (/ "Longitude of grid cell center", "Latitude of grid cell center", "Longitudes of grid cell vertices", "Latitudes of grid cell vertices", "Mask of active cells: 1=active", "Fraction of grid cell that is active", "Area of grid cell", "CO2 concentration" /); print( "Define variables: "+vars ); do i= 0, dimsizes(vars)-1 if ( vars(i) .eq. "lonv" .or. vars(i) .eq. "latv" )then @@ -95,6 +101,9 @@ begin nco->$vars(i)$@units = units(i); nco->$vars(i)$@lname = lname(i); end do + filevardef ( nco, "time", "float", (/ "time" /) ); + filevardef ( nco, "time_bnds", "float", (/ "time", "bounds" /) ); + filevardef ( nco, "date", "integer", (/ "time" /) ); varstatic = (/ "mask", "frac", "area" /); do i = 0, dimsizes(varstatic)-1 nco->$varstatic(i)$@coordinate = "latc lonc"; @@ -108,30 +117,77 @@ begin nco@history = ldate+": Convert by getco2_historical.ncl"; nco@source = "Convert from:"+ghgfile; nco@Version = systemfunc( "git describe" ); + filevarattdef( nco, "time", ncg->time ); + filevarattdef( nco, "date", ncg->date ); + nco->time_bnds@long_name = nco->time@long_name; + nco->time_bnds@units = nco->time@units; + nco->time_bnds@calendar = nco->time@calendar; ; ; Set static variables ; - pi = 3.14159265358979323846d00; - nco->area = 4.0*pi; - nco->mask = 1; - nco->frac = 1.0; - nco->latv(0,0,0:1) = 90.0; - nco->latc = 0.0; - nco->latv(0,0,2:3) = -90.0; - nco->lonv(0,0,0:3:3) = 0.0; - nco->lonc = 180.0; - nco->lonv(0,0,1:2) = 360.0; + pi = 3.14159265358979323846d00; + nco->mask = 1; + nco->frac = 1.0; + if ( nlat .gt. 1 )then + nco->latc = (/ ncg->lat/); + else + nco->latc = (/ 0.0d00 /); + end if + nco->latv(nlat-1,0,0) = 90.0d00; + nco->latv(nlat-1,0,3) = 90.0d00; + if ( nlat .gt. 1 )then + nco->latv(0:nlat-2,0,0) = ( (/ ncg->lat(0:nlat-2) /) + (/ncg->lat(1:nlat-1) /) )*0.5d00 + nco->latv(0:nlat-2,0,3) = (/ nco->latv(0:nlat-2,0,0) /); + nco->latv(1:nlat-1,0,1) = (/ nco->latv(0:nlat-2,0,0) /); + nco->latv(1:nlat-1,0,2) = (/ nco->latv(1:nlat-1,0,1) /); + end if + nco->latv(0,0,1) = -90.0d00; + nco->latv(0,0,2) = -90.0d00; + nco->lonv(:,0,0) = 0.0d00; + nco->lonv(:,0,3) = 0.0d00; + nco->lonc = 180.0d00; + nco->lonv(:,0,1) = 360.0d00; + nco->lonv(:,0,2) = 360.0d00; + clkws = gc_clkwise( nco->latv, nco->lonv ); + if ( any(clkws .eq. False) )then + print( "Some varticies are NOT clockwise" ); + exit + end if + ; EBK -- NOTE The NCL function wasn't giving me the correct answer so I used the mathmatical expression + ;nco->area = dble2flt( gc_qarea( nco->latv, nco->lonv ) ); + conv2rad = pi/180.0d00 + nco->area(:,0) = 2.0d00*pi*abs( sin((/nco->latv(:,0,0)/)*conv2rad) - sin((/nco->latv(:,0,1)/)*conv2rad) ); + if ( abs(sum(nco->area) - 4.0d00*pi) .gt. 1.d-14 )then + print( "Area of globe does not sum to 4*pi as expected" ); + exit + end if ; ; Time and date ; - nco->time = ncg->time; - nco->date = ncg->date; + nco->date = (/ ncg->date /); + nco->time = (/ ncg->time /); nco->date@comment = "This variable is NOT used when read by datm, the time coordinate is used"; ; ; CO2 ; - print( "Copy CO2 for "+nyrs+" years of data" ); - nco->CO2(:,0,0) = (/ ncg->CO2(:) /); + print( "Copy CO2 for "+ntime+" time samples of data" ); + if ( nlat .gt. 1 )then + do y = 0, nlat-1 + print( "latitude: "+ nco->latc(y,0) ); + nco->CO2(:,y,0) = (/ ncg->CO2_LBC(:,y) /) * 1.e6; + end do + else + ; make sure all latitudes on file are the same for each time + do itime = 0, ntime-1 + if ( max(ncg->CO2_LBC(itime,:)) .ne. min(ncg->CO2_LBC(itime,:)) )then + print( "Global average, but latitudes are NOT constant" ); + exit + end if + end do + nco->CO2(:,0,0) = (/ ncg->CO2_LBC(:,0) /) * 1.e6; + end if + print( "Average Global First CO2 ppmv value: Date="+nco->date(0)+" CO2="+avg(nco->CO2(0,:,0) ) ); + print( "Average Global Last CO2 ppmv value: Date="+nco->date(ntime-1)+" CO2="+avg(nco->CO2(ntime-1,:,0)) ); print( "================================================================================================" ); print( "Successfully created output historical CO2 file: "+outco2filename);