From de738e3a57589be5b005a6d0696b3259fd6e02f8 Mon Sep 17 00:00:00 2001 From: dwongepa Date: Sat, 13 Mar 2021 16:48:24 -0500 Subject: [PATCH 01/15] TYPE: new feature KEYWORDS: CMAQ, coupled model, constructing, WRF-CMAQ, executable SOURCE: David Wong, US EPA (wong.david-c@epa.gov) DESCRIPTION OF CHANGES: These modifications are for facilitating constructing a WRF-CMAQ coupled model without compromising creating a regular WRF executable. Problem: A WRF-CMAQ coupled model was developed in US EPA in around 2008. The work has been published in a journal article [1]. Since then, the coupled model has been applied to various studies and region to validate the effectiveness of the model [2-12]. Throughout the WRF released cycle, the coupled has been updated manually into a subset of WRF version. The update process was tedious and error prone. With this PR, the manual update process can be eliminated. LIST OF MODIFIED FILES: list of changed files (use git diff --name-status master to get formatted list) newly added file: Registry/registry.WRF-CMAQ-twoway phys/complex_number_module.F phys/module_twoway_ra_rrtmg_sw.F phys/module_twoway_rrtmg_aero_optical_util.F modified file: Makefile clean configure Registry/Registry.EM Registry/registry.em_shared_collection arch/Config.pl dyn_em/module_first_rk_step_part1.F dyn_em/solve_em.F external/io_netcdf/makefile main/Makefile main/depend.common phys/Makefile phys/module_ra_rrtmg_sw.F phys/module_radiation_driver.F AGAIN cherry pick to try to fix old WRF with CMAQ modified: Makefile modified: Registry/Registry.EM new file: Registry/registry.WRF-CMAQ-twoway modified: Registry/registry.em_shared_collection modified: arch/Config.pl modified: clean modified: configure modified: dyn_em/module_first_rk_step_part1.F modified: dyn_em/solve_em.F modified: external/io_netcdf/makefile modified: main/Makefile modified: main/depend.common modified: phys/Makefile new file: phys/complex_number_module.F modified: phys/module_ra_rrtmg_sw.F modified: phys/module_radiation_driver.F new file: phys/module_twoway_ra_rrtmg_sw.F new file: phys/module_twoway_rrtmg_aero_optical_util.F TESTS CONDUCTED: Three different models were built with given lines of instructions: * model 1 (used the original WRF source code without any changes in this PR) + went through a regular configure process to build a regular WRF model with the configure command + compiled the source code with usual command: compile em_real to generate a WRF executable (for discussion purposes, the WRF executable is named p.exe) * model 2 (used the new WRF source code with all the changes in this PR) + went through a regular configure process to build a regular WRF model with the configure command + compiled the source code with usual command: compile em_real to generate a WRF executable. (for discussion purposes, the WRF executable is named np.exe) * model 3 (used the new WRF source code with all the changes in this PR) + issue the following two environment setting commands: setenv WRF_CMAQ 1 setenv IOAPI /home/wdx/lib/x86_64/ifc-18.0/ioapi-3.2-2020220 (Note: ioapi_3.2-2020220 is another third-party library) + went through a regular configure process to build a regular WRF model with the configure command + compiled the source code with usual command: compile em_real to generate a WRF executable. (for discussion purposes, the WRF executable is named nc.exe) Case p was based on the executable of model 1, case np was based on the executable of model 2, and case ncX, where X represents various available options in running the WRF-CMAQ coupled model and detailed information is provided below, was based on the executable of model 3. All cases were tested with two days of simulation, 2016/7/1 - 2016/7/2. To verify the correctness of this new PR, a bit-by-bit comparison of every variable in all time steps on the second day of the wrfout file between two selected cases. One of the WRF-CMAQ coupled model running setting, wrf_cmaq_option: 0 = run WRF only 1 = run WRF only and w producing MCIP like GRID and MET files 2 = run WRF-CMAQ coupled model w/o producing MCIP like GRID and MET files 3 = run WRF-CMAQ coupled model w producing MCIP like GRID and MET files In these scenarios, X was set to the wrf_cmaq_option and another coupled model setting, direct_sw_feedback was set to .false., where direct_sw_feedback is a switch to turn on short wave aerosol direct effect or not. For an additional case, X = s means wrf_cmaq_option = 3, and direct_sw_feedback was set to .true. * comparing case p and case np Identical result was rendered. This means that this new PR does not alter anything in terms of constructing a regular WRF executable. * comparing case np and nc0 Identical result was rendered. This means that this new PR can make the WRF-CMAQ coupled behave as a regular WRF model * comparing case nc0, nc1, nc2 and nc3 Identical result was rendered. This means that this new PR with or without the chemistry portion of the coupled model does not change the WRF result * comparing case nc0 and ncs Different result was rendered with respect to any variables can be affect by radiation calculation. This means the aerosol information from CMAQ was transferred to WRF and radiation calculation through the short-wave aerosol direct effect. Here are some of the research work based on the WRF-CMAQ coupled model: 12. Jia Xing, Jiandong Wang, Rohit Mathur, Shuxiao Wang, Golam Sarwar, Jonathan Pleim, Christian Hogrefe, Yuqiang Zhang, Jingkun Jiang, David C. Wong, Jimin Hao, "Impacts of aerosol direct effects on tropospheric ozone through changes in atmospheric dynamics and photolysis rates", Atmos. Chem. Phys., 17, 9869-9883, 2017. 11. Wang, J., Xing, J., Mathur, R., Pleim, J., Wang, S. Hogrefe, C., Gan, C-M, Wong. D. C., and Hao, J., "Historical trends in PM2.5 related premature mortality during 1990-2010 across the northern hemisphere", Environmental Health Perspectives, Vol 125, Number 3, p400-408, March 2017. 10. Xing, J., Wang, J. Mathur, R., Pleim, J., Wang, S. Hogrefe, C., Gan, C-M, Wong. D. C., and Hao, J., "Unexpected benefits of reducing aerosol cooling effects", Environmental Science and Technology, 2016 Jul 19, 50(14):7527-34. 9. Gan, C-M., Hogrefe, C., Pleim, J., Mathur, R., Xing, J., Wong, D. C., Gilliam, R., Pouliot, G., and Wei, C., "Assessment of the Effects of Horizontal Grid Resolution on Long-Term Air Quality Trends using Coupled WRF-CMAQ Simulations", Atmospheric Environment, Volume 132, May 2016, Pages 207-216. 8. Xing, J., Mathur, R., Pleim, J., Hogrefe, C., Gan, C.-M., Wong, D. C., Wei C., and Wang J., "Air pollution and climate response to aerosol direct radiative effects: a modeling study of decadal trends across the northern hemisphere", Journal of Geophysical Research: Atmospheres, Vol. 120, Issue 23, Pages 12221 - 12236, Dec 2015. 7. Xing, J., Mathur, R., Pleim, J., Hogrefe, C., Gan, C.-M., Wong, D. C., and Wei C., "Can a coupled meteorology-chemistry model reproduce the historical trend in aerosol direct radiative effects over the Northern Hemisphere?", Atmos. Chem. Phys., 15, 9997-10018, 2015. 6. Gan, C-M., Pleim, J., Mathur, R., Hogrefe, C., Long, C., Xing, J., Wong, D. C., Gilliam, R., and Wei, C., "Assessment of long-term WRF-CMAQ simulations for understanding direct aerosol effects on radiation "brightening" in the United States", Atmos. Chem. Phys., 15, 12193-12209, 2015 5. Gan, C-M., Binkowski, F., Pleim, J., Xing, J., Wong, D. C., Mathur, R., and Gilliam, R., "Assessment of the aerosol optics component of the coupled WRF-CMAQ model using CARES field campaign data and a single column model", Atmospheric Environment, Volume 115, Pages 670-682, August 2015. 4. Hogrefe, C., Pouliot, G., Wong, D. C., Torian, A., Roselle, S., Pleim, J., and Mathur, R., "Annual application and evaluation of the online coupled WRF-CMAQ system over North America under AQMEII phase 2", Atmospheric Environment, Volume 115, Pages 1-756, August 2015. 3. Xing, J., Mathur, R., Pleim, J., Hogrefe, C., Gan, C.-M., Wong, D. C., Wei, C., Gilliam, R., and Pouliot, G., "Observations and modeling of air quality trends over 1990-2010 across the Northern Hemisphere: China, the United States and Europe", Atmos. Chem. Phys., 15, 2723-2747, 2015. 2. Wang, J., Wang, S., Jiang, J., Ding, A., Zheng, M., Zhao, B., Wong, D. C., Zhou, W., Zheng, G., Wang, L., Pleim, J. E., and Hao, J., "Impact of aerosol-meteorology interactions on fine particle pollution during China's severe haze episode in January 2013", Environ. Res. Lett. 9, September 2014. 1. Wong, D. C., Pleim, J., Mathur, R., Binkowski, F., Otte, T., Gilliam, R., Pouliot, G., Xiu, A., and Kang, D., "WRF-CMAQ two-way coupled system with aerosol feedback: software development and preliminary results", Geosci. Model Dev., 5, 299-312, 2012. --- Makefile | 4 + Registry/Registry.EM | 8 + Registry/registry.WRF-CMAQ-twoway | 51 + Registry/registry.em_shared_collection | 1 + arch/Config.pl | 155 +- clean | 8 +- configure | 7 + dyn_em/module_first_rk_step_part1.F | 55 +- dyn_em/solve_em.F | 144 +- external/io_netcdf/makefile | 2 +- main/Makefile | 2 +- main/depend.common | 3 + phys/Makefile | 2 + phys/complex_number_module.F | 253 ++ phys/module_ra_rrtmg_sw.F | 219 +- phys/module_radiation_driver.F | 201 +- phys/module_twoway_ra_rrtmg_sw.F | 327 +++ phys/module_twoway_rrtmg_aero_optical_util.F | 2768 ++++++++++++++++++ 18 files changed, 4199 insertions(+), 11 deletions(-) create mode 100644 Registry/registry.WRF-CMAQ-twoway create mode 100644 phys/complex_number_module.F create mode 100644 phys/module_twoway_ra_rrtmg_sw.F create mode 100644 phys/module_twoway_rrtmg_aero_optical_util.F diff --git a/Makefile b/Makefile index 1a9d628130..55634e35ef 100644 --- a/Makefile +++ b/Makefile @@ -988,6 +988,10 @@ physics : @ echo '--------------------------------------' if [ $(WRF_CHEM) -eq 0 ] ; then \ ( cd phys ; $(MAKE) CF2=" " ) ; \ + if [ $(WRF_CMAQ) -eq 1 ] ; then \ + @ echo '----------- make cmaq ----------------' ; \ + ( rm -f main/libcmaqlib.a; cd cmaq ; $(MAKE) -f Makefile.twoway ) ; \ + fi \ else \ ( cd phys ; $(MAKE) CF2="$(CHEM_FILES2)" ) ; \ fi diff --git a/Registry/Registry.EM b/Registry/Registry.EM index abfc5030bd..ee7a57d806 100644 --- a/Registry/Registry.EM +++ b/Registry/Registry.EM @@ -4,6 +4,14 @@ include registry.dimspec include registry.em_shared_collection +# WRF-CMAQ coupled model +rconfig integer wrf_cmaq_option namelist,wrf_cmaq 1 0 +rconfig integer wrf_cmaq_freq namelist,wrf_cmaq 1 1 +rconfig integer met_file_tstep namelist,wrf_cmaq 1 10000 + +rconfig logical direct_sw_feedback namelist,wrf_cmaq 1 .false. +rconfig logical feedback_restart namelist,wrf_cmaq 1 .false. + # added to output 5 for ESMF state real landmask ij misc 1 - i0125rh056d=(interp_fcnm_imask)u=(copy_fcnm) "LANDMASK" "LAND MASK (1 FOR LAND, 0 FOR WATER)" "" state real lakemask ij misc 1 - i012rhd=(interp_fcnm_imask)u=(copy_fcnm) "LAKEMASK" "LAKE MASK (1 FOR LAKE, 0 FOR NON-LAKE)" "" diff --git a/Registry/registry.WRF-CMAQ-twoway b/Registry/registry.WRF-CMAQ-twoway new file mode 100644 index 0000000000..fcbd101e72 --- /dev/null +++ b/Registry/registry.WRF-CMAQ-twoway @@ -0,0 +1,51 @@ +#state real mass_ws_i ikj twoway_feedback_data 1 - r "mass_ws_i" "Water soluble i mode" "ug/m**3" +#state real mass_ws_j ikj twoway_feedback_data 1 - r "mass_ws_j" "Water soluble j mode" "ug/m**3" +#state real mass_ws_k ikj twoway_feedback_data 1 - r "mass_ws_k" "Water soluble k mode" "ug/m**3" +#state real mass_in_i ikj twoway_feedback_data 1 - r "mass_in_i" "Water insoluble i mode" "ug/m**3" +#state real mass_in_j ikj twoway_feedback_data 1 - r "mass_in_j" "Water insoluble j mode" "ug/m**3" +#state real mass_in_k ikj twoway_feedback_data 1 - r "mass_in_k" "Water insoluble k mode" "ug/m**3" +#state real mass_ec_i ikj twoway_feedback_data 1 - r "mass_ec_i" "Elemental carbon i mode" "ug/m**3" +#state real mass_ec_j ikj twoway_feedback_data 1 - r "mass_ec_j" "Elemental carbon j mode" "ug/m**3" +#state real mass_ec_k ikj twoway_feedback_data 1 - r "mass_ec_k" "Elemental carbon k mode" "ug/m**3" +#state real mass_ss_i ikj twoway_feedback_data 1 - r "mass_ss_i" "Seasalt i mode" "ug/m**3" +#state real mass_ss_j ikj twoway_feedback_data 1 - r "mass_ss_j" "Seasalt j mode" "ug/m**3" +#state real mass_ss_k ikj twoway_feedback_data 1 - r "mass_ss_k" "Seasalt k mode" "ug/m**3" +#state real mass_h2o_i ikj twoway_feedback_data 1 - r "mass_h2o_i" "water i mode" "ug/m**3" +#state real mass_h2o_j ikj twoway_feedback_data 1 - r "mass_h2o_j" "water j mode" "ug/m**3" +#state real mass_h2o_k ikj twoway_feedback_data 1 - r "mass_h2o_k" "water k mode" "ug/m**3" +#state real dgn_i ikj twoway_feedback_data 1 - r "dgn_i" "diameter i mode" "m" +#state real dgn_j ikj twoway_feedback_data 1 - r "dgn_j" "diameter j mode" "m" +#state real dgn_k ikj twoway_feedback_data 1 - r "dgn_k" "diameter k mode" "m" +#state real sig_i ikj twoway_feedback_data 1 - r "sig_i" "standard deviations i mode" "" +#state real sig_j ikj twoway_feedback_data 1 - r "sig_j" "standard deviations j mode" "" +#state real sig_k ikj twoway_feedback_data 1 - r "sig_k" "standard deviations k mode" "" + +#state real prev_rainnc ij twoway_feedback_data 1 - r "prev_rainnc" "previous accumlated rainnc" "mm" +#state real prev_rainc ij twoway_feedback_data 1 - r "prev_rainc" "previous accumlated rainc" "mm" + +#state real ozone ikj twoway_feedback_data 1 - r "OZONE" "ozone value" "" + +#state real sw_gtauxar_01 ikj misc 1 - hr "SW_GTAUXAR_01" "SW Optical depth of lamda = 0.388 um" "" +#state real sw_gtauxar_02 ikj misc 1 - hr "SW_GTAUXAR_02" "SW Optical depth of lamda = 0.533 um" "" +#state real sw_gtauxar_03 ikj misc 1 - hr "SW_GTAUXAR_03" "SW Optical depth of lamda = 0.702 um" "" +#state real sw_gtauxar_04 ikj misc 1 - hr "SW_GTAUXAR_04" "SW Optical depth of lamda = 1.010 um" "" +#state real sw_gtauxar_05 ikj misc 1 - hr "SW_GTAUXAR_05" "SW Optical depth of lamda = 1.271 um" "" +#state real sw_ttauxar_01 ij misc 1 - hr "SW_TTAUXAR_01" "Optical depth sum of SW_GTAUXAR_01" "" +#state real sw_ttauxar_02 ij misc 1 - hr "SW_TTAUXAR_02" "Optical depth sum of SW_GTAUXAR_02" "" +#state real sw_ttauxar_03 ij misc 1 - hr "SW_TTAUXAR_03" "Optical depth sum of SW_GTAUXAR_03" "" +#state real sw_ttauxar_04 ij misc 1 - hr "SW_TTAUXAR_04" "Optical depth sum of SW_GTAUXAR_04" "" +#state real sw_ttauxar_05 ij misc 1 - hr "SW_TTAUXAR_05" "Optical depth sum of SW_GTAUXAR_05" "" +#state real sw_asy_fac_01 ikj misc 1 - hr "SW_ASY_FAC_01" "Corresponding Asymmetry Factor of SW_GTAUXAR_01" "" +#state real sw_asy_fac_02 ikj misc 1 - hr "SW_ASY_FAC_02" "Corresponding Asymmetry Factor of SW_GTAUXAR_02" "" +#state real sw_asy_fac_03 ikj misc 1 - hr "SW_ASY_FAC_03" "Corresponding Asymmetry Factor of SW_GTAUXAR_03" "" +#state real sw_asy_fac_04 ikj misc 1 - hr "SW_ASY_FAC_04" "Corresponding Asymmetry Factor of SW_GTAUXAR_04" "" +#state real sw_asy_fac_05 ikj misc 1 - hr "SW_ASY_FAC_05" "Corresponding Asymmetry Factor of SW_GTAUXAR_05" "" +#state real sw_ssa_01 ikj misc 1 - hr "SW_SSA_01" "Corresponding Single scattering Albedo of SW_GTAUXAR_01" "" +#state real sw_ssa_02 ikj misc 1 - hr "SW_SSA_02" "Corresponding Single scattering Albedo of SW_GTAUXAR_02" "" +#state real sw_ssa_03 ikj misc 1 - hr "SW_SSA_03" "Corresponding Single scattering Albedo of SW_GTAUXAR_03" "" +#state real sw_ssa_04 ikj misc 1 - hr "SW_SSA_04" "Corresponding Single scattering Albedo of SW_GTAUXAR_04" "" +#state real sw_ssa_05 ikj misc 1 - hr "SW_SSA_05" "Corresponding Single scattering Albedo of SW_GTAUXAR_05" "" + +#state real sw_zbbcddir ij misc 1 - hr "SW_ZBBCDDIR" "Clear sky downward direct shortwave flux" "W m-2" +#state real sw_dirdflux ij misc 1 - hr "SW_DIRDFLUX" "Direct downward shortwave surface flux" "W m-2" +#state real sw_difdflux ij misc 1 - hr "SW_DIFDFLUX" "Diffuse downward shortwave surface flux" "W m-2" diff --git a/Registry/registry.em_shared_collection b/Registry/registry.em_shared_collection index 3a633e8b26..c58105e805 100644 --- a/Registry/registry.em_shared_collection +++ b/Registry/registry.em_shared_collection @@ -28,3 +28,4 @@ include registry.new3d_wif include registry.trad_fields include registry.solar_fields include registry.diags +include registry.WRF-CMAQ-twoway diff --git a/arch/Config.pl b/arch/Config.pl index ecccb670cd..8bb169ebe8 100644 --- a/arch/Config.pl +++ b/arch/Config.pl @@ -5,6 +5,9 @@ # Be sure to run as ./configure (to avoid getting a system configure command by mistake) # +use Cwd qw(getcwd); +$wrf_cmaq_option = $ENV{'WRF_CMAQ'}; # determine building WRF-CMAQ coupled model or not + select((select(STDOUT), $|=1)[0]); $sw_perl_path = perl ; $sw_netcdf_path = "" ; @@ -444,11 +447,157 @@ $latchon = 0 ; while ( ) { + if ( $_ =~ /ifort compiler/ ) + { $lioapi = 'Linux2_x86_64ifort'; + } + elsif ( $_ =~ /PGI compiler/ ) + { $lioapi = 'Linux2_x86_64pg'; + } + elsif ( $_ =~ /gfortran compiler/ ) + { $lioapi = 'Linux2_x86_64gfort'; + } + if ( substr( $_, 0, 5 ) eq "#ARCH" && $latchon == 1 ) { close CONFIGURE_DEFAULTS ; if ( $sw_opt_level eq "-f" ) { - open CONFIGURE_DEFAULTS, "cat ./arch/postamble ./arch/noopt_exceptions_f |" or die "horribly" ; + + # determine whether variable EM_MODULE_DIR contains -I../cmaq or not + my $file = "Makefile"; + open(FH, $file) or die("File $file not found"); + + $lib_path_wo_cmaq = 1; + while ( my $String = ) + { if($String =~ /-I..\/dyn_em -I..\/cmaq/) + { $lib_path_wo_cmaq = 0; + } + } + close (FH); + + # determine whether declarations in Registry/registry.WRF-CMAQ-twoway is commented out or not + my $file = "Registry/registry.WRF-CMAQ-twoway"; + open (FH, $file) or die("File $file not found"); + + $registry_wo_cmaq = 0; + while ( my $String = ) + { if($String =~ /#state/) + { $registry_wo_cmaq = 1; + } + } + close (FH); + + # determine whether express #NOWIN LIB_BUNDLED contains $(IOAPI_LIB) or not + my $file = "arch/preamble"; + open (FH, $file) or die("File $file not found"); + + $bundle_wo_ioapi = 1; + while ( my $String = ) + { if ( $String =~ /IOAPI_LIB/ ) + { $bundle_wo_ioapi = 0; + } + } + close (FH); + + if ( $wrf_cmaq_option eq 1 ) # build WRF-CMAQ coupled model + { system("cp Registry/registry.WRF-CMAQ-twoway_full Registry/registry.WRF-CMAQ-twoway"); + + if ( $lib_path_wo_cmaq == 1 ) + { open (FILE, "; + close (FILE); + + foreach ( @lines ) + { $_ =~ s/ -I..\/dyn_em/ -I..\/dyn_em -I..\/cmaq/g; + } + + open (FILE, ">Makefile") || die "File not found"; + print FILE @lines; + close (FILE); + } + + if ( $registry_wo_cmaq == 1 ) + { open (FILE, "; + close (FILE); + + foreach (@lines) + { $_ =~ s/#state/state/g; + } + + open (FILE, ">Registry/registry.WRF-CMAQ-twoway") || die "File not found"; + print FILE @lines; + close (FILE); + } + + if ( $bundle_wo_ioapi == 1 ) + { open (FILE, "; + close (FILE); + + foreach (@lines) + { $_ =~ s/#NOWIN LIB_BUNDLED = \\/#NOWIN LIB_BUNDLED = \$(IOAPI_LIB) \\/g; + } + + open (FILE, ">arch/preamble") || die "File not found"; + print FILE @lines; + close (FILE); + } + + open (FH, '>', wrf_cmaq_path) or die $! ; + $ioapi_path = $ENV{'IOAPI'} ; + print FH "CMAQLIB = libcmaqlib.a \n" ; + print FH "IOAPI = $ioapi_path\n" ; + print FH "LIOAPI = $lioapi\n" ; + print FH "IOAPI_LIB = -L$ioapi_path/$lioapi -lioapi \n" ; + close (FH) ; + open CONFIGURE_DEFAULTS, "cat wrf_cmaq_path ./arch/postamble ./arch/noopt_exceptions_f |" or die "horribly" ; + } + else + { if ( $lib_path_wo_cmaq == 0 ) + { open (FILE, "; + close (FILE); + + foreach (@lines) + { $_ =~ s/ -I..\/cmaq//g; + } + + open (FILE, ">Makefile") || die "File not found"; + print FILE @lines; + close (FILE); + + } + + if ( $registry_wo_cmaq == 0 ) + { open (FILE, "; + close (FILE); + + foreach(@lines) + { $_ =~ s/state/#state/g; + } + + open (FILE, ">Registry/registry.WRF-CMAQ-twoway") || die "File not found"; + print FILE @lines; + close (FILE); + } + + if ( $bundle_wo_ioapi == 0 ) + { open (FILE, "; + close (FILE); + + foreach (@lines) + { $_ =~ s/ \$\(IOAPI_LIB\)//g; + } + + open (FILE, ">arch/preamble") || die "File not found"; + print FILE @lines; + close (FILE); + } + + open CONFIGURE_DEFAULTS, "cat ./arch/postamble ./arch/noopt_exceptions_f |" or die "horribly" ; + } } else { open CONFIGURE_DEFAULTS, "cat ./arch/postamble ./arch/noopt_exceptions |" or die "horribly" ; } @@ -777,6 +926,10 @@ close POSTAMBLE ; close ARCH_NOOPT_EXCEPTIONS ; +if ( $wrf_cmaq_option eq 1 ) + { unlink "wrf_cmaq_path"; + } + open CONFIGURE_WRF, "> configure.wrf" or die "cannot append configure.wrf" ; open ARCH_PREAMBLE, "< arch/preamble" or die "cannot open arch/preamble" ; my @preamble; diff --git a/clean b/clean index 007447603e..f0f55824f2 100755 --- a/clean +++ b/clean @@ -3,11 +3,15 @@ set nonomatch -foreach dir ( frame chem share dyn_em phys main tools wrftladj ) +foreach dir ( frame chem share dyn_em phys cmaq main tools wrftladj ) if ( -d $dir ) then - ( cd $dir ; echo $dir ; /bin/rm -f core wrf *.f90 *.exe *.kmo *.mod *.o *.obj *.inc *.F90 *.a \ + if ( $dir == cmaq ) then + ( cd $dir ; echo $dir ; /bin/rm -f *.o *.mod ) >& /dev/null + else + ( cd $dir ; echo $dir ; /bin/rm -f core wrf *.f90 *.exe *.kmo *.mod *.o *.obj *.inc *.F90 *.a \ db_* Warnings module_state_description.F module_dm.F gmeta \ wrfdata whatiread rsl.* show_domain* ) >& /dev/null + endif endif end diff --git a/configure b/configure index 5e3df42d8c..4010981de6 100755 --- a/configure +++ b/configure @@ -25,6 +25,7 @@ while [ $# -ge 1 ]; do -r8) rword="-r8" ;; -time) shift ; FORTRAN_COMPILER_TIMER=$1 ;; chem) WRF_CHEM=1 ;; + cmaq) WRF_CMAQ=1 ;; kpp) WRF_KPP=1 ;; radardfi) WRF_DFI_RADAR=1 ;; wrfda) wrf_core=DA_CORE ;; @@ -460,6 +461,12 @@ if [ -n "$WRF_DFI_RADAR" ] ; then compileflags="${compileflags}!-DWRF_DFI_RADAR=1" fi fi + +if [ -n "$WRF_CMAQ" ] ; then + echo building WRF with CMAQ option + compileflags="${compileflags}!-DWRF_CMAQ" +fi + if [ -n "$WRF_CHEM" ] ; then if [ $WRF_CHEM = 1 ] ; then echo building WRF with chemistry option diff --git a/dyn_em/module_first_rk_step_part1.F b/dyn_em/module_first_rk_step_part1.F index 456126669a..edac27e779 100644 --- a/dyn_em/module_first_rk_step_part1.F +++ b/dyn_em/module_first_rk_step_part1.F @@ -34,6 +34,7 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & , f_flux & , aerocu & , restart_flag & + , feedback_is_ready & ) USE module_state_description USE module_model_constants @@ -118,6 +119,9 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & LOGICAL, INTENT(IN), OPTIONAL :: f_flux LOGICAL, INTENT(IN) :: restart_flag + + LOGICAL, INTENT(IN) :: feedback_is_ready ! For WRF-CMAQ coupled model, indicates feedback information is available + ! Local real :: HYDRO_dt REAL, DIMENSION( ims:ime, jms:jme ) :: exch_temf ! 1/7/09 WA @@ -450,7 +454,56 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & ,MP_PHYSICS=CONFIG_FLAGS%MP_PHYSICS & & ,EFCG=grid%EFCG,EFCS=grid%EFCS,EFIG=grid%EFIG & & ,EFIS=grid%EFIS,EFSG=grid%EFSG,aercu_opt=config_flags%aercu_opt & - & ,EFSS=grid%EFSS,QS_CU=grid%QS_CU) + & ,EFSS=grid%EFSS,QS_CU=grid%QS_CU & + & ,feedback_is_ready=feedback_is_ready & +#if ( WRF_CMAQ == 1) + & ,mass_ws_i=grid%mass_ws_i & + & ,mass_ws_j=grid%mass_ws_j & + & ,mass_ws_k=grid%mass_ws_k & + & ,mass_in_i=grid%mass_in_i & + & ,mass_in_j=grid%mass_in_j & + & ,mass_in_k=grid%mass_in_k & + & ,mass_ec_i=grid%mass_ec_i & + & ,mass_ec_j=grid%mass_ec_j & + & ,mass_ec_k=grid%mass_ec_k & + & ,mass_ss_i=grid%mass_ss_i & + & ,mass_ss_j=grid%mass_ss_j & + & ,mass_ss_k=grid%mass_ss_k & + & ,mass_h2o_i=grid%mass_h2o_i & + & ,mass_h2o_j=grid%mass_h2o_j & + & ,mass_h2o_k=grid%mass_h2o_k & + & ,dgn_i=grid%dgn_i & + & ,dgn_j=grid%dgn_j & + & ,dgn_k=grid%dgn_k & + & ,sig_i=grid%sig_i & + & ,sig_j=grid%sig_j & + & ,sig_k=grid%sig_k & + & ,sw_gtauxar_01=grid%sw_gtauxar_01 & + & ,sw_gtauxar_02=grid%sw_gtauxar_02 & + & ,sw_gtauxar_03=grid%sw_gtauxar_03 & + & ,sw_gtauxar_04=grid%sw_gtauxar_04 & + & ,sw_gtauxar_05=grid%sw_gtauxar_05 & + & ,sw_ttauxar_01=grid%sw_ttauxar_01 & + & ,sw_ttauxar_02=grid%sw_ttauxar_02 & + & ,sw_ttauxar_03=grid%sw_ttauxar_03 & + & ,sw_ttauxar_04=grid%sw_ttauxar_04 & + & ,sw_ttauxar_05=grid%sw_ttauxar_05 & + & ,sw_asy_fac_01=grid%sw_asy_fac_01 & + & ,sw_asy_fac_02=grid%sw_asy_fac_02 & + & ,sw_asy_fac_03=grid%sw_asy_fac_03 & + & ,sw_asy_fac_04=grid%sw_asy_fac_04 & + & ,sw_asy_fac_05=grid%sw_asy_fac_05 & + & ,sw_ssa_01=grid%sw_ssa_01 & + & ,sw_ssa_02=grid%sw_ssa_02 & + & ,sw_ssa_03=grid%sw_ssa_03 & + & ,sw_ssa_04=grid%sw_ssa_04 & + & ,sw_ssa_05=grid%sw_ssa_05 & + & ,ozone=grid%ozone & + & ,sw_zbbcddir=grid%sw_zbbcddir & + & ,sw_dirdflux=grid%sw_dirdflux & + & ,sw_difdflux=grid%sw_difdflux & +#endif + ) BENCH_END(rad_driver_tim) diff --git a/dyn_em/solve_em.F b/dyn_em/solve_em.F index d6fc975b41..47a38f9af1 100644 --- a/dyn_em/solve_em.F +++ b/dyn_em/solve_em.F @@ -78,6 +78,9 @@ SUBROUTINE solve_em ( grid , config_flags & USE module_llxy, ONLY : proj_cassini USE module_avgflx_em, ONLY : zero_avgflx, upd_avgflx USE module_cpl, ONLY : coupler_on, cpl_settime, cpl_store_input +#if (WRF_CMAQ == 1) + use twoway_data_module +#endif IMPLICIT NONE @@ -118,6 +121,29 @@ SUBROUTINE solve_em ( grid , config_flags & REAL :: t_new, time_duration_of_lbcs +! begin WRF-CMAQ twoway coupled model block + integer :: twoway_jdate, & ! CMAQ current job date + twoway_jtime, & ! CMAQ current job time + met_file_tstep ! MCIP like MET file time step + + integer, save :: cmaq_nstep, & ! total number of CMAQ steps + wrf_end_step, & ! WRF ending step # + counter = -1, & ! step counter + wrf_cmaq_freq, & ! call frequency between WRF and CMAQ + wrf_cmaq_option ! WRF-CMAQ coupled model option + ! 0 = run WRF only + ! 1 = run WRF only w producing MCIP like GRID and MET files + ! 2 = run WRF-CMAQ coupled model w/o producing MCIP like GRID and MET files + ! 3 = run WRF-CMAQ coupled model w producing MCIP like GRID and MET files + + logical :: cmaq_step ! CMAQ step number + + logical, save :: firstime = .true., & ! logical variable indicating first time + feedback_is_ready, & ! logical variable indicating feedback process can proceed + feedback_restart, & ! logical variable indicating feedback information is available + direct_sw_feedback ! logical variable indicating direct aerosol sw feedback is on or not +! end WRF-CMAQ twoway coupled model block + ! Changes in tendency at this timestep real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,grid%sm33:grid%em33) :: h_tendency, & z_tendency @@ -181,6 +207,20 @@ SUBROUTINE solve_em ( grid , config_flags & INTEGER :: num, den TYPE(WRFU_TimeInterval) :: dtInterval, intervaltime,restartinterval +#if (WRF_CMAQ == 1) + interface + SUBROUTINE CMAQ_DRIVER ( MODEL_STDATE, MODEL_STTIME, MODEL_TSTEP, & + MODEL_JDATE, MODEL_JTIME, LAST_STEP, & + COUPLE_TSTEP, NCOLS_IN, NLAYS_IN) + INTEGER, INTENT( IN ) :: MODEL_STDATE, MODEL_STTIME, MODEL_TSTEP + INTEGER, INTENT( OUT ) :: MODEL_JDATE, MODEL_JTIME + LOGICAL, INTENT( IN ) :: LAST_STEP + INTEGER, INTENT( IN ), OPTIONAL :: COUPLE_TSTEP + INTEGER, INTENT( IN ), OPTIONAL :: NCOLS_IN, NLAYS_IN + END SUBROUTINE CMAQ_DRIVER + end interface +#endif + ! Define benchmarking timers if -DBENCH is compiled #include "bench_solve_em_def.h" @@ -222,6 +262,19 @@ SUBROUTINE solve_em ( grid , config_flags & ! Initialize timers if compiled with -DBENCH #include "bench_solve_em_init.h" +#if (WRF_CMAQ == 1) + if (firstime) then + CALL nl_get_feedback_restart ( .false., feedback_restart ) + if (feedback_restart) then + feedback_is_ready = .true. + else + feedback_is_ready = .false. + end if + end if +#else + feedback_is_ready = .false. +#endif + ! set runge-kutta solver (2nd or 3rd order) dynamics_option = config_flags%rk_ord @@ -773,6 +826,7 @@ SUBROUTINE solve_em ( grid , config_flags & , f_flux & , aerocu & , restart_flag & + , feedback_is_ready=feedback_is_ready & ) #ifdef DM_PARALLEL @@ -2823,7 +2877,6 @@ SUBROUTINE solve_em ( grid , config_flags & grid%j_start(ij),grid%j_end(ij), & k_start,k_end) ELSEIF (is.EQ.P_QNN) THEN ! for ntu3m -! IF ( is .EQ. P_QNN ) THEN ! for ntu3m CALL flow_dep_bdy_qnn ( scalar(ims,kms,jms,is), & grid%ru_m, grid%rv_m, config_flags, & grid%spec_zone, & @@ -4657,7 +4710,7 @@ SUBROUTINE solve_em ( grid , config_flags & CALL after_all_rk_steps ( grid, config_flags, & moist, chem, tracer, scalar, & - th_phy, pi_phy, p_phy, & + th_phy, pi_phy, p_phy, rho_phy, & p8w, t8w, dz8w, & REAL(curr_secs,8), curr_secs2, & diag_flag, & @@ -4763,6 +4816,93 @@ SUBROUTINE solve_em ( grid , config_flags & ! Finish timers if compiled with -DBENCH. #include "bench_solve_em_end.h" +#if (WRF_CMAQ == 1) + if (firstime) then + CALL nl_get_wrf_cmaq_option ( 1, wrf_cmaq_option ) + CALL nl_get_wrf_cmaq_freq ( 1, wrf_cmaq_freq ) + CALL nl_get_direct_sw_feedback ( .false., direct_sw_feedback ) + CALL nl_get_met_file_tstep ( 1, met_file_tstep ) + + cmaq_wrf_feedback = direct_sw_feedback + + if (wrf_cmaq_option .gt. 0) then + cmaq_nstep = ((grid%run_days * 24 + grid%run_hours) * 3600 + grid%run_minutes * 60 + grid%run_seconds) / & + (grid%time_step * WRF_CMAQ_FREQ) + + wrf_end_step = cmaq_nstep * WRF_CMAQ_FREQ - 1 + end if + end if + + if (wrf_cmaq_option .gt. 0) then + COUNTER = COUNTER + 1 + + if ( .not. cmaq_wrf_feedback .and. firstime) then + grid%prev_rainnc = 0.0 + grid%prev_rainc = 0.0 + end if + + CMAQ_STEP = (mod(COUNTER, WRF_CMAQ_FREQ) .EQ. 0) + + if (CMAQ_STEP) then + CALL aqprep (grid, config_flags, grid%t_phy, p_phy, & + grid%rho, grid%z_at_w, dz8w, p8w, t8w, & + model_config_rec%num_land_cat, 'V4.1.1', & + wrf_cmaq_option, wrf_cmaq_freq, & + ids, ide, jds, jde, kds, kde, & + ims, ime, jms, jme, kms, kme, & + ips, ipe, jps, jpe, kps, kpe, & + moist(:,:,:,p_qv), & ! optional + moist(:,:,:,p_qc), & ! optional + moist(:,:,:,p_qr), & ! optional + moist(:,:,:,p_qi), & ! optional + moist(:,:,:,p_qs), & ! optional + moist(:,:,:,p_qg) & ! optional + ) + grid%prev_rainnc = grid%rainnc + grid%prev_rainc = grid%rainc + end if + + if ((counter >= 1) .and. (CMAQ_STEP) .and. (wrf_cmaq_option .gt. 1)) then + + CALL CMAQ_DRIVER (cmaq_sdate, cmaq_stime, grid%time_step*WRF_CMAQ_FREQ, & + twoway_jdate, twoway_jtime, .false.) + + if (direct_sw_feedback) then + CALL FEEDBACK_READ (grid, twoway_jdate, twoway_jtime) + feedback_is_ready = .true. + end if + + end if + +! call aqprep and cmaq one last time before the entire twoway model ends + if (wrf_end_step == counter) then + CALL aqprep (grid, config_flags, grid%t_phy, p_phy, & + grid%rho, grid%z_at_w, dz8w, p8w, t8w, & + model_config_rec%num_land_cat, 'V4.1.1', & + wrf_cmaq_option, wrf_cmaq_freq, & + ids, ide, jds, jde, kds, kde, & + ims, ime, jms, jme, kms, kme, & + ips, ipe, jps, jpe, kps, kpe, & + moist(:,:,:,p_qv), & ! optional + moist(:,:,:,p_qc), & ! optional + moist(:,:,:,p_qr), & ! optional + moist(:,:,:,p_qi), & ! optional + moist(:,:,:,p_qs), & ! optional + moist(:,:,:,p_qg) & ! optional + ) + + if (wrf_cmaq_option .gt. 1) then + + CALL CMAQ_DRIVER (cmaq_sdate, cmaq_stime, grid%time_step*WRF_CMAQ_FREQ, & + twoway_jdate, twoway_jtime, .true.) + end if + + end if + + end if + firstime = .false. +#endif + RETURN END SUBROUTINE solve_em diff --git a/external/io_netcdf/makefile b/external/io_netcdf/makefile index fa638d455f..9ac0d94eab 100644 --- a/external/io_netcdf/makefile +++ b/external/io_netcdf/makefile @@ -29,7 +29,7 @@ wrf_io.o: wrf_io.F90 $(CODE) if [ $$a -a "$$WRFIO_NCD_LARGE_FILE_SUPPORT" = "1" ] ; then \ $(CPP1) -DWRFIO_NCD_LARGE_FILE_SUPPORT -I../ioapi_share wrf_io.F90 | $(M4) - > wrf_io.f ; \ else \ - $(CPP1) -I../ioapi_share wrf_io.F90 | $(M4) - > wrf_io.f ; \ + $(CPP1) -DWRFIO_NCD_LARGE_FILE_SUPPORT -I../ioapi_share wrf_io.F90 | $(M4) - > wrf_io.f ; \ fi $(FC) -o $@ $(FFLAGS) -c wrf_io.f diff --git a/main/Makefile b/main/Makefile index 5e8e1c5cfe..9efb8aed43 100644 --- a/main/Makefile +++ b/main/Makefile @@ -14,7 +14,7 @@ include ../configure.wrf $(SOLVER)_wrf : wrf.o ../main/module_wrf_top.o $(RANLIB) $(RLFLAGS) $(LIBWRFLIB) - $(LD) -o wrf.exe $(LDFLAGS) wrf.o ../main/module_wrf_top.o $(LIBWRFLIB) $(LIB) + $(LD) -o wrf.exe $(LDFLAGS) wrf.o ../main/module_wrf_top.o $(LIBWRFLIB) $(CMAQLIB) $(LIB) $(SOLVER)_wrfplus : wrf.o ../main/module_wrf_top.o $(RANLIB) $(RLFLAGS) $(LIBWRFLIB) diff --git a/main/depend.common b/main/depend.common index 831a28381e..e7eec48c15 100644 --- a/main/depend.common +++ b/main/depend.common @@ -487,6 +487,8 @@ module_sf_ruclsm.o: ../frame/module_wrf_error.o module_data_gocart_dust.o module_sf_pxlsm.o: ../share/module_model_constants.o module_sf_pxlsm_data.o +module_twoway_ra_rrtmg_sw.o: module_twoway_rrtmg_aero_optical_util.o + module_ra_rrtmg_sw.o: module_ra_rrtmg_lw.o module_ra_rrtmg_swf.o: module_ra_rrtmg_lwf.o module_ra_rrtmg_swk.o: module_ra_rrtmg_lwk.o module_ra_effective_radius.o @@ -703,6 +705,7 @@ module_radiation_driver.o: \ module_ra_rrtm.o \ module_ra_rrtmg_lw.o \ module_ra_rrtmg_sw.o \ + module_twoway_ra_rrtmg_sw.o \ module_ra_rrtmg_lwf.o \ module_ra_rrtmg_swf.o \ module_ra_rrtmg_lwk.o \ diff --git a/phys/Makefile b/phys/Makefile index 77da6da046..3bdf54d7ff 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -7,7 +7,9 @@ RM = rm -f MODULES = \ module_bep_bem_helper.o \ + complex_number_module.o \ module_cam_shr_kind_mod.o \ + module_twoway_rrtmg_aero_optical_util.o \ module_cam_support.o \ module_cam_shr_const_mod.o \ module_cam_physconst.o \ diff --git a/phys/complex_number_module.F b/phys/complex_number_module.F new file mode 100644 index 0000000000..3c4e1da9fc --- /dev/null +++ b/phys/complex_number_module.F @@ -0,0 +1,253 @@ +!------------------------------------------------------------------------! +! The Community Multiscale Air Quality (CMAQ) system software is in ! +! continuous development by various groups and is based on information ! +! from these groups: Federal Government employees, contractors working ! +! within a United States Government contract, and non-Federal sources ! +! including research institutions. These groups give the Government ! +! permission to use, prepare derivative works of, and distribute copies ! +! of their work in the CMAQ system to the public and to permit others ! +! to do so. The United States Environmental Protection Agency ! +! therefore grants similar permission to use the CMAQ system software, ! +! but users are requested to provide copies of derivative works or ! +! products designed to operate in the CMAQ system to the United States ! +! Government without restrictions as to use by others. Software ! +! that is used with the CMAQ system but distributed under the GNU ! +! General Public License or the GNU Lesser General Public License is ! +! subject to their copyright restrictions. ! +!------------------------------------------------------------------------! + +! complex number general function + +! Revision History: +! 2012/07/31 David Wong Original version +! 2012/10/22 David Wong Added treatment to avoid division by 0 +! 2013/11/20 David Wong modified the way to compute +! max(abs(c_div_cc%real_part), min_val) in subroutine +! c_div_cc to satisfy absoft compiler requirement +! 2015/12/16 David Wong renamed argument list for c_add_cr, c_add_rc, +! c_sub_cr and c_sub_rc to avoid gfortran not able to +! distinguish those routines in the interface block issue +! 2016/02/23 David Wong extracted the entire module and put it in a +! file alone. + + module complex_number_module + + implicit none + +! integer, parameter :: loc_real_precision = selected_real_kind(p=16, r=60) + integer, parameter :: loc_real_precision = 8 + + real (kind=loc_real_precision), parameter, private :: min_val = 1.0e-30_loc_real_precision + + type complex_number + real(kind=loc_real_precision) :: real_part, imag_part + end type complex_number + + interface c_add + module procedure c_add_cc, & ! z1 + z2 + c_add_cr, & ! z1 + num, where num is a real number + c_add_rc ! num + z1, where num is a real number + end interface + + interface c_sub + module procedure c_sub_cc, & ! z1 - z2 + c_sub_cr, & ! z1 - num, where num is a real number + c_sub_rc ! num - z1, where num is a real number + end interface + + interface c_mul + module procedure c_mul_cc, & ! z1 * z2 + c_mul_rc ! num * z1, where num is a real number + end interface + + interface c_div + module procedure c_div_cc, & ! z1 / z2 + c_div_rc ! num / z1, where num is a real number + end interface + + contains + +! -------------------------------------------------------------------------- + type (complex_number) function c_set (x, y) + +! initialize a complex number + + real(kind=loc_real_precision), intent(in) :: x, y + + character (len = 80) :: str + + write (str, *) x + read(str, *) c_set%real_part + write (str, *) y + read(str, *) c_set%imag_part + + end function c_set + +! -------------------------------------------------------------------------- + type (complex_number) function c_add_cc (z1, z2) + + type (complex_number), intent(in) :: z1, z2 + + c_add_cc%real_part = z1%real_part + z2%real_part + c_add_cc%imag_part = z1%imag_part + z2%imag_part + + end function c_add_cc + +! -------------------------------------------------------------------------- + type (complex_number) function c_add_cr (z3, num1) + + type (complex_number), intent(in) :: z3 + real(kind=loc_real_precision), intent(in) :: num1 + + c_add_cr%real_part = z3%real_part + num1 + c_add_cr%imag_part = z3%imag_part + + end function c_add_cr + +! -------------------------------------------------------------------------- + type (complex_number) function c_add_rc (num2, z4) + + type (complex_number), intent(in) :: z4 + real(kind=loc_real_precision), intent(in) :: num2 + + c_add_rc%real_part = z4%real_part + num2 + c_add_rc%imag_part = z4%imag_part + + end function c_add_rc + +! -------------------------------------------------------------------------- + type (complex_number) function c_sub_cc (z1, z2) + + type (complex_number), intent(in) :: z1, z2 + + c_sub_cc%real_part = z1%real_part - z2%real_part + c_sub_cc%imag_part = z1%imag_part - z2%imag_part + + end function c_sub_cc + +! -------------------------------------------------------------------------- + type (complex_number) function c_sub_cr (z3, num1) + + type (complex_number), intent(in) :: z3 + real(kind=loc_real_precision), intent(in) :: num1 + + c_sub_cr%real_part = z3%real_part - num1 + c_sub_cr%imag_part = z3%imag_part + + end function c_sub_cr + +! -------------------------------------------------------------------------- + type (complex_number) function c_sub_rc (num2, z4) + + type (complex_number), intent(in) :: z4 + real(kind=loc_real_precision), intent(in) :: num2 + + c_sub_rc%real_part = num2 - z4%real_part + c_sub_rc%imag_part = - z4%imag_part + + end function c_sub_rc + +! -------------------------------------------------------------------------- + type (complex_number) function c_mul_cc (z1, z2) + + type (complex_number), intent(in) :: z1, z2 + + c_mul_cc%real_part = z1%real_part * z2%real_part & + - z1%imag_part * z2%imag_part + c_mul_cc%imag_part = z1%real_part * z2%imag_part & + + z1%imag_part * z2%real_part + + end function c_mul_cc + +! -------------------------------------------------------------------------- + type (complex_number) function c_mul_rc (x, z1) + + type (complex_number), intent(in) :: z1 + real(kind=loc_real_precision), intent(in) :: x + + c_mul_rc%real_part = z1%real_part * x + c_mul_rc%imag_part = z1%imag_part * x + + end function c_mul_rc + +! -------------------------------------------------------------------------- + type (complex_number) function c_div_cc (z1, z2) + + type (complex_number), intent(in) :: z1, z2 + + real(kind=loc_real_precision) :: denom + real(kind=loc_real_precision) :: temp(2) + + denom = 1.0 / ( z2%real_part * z2%real_part & + + z2%imag_part * z2%imag_part) + + c_div_cc%real_part = ( z1%real_part * z2%real_part & + + z1%imag_part * z2%imag_part) * denom + + temp(1) = abs(c_div_cc%real_part) + temp(2) = min_val + + c_div_cc%real_part = sign(maxval(temp), c_div_cc%real_part) + + c_div_cc%imag_part = ( z1%imag_part * z2%real_part & + - z1%real_part * z2%imag_part) * denom + + temp(1) = abs(c_div_cc%imag_part) + + c_div_cc%imag_part = sign(maxval(temp), c_div_cc%imag_part) + + end function c_div_cc + +! -------------------------------------------------------------------------- + type (complex_number) function c_div_rc (num, z1) + +! compute 1 / z1 + + real(kind=loc_real_precision), intent(in) :: num + type (complex_number), intent(in) :: z1 + + real(kind=loc_real_precision) :: denom, temp + + temp = z1%real_part * z1%real_part + z1%imag_part * z1%imag_part + temp = sign(max(abs(temp), min_val), temp) + + denom = num / temp + c_div_rc%real_part = z1%real_part * denom + c_div_rc%imag_part = -1.0 * z1%imag_part * denom + + end function c_div_rc + +! -------------------------------------------------------------------------- + type (complex_number) function c_sin (z1) + +! compute sin of a complex number + + type (complex_number), intent(in) :: z1 + + c_sin%real_part = sin(z1%real_part) * cosh(z1%imag_part) + c_sin%imag_part = cos(z1%real_part) * sinh(z1%imag_part) + + end function c_sin + +! -------------------------------------------------------------------------- + type (complex_number) function c_cos (z1) + + type (complex_number), intent(in) :: z1 + + c_cos%real_part = cos(z1%real_part) * cosh(z1%imag_part) + c_cos%imag_part = -1.0 * sin(z1%real_part) * sinh(z1%imag_part) + + end function c_cos + +! -------------------------------------------------------------------------- + real(kind=loc_real_precision) function c_abs (z1) + +! computer absolute value of a complex number + + type (complex_number), intent(in) :: z1 + + c_abs = sqrt(z1%real_part**2 + z1%imag_part**2) + + end function c_abs + + end module complex_number_module diff --git a/phys/module_ra_rrtmg_sw.F b/phys/module_ra_rrtmg_sw.F index 95905fb53c..0bf899c077 100644 --- a/phys/module_ra_rrtmg_sw.F +++ b/phys/module_ra_rrtmg_sw.F @@ -8857,6 +8857,7 @@ subroutine rrtmg_sw & swdkdir,swdkdif, & ! jararias, 2013/08/10 swdkdirc & ! PAJ ,calc_clean_atm_diag & + ,sw_zbbcddir, sw_dirdflux, sw_difdflux & ! WRF-CMAQ twoway coupled model ) @@ -9090,6 +9091,10 @@ subroutine rrtmg_sw & swdkdif(:,:), & ! Total shortwave downward diffuse flux (W/m2), Dimensions: (ncol,nlay) jararias, 2013/08/10 swdkdirc(:,:) ! Total shortwave downward direct flux clear sky (W/m2), Dimensions: (ncol,nlay) + real, intent(out) :: sw_zbbcddir, & ! WRF-CMAQ twoway coupled model + sw_dirdflux, & ! WRF-CMAQ twoway coupled model + sw_difdflux ! WRF-CMAQ twoway coupled model + @@ -9594,6 +9599,12 @@ subroutine rrtmg_sw & ! End longitude loop enddo +! begin WRF-CMAQ twoway coupled model block + sw_zbbcddir = zbbcddir(1) + sw_dirdflux = dirdflux(1) + sw_difdflux = difdflux(1) +! end WRF-CMAQ twoway coupled model block + end subroutine rrtmg_sw !************************************************************************* @@ -10001,6 +10012,7 @@ MODULE module_ra_rrtmg_sw use module_ra_rrtmg_lw, only : inirad, o3data, relcalc, reicalc, retab ! mcica_random_numbers, randomNumberSequence, & ! new_RandomNumberSequence, getRandomReal +use module_twoway_ra_rrtmg_sw CONTAINS @@ -10051,8 +10063,26 @@ SUBROUTINE RRTMG_SWRAD( & tauaer3d_sw,ssaaer3d_sw,asyaer3d_sw, & ! jararias 2013/11 swddir, swddni, swddif, & ! jararias 2013/08 swdownc, swddnic, swddirc, & ! PAJ - xcoszen,yr,julian, & ! jararias 2013/08 - obscur & ! amontornes-bcodina 2015/09 solar eclipses + xcoszen,yr,julian, & ! jararias 2013/08 + obscur, & ! amontornes-bcodina 2015/09 solar eclipses + proceed_twoway_sw, & ! WRF-CMAQ twoway coupled model + nmode, & ! WRF-CMAQ twoway coupled model + mass_ws_i, mass_ws_j, mass_ws_k, & ! WRF-CMAQ twoway coupled model + mass_in_i, mass_in_j, mass_in_k, & ! WRF-CMAQ twoway coupled model + mass_ec_i, mass_ec_j, mass_ec_k, & ! WRF-CMAQ twoway coupled model + mass_ss_i, mass_ss_j, mass_ss_k, & ! WRF-CMAQ twoway coupled model + mass_h2o_i, mass_h2o_j, mass_h2o_k, & ! WRF-CMAQ twoway coupled model + dgn_i, dgn_j, dgn_k, & ! WRF-CMAQ twoway coupled model + sig_i, sig_j, sig_k, & ! WRF-CMAQ twoway coupled model + gtauxar_01, gtauxar_02, gtauxar_03, & ! WRF-CMAQ twoway coupled model + gtauxar_04, gtauxar_05, & ! WRF-CMAQ twoway coupled model + asy_fac_01, asy_fac_02, asy_fac_03, & ! WRF-CMAQ twoway coupled model + asy_fac_04, asy_fac_05, & ! WRF-CMAQ twoway coupled model + ssa_01, ssa_02, ssa_03, & ! WRF-CMAQ twoway coupled model + ssa_04, ssa_05 & ! WRF-CMAQ twoway coupled model + ,sw_zbbcddir & ! WRF-CMAQ twoway coupled model + ,sw_dirdflux & ! WRF-CMAQ twoway coupled model + ,sw_difdflux & ! WRF-CMAQ twoway coupled model ) !------------------------------------------------------------------ USE MODULE_RA_CLWRF_SUPPORT, ONLY : read_CAMgases @@ -10223,6 +10253,49 @@ SUBROUTINE RRTMG_SWRAD( & ! obscur --> degree of obscuration for solar eclipses prediction (2D) REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: obscur +! begin WRF-CMAQ twoway coupled model block + LOGICAL, INTENT(IN) :: proceed_twoway_sw + +! ** FSB items needed for new aerosol code from CMAQ + integer, optional, intent(in) :: nmode ! number of log-normal modes + + real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(in) :: & + mass_ws_i, mass_ws_j, mass_ws_k, & ! mass cocentrations in [ ug/m**3 ] for water + ! soluble species in each mode + mass_in_i, mass_in_j, mass_in_k, & ! mass cocentrations in [ ug/m**3 ] for water + ! insoluble species in each mode + mass_ec_i, mass_ec_j, mass_ec_k, & ! mass cocentrations in [ ug/m**3 ] for elemental + ! carbon species in each mode + mass_ss_i, mass_ss_j, mass_ss_k, & ! mass cocentrations in [ ug/m**3 ] for aerosol + ! water species in each mode + mass_h2o_i, mass_h2o_j, mass_h2o_k, & ! mass cocentrations in [ ug/m**3 ] for sea + ! salt species in each mode + dgn_i, dgn_j, dgn_k, & ! geometric mean diameter of each mode [ m ] + sig_i, sig_j, sig_k ! geometric standard deviation of each mode + + real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(out) :: & + gtauxar_01, & ! Aerosol optical depth of RRTMG SW band 11 + gtauxar_02, & ! Aerosol optical depth of RRTMG SW band 10 + gtauxar_03, & ! Aerosol optical depth of RRTMG SW band 9 + gtauxar_04, & ! Aerosol optical depth of RRTMG SW band 8 + gtauxar_05, & ! Aerosol optical depth of RRTMG SW band 7 + asy_fac_01, & ! asymmetry factor of RRTMG SW band 11 + asy_fac_02, & ! asymmetry factor of RRTMG SW band 10 + asy_fac_03, & ! asymmetry factor of RRTMG SW band 9 + asy_fac_04, & ! asymmetry factor of RRTMG SW band 8 + asy_fac_05, & ! asymmetry factor of RRTMG SW band 7 + ssa_01, & ! single scattering albedo of RRTMG SW band 11 + ssa_02, & ! single scattering albedo of RRTMG SW band 10 + ssa_03, & ! single scattering albedo of RRTMG SW band 9 + ssa_04, & ! single scattering albedo of RRTMG SW band 8 + ssa_05 ! single scattering albedo of RRTMG SW band 7 + + REAL, DIMENSION( ims:ime, jms:jme ), optional, INTENT(OUT) :: & + sw_zbbcddir, & ! clearsky downward direct shortwave flux (w/m2) + sw_dirdflux, & ! Direct downward shortwave surface flux + sw_difdflux ! Diffuse downward shortwave surface flux +! end WRF-CMAQ twoway coupled model block + ! LOCAL VARS REAL, DIMENSION( kts:kte+1 ) :: Pw1D, & @@ -10421,6 +10494,34 @@ SUBROUTINE RRTMG_SWRAD( & REAL :: da, eot ! jararias, 14/08/2013 +! begin WRF-CMAQ twoway coupled model block +#if ( WRF_CMAQ == 1 ) + REAL, DIMENSION (3) :: INMASS_ws, & ! holds mass cocentrations in [ ug/m**3 ] for + ! water soluble species in all three modes + INMASS_in, & ! holds mass cocentrations in [ ug/m**3 ] for + ! water insoluble species in all three modes + INMASS_ec, & ! holds mass cocentrations in [ ug/m**3 ] for + ! elemental carbon species in all three modes + INMASS_ss, & ! holds mass cocentrations in [ ug/m**3 ] for + ! aerosol water species in all three modes + INMASS_h2o, & ! holds mass cocentrations in [ ug/m**3 ] for + ! sea salt species in all three modes + INDGN, & ! holds geometric mean diameter in all three modes + INSIG ! holds geometric standard deviation in all three modes +#endif + REAL :: xtauaer, & ! temporary variable for Aerosol Optical Depth + waer, & ! temporary variable for single scattering albedo + gaer, & ! temporary variable for symmetry factor + delta_z, & ! layer thickness + loc_sw_zbbcddir, & ! clearsky downward direct shortwave flux (w/m2) + loc_sw_dirdflux, & ! Direct downward shortwave surface flux + loc_sw_difdflux ! Diffuse downward shortwave surface flux + + INTEGER :: modes ! number of modes + + character (len = 50) :: mystr ! temporary character string +! end WRF-CMAQ twoway coupled model block + !------------------------------------------------------------------ ! Annual function for co2 in WRF v4.2 co2 = (280. + 90.*exp(0.02*(yr-2000)))*1.e-6 @@ -11131,6 +11232,91 @@ SUBROUTINE RRTMG_SWRAD( & ! ... Aerosol effects. Added aerosol feedbacks from Chem , 03/2010 -czhao ! +#if ( WRF_CMAQ == 1 ) + do nb = 1, nbndsw + do k = kts,kte+1 + tauaer(ncol,k,nb) = 0. + ssaaer(ncol,k,nb) = 1. + asmaer(ncol,k,nb) = 0. + + if (proceed_twoway_sw) then + INMASS_ws(1) = mass_ws_i(i,k,j) + INMASS_ws(2) = mass_ws_j(i,k,j) + INMASS_ws(3) = mass_ws_k(i,k,j) + INMASS_in(1) = mass_in_i(i,k,j) + INMASS_in(2) = mass_in_j(i,k,j) + INMASS_in(3) = mass_in_k(i,k,j) + INMASS_ec(1) = mass_ec_i(i,k,j) + INMASS_ec(2) = mass_ec_j(i,k,j) + INMASS_ec(3) = mass_ec_k(i,k,j) + INMASS_ss(1) = mass_ss_i(i,k,j) + INMASS_ss(2) = mass_ss_j(i,k,j) + INMASS_ss(3) = mass_ss_k(i,k,j) + INMASS_h2o(1) = mass_h2o_i(i,k,j) + INMASS_h2o(2) = mass_h2o_j(i,k,j) + INMASS_h2o(3) = mass_h2o_k(i,k,j) + INDGN(1) = dgn_i(i,k,j) + INDGN(2) = dgn_j(i,k,j) + INDGN(3) = dgn_k(i,k,j) + INSIG(1) = sig_i(i,k,j) + INSIG(2) = sig_j(i,k,j) + INSIG(3) = sig_k(i,k,j) + + delta_z = dz8w(i,k,j) + + call get_aerosol_Optics_RRTMG_SW( nb,nmode,delta_z, & + INMASS_ws, INMASS_in, INMASS_ec, INMASS_ss, & + INMASS_h2o, INDGN, INSIG, & + xtauaer, waer, gaer ) + + write (mystr, *) xtauaer + if (trim(mystr) == ' NaN') then + write (6, '(a13, 2i5)') ' ==d== ', nb, nmode + write (6, '(a13, 5e18.10)') ' ==d== delta ', delta_z + write (6, '(a13, 5e18.10)') ' ==d== ws ', INMASS_ws + write (6, '(a13, 5e18.10)') ' ==d== in ', INMASS_in + write (6, '(a13, 5e18.10)') ' ==d== ec ', INMASS_ec + write (6, '(a13, 5e18.10)') ' ==d== ss ', INMASS_ss + write (6, '(a13, 5e18.10)') ' ==d== h2o ', INMASS_h2o + write (6, '(a13, 5e18.10)') ' ==d== indgn ', INDGN + write (6, '(a13, 5e18.10)') ' ==d== insig ', INSIG + end if + + if (nb == 11) then + gtauxar_01 (i,k,j) = xtauaer + asy_fac_01 (i,k,j) = gaer + ssa_01 (i,k,j) = waer + else if (nb == 10) then + gtauxar_02 (i,k,j) = xtauaer + asy_fac_02 (i,k,j) = gaer + ssa_02 (i,k,j) = waer + else if (nb == 9) then + gtauxar_03 (i,k,j) = xtauaer + asy_fac_03 (i,k,j) = gaer + ssa_03 (i,k,j) = waer + else if (nb == 8) then + gtauxar_04 (i,k,j) = xtauaer + asy_fac_04 (i,k,j) = gaer + ssa_04 (i,k,j) = waer + else if (nb == 7) then + gtauxar_05 (i,k,j) = xtauaer + asy_fac_05 (i,k,j) = gaer + ssa_05 (i,k,j) = waer + end if + + tauaer(ncol,k,nb) = xtauaer + ssaaer(ncol,k,nb) = waer + asmaer(ncol,k,nb) = gaer + end if + enddo ! loop over layers + if (proceed_twoway_sw) then +! No aerosols in top layer above model top (kte+1). + tauaer(ncol, kte+1 ,nb) = 0. + ssaaer(ncol, kte+1 ,nb) = 1. + asmaer(ncol, kte+1 ,nb) = 0. + end if + enddo ! loop over wavelengths +#else do nb = 1, nbndsw do k = kts,kte+1 tauaer(ncol,k,nb) = 0. @@ -11149,6 +11335,7 @@ SUBROUTINE RRTMG_SWRAD( & end do end do end if +#endif #if ( WRF_CHEM == 1 ) IF ( AER_RA_FEEDBACK == 1) then @@ -11278,8 +11465,17 @@ SUBROUTINE RRTMG_SWRAD( & swdkdir, swdkdif, & ! jararias, 2012/08/10 swdkdirc & ! PAJ ,calc_clean_atm_diag & + ,loc_sw_zbbcddir & ! WRF-CMAQ twoway coupled model + ,loc_sw_dirdflux & ! WRF-CMAQ twoway coupled model + ,loc_sw_difdflux & ! WRF-CMAQ twoway coupled model ) + ! WRF-CMAQ twoway coupled model + if (present(sw_zbbcddir)) then + sw_zbbcddir(i,j) = loc_sw_zbbcddir + sw_dirdflux(i,j) = loc_sw_dirdflux + sw_difdflux(i,j) = loc_sw_difdflux + end if ! Output net absorbed shortwave surface flux and shortwave cloud forcing ! at the top of atmosphere (W/m2) @@ -11337,6 +11533,25 @@ SUBROUTINE RRTMG_SWRAD( & rthratenswc(i,k,j) = tten1d(k)/pi3d(i,k,j) enddo else + + if (proceed_twoway_sw) then ! this is for WRF-CMAQ twoway coupled model + gtauxar_01 (i,:,j) = 0.0 + gtauxar_02 (i,:,j) = 0.0 + gtauxar_03 (i,:,j) = 0.0 + gtauxar_04 (i,:,j) = 0.0 + gtauxar_05 (i,:,j) = 0.0 + asy_fac_01 (i,:,j) = 0.0 + asy_fac_02 (i,:,j) = 0.0 + asy_fac_03 (i,:,j) = 0.0 + asy_fac_04 (i,:,j) = 0.0 + asy_fac_05 (i,:,j) = 0.0 + ssa_01 (i,:,j) = 0.0 + ssa_02 (i,:,j) = 0.0 + ssa_04 (i,:,j) = 0.0 + ssa_04 (i,:,j) = 0.0 + ssa_05 (i,:,j) = 0.0 + end if + if (present(swupt)) then ! Output up and down toa fluxes for total and clear sky swupt(i,j) = 0. diff --git a/phys/module_radiation_driver.F b/phys/module_radiation_driver.F index 3b11385c96..3187b9b09e 100644 --- a/phys/module_radiation_driver.F +++ b/phys/module_radiation_driver.F @@ -181,6 +181,52 @@ SUBROUTINE radiation_driver ( & ,chem_opt & ,gsfcrad_gocart_coupling & #endif + ,feedback_is_ready & ! WRF-CMAQ twoway coupled model + ,mass_ws_i & ! WRF-CMAQ twoway coupled model + ,mass_ws_j & ! WRF-CMAQ twoway coupled model + ,mass_ws_k & ! WRF-CMAQ twoway coupled model + ,mass_in_i & ! WRF-CMAQ twoway coupled model + ,mass_in_j & ! WRF-CMAQ twoway coupled model + ,mass_in_k & ! WRF-CMAQ twoway coupled model + ,mass_ec_i & ! WRF-CMAQ twoway coupled model + ,mass_ec_j & ! WRF-CMAQ twoway coupled model + ,mass_ec_k & ! WRF-CMAQ twoway coupled model + ,mass_ss_i & ! WRF-CMAQ twoway coupled model + ,mass_ss_j & ! WRF-CMAQ twoway coupled model + ,mass_ss_k & ! WRF-CMAQ twoway coupled model + ,mass_h2o_i & ! WRF-CMAQ twoway coupled model + ,mass_h2o_j & ! WRF-CMAQ twoway coupled model + ,mass_h2o_k & ! WRF-CMAQ twoway coupled model + ,dgn_i & ! WRF-CMAQ twoway coupled model + ,dgn_j & ! WRF-CMAQ twoway coupled model + ,dgn_k & ! WRF-CMAQ twoway coupled model + ,sig_i & ! WRF-CMAQ twoway coupled model + ,sig_j & ! WRF-CMAQ twoway coupled model + ,sig_k & ! WRF-CMAQ twoway coupled model + ,sw_gtauxar_01 & ! WRF-CMAQ twoway coupled model + ,sw_gtauxar_02 & ! WRF-CMAQ twoway coupled model + ,sw_gtauxar_03 & ! WRF-CMAQ twoway coupled model + ,sw_gtauxar_04 & ! WRF-CMAQ twoway coupled model + ,sw_gtauxar_05 & ! WRF-CMAQ twoway coupled model + ,sw_ttauxar_01 & ! WRF-CMAQ twoway coupled model + ,sw_ttauxar_02 & ! WRF-CMAQ twoway coupled model + ,sw_ttauxar_03 & ! WRF-CMAQ twoway coupled model + ,sw_ttauxar_04 & ! WRF-CMAQ twoway coupled model + ,sw_ttauxar_05 & ! WRF-CMAQ twoway coupled model + ,sw_asy_fac_01 & ! WRF-CMAQ twoway coupled model + ,sw_asy_fac_02 & ! WRF-CMAQ twoway coupled model + ,sw_asy_fac_03 & ! WRF-CMAQ twoway coupled model + ,sw_asy_fac_04 & ! WRF-CMAQ twoway coupled model + ,sw_asy_fac_05 & ! WRF-CMAQ twoway coupled model + ,sw_ssa_01 & ! WRF-CMAQ twoway coupled model + ,sw_ssa_02 & ! WRF-CMAQ twoway coupled model + ,sw_ssa_03 & ! WRF-CMAQ twoway coupled model + ,sw_ssa_04 & ! WRF-CMAQ twoway coupled model + ,sw_ssa_05 & ! WRF-CMAQ twoway coupled model + ,ozone & ! WRF-CMAQ twoway coupled model + ,sw_zbbcddir & ! WRF-CMAQ twoway coupled model + ,sw_dirdflux & ! WRF-CMAQ twoway coupled model + ,sw_difdflux & ! WRF-CMAQ twoway coupled model ) @@ -222,6 +268,10 @@ SUBROUTINE radiation_driver ( & ! *** add new modules of schemes here +#if ( WRF_CMAQ == 1) + USE module_twoway_ra_rrtmg_sw +#endif + USE module_ra_sw , ONLY : swrad USE module_ra_gsfcsw , ONLY : gsfcswrad USE module_ra_rrtm , ONLY : rrtmlwrad @@ -244,6 +294,7 @@ SUBROUTINE radiation_driver ( & USE module_ra_aerosol , ONLY : calc_aerosol_goddard_sw, & calc_aerosol_rrtmg_sw USE module_ra_farms , ONLY : farms_driver + ! amontornes-bcodina 2015/09 solar eclipses USE module_ra_eclipse , ONLY : solar_eclipse @@ -827,6 +878,46 @@ SUBROUTINE radiation_driver ( & REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & INTENT(INOUT) :: o3rad + ! begin WRF-CMAQ coupled model block + REAL, DIMENSION (ims:ime, kms:kme, jms:jme ), optional, & + INTENT(INOUT) :: mass_ws_i, mass_ws_j, mass_ws_k, & + mass_in_i, mass_in_j, mass_in_k, & + mass_ec_i, mass_ec_j, mass_ec_k, & + mass_ss_i, mass_ss_j, mass_ss_k, & + mass_h2o_i, mass_h2o_j, mass_h2o_k, & + dgn_i, dgn_j, dgn_k, & + sig_i, sig_j, sig_k + + REAL, DIMENSION (ims:ime, kms:kme, jms:jme ), optional, INTENT(OUT) :: sw_gtauxar_01, & + sw_gtauxar_02, & + sw_gtauxar_03, & + sw_gtauxar_04, & + sw_gtauxar_05, & + sw_asy_fac_01, & + sw_asy_fac_02, & + sw_asy_fac_03, & + sw_asy_fac_04, & + sw_asy_fac_05, & + sw_ssa_01, & + sw_ssa_02, & + sw_ssa_03, & + sw_ssa_04, & + sw_ssa_05 + + REAL, DIMENSION( ims:ime, jms:jme ), optional, INTENT(OUT) :: sw_ttauxar_01, & + sw_ttauxar_02, & + sw_ttauxar_03, & + sw_ttauxar_04, & + sw_ttauxar_05 + + REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), optional, INTENT(IN) :: ozone + REAL, DIMENSION( ims:ime, jms:jme ), optional, INTENT(OUT) :: sw_zbbcddir, & + sw_dirdflux, & + sw_difdflux + + LOGICAL, INTENT(IN) :: feedback_is_ready + ! end WRF-CMAQ coupled model block + ! vert nesting REAL, OPTIONAL , INTENT(IN ) :: p_top REAL :: p_top_dummy @@ -903,6 +994,22 @@ SUBROUTINE radiation_driver ( & !---------- local test vars ! real, dimension(ims:ime, kms:kme, jms:jme) :: hlw, hsw + LOGICAL :: proceed_twoway_sw + + logical, save :: firstime = .true. + logical, save :: feedback_restart, direct_sw_feedback + +#if ( WRF_CMAQ == 1 ) + if (firstime) then + firstime = .false. + CALL nl_get_direct_sw_feedback ( .false. , direct_sw_feedback ) + CALL nl_get_feedback_restart ( .false. , feedback_restart ) + end if +#else + direct_sw_feedback = .false. + feedback_restart = .false. +#endif + ! This allows radiation schemes to correctly ! interface with the convection scheme when convection is only ! enabled in some domains: @@ -2391,6 +2498,43 @@ SUBROUTINE radiation_driver ( & CALL wrf_debug(100, 'CALL rrtmg_sw') + if ( direct_sw_feedback .and. feedback_is_ready ) then + proceed_twoway_sw = .true. + + if (present(mass_ws_i)) then + + mass_ws_i(:,kte+1,:) = mass_ws_i(:,kte,:) + mass_ws_j(:,kte+1,:) = mass_ws_j(:,kte,:) + mass_ws_k(:,kte+1,:) = mass_ws_k(:,kte,:) + + mass_in_i(:,kte+1,:) = mass_in_i(:,kte,:) + mass_in_j(:,kte+1,:) = mass_in_j(:,kte,:) + mass_in_k(:,kte+1,:) = mass_in_k(:,kte,:) + + mass_ec_i(:,kte+1,:) = mass_ec_i(:,kte,:) + mass_ec_j(:,kte+1,:) = mass_ec_j(:,kte,:) + mass_ec_k(:,kte+1,:) = mass_ec_k(:,kte,:) + + mass_ss_i(:,kte+1,:) = mass_ss_i(:,kte,:) + mass_ss_j(:,kte+1,:) = mass_ss_j(:,kte,:) + mass_ss_k(:,kte+1,:) = mass_ss_k(:,kte,:) + + mass_h2o_i(:,kte+1,:) = mass_h2o_i(:,kte,:) + mass_h2o_j(:,kte+1,:) = mass_h2o_j(:,kte,:) + mass_h2o_k(:,kte+1,:) = mass_h2o_k(:,kte,:) + + dgn_i(:,kte+1,:) = dgn_i(:,kte,:) + dgn_j(:,kte+1,:) = dgn_j(:,kte,:) + dgn_k(:,kte+1,:) = dgn_k(:,kte,:) + + sig_i(:,kte+1,:) = sig_i(:,kte,:) + sig_j(:,kte+1,:) = sig_j(:,kte,:) + sig_k(:,kte+1,:) = sig_k(:,kte,:) + end if + else + proceed_twoway_sw = .false. + end if + CALL RRTMG_SWRAD( & RTHRATENSW=RTHRATENSW, & RTHRATENSWC=RTHRATENSWC, & @@ -2453,7 +2597,62 @@ SUBROUTINE radiation_driver ( & swddir=swddir,swddni=swddni,swddif=swddif, & ! jararias 2013/08/10 swdownc=swdownc, swddnic=swddnic, swddirc=swddirc, & ! PAJ obscur=obscur_loc, & - xcoszen=coszen,yr=yr,julian=julian,mp_physics=mp_physics ) ! jararias 2013/08/14 + xcoszen=coszen,yr=yr,julian=julian,mp_physics=mp_physics, & ! jararias 2013/08/14 + proceed_twoway_sw=proceed_twoway_sw, & ! WRF-CMAQ coupled model + nmode=3, & ! WRF-CMAQ coupled model + mass_ws_i=mass_ws_i, & ! WRF-CMAQ coupled model + mass_ws_j=mass_ws_j, & ! WRF-CMAQ coupled model + mass_ws_k=mass_ws_k, & ! WRF-CMAQ coupled model + mass_in_i=mass_in_i, & ! WRF-CMAQ coupled model + mass_in_j=mass_in_j, & ! WRF-CMAQ coupled model + mass_in_k=mass_in_k, & ! WRF-CMAQ coupled model + mass_ec_i=mass_ec_i, & ! WRF-CMAQ coupled model + mass_ec_j=mass_ec_j, & ! WRF-CMAQ coupled model + mass_ec_k=mass_ec_k, & ! WRF-CMAQ coupled model + mass_ss_i=mass_ss_i, & ! WRF-CMAQ coupled model + mass_ss_j=mass_ss_j, & ! WRF-CMAQ coupled model + mass_ss_k=mass_ss_k, & ! WRF-CMAQ coupled model + mass_h2o_i=mass_h2o_i, & ! WRF-CMAQ coupled model + mass_h2o_j=mass_h2o_j, & ! WRF-CMAQ coupled model + mass_h2o_k=mass_h2o_k, & ! WRF-CMAQ coupled model + dgn_i=dgn_i, & ! WRF-CMAQ coupled model + dgn_j=dgn_j, & ! WRF-CMAQ coupled model + dgn_k=dgn_k, & ! WRF-CMAQ coupled model + sig_i=sig_i, & ! WRF-CMAQ coupled model + sig_j=sig_j, & ! WRF-CMAQ coupled model + sig_k=sig_k, & ! WRF-CMAQ coupled model + gtauxar_01=sw_gtauxar_01, & ! WRF-CMAQ coupled model + gtauxar_02=sw_gtauxar_02, & ! WRF-CMAQ coupled model + gtauxar_03=sw_gtauxar_03, & ! WRF-CMAQ coupled model + gtauxar_04=sw_gtauxar_04, & ! WRF-CMAQ coupled model + gtauxar_05=sw_gtauxar_05, & ! WRF-CMAQ coupled model + asy_fac_01=sw_asy_fac_01, & ! WRF-CMAQ coupled model + asy_fac_02=sw_asy_fac_02, & ! WRF-CMAQ coupled model + asy_fac_03=sw_asy_fac_03, & ! WRF-CMAQ coupled model + asy_fac_04=sw_asy_fac_04, & ! WRF-CMAQ coupled model + asy_fac_05=sw_asy_fac_05, & ! WRF-CMAQ coupled model + ssa_01=sw_ssa_01, & ! WRF-CMAQ coupled model + ssa_02=sw_ssa_02, & ! WRF-CMAQ coupled model + ssa_03=sw_ssa_03, & ! WRF-CMAQ coupled model + ssa_04=sw_ssa_04, & ! WRF-CMAQ coupled model + ssa_05=sw_ssa_05, & ! WRF-CMAQ coupled model + sw_zbbcddir=sw_zbbcddir, & ! WRF-CMAQ coupled model + sw_dirdflux=sw_dirdflux, & ! WRF-CMAQ coupled model + sw_difdflux=sw_difdflux & ! WRF-CMAQ coupled model + ) + + ! = WRF-CMAQ twoway coupled model + if (proceed_twoway_sw) then + DO j=jts,jte + DO i=its,ite + sw_ttauxar_01(i,j) = sum(sw_gtauxar_01(i,:,j)) + sw_ttauxar_02(i,j) = sum(sw_gtauxar_02(i,:,j)) + sw_ttauxar_03(i,j) = sum(sw_gtauxar_03(i,:,j)) + sw_ttauxar_04(i,j) = sum(sw_gtauxar_04(i,:,j)) + sw_ttauxar_05(i,j) = sum(sw_gtauxar_05(i,:,j)) + END DO + END DO + end if DO j=jts,jte DO k=kts,kte diff --git a/phys/module_twoway_ra_rrtmg_sw.F b/phys/module_twoway_ra_rrtmg_sw.F new file mode 100644 index 0000000000..3a397e49a8 --- /dev/null +++ b/phys/module_twoway_ra_rrtmg_sw.F @@ -0,0 +1,327 @@ +MODULE module_twoway_ra_rrtmg_sw + +contains + +!------------------------------------------------------------------ + Subroutine get_aerosol_Optics_RRTMG_SW ( ns, nmode,delta_z, INMASS_ws, & + INMASS_in, INMASS_ec, INMASS_ss, & + INMASS_h2o, INDGN, INSIG, & + tauaer, waer, gaer ) + +!FSB This version switches between BHCOAT to BHMIE depending upon whether +! EC is present or not. 04/15/2012. + +!FSB this version does a core-shell calculation with BHCOAT 04/11/2012 +! This version is set up to be used with RRTMG_SW <<<<<<<< +! wavelenght is calculated internally +! FSB This routine calculates the aerosol information ( tauaer, waer, +! gaer, needed to calculate the solar radiation) The calling +! program specifies the location ( row, column, layer, +! layer thicknes, and wave length for the calculation. +! FSB 02/09/2011 Modifications made to subroutine ghintBH. +! FSB 04/14/2012 REmoved MODULUS, made changes to ghintBH. +! Put in option for core-shell (coated-sphere). 2 + +! FSB Input variables: + + use rrtmg_aero_optical_util_module + + implicit none + + integer,intent(in) :: ns ! index for wavelength should be + ! between 1 and 14. <<< RRTMG_SW + integer,intent(in) :: nmode ! should be 3 for WRF/CMAQ calculation + real,intent(in) :: delta_z ! layer thickness [m] +! FSB mode types for WRF/CMAQ +! nmode = 1 Aitken +! nmode = 2 accumulation +! nmode = 3 coarse +! FSB modal mass concentration by species [ ug / m**3] NOTE: MKS + real, intent(in) :: INMASS_ws(nmode) ! water soluble + real, intent(in) :: INMASS_in(nmode) ! insolugle + real, intent(in) :: INMASS_ec(nmode) ! elemental carbon or soot like + real, intent(in) :: INMASS_ss(nmode) ! sea salt + real, intent(in) :: INMASS_h2o(nmode) ! water +! FSB particle size-distribution information + real, intent(in) :: INDGN( nmode) ! geometric mean diameter [ m ] NOTE: MKS + real, intent(in) :: INSIG( nmode) ! geometric standard deviation + +!FSB output aerosol radiative properties [dimensionless] + real, intent(out) :: tauaer ! aerosol extinction optical depth + real, intent(out) :: waer ! aerosol single scattering albedo + real, intent(out) :: gaer ! aerosol assymetry parameter + +! FSB Internal variables + + real :: NR(nmode), NI(nmode) ! refractive indices + complex :: refcor(nmode), refshell(nmode) ! complex refracive indices + complex :: crefin(nmode) ! complex refractive index + +! FSB special values for EC CORE-shell calculation + real :: DGNSHELL(nmode) ! modal geometric mean diameter [m] + real :: DGNCORE (nmode) ! modal geometric mean diameter [m] + +! FSB Modal volumes [ m**3 / m**3 ] + real :: MVOL_ws(nmode) ! water soluble + real :: MVOL_in(nmode) ! insolugle + real :: MVOL_ec(nmode) ! soot like + real :: MVOL_ss(nmode) !sea salt + real :: MVOL_h2o(nmode) ! water +! real :: VOL(nmode) ! total modal volume [m** 3 / m**3] +! FSB special values for EC CORE-shell calculation + real :: VOLCOR(nmode) ! volume of EC core [m** 3 / m**3] + real :: VOLSHELL(nmode) ! volume of shell [m** 3 / m**3] + + integer :: m ! loop index + real :: bext ! extinction coefficient [1 / m] + real :: bscat ! scattering coefficient [1 / m] + real :: gfac ! asymmetry factor + + real :: bextsum, bscatsum, bsgsum + +! FSB History variables by wavelength and mode +! real :: bext_wm(ns,nmode) +! real :: bscat_wm(ns,nmode) +! real :: gfac_wm(ns,nmode) + + real, parameter :: one3rd = 1.0 / 3.0 + real :: dfac ! ratio of (volcor/vol) ** one3rd + ! used for calculating the diameter + ! of the EC core + + logical :: succesS + +!...component densities [ g/ cm**3 ] <<<<< cgs + + real, parameter :: rhows = 1.8 ! bulk density of water soluble aerosol + + real, parameter :: rhoin = 2.2 ! bulk density forinsoluble aerosol + +! real, parameter :: rhoec = 1.7 ! bulk density for soot aerosol + real, parameter :: rhoec = 1.8 ! new value + + real, parameter :: rhoh2o = 1.0 ! bulk density of aerosol water + + real, parameter :: rhoss = 2.2 ! bulk density of seasalt + +! FSB scale factor for volume calculation +! 1.0d-12 * [ cm**3 / g] -> [ m** 3 / ug ] + real, parameter :: scalefactor = 1.0e-12 + +! FSB scale factor for [1/g] to [1/ug] + real, parameter :: cug2g = 1.0e-06 + +! FSB reciprocal component densities[ m ** 3 / ug ] + + real, parameter :: rhows1 = scalefactor / rhows ! water soluble aerosol + + real, parameter :: rhoin1 = scalefactor / rhoin ! insoluble aerosol + + real, parameter :: rhoec1 = scalefactor / rhoec ! soot aerosol + + real, parameter :: rhoh2o1 = scalefactor / rhoh2o ! aerosol water + + real, parameter :: rhoss1 = scalefactor / rhoss ! seasalt + + integer,parameter :: nspint_sw = 14 ! number of spectral intervals for RRTMG_SW + +! FSB Band numbers and wavelengths for RRTMG_SW + integer, parameter :: Band(nspint_sw) = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 /) + + real, parameter :: LAMDA_SW(nspint_sw) = (/ 3.4615, 2.7885, 2.325, 2.046, 1.784, & + 1.4625, 1.2705, 1.0101, 0.7016, 0.53325, & + 0.38815, 0.299, 0.2316, 8.24 /) ! wavelength [ um ] + +! *** refractive indices + +! *** Except as otherwise noted reference values of refractive +! indices for aerosol particles are from the OPAC Data base. +! Hess, Koepke, and Schult, Optical properties of +! aerosols and clouds: The software package OPAC, Bulletan of +! the American Meteorological Society, Vol 79, No 5, +! pp 831 - 844, May 1998. +! OPAC is a downloadable data set of optical properties of +! 10 aerosol components, 6 water clouds and 3 cirrus clouds +! at UV, visible and IR wavelengths +! www.lrz-muenchen.de/~uh234an/www/radaer/opac.htm + + +! FSB water soluble + real, parameter :: xnreal_ws(nspint_sw) = (/ 1.443, 1.420, 1.420, 1.420, 1.463, 1.510, 1.510, & + 1.520, 1.530, 1.530, 1.530, 1.530, 1.530, 1.710 /) + real, parameter :: xnimag_ws(nspint_sw) = (/ 5.718E-3, 1.777E-2, 1.060E-2, 8.368E-3, 1.621E-2, & + 2.198E-2, 1.929E-2, 1.564E-2, 7.000E-3, 5.666E-3, & + 5.000E-3, 8.440E-3, 3.000E-2, 1.100E-1 /) + +! FSB sea salt + real, parameter :: xnreal_ss(nspint_sw) = (/ 1.480, 1.534, 1.437, 1.448, 1.450, 1.462, 1.469, & + 1.470, 1.490, 1.500, 1.502, 1.510, 1.510, 1.510 /) + + real, parameter :: xnimag_ss(nspint_sw) = (/ 1.758E-3, 7.462E-3, 2.950E-3, 1.276E-3, 7.944E-4, & + 5.382E-4, 3.754E-4, 1.498E-4, 2.050E-7, 1.184E-8, & + 9.938E-8, 2.060E-6, 5.000E-6, 1.000E-2 /) + +! FSB insoluble + real, parameter :: xnreal_in(nspint_sw) = (/ 1.272, 1.168, 1.208, 1.253, 1.329, 1.418, 1.456, & + 1.518, 1.530, 1.530, 1.530, 1.530, 1.530, 1.470 /) + real, parameter :: xnimag_in(nspint_sw) = (/ 1.165E-2, 1.073E-2, 8.650E-3, 8.092E-3, 8.000E-3, & + 8.000E-3, 8.000E-3, 8.000E-3, 8.000E-3, 8.000E-3, & + 8.000E-3, 8.440E-3, 3.000E-2,9.000E-2 /) + +! FSB 02/11/2012 These values are replaced. +! data xnreal_ec /1.877, 1.832, 1.813, 1.802, 1.791, 1.769, 1.761, & +! 1.760, 1.750, 1.740, 1.750, 1.738, 1.620, 2.120/ +! data xnimag_ec/ 5.563E-1, 5.273E-1, 5.030E-1, 4.918E-1, 4.814E-1, & +! 4.585E-1, 4.508E-1, 4.404E-1, 4.300E-1, 4.400E-1, & +! 4.600E-1, 4.696E-1, 4.500E-1, 5.700E-1/ + +! New Refractive indices for EC at RRTMG Wavelengths +! Source lamda xnreal_ec xnimag_ec +! C&C Values +! 3.4615 2.089 1.070 +! 2.7885 2.014 0.939 +! 2.325 1.962 0.843 +! 2.046 1.950 0.784 +! Bond values +! 1.784 1.940 0.760 +! 1.4625 1.930 0.749 +! 1.2705 1.905 0.737 +! 1.0101 1.870 0.726 +! B&B Values +! 0.7016 1.85 0.71 +! 0.53325 1.85 0.71 +! 0.38815 1.85 0.71 +! 0.299 1.85 0.71 +! 0.2316 1.85 0.71 +! C & C values +! 8.24 2.589 1.771 +!References: +! Bond Personal Communication from Tami Bond +! B&B Bond, T.C. & R.W. Bergstrom (2006) Light absorption by +! Carbonaceous Particles: An investigative review, +! Aerosol Science and Technology. Vol. 40. pp 27-67 +! +! C&C Chang,H and T.T. Charalmpopoulos (1990) Determination of the +! wavelength dependence of refractive indices of flame soot, +! Proceeding of the Royal Society of London A, Vol. 430, pp 577-591. +! FSB new values + +! FSB elemental carbon - soot like + + real, parameter :: xnreal_ec(nspint_sw) = (/ 2.089, 2.014, 1.962, 1.950, 1.940, 1.930, 1.905, & + 1.870, 1.85, 1.85, 1.85, 1.85, 1.85, 2.589 /) + real, parameter :: xnimag_ec(nspint_sw) = (/ 1.070, 0.939, 0.843, 0.784, 0.760, 0.749, 0.737, & + 0.726, 0.71, 0.71, 0.71, 0.71, 0.71, 1.771 /) + +! FSB water + real, save :: xnreal_h2o(nspint_sw) = (/ 1.408, 1.324, 1.277, 1.302, 1.312, 1.321, 1.323, & + 1.327, 1.331, 1.334, 1.340, 1.349, 1.362, 1.260 /) + real, save :: xnimag_h2o(nspint_sw) = (/ 1.420E-2, 1.577E-1, 1.516E-3, 1.159E-3, 2.360E-4, & + 1.713E-4, 2.425E-5, 3.125E-6, 3.405E-8, 1.639E-9, & + 2.955E-9, 1.635E-8, 3.350E-8, 6.220E-2 /) + + +! FSB Begin code ====================================================== + + bextsum = 0.0 + bscatsum = 0.0 + bsgsum = 0.0 + do m = 1, nmode +! FSB calculate volumes [ m**3 / m**3 ] +! FSB the reciprocal densities have been scaled to [ m**3 / ug ] + + MVOL_ws(m) = rhows1 * INMASS_ws(m) + MVOL_in(m) = rhoin1 * INMASS_in(m) + MVOL_ec(m) = rhoec1 * INMASS_ec(m) + MVOL_ss(m) = rhoss1 * INMASS_ss(m) + MVOL_h2o(m) = rhoh2o1 * INMASS_h2o(m) + + VOLSHELL(m) = MVOL_ws(m) + MVOL_in(m) + MVOL_ss(m) + MVOL_h2o(m) + VOLCOR(m) = MVOL_ec(m) +! VOL(m) = VOLSHELL(m) + VOLCOR(m) ! VOL is total volume + + if ( VOLCOR(m) .gt. 0.0 ) then +! FSB EC is present +! calculate the ratio of core to shell volume +! take cube root for scaling the diameter of +! the core to that of the shell. + +! dfac = ( VOLCOR(m) / VOL(m) ) ** one3rd + dfac = ( VOLCOR(m) / ( VOLSHELL(m) + VOLCOR(m) ) ) ** one3rd +! dfac = ( VOLCOR(m) / ( VOLSHELL(m) + VOLCOR(m) ) ) +! FSB Set shell and core diameters + DGNSHELL(m) = INDGN(m) + DGNCORE(M) = dfac * INDGN(m) +! FSB note that VOLSHELL(m) is the total volume when EC is not present + end if + +! internal mixture of non-EC species. + +! modal real refractive index No EC + nr(m) = (MVOL_ws(m) * xnreal_ws(ns) + & + MVOL_in(m) * xnreal_in(ns) + & + MVOL_ss(m) * xnreal_ss(ns) + & +! MVOL_h2o(m) * xnreal_h2o(ns)) / VOL(m) + MVOL_h2o(m) * xnreal_h2o(ns)) / VOLSHELL(m) + +! modal imaginary refractive index no EC + ni(m) = (MVOL_ws(m) * xnimag_ws(ns) + & + MVOL_in(m) * xnimag_in(ns) + & + MVOL_ss(m) * xnimag_ss(ns) + & +! MVOL_h2o(m) * xnimag_h2o(ns)) / VOL(m) + MVOL_h2o(m) * xnimag_h2o(ns)) / VOLSHELL(m) + + if ( VOLCOR(m) .gt. 0.0) then + +! FSB calculate the complex refractive indices for the CORE and +! the SHELL for case when and EC core is assumed to exist. + + refcor(m) = cmplx( xnreal_ec(ns), xnimag_ec(ns) ) + refshell(m) = cmplx(nr(m), ni(m) ) +! FSB do BHCOAT case + CALL aero_optical_CS( LAMDA_SW(ns), refcor(m), refshell(m), & + VOLCOR(m),VOLSHELL(m), DGNCORE(m), & + DGNSHELL(m), INSIG(m), & + bext, bscat, gfac, succesS ) +! else if ( VOL(m) .gt. 0.0) then + else if ( VOLSHELL(m) .gt. 0.0) then +! FSB do BHMIE case for the case when EC is not present. + crefin(m) = cmplx(nr(m), ni(m) ) +! CALL aero_optical2( LAMDA_SW(ns), crefin(m), VOL(m), & + CALL aero_optical2( LAMDA_SW(ns), crefin(m), VOLSHELL(m), & + INDGN(m), INSIG(m), & + bext, bscat, gfac, success ) + else + bext = 0.0 + bscat = 0.0 + gfac = 0.0 + end if + +! FSB sum for total values + bextsum = bextsum + bext + bscatsum = bscatsum +bscat + bsgsum = bsgsum + bscat * gfac +! FSB get history +! bext_wm(ns,m) = bext +! bscat_wm(ns,m) = bscat +! gfac_wm(ns,m) = gfac + end do ! loop over modes + +! FSB construct output variables + tauaer = bextsum * delta_z + waer = bscatsum / bextsum + gaer = bsgsum / bscatsum + +! Write out modal values. + +! write(100,*) ' lamda mode bext bscat gaer ' +! m=1 +! write(100,*) lamda(ns), m, bext_wm(ns,m),bscat_wm(ns,m), gfac_wm(ns,m) +! m=2 +! write(100,*) lamda(ns), m, bext_wm(ns,m),bscat_wm(ns,m), gfac_wm(ns,m) +! m=3 +! write(100,*) lamda(ns), m, bext_wm(ns,m),bscat_wm(ns,m), gfac_wm(ns,m) + + end subroutine get_aerosol_Optics_RRTMG_SW + +END MODULE module_twoway_ra_rrtmg_sw diff --git a/phys/module_twoway_rrtmg_aero_optical_util.F b/phys/module_twoway_rrtmg_aero_optical_util.F new file mode 100644 index 0000000000..e074d34fa5 --- /dev/null +++ b/phys/module_twoway_rrtmg_aero_optical_util.F @@ -0,0 +1,2768 @@ +! Revision History: +! 2016/02/23 David Wong extracted the complex number module and put it in a file +! 2016/05/23 David Wong - replaced rrtmg_aero_optical_util_module with +! cmaq_rrtmg_aero_optical_util_module to avoid duplication +! of the same module name on WRF side of the two-way model + +MODULE rrtmg_aero_optical_util_module + + Integer :: AERO_UTIL_LOG = 0 + + private + public :: aero_optical, aero_optical2, aero_optical_CS, AERO_UTIL_LOG + + interface ghintBH + module procedure ghintBH_1, ghintBH_2, ghintBH_Odd + end interface + + interface ghintBH_CS + module procedure ghintBH_CS_even, ghintBH_CS_odd + end interface + + Logical, Parameter :: Use_Odd_Quadrature = .True. + Integer, Parameter :: Quadrature_Points = 3 + +!B.Hutzell One point quadature IGH = 1 + + real, parameter :: ghxi_1(1) = 0.00000000000 + real, parameter :: ghwi_1(1) = 1.77245385091 + +!B.Hutzell Three point quadature IGH = 3 + real, parameter :: ghxi_3(3) = (/ -1.22474487139, & + 0.00000000000, & + 1.22474487139 /) + + real, parameter :: ghwi_3(3) = (/ 0.295408975151, & + 1.181635900000, & + 0.295408975151 /) + +!B.Hutzell Five point quadature IGH = 5 + real(8), parameter :: ghxi_5(5) = (/ -2.02018287046d0, & + -0.958572464614d0, & + 0.00000000000d0, & + 0.958572464614d0, & + 2.02018287046d0 /) + + real(8), parameter :: ghwi_5(5) = (/ 0.019953242059d0, & + 0.393619323152d0, & + 0.945308720483d0, & + 0.393619323152d0, & + 0.019953242059d0 /) + +!B.Hutzell Nine point quadature IGH = 9 points +!No. Abscissas Weight Total Weight + real, parameter :: ghxi_9(9) = (/ -3.19099320178, & + -2.26658058453, & + -1.46855328922, & + -0.72355101875, & + 0.00000000000, & + 0.72355101875, & + 1.46855328922, & + 2.26658058453, & + 3.19099320178 /) + + real, parameter :: ghwi_9(9) = (/ 3.96069772633E-5, & + 0.00494362428, & + 0.08847452739, & + 0.43265155900, & + 0.72023521561, & + 0.43265155900, & + 0.08847452739, & + 0.004943624275, & + 3.96069772633E-5 /) + + contains + +! ------------------------------------------------------------------ + subroutine getqext_BH (xx, crefin, qextalf, qscatalf, gscatalfg,SUCCESS) + + implicit none + + real, intent(in) :: XX + real, intent(out) :: qextalf, qscatalf, gscatalfg + complex, intent(in) :: CREFIN + logical, intent(out) :: success + +! local variables + real( 8 ), parameter :: one_third = 1.0d0 / 3.0d0 + + integer :: NXX + integer :: nstop, modulus + + real :: QEXT, QSCA, QBACK, G_MIE, xx1 + + real( 8 ) :: x + complex( 8 ) :: refractive_index + + x = real( XX, 8 ) + refractive_index = dcmplx( real( CREFIN ), imag( CREFIN ) ) + + modulus = int( abs( x * refractive_index ) ) + nstop = int( x + 4.0d0 * x**one_third + 2.0d0 ) + + nxx = max( modulus, nstop ) + 15 + + xx1 = 1.0 / XX + + CALL BHMIE_FLEXI (XX, NXX, NSTOP, CREFIN,QEXT,QSCA,QBACK,G_MIE, SUCCESS) + + qextalf = QEXT * xx1 + qscatalf = QSCA * xx1 + gscatalfg = qscatalf * G_MIE + + end subroutine getqext_bh + +! ------------------------------------------------------------------ + SUBROUTINE BHMIE (X, REFREL, QQEXT, QQSCA, QBACK, GSCA, SUCCESS) + +! FSB Changed the call vector to return only QEXT, QSCAT QBACK GSCA +! and ignore NANG, S1 and S2 and all calculations for them + + implicit none + +! Arguments: + real, intent(in) :: X ! X = pi*particle_diameter / Wavelength + complex, intent(in) :: REFREL + +! REFREL = (complex refr. index of sphere)/(real index of medium) +! in the current use the index of refraction of the the medium +! i taken at 1.0 real. +! +! Output + + real, intent(out) :: QQEXT, QQSCA, QBACK, GSCA + logical, intent(out) :: SUCCESS + +! QQEXT Efficiency factor for extinction +! QQSCA Efficiency factor for scattering +! QQBACK Efficiency factor for back scatter +! GSCA asymmetry factor +! SUCCESS flag for successful calculation +! REFERENCE: +! Bohren, Craig F. and Donald R. Huffman, Absorption and +! Scattering of Light by Small Particles, Wiley-Interscience +! copyright 1983. Paperback Published 1998. +! FSB +! This code was originally listed in Appendix A. pp 477-482. +! As noted below, the original code was subsequently +! modified by Prof. Bruce T. Drain of Princetion University. +! The code was further modified for a specific application +! in a large three-dimensional code requiring as much +! computational efficiency as possible. +! Prof. Francis S. Binkowski of The University of North +! Carolina at Chapel Hill. + +! Declare parameters: +! Note: important that MXNANG be consistent with dimension of S1 and S2 +! in calling routine! + + integer, parameter :: MXNANG=10, NMXX=600000 ! FSB new limits + real*8, parameter :: PII = 3.1415916536D0 + real*8, parameter :: ONE = 1.0D0, TWO = 2.0D0 + +! Local variables: + integer :: NANG + integer :: N,NSTOP,NMX,NN + real*8 :: QSCA, QEXT, DX1, DXX1 + real*8 :: CHI,CHI0,CHI1,DX,EN,P,PSI,PSI0,PSI1,XSTOP,YMOD + real*8 :: TWO_N_M_ONE, TWO_N_P_ONE, EN1, FACTOR + complex*16 :: AN,AN1,BN,BN1,DREFRL,XI,XI1,Y, Y1, DREFRL1 + complex*16 :: D(NMXX), FAC1, FAC2 + complex*16 :: XBACK + +!*********************************************************************** +! Subroutine BHMIE is the Bohren-Huffman Mie scattering subroutine +! to calculate scattering and absorption by a homogenous isotropic +! sphere. +! Given: +! X = 2*pi*a/lambda +! REFREL = (complex refr. index of sphere)/(real index of medium) +! real refractive index of medium taken as 1.0 +! Returns: +! QEXT = efficiency factor for extinction +! QSCA = efficiency factor for scattering +! QBACK = efficiency factor for backscatter +! see Bohren & Huffman 1983 p. 122 +! GSCA = asymmetry for scattering +! +! Original program taken from Bohren and Huffman (1983), Appendix A +! Modified by Prof. Bruce T.Draine, Princeton Univ. Obs., 90/10/26 +! in order to compute +! 91/05/07 (BTD): Modified to allow NANG=1 +! 91/08/15 (BTD): Corrected error (failure to initialize P) +! 91/08/15 (BTD): Modified to enhance vectorizability. +! 91/08/15 (BTD): Modified to make NANG=2 if called with NANG=1 +! 91/08/15 (BTD): Changed definition of QBACK. +! 92/01/08 (BTD): Converted to full double precision and double complex +! eliminated 2 unneed lines of code +! eliminated redundant variables (e.g. APSI,APSI0) +! renamed RN -> EN = double precision N +! Note that DOUBLE COMPLEX and DCMPLX are not part +! of f77 standard, so this version may not be fully +! portable. In event that portable version is +! needed, use src/bhmie_f77.f +! 93/06/01 (BTD): Changed AMAX1 to generic function MAX +! FSB April 09,2012 This code was modified by: +! Prof. Francis S. Binkowski University of North Carolina at +! Chapel Hill, Institue for the Environment. +! +! The modifications were made to enhance computation speed +! for use in a three-dimensional code. This was done by +! removing code that calculated angular scattering. The method +! of calculating QEXT, QBACK was also changed. + +!*********************************************************************** +!*** Safety checks + + SUCCESS = .TRUE. + NANG = 2 ! FSB only this value +! IF(NANG.GT.MXNANG)STOP'***Error: NANG > MXNANG in bhmie' +! IF (NANG .LT. 2) NANG = 2 + + DX = REAL( X, 8 ) +! FSB Define reciprocals so that divisions can be replaced by multiplications. + DX1 = ONE / DX + DXX1 = DX1 * DX1 + DREFRL = DCMPLX( REFREL ) + DREFRL1 = ONE / DREFRL + Y = DX * DREFRL + Y1 = ONE / Y + YMOD = ABS(Y) + +!*** Series expansion terminated after NSTOP terms +! Logarithmic derivatives calculated from NMX on down + XSTOP = REAL( X + 4.0 * X**0.3333 + 2.0, 8) + NMX = INT( MAX(XSTOP,YMOD) ) + 15 + +! BTD experiment 91/1/15: add one more term to series and compare results +! NMX=AMAX1(XSTOP,YMOD)+16 +! test: compute 7001 wavelengths between .0001 and 1000 micron +! for a=1.0micron SiC grain. When NMX increased by 1, only a single +! computed number changed (out of 4*7001) and it only changed by 1/8387 +! conclusion: we are indeed retaining enough terms in series! + NSTOP = INT( XSTOP ) + FACTOR = 1.0D0 + + IF (NMX .GT. NMXX) THEN + WRITE(6,*)'Error: NMX > NMXX=',NMXX,' for |m|x=',YMOD + SUCCESS = .FALSE. + RETURN + END IF + +! FSB all code relating to scattering angles is removed out for +! reasons of efficiency when running in a three-dimensional +! code. We only need QQSCA, QQEXT, GSCA AND QBACK + + +!*** Logarithmic derivative D(J) calculated by downward recurrence +! beginning with initial value (0.,0.) + + D(NMX) = DCMPLX(0.0D0,0.0D0) + NN = NMX - 1 + DO N = 1,NN + EN = REAL(NMX - N + 1, 8 ) +! FSB In the following division by Y has been replaced by +! multiplication by Y1, the reciprocal of Y. + D(NMX-N) = ( EN * Y1 ) - (ONE / ( D(NMX-N+1) + EN * Y1)) + END DO + +!*** Riccati-Bessel functions with real argument X +! calculated by upward recurrence + + PSI0 = COS(DX) + PSI1 = SIN(DX) + CHI0 = -SIN(DX) + CHI1 = PSI0 + XI1 = DCMPLX(PSI1,-CHI1) + QSCA = 0.0D0 + GSCA = 0.0D0 + QEXT = 0.0D0 + P = -ONE + XBACK = (0.0d0,0.0d0) + +! FSB Start main loop + DO N = 1,NSTOP + EN = REAL( N, 8) + EN1 = ONE / EN + TWO_N_M_ONE = TWO * EN - ONE +! for given N, PSI = psi_n CHI = chi_n +! PSI1 = psi_{n-1} CHI1 = chi_{n-1} +! PSI0 = psi_{n-2} CHI0 = chi_{n-2} +! Calculate psi_n and chi_n + PSI = TWO_N_M_ONE * PSI1 * DX1 - PSI0 + CHI = TWO_N_M_ONE * CHI1 * DX1 - CHI0 + XI = DCMPLX(PSI,-CHI) + +!*** Compute AN and BN: +! FSB Rearrange to get common terms + FAC1 = D(N) * DREFRL1 + EN * DX1 + AN = (FAC1) * PSI - PSI1 + AN = AN / ( (FAC1 )* XI - XI1 ) + FAC2 = ( DREFRL * D(N) + EN * DX1) + BN = ( FAC2) * PSI -PSI1 + BN = BN / ((FAC2) * XI - XI1 ) + +! FSB calculate sum for QEXT as done by Wiscombe +! get common factor + TWO_N_P_ONE = (TWO * EN + ONE) + QEXT = QEXT + (TWO_N_P_ONE) * (REAL(AN) + REAL(BN) ) + QSCA = QSCA + (TWO_N_P_ONE) * ( ABS(AN)**2+ ABS(BN)**2 ) + +! FSB calculate XBACK from B & H Page 122 + FACTOR = -1.0d0 * FACTOR ! calculate (-1.0 ** N) + XBACK = XBACK + (TWO_N_P_ONE) * factor * (AN - BN) + +! FSB calculate asymmetry factor + GSCA = GSCA + REAL( ((TWO_N_P_ONE)/(EN * (EN + ONE))) * & + (REAL(AN)*REAL(BN)+IMAG(AN)*IMAG(BN))) + + IF (N .GT. 1)THEN + GSCA = GSCA + REAL( (EN - EN1) * & + (REAL(AN1)*REAL(AN) + IMAG(AN1)*IMAG(AN) + & + REAL(BN1)*REAL(BN) + IMAG(BN1)*IMAG(BN))) + ENDIF + +!*** Store previous values of AN and BN for use in computation of g= + AN1 = AN + BN1 = BN + +! FSB set up for next iteration + PSI0 = PSI1 + PSI1 = PSI + CHI0 = CHI1 + CHI1 = CHI + XI1 = DCMPLX(PSI1,-CHI1) + + END DO ! main loop on n + +!*** Have summed sufficient terms. + +! Now compute QQSCA,QQEXT,QBACK,and GSCA + GSCA = REAL( TWO / QSCA ) * GSCA + +! FSB in the following, divisions by DX * DX has been replaced by +! multiplication by DXX1 the reciprocal of 1.0 / (DX *DX) + QQSCA = REAL( TWO * QSCA * DXX1 ) + QQEXT = REAL( TWO * QEXT * DXX1 ) + QBACK = REAL( REAL ( 0.5d0 * XBACK * CONJG(XBACK), 8 ) * DXX1 ) ! B&H Page 122 + + END subroutine BHMIE + +! ------------------------------------------------------------------ + subroutine aero_optical ( lamda_in, nmode, nr, ni, Vol, & + dgn, sig, bext, bscat, g_bar, & + success, modulus ) + +! *** calculate the extinction and scattering coefficients and +! assymetry factors for each wavelength as a sum over the +! individual lognormal modes. Each mode may have a different +! set of refractive indices. + + IMPLICIT NONE +! *** input variables + real, intent(in) :: lamda_in ! wavelengths [micro-m] + INTEGER, intent(in) :: nmode ! number of lognormal modes + real, intent(in) :: nr( nmode), ni(nmode) ! real and imaginary + ! refractive indices + real, intent(in) :: Vol(nmode) ! modal aerosol volumes [m**3 /m**3] + real, intent(in) :: dgn(nmode) ! geometric mean diameters + ! for number distribution [ m] + real, intent(in) :: sig(nmode) ! geometric standard deviation + + real, intent(in), optional :: modulus(nmode) ! modulus of refracive index + +! *** output variables + real, intent(out) :: bext ! extinction coefficient [ 1 / m ] + real, intent(out) :: bscat ! scattering coefficient [ 1 / m ] + real, intent(out) :: g_bar ! assymetry factor for Mie and molecular scattering + logical, intent(out) :: success ! flag for successful calculation +! *** internal variables + INTEGER :: j ! loop index +! real :: xlnsig(nmode) ! natural log of geometric standard deviations + real :: beta_Sc, bsc !aerosol scattering coefficient + + real :: beta_Ex ! aerosol extinction coefficients + real :: G ! modal aerosol assymetry factors + real :: sum_g + real :: LSIGX + real :: lamdam1 ! 1/ lamda + real :: alphav ! Mie size parameter + real :: vfac + real :: modalph + + real, parameter :: pi = 3.14159265359 + + Logical, Save :: Initialize = .True. + +! *** coded 09/08/2004 by Dr. Francis S. Binkowski +! FSB Modified for RRTMG version December 2009. +! FSB modified 10/06/2004, 10/12/2004, 10/18/2005 +! FSB 01/12/2006 +! Formerly Carolina Environmental Program +! FSB now the Institute for the Environment +! University of North Carolina at Chapel Hill +! email: frank_binkowski@unc.edu + + +! *** initialize variables + lamdam1 = 1.0e6 / lamda_in ! lamda now in [ m ] + bext = 0.0 + bscat = 0.0 + sum_g = 0.0 + + DO j = 1, nmode +! calculate the extinction and scattering coefficients +! for each mode + LSIGX = log(sig(j)) + +! calculate Mie size parameter for volume distribution +! exp(3.0 * xlnsig*xlnsig) converts dgn to dgv (volume diameter) + alphav = pi * dgn(j) * exp(3.0 * LSIGX * LSIGX) * lamdam1 + + if (present(modulus)) then + modalph = alphav * modulus(j) + end if + + CALL ghintBH (nr(j), ni(j), alphav, LSIGX, beta_EX, beta_Sc, G, success) + +! *** ghintBH returns the normalized values +! Calculate the actual extinction and scattering coefficients +! by multplying by the modal volume and dividing by the wavelength + + vfac = Vol(j) * lamdam1 + +! *** sum to get total extinction and scattering +! and contribution to the overal assymetry factor + + bext = bext + vfac * beta_Ex ! [ 1 / m ] + bsc = vfac * beta_Sc + bscat = bscat + bsc + sum_g = sum_g + bsc * G + + END DO ! loop on modes + +! *** calculate combined assymetry factor for all modes + + g_bar = sum_g / bscat ! changed to divide by bscat + + END SUBROUTINE aero_optical + +! ------------------------------------------------------------------ + subroutine ghintBH_1 (nr, ni, alfv, xlnsig, Qext_GH, Qscat_GH, g_gh, success) + +! FSB *********** This is the newest (05_30_2012) version of GhintBH +! this version does the Mie method and calculates the optimum set of +! set of Gauss-Hermite abscissas and weights. +! FSB Calls Penndorf codes for alfv .le. 0.3 + +! Dr. Francis S. Binkowski, The University of North Carolina +! at Chapel Hill +! FSB this code file now contains all of the necessary subroutines that +! are called to perform an integral of the Bohren and Huffman +! Mie codes ( as updated by Prof. Bruce C. Drain of Princeton) +! calculates the extinction and scattering coefficients +! normalized by wavelength and total particle volume +! concentration for a log normal particle distribution +! with the logarithm of the geometric standard deviation +! given by xlnsig. The integral of the +! asymmetry factor g is also calculated. +! FSB Change 12/20/2011 This code now has a choice of IGH based +! upon alfv and nr. +! *** Does Gauss-Hermite quadrature of Qext / alfa & Qscat / alfa +! and asymmetry factor over log normal distribution using +! symmetric points. + + implicit none + + real, intent(in) :: nr, ni ! refractive indices + real, intent(in) :: alfv ! Mie parameter for dgv + real, intent(in) :: xlnsig ! log of geometric standard deviation + real, intent(out) :: Qext_GH ! normalized extinction efficiency + real, intent(out) :: Qscat_GH ! normalized scattering efficiency + real, intent(out) :: g_GH ! asymmetry factor + logical, intent(out) :: success ! flag for successful calculation + + real :: bext_P, bscat_P, babs_P, g_PCS, xlnsg2 ! see below for definition + + real :: aa1 ! see below for definition + real :: alfaip, alfaim ! Mie parameters at abscissas + +! *** these are Qext/alfa and Qscat/alfv at the abscissas + real :: qalfip_e, qalfim_e ! extinction + real :: qalfip_s, qalfim_s ! scattering + real :: gsalfp, gsalfm ! scattering times asymmetry factor + integer :: IGH ! index for GH quadrature + +! FSB define parameters + real, parameter :: pi = 3.14159265 + real, parameter :: sqrtpi = 1.772454 + real, parameter :: sqrtpi1 = 1.0 / sqrtpi + real, parameter :: sqrt2 = 1.414214 + real, parameter :: three_pi_two = 3.0 * pi / 2.0 + real, parameter :: const = three_pi_two * sqrtpi1 + + integer :: i + complex :: crefin ! complex index of refraction + real :: sum_e,sum_s, xi,wxi,xf + real :: sum_sg + +! Gauss-Hermite abscissas and weights +! *** the following weights and abscissas are from Abramowitz +! Stegun, Table 25.10 page 924 +! FSB full precision from Table 25.10 + +! FSB ten-point - IGH = 5 + real, parameter :: ghxi_10(5) = (/ 0.342901327223705, & + 1.036610829789514, & + 1.756683649299882, & + 2.532731674232790, & + 3.436159118837738 /) + + real, parameter :: ghwi_10(5) = (/ 6.108626337353e-01, & + 2.401386110823e-01, & + 3.387439445548e-02, & + 1.343645746781e-03, & + 7.640432855233e-06 /) + +! FSB six-point - IGH = 3 + real, parameter :: ghxi_6(3) = (/ 0.436077411927617, & + 1.335849074013597, & + 2.350604973674492 /) + + real, parameter :: ghwi_6(3) = (/ 7.246295952244e-01, & + 1.570673203229e-01, & + 4.530009905509e-03 /) + +! FSB two-point - IGH = 1 + real, parameter :: ghxi_2(1) = (/ 0.707106781186548 /) + + real, parameter :: ghwi_2(1) = (/ 8.862269254528e-01 /) + + real :: GHXI(5), GHWI(5) ! weight and abscissas + integer :: NMAX ! number of weights and abscissa + +! FSB Check for valid range of Penndorf application. + if ( alfv .le. 0.3) then + xlnsg2 = xlnsig*xlnsig + call pennfsb (nr,ni,alfv,xlnsg2,bext_P,bscat_P,babs_P,g_PCS) + Qext_GH = bext_P + Qscat_GH = bscat_p + g_GH = g_PCS * exp(4.0 * xlnsg2) ! match GH integral + else + +! FSB We need to do a full Mie calculation now +! Choose IGH. These choices are designed to improve +! the computational efficiency without sacrificing accuracy. + + IGH=3 ! default value; six_point is sufficient generally +! six point + NMAX = 3 + + if (nr .ge. 1.7) then +! 10 point + IGH = 5 ! more points needed here + NMAX = 5 + end if + + if ( alfv .gt. 20.0 .or. alfv .lt. 0.5 ) then + IGH = 1 ! in this range fewer points are needed + NMAX = 1 + end if + + if (IGH == 1) then + GHXI(1) = ghxi_2(1) + GHWI(1) = ghwi_2(1) + else if (IGH == 3) then + do i = 1, NMAX + GHXI(i) = ghxi_6(i) + GHWI(i) = ghwi_6(i) + end do + else + do i = 1,NMAX + GHXI(i) = ghxi_10(i) + GHWI(i) = ghwi_10(i) + end do + end if ! set up number of abscissas and weights + +! FSB set complex refractive index. + crefin= cmplx(nr,ni) + +! FSB now start the integration code + aa1 = sqrt2 * xlnsig ! This 1.0 / Sqrt( A ) in derivation of the integral + ! where A = 1.0 / ( 2.0 * xlnsg**2 ) + +! Then alpha = alfv * exp[ u / sqrt(A) ] +! For Gauss-Hermite Quadrature u = xi +! Therefore, xf = exp( xi / sqrt(A) ), +! or xf = exp( xi * aa1 ) + sum_e = 0.0 + sum_s = 0.0 + sum_sg = 0.0 +! FSB do NMAX calls to the MIE codes + do i = 1,NMAX + xi = GHXI(i) + wxi = GHWI(i) + xf = exp( xi * aa1 ) + alfaip = alfv * xf + alfaim = alfv / xf ! division cheaper than another exp() +! *** call subroutine to fetch the effficiencies + + call getqext_BH (alfaip, crefin, qalfip_e, qalfip_s, gsalfp, success) + call getqext_BH (alfaim, crefin, qalfim_e, qalfim_s, gsalfm, success) + + sum_e = sum_e + wxi * ( qalfip_e + qalfim_e ) + sum_s = sum_s + wxi * ( qalfip_s + qalfim_s ) + sum_sg = sum_sg + wxi * ( gsalfp + gsalfm ) + end do + + g_GH = sum_sg / sum_s ! this is + Qext_GH = const * sum_e ! + Qscat_GH = const * sum_s + end if + + end subroutine ghintBH_1 + +! ------------------------------------------------------------------ + subroutine pennfsb (n, k, xx, lnsg2, bext, bscat, babs, g) + +! FSB a new version of Penndorf's equations. This version does +! analytical integration for Qext, Qscat, Qabs to generate +! bext, bscat, babs. Note that the expressions for Qext & Qscat +! hve been divide through by xx. +! +! Reference: +! Caldas, M., V. Semiao, 2001, Radiative properties of small +! particles: and extension of the Penndorff Model. Journal +! of the Optical Society of America A, Vol. 18, No. 4, +! pp 831-838. + +! Penndorf, R., 1962a,Scattering and extinction coefficients for small +! absorbing and nonabsorbing aerosols, +! J. Optical Society of America, 52, 896-904. + +! Penndorf, P., 1962b,Scattering and extinction coefficients for +! small Spherical aerosols, J. Atmos. Sci., 19, p 193 + +! FSB Coded by Dr. Francis S. Binkowski on October 25, 2011 by combining +! two previous versions to get a common code for the Penndorf and +! and Caldas & Semiao approaches. The Penndorf Qext, Qscat are much +! better than the versions from Caldas & Semiao despite claims to +! the contrary. The values of the asymmetry factor from Caldas & Semiao +! are better than can be obtained from Penndorf. + +! FSB This version does the analytical integral ove a lognormal +! size distribution. + + implicit none +! input variables + real, intent(in) :: n, k ! refractive index + real, intent(in) :: xx ! pi * diameter / wavelength + real, intent(in) :: lnsg2 ! log(sigma_g)**2 + real, intent(out) :: bext ! extinction coefficient + real, intent(out) :: bscat ! scattering coefficient + real, intent(out) :: babs ! absorption coefficient + real, intent(out) :: g ! asmmetry factor + +! internal variables + complex*16 :: m, m2,m4,m6,m21,m22 + complex*16 :: P,Q,R,S,T,U,V,W + complex*16 :: Qprime, Rprime,Sprime,Tprime + complex*16 :: Uprime, Vprime, Wprime + real*8 :: Qs, gQs, gpennCS + real*8 :: P1,P2, Q1, Q2 , S2,V1, V2 ! see usage + real*8 :: P1SQ, P2SQ ! see usage + real*8 :: y, y2, y3, y4, y6, y7, y8, y9 + real*8 :: x, x2, x3, x4, x6, x7, x8, x9 + real :: mag, modalf +! FSB define useful numbers and fractions + real, parameter :: pi = 3.14159265358979324d0 + real, parameter :: three_pi_two = 1.5d0 * pi + + real*8, parameter :: one = 1.0d0 + real*8, parameter :: two = 2.0d0 + real*8, parameter :: three = 3.0d0 + real*8, parameter :: four = 4.0d0 + real*8, parameter :: five = 5.0d0 + real*8, parameter :: six = 6.0d0 + real*8, parameter :: eight = 8.0d0 + real*8, parameter :: nine = 9.0d0 + real*8, parameter :: fifteen = 15.0d0 + real*8, parameter :: fortyfive = 45.0d0 +! real*8, parameter :: two5ths = two / five + real*8, parameter :: twothrds = two / three + real*8, parameter :: fourthirds = four / three + real*8, parameter :: onefifteenth = one / fifteen + real*8, parameter :: twofifteenths = two * onefifteenth +! real*8, parameter :: fourninths = four / nine + real*8, parameter :: eightthirds = two * fourthirds + real*8, parameter :: one_big = one / 31500.0d0 + real*8, parameter :: two_fortyfive = two / fortyfive + real*8, parameter :: four_225 = four / 225.0d0 + real*8, parameter :: one_210 = one / 210.0d0 +! real*8, parameter :: one_half = one / two +! real*8, parameter :: four_two = two + real*8, parameter :: nine_two = 4.5d0 +! real*8, parameter :: sixteen_two = eight +! real*8, parameter :: thirtysix_two = 36.0 / two +! real*8, parameter :: twentyfive_two = 25.0d0 / two +! real*8, parameter :: sixtyfour_two = 64.0d0 / two +! real*8, parameter :: fortynine_two = 49.0d0 / two +! real*8, parameter :: eightyone_two = 81.0d0 / two + real*8 :: A,B,C,D,E, AA,BB,CC + +! FSB start code + mag = sqrt( n * n + k * k ) + modalf = mag * xx + y = REAL( xx, 8 ) ! convert to real*8 +! FSB get powers of y + y2 = y * y + y3 = y2 * y + y4 = y3 * y + y6 = y3 * y3 + y7 = y3 * y4 + y8 = y4 * y4 + y9 = y6 * y3 + +! FSB Calculate integrals ove the lognormal distribution +! this is done term by term and the form is +! xn = yn * exp( (n**2) * lnsig2 /2.0d0) + + x = y + x2 = y2 * exp( two * lnsg2) + x3 = y3 * exp( nine_two * lnsg2) + x4 = y4 ! * exp( eight * lnsg2) + x6 = y6 ! * exp( thirtysix_two * lnsg2) + x7 = y7 ! * exp( fortynine_two * lnsg2) + x8 = y8 ! * exp( fortynine_two * lnsg2) + x9 = y9 ! * exp( eightyone_two * lnsg2) + + +! FSB explicitly calculate complex refrative index m + m = dcmplx(n,-k) +! FSB get powers and functions of m + m2 = m * m + m4 = m2 * m2 + m6 = m2 * m4 + m21 = m2 - one + m22 = m2 + two + +! FSB calculate Penndorf's definitions from Table II of Penndorf (1962a) + P = m21 / m22 + Q = (m2 - two ) / m22 + S = m21 / ( two * m2 + three) + V = m21 +! FSB get real & imaginary parts following Penndorf's mptation + P1 = real(P) + P2 = -aimag(P) + P1SQ = P1 * P1 + P2SQ = P2 * P2 + + Q1 = real(Q) + Q2 = -aimag(Q) + S2 = -aimag(S) + V1 = real(V) + v2 = -aimag(V) + + +! FSB Get bext from Penndorf (1962a) Equation (7) up to x4 +! consistent with equation (8) +! We have then divided through by x and integrated analytically + bext = REAL( four * P2 + ( 2.4d0 * (P1 * Q2 + P2 * Q1 ) + twothrds * S2 & + + twofifteenths * V2 ) * x2 + ( eightthirds * ( P1SQ - P2SQ ) ) * x3, 4 ) + +! FSB get bscat from Penndorf Equation (9) up to x4 +! we have divided through by x and integrated analytically + bscat = REAL( eightthirds * ( P1SQ + P2SQ ) * x3 ) +! FSB calculate babs +! babs = bext - bscat + +! FSB now get asymmetry factor from Caldas & Semiao (2001) +! +! *** The following additional variables from Caldas & Semiao (2001) +! are defined in Equations 10a to 10h. + + R = (m6 + 20.0d0*m4 -200.0d0*m2 + 200.0d0) / m22**2 + T = m21 / ( ( 2.0d0 * M2 + 3.0d0) **2 ) + U = m21 / (3.0d0 * M2 + 4.0d0 ) + W = m21 * ( 2.0d0 * m2 - 5.0d0) + +! *** further definitions from Caldas & Semiao (2001) + Qprime = Q + Rprime = 18.0d0 * R + Sprime = 5.0d0 * S / P + Tprime = 375.0d0 * T / P +! Uprime = 28.0d0 * U / P + Vprime = V / P + Wprime = 5.0d0 * W / P + +! FSB calculate gQs and Qs from Caldas & Semiao (2001) +! *** calculate Qs equation 13 +! Qs = eightthirds * abs(P)**2 & +! * (x4 + onefifteenth * real(Qprime) * x6 & +! + fourthirds * aimag(P) * x7 & +! + one_big * ( 35.0d0 * abs(Qprime)**2 & +! + 20.0d0 * real(Rprime) + 35.0d0 * abs(Vprime)**2 & +! + 21.0d0 * abs(Sprime)**2 ) * x8 & +! + two_fortyfive * aimag( Qprime * ( P - conjg(P) )) * x9 ) + +! *** calculate gQs equation 15 + +! gQs = four_225 * abs(P)**2 * ( & +! (5.0d0 * Real(Vprime) + 3.0d0 * real(Sprime) ) * x6 & +! + one_210 * ( 35.0d0 * real(Vprime*conjg(Qprime) ) & +! + 21.0d0 * real(Sprime * conjg(Qprime) ) & +! + 10.0d0 * real(Wprime)- 6.0d0 * real(Tprime) ) * x8 & +! - twothrds * ( 5.0d0 * aimag(Vprime * conjg(P) ) & +! + 3.0d0 * aimag(Sprime * conjg(P) ) ) * x9 ) + +! FSB recast into specific terms + A = 1.0D0 * x4 + B = onefifteenth * real(Qprime) * x6 + C = fourthirds * aimag(P) * x7 + D = one_big * ( 35.0d0 * abs(Qprime)**2 & + + 20.0d0 * real(Rprime) + 35.0d0 * abs(Vprime)**2 & + + 21.0d0 * abs(Sprime)**2 ) * x8 + E = two_fortyfive * aimag( Qprime * ( P - conjg(P) )) * x9 + + Qs = eightthirds * abs(P)**2 *( A + B + C + D + E ) + + AA = (5.0d0 * Real(Vprime) + 3.0d0 * real(Sprime) ) * x6 + BB = one_210 * ( 35.0d0 * real(Vprime*conjg(Qprime) ) & + + 21.0d0 * real(Sprime * conjg(Qprime) ) & + + 10.0d0 * real(Wprime)- 6.0d0 * real(Tprime) ) * x8 + CC = twothrds * ( 5.0d0 * aimag(Vprime * conjg(P) ) & + + 3.0d0 * aimag(Sprime * conjg(P) ) ) * x9 + + gQs = four_225 * abs(P)**2 * ( AA + BB + CC ) + +! FSB calculate asymmetry factor and adjust with empirical term. + g = REAL(gQs / Qs) +! FSB now multiply by three_pi_two get output values + bext = three_pi_two * bext + bscat = three_pi_two * bscat +! FSB calculate babs + babs = bext - bscat + + end subroutine pennfsb + +! ------------------------------------------------------------------ + +! FSB a new version of Penndorf's equations. This version does +! analytical integration for Qext, Qscat, Qabs to generate +! bext, bscat, babs. Note that the expressions for Qext & Qscat +! hve been divide through by xx. +! +! References: +! Penndorf, R., 1962a,Scattering and extinction coefficients for small +! absorbing and nonabsorbing aerosols, +! J. Optical Society of America, 52, 896-904. + +! Penndorf, P., 1962b,Scattering and extinction coefficients for +! small Spherical aerosols, J. Atmos. Sci., 19, p 193 + +! FSB Coded by Dr. Francis S. Binkowski on October 25, 2011 by combining +! two previous versions to get a common code for the Penndorf and +! and Caldas & Semiao approaches. The Penndorf Qext, Qscat are much +! better than the versions from Caldas & Semiao despite claims to +! the contrary. The values of the asymmetry factor from Caldas & Semiao +! are better than can be obtained from Penndorf. +! FSB Modified by FSB on 12/02/2013 to remove the Caldas & Semiao parts. +! This version is set up to calculate LW bext and bscat so that absorption +! in the LW can be calculated. This will work with RRTMG LW. +! FSB This version does the analytical integral ove a lognormal +! size distribution. + + subroutine pennfsbLW (crefin, xx, lnsig, bext, bscat) + + implicit none +! input variables + complex, intent(in) :: crefin !complex refractive index + real, intent(in) :: xx ! pi * diameter / wavelength + real, intent(in) :: lnsig ! log(sigma_g) + real, intent(out) :: bext ! extinction coefficient + real, intent(out) :: bscat ! scattering coefficient +! internal variables + real*8 :: Qext, Qscat + real*8 :: lnsg2 +! internal variables + complex*16 :: m, m2,m21,m22 + complex*16 :: P,Q,S,V + real*8 :: P1,P2, Q1, Q2 , S2,V1, V2 ! see usage + real*8 :: P1SQ, P2SQ ! see usage + real*8 :: y, y2, y3 + real*8 :: x2, x3 +! FSB define useful numbers and fractions + real*8, parameter :: one= 1.0d0 + real*8, parameter :: two = 2.0d0 + real*8, parameter :: three = 3.0d0 + real*8, parameter :: four = 4.0d0 + real*8, parameter :: eight = 8.0d0 + real*8, parameter :: nine = 9.0d0 + real*8, parameter :: fifteen = 15.0d0 + real*8, parameter :: twothrds = two/three + real*8, parameter :: twofifteenths = two / fifteen + real*8, parameter :: eightthirds = eight / three + real*8, parameter :: nine_two = nine / two + +! FSB define useful numbers and fractions + real, parameter :: pi = four*atan(one) + real, parameter :: three_pi_two = 1.5d0 * pi + +! FSB start code + +! FSB Calculate integrals ove the lognormal distribution +! this is done term by term; the form is +! xn = yn * exp( (n**2) * lnsig2 /2.0d0) + lnsg2 = lnsig * lnsig + y = xx + y2 = y * y + y3 = y * y2 + x2 = y2 * exp( two * lnsg2) + x3 = y3 * exp( nine_two * lnsg2) + + m= conjg(crefin) ! Penndorf asuumes k is negative. + m2 = m * m + m21 = m2 - one + m22 = m2 + two +! FSB calculate Penndorf's definitions from Table II of Penndorf (1962a) + P = m21 / m22 + Q = (m2 - two ) / m22 + S = m21 / ( two * m2 + three) + V = m21 +! FSB get real & imaginary parts following Penndorf's mptation + P1 = real(P) + P2 = -aimag(P) + P1SQ = P1 * P1 + P2SQ = P2 * P2 + + Q1 = real(Q) + Q2 = -aimag(Q) + S2 = -aimag(S) + V1 = real(V) + v2 = -aimag(V) + +! FSB Get bext from Penndorf (1962a) Equation (7) up to x4 +! consistent with equation (8) +! We have then divided through by x and integrated analytically + Qext = four * P2 & + + ( 2.4d0 * (P1 * Q2 + P2 * Q1 ) + twothrds * S2 & + + twofifteenths * V2 ) * x2 & + + ( eightthirds * ( P1SQ - P2SQ ) ) * x3 + +! FSB get bscat from Penndorf Equation (9) up to x4 +! we have divided through by x and integrated analytically + Qscat = eightthirds * ( P1SQ + P2SQ ) * x3 + +! FSB now multiply by three_pi_two get output values + bext = three_pi_two * Qext + bscat = three_pi_two * Qscat + + end subroutine pennfsbLW + +! ------------------------------------------------------------------ + subroutine aero_optical2( lamda_in, crefin, Vol, dgn, & + sig, bext, bscat, gfac, success ) + +! FSB NOTE: this subroutine calculates for single mode + +! *** calculate the extinction and scattering coefficients and +! assymetry factors for each wavelength as a sum over the +! individual lognormal modes. Each mode may have a different +! set of refractive indices. + + IMPLICIT NONE + +! *** input variables + real, intent(in) :: lamda_in ! wavelengths [micro-m] + complex, intent(in) :: crefin ! Complex refractive index + real, intent(in) :: Vol ! modal aerosol volumes [m**3 /m**3] + real, intent(in) :: dgn ! geometric mean diameters + ! for number distribution [ m] + real, intent(in) :: sig ! geometric standard deviation + +! *** output variables + real, intent(out) :: bext ! extinction coefficient [ 1 / m ] + real, intent(out) :: bscat ! scattering coefficient [ 1 / m ] + real, intent(out) :: gfac ! assymetry factor for Mie and molecular scattering + logical, intent(out) :: success ! flag for successful calculation + +! *** internal variables +! real :: xlnsig(nmode) ! natural log of geometric standard deviations + real :: beta_Sc ! aerosol scattering coefficient + + real :: beta_Ex ! aerosol extinction coefficients + real :: G ! modal aerosol assymetry factors + real :: sum_g + real :: LSIGX + real :: lamdam1 ! 1/ lamda + real :: alphav ! Mie size parameter + real :: vfac + real, parameter :: pi = 3.14159265359 + + Logical, Save :: Initialize = .True. + +! FSB coded 04/15/2012 by Dr. Francis S. Binkowski +! modified from an earlier version +! Center for Environmental Modeling for PolicyDevelopment +! Institute for the Environment +! University of North Carolina at Chapel Hill +! email: frank_binkowski@unc.edu + +! *** initialize variables + lamdam1 = 1.0e6 / lamda_in ! lamda now in [ m ] + bext = 0.0 + bscat = 0.0 + sum_g = 0.0 + LSIGX = log(sig) + +! calculate Mie size parameter for volume distribution +! exp(3.0 * xlnsig*xlnsig) converts dgn to dgv (volume diameter) + alphav = pi * dgn * exp(3.0 * LSIGX * LSIGX) * lamdam1 + + If(Initialize .And. AERO_UTIL_LOG .GT. 0 )Then + If( Use_Odd_Quadrature )then + write(AERO_UTIL_LOG,99501)Quadrature_Points + else + write(AERO_UTIL_LOG,99504) + Initialize = .False. + End If + End If + + If( Use_Odd_Quadrature )then + CALL ghintBH (Initialize, crefin, alphav, LSIGX, beta_EX, beta_Sc, G, success) + Else + CALL ghintBH (crefin, alphav, LSIGX, beta_EX, beta_Sc, G, success) + End If + +! *** ghintBH returns the normalized values +! Calculate the actual extinction and scattering coefficients +! by multplying by the modal volume and dividing by the wavelength + + vfac = Vol * lamdam1 + bext = vfac * beta_Ex ! [ 1 / m ] + bscat = vfac * beta_Sc ! [ 1 / m ] + gfac = G +99501 Format(I2,' Quadrature Points for Volume Averaged Aerosol Optics') +99504 Format('Even Number Quadrature Points for Volume Averaged Aerosol Optics') + + END SUBROUTINE aero_optical2 + +! ------------------------------------------------------------------ + subroutine aero_optical_CS ( lamda_in, refcor,refshell, VOLCOR, & + VOLSHELL, DGNCOR, DGNSHELL, SIG, & + bext, bscat, gfac, succesS ) + +! FSB NOTE: values for one mode are returend +! *** calculate the extinction and scattering coefficients and +! assymetry factors for each wavelength as a sum over the +! individual lognormal modes. Each mode may have a different +! set of refractive indices. + + IMPLICIT NONE +! *** input variables + real,intent(in) :: lamda_in ! wavelengths [micro-m] + complex,intent(in) :: refcor ! Complex refractive index -core + complex,intent(in) :: refshell ! Complex refractive index -shell + real,intent(in) :: VOLCOR ! volume of core + real,intent(in) :: VOLSHELL ! volume of shell + real,intent(in) :: DGNCOR ! geometric mean diameters + ! for number distribution [m] + real,intent(in) :: DGNSHELL ! geometric mean diameters + ! for number distribution [m] + real,intent(in) :: SIG ! geometric standard deviation + +! *** output variables + real,intent(out) :: bext ! extinction coefficient [ 1 / m ] + real,intent(out) :: bscat ! scattering coefficient [ 1 / m ] + real,intent(out) :: gfac ! assymetry factor + logical, intent(OUT) :: success ! flag for successful calculation + +! *** internal variables +! real :: xlnsig(nmode) ! natural log of geometric standard deviations + real :: beta_Sc ! aerosol scattering coefficient + + real :: beta_Ex ! aerosol extinction coefficients + real :: G ! modal aerosol assymetry factors + real :: LSIGX + real :: XX, YY ! Mie size parameter + real :: expfac + real :: lamdam1 ! 1/ lamda + real :: vfac + + Logical, Save :: Initialize = .True. + + real, parameter :: pi = 3.14159265359 + +! FSB coded 04/15/2012 by Dr. Francis S. Binkowski +! modified from an earlier version +! Center for Environmental Modeling for PolicyDevelopment +! Institute for the Environment +! University of North Carolina at Chapel Hill +! email: frank_binkowski@unc.edu + + +! *** initialize variables + lamdam1 = 1.0e6 / lamda_in ! lamda now in [ m ] + +! calculate the extinction and scattering coefficients + LSIGX = log(SIG) + expfac = pi * exp(3.0 * LSIGX * LSIGX) * lamdam1 + +! calculate Mie size parameter for volume distribution +! exp(3.0 * xlnsig*xlnsig) converts dgn to dgv (volume diameter) + XX = DGNCOR * expfac + YY = DGNSHELL * expfac + + If(Initialize .And. AERO_UTIL_LOG .GT. 0 )Then + If( Use_Odd_Quadrature )then + write(AERO_UTIL_LOG,99500)Quadrature_Points + else + write(AERO_UTIL_LOG,99502) + Initialize = .False. + End If + End If + + If( Use_Odd_Quadrature )then + CALL ghintBH_CS(Initialize,refcor,refshell,XX,YY,LSIGX,beta_EX,beta_Sc,G, success) + Else + CALL ghintBH_CS(refcor,refshell,XX,YY,LSIGX,beta_EX,beta_Sc,G, success) + End If + +! FSB ghintBH_CS returns the normalized values +! Calculate the actual extinction and scattering coefficients +! by multplying by the modal volume and dividing by the wavelength. +! For the coated-sphere (core-shell) calculation use the combined +! volume + + vfac = (VOLCOR + VOLSHELL) * lamdam1 + bext = vfac * beta_Ex ! [ 1 / m ] + bscat = vfac * beta_Sc ! [ 1 / m ] + gfac = G +99500 Format(I2,' Quadrature Points for Core-Shell Aerosol Optics') +99502 Format('Even Number Quadrature Points for Core-Shell Aerosol Optics') + + END SUBROUTINE aero_optical_CS + +! ------------------------------------------------------------------ + + subroutine aero_optical_LW (lamda_in, crefin, Vol, & + dgn, sig, bext, bscat ) + +! *** calculate the extinction and scattering coefficients and +! asymmetry factor for LW radiation. +! The calling code only uses + + IMPLICIT NONE +! *** input variables + real, intent(in) :: lamda_in ! wavelengths [micro-m] + + complex, intent(in) :: crefin ! complex refractive index + real, intent(in) :: Vol ! modal aerosol volumes [m**3 /m**3] + real, intent(in) :: dgn ! geometric mean diameter + ! for number distribution [ m] + real, intent(in) :: sig ! geometric standard deviation + + +! *** output variables + real, intent(out) :: bext ! extinction coefficient [ 1 / m ] + real, intent(out) :: bscat ! scattering coefficient [ 1 / m ] + + +! *** internal variables + real, parameter :: pi = 3.14159265359 + + real :: lamda ! wavelength [ m] + real :: beta_Sc ! aerosol scattering coefficient + + real :: beta_Ex ! aerosol extinction coefficients + real :: G ! modal aerosol assymetry factors + real :: VLX, DGX, SIGX, LSIGX + real :: lamdam1 ! 1/ lamda + real :: alphav ! Mie size parameter + real vfac + + logical :: success + +! *** coded 09/08/2004 by Dr. Francis S. Binkowski +! FSB Modified for RRTMG version December 2009. +! FSB modified 10/06/2004, 10/12/2004, 10/18/2005 +! FSB 01/12/2006 +! FSB 12/02/2013 +! FSB Institute for the Environment +! University of North Carolina at Chapel Hill +! email: frank_binkowski@unc.edu + +! FSB Start Code: +! *** initialize variables +! lamda = 1.0e-6 * lamda_in ! lamda now in [ m ] + bext = 0.0 + bscat = 0.0 + +! lamdam1 = 1.0 / lamda ! 1 / [m] + lamdam1 = 1.0e6 / lamda_in ! 1 / [m] + + VLX = Vol + DGX = dgn + SIGX = sig + LSIGX = log(SIGX) + +! calculate Mie size parameter for volume distribution +! exp(3.0 * xlnsig*xlnsig) converts dgn to dgv (volume diameter) + alphav = pi * DGX * exp(3.0 * LSIGX * LSIGX) * lamdam1 + vfac = VLX * lamdam1 + + ! FSB Check for valid range of Penndorf application. + if ( alphav .le. 0.3) then + call pennfsbLW(crefin, alphav, LSIGX, beta_EX, beta_Sc) + G = 0.0 + else + CALL ghintBH(crefin, alphav, LSIGX, beta_EX, beta_Sc, G, success) + end if + + bext = vfac * beta_Ex ! [ 1 / m ] + bscat = vfac * beta_Sc + + END SUBROUTINE aero_optical_LW + +! ------------------------------------------------------------------ + subroutine ghintBH_2 (crefin,alfv,xlnsig,Qext_GH,Qscat_GH,g_gh, success) + +! *************** REVISED VERSION < NOTE +! FSB *********** This is the newest (04_14_2012) version of GhintBH +! this version does the Mie method and calculates the optimum set of +! set of Gauss-Hermite abscissas and weights. +! Dr. Francis S. Binkowski, The University of North Carolina +! at Chapel Hill +! FSB this code file now contains all of the necessary subroutines that +! are called to perform an integral of the Bohren and Huffman +! Mie codes ( as updated by Prof. Bruce C. Drain of Princeton) +! calculates the extinction and scattering coefficients +! normalized by wavelength and total particle volume +! concentration for a log normal particle distribution +! with the logarithm of the geometric standard deviation +! given by xlnsig. The integral of the +! asymmetry factor g is also calculated. +! FSB Change 12/20/2011 This code now has a choice of IGH based +! upon alfv and nr. +! FBB Changes Simplified code. Eliminated Penndorf code +! *** Does Gauss-Hermite quadrature of Qext / alfa & Qscat / alfa +! and asymmetry factor over log normal distribution using +! symmetric points. +! + implicit none + + complex, intent(in) :: crefin ! complex index of refraction + real, intent(in) :: alfv ! Mie parameter for dgv + real, intent(in) :: xlnsig ! log of geometric standard deviation + real, intent(out) :: Qext_GH ! normalized extinction efficiency + real, intent(out) :: Qscat_GH ! normalized scattering efficiency + real, intent(out) :: g_GH ! asymmetry factor + logical, intent(out) :: success ! flag for successful calculation + + real :: nr ! real part of refractive index + real :: aa1 ! see below for definition + real :: alfaip, alfaim ! Mie parameters at abscissas + +! *** these are Qext/alfa and Qscat/alfv at the abscissas + real :: qalfip_e, qalfim_e ! extinction + real :: qalfip_s, qalfim_s ! scattering + real :: gsalfp, gsalfm ! scattering times asymmetry factor + integer :: IGH ! index for GH quadrature + +! FSB define parameters + real, parameter :: pi = 3.14159265 + real, parameter :: sqrtpi = 1.772454 + real, parameter :: sqrtpi1 = 1.0 / sqrtpi + real, parameter :: sqrt2 = 1.414214 + real, parameter :: three_pi_two = 3.0 * pi / 2.0 + real, parameter :: const = three_pi_two * sqrtpi1 + + integer :: i + real :: sum_e,sum_s, xi,wxi,xf + real :: sum_sg + +! Gauss-Hermite abscissas and weights +! *** the following weights and abscissas are from Abramowitz +! Stegun, Table 25.10 page 924 +! FSB full precision from Table 25.10 + +! FSB ten-point - IGH = 5 + real, parameter :: ghxi_10(5) = (/ 0.342901327223705, & + 1.036610829789514, & + 1.756683649299882, & + 2.532731674232790, & + 3.436159118837738 /) + + real, parameter :: ghwi_10(5) = (/ 6.108626337353e-01, & + 2.401386110823e-01, & + 3.387439445548e-02, & + 1.343645746781e-03, & + 7.640432855233e-06 /) + +! FSB six-point - IGH = 3 + real, parameter :: ghxi_6(3) = (/ 0.436077411927617, & + 1.335849074013597, & + 2.350604973674492 /) + + real, parameter :: ghwi_6(3) = (/ 7.246295952244e-01, & + 1.570673203229e-01, & + 4.530009905509e-03 /) + +! FSB two-point - IGH = 1 + real, parameter :: ghxi_2(1) = (/ 0.707106781186548 /) + + real, parameter :: ghwi_2(1) = (/ 8.862269254528e-01 /) + + real :: GHXI(5), GHWI(5) ! weight and abscissas + integer :: NMAX ! number of weights and abscissa + + +! start code +! FSB now choose IGH. These choices are designed to improve +! the computational efficiency without sacrificing accuracy. + + nr = real(crefin) + + IGH=3 ! default value; six_point is sufficient generally +! six point + NMAX = 3 + + if (nr .ge. 1.7) then +! 10 point + IGH = 5 ! more points needed here + NMAX = 5 + end if + + if( alfv .gt. 20.0 .or. alfv .lt. 0.5 ) then + IGH = 1 ! in this range fewer points are needed + NMAX = 1 + end if + + if (IGH == 1) then +! two point + GHXI(1) = ghxi_2(1) + GHWI(1) = ghwi_2(1) + else if (IGH == 3) then + do i = 1, NMAX + GHXI(i) = ghxi_6(i) + GHWI(i) = ghwi_6(i) + end do + else + do i = 1,NMAX + GHXI(i) = ghxi_10(i) + GHWI(i) = ghwi_10(i) + end do + end if ! set up number of abscissas and weights + +! FSB now start the integration code + aa1 = sqrt2 * xlnsig ! This 1.0 / Sqrt( A ) in derivation of the integral + ! where A = 1.0 / ( 2.0 * xlnsg**2 ) + +! Then alpha = alfv * exp[ u / sqrt(A) ] +! For Gauss-Hermite Quadrature u = xi +! Therefore, xf = exp( xi / sqrt(A) ), +! or xf = exp( xi * aa1 ) + + sum_e = 0.0 + sum_s = 0.0 + sum_sg = 0.0 +! FSB do NMAX calls to the MIE codes + do i = 1,NMAX + xi = GHXI(i) + wxi = GHWI(i) + xf = exp( xi * aa1 ) + alfaip = alfv * xf + alfaim = alfv / xf ! division cheaper than another exp() +! *** call subroutine to fetch the effficiencies + + call getqext_BH(alfaip,crefin,qalfip_e,qalfip_s, gsalfp, success) + call getqext_BH(alfaim,crefin,qalfim_e,qalfim_s, gsalfm, success) + + sum_e = sum_e + wxi * ( qalfip_e + qalfim_e ) + sum_s = sum_s + wxi * ( qalfip_s + qalfim_s ) + sum_sg = sum_sg + wxi * ( gsalfp + gsalfm ) + + end do + + g_GH = sum_sg / sum_s ! this is + Qext_GH = const * sum_e ! + Qscat_GH = const * sum_s + + end subroutine ghintBH_2 + +! ------------------------------------------------------------------ + subroutine ghintBH_CS_even (RCORE, RSHELL , XX, YY, xlnsig, & + Qext_GH,Qscat_GH, g_gh, success) + +! FSB code for coated-sphere (core-shell) version + +! *************** REVISED VERSION < NOTE +! FSB *********** This is the newest (04_14_2012) version of ghintBH_CS +! for the coated-sphere (core-shell) method using BHCOAT +! this version does the Mie method and calculates the optimum set of +! set of Gauss-Hermite abscissas and weights. +! Dr. Francis S. Binkowski, The University of North Carolina +! at Chapel Hill + +! FSB this code file now contains all of the necessary subroutines that +! are called to perform an integral of the Bohren and Huffman +! Mie codes ( as updated by Prof. Bruce C. Drain of Princeton) +! calculates the extinction and scattering coefficients +! normalized by wavelength and total particle volume +! concentration for a log normal particle distribution +! with the logarithm of the geometric standard deviation +! given by xlnsig. The integral of the +! asymmetry factor g is also calculated. +! FSB Change 12/20/2011 This code now has a choice of IGH based +! upon alfv and nr. +! FBB Changes Simplified code. Eliminated Penndorf code +! *** Does Gauss-Hermite quadrature of Qext / alfa & Qscat / alfa +! and asymmetry factor over log normal distribution using +! symmetric points. +! + implicit none + complex, intent(in) :: RCORE ! refractive index of core + complex, intent(in) :: RSHELL ! refractive index of shell + real, intent(in) :: XX ! Mie parameter for core + real, intent(in) :: YY ! Mie parameter for shell + real, intent(in) :: xlnsig ! log of geometric standard deviation + real, intent(out) :: Qext_GH ! normalized extinction efficiency + real, intent(out) :: Qscat_GH ! normalized scattering efficiency + real, intent(out) :: g_GH ! asymmetry factor + logical, intent(out) :: success ! flag for successful calculation + + real :: nr ! real part of refractive index + real :: aa1 ! see below for definition + real :: XXP, XXM ! Mie parameters at abscissas - CORE + real :: YYP, YYM ! Mie parameters at abscissas - SHELL + +! FSB define parameters + real, parameter :: pi = 3.14159265 + real, parameter :: sqrtpi = 1.772454 + real, parameter :: sqrtpi1 = 1.0 / sqrtpi + real, parameter :: sqrt2 = 1.414214 + real, parameter :: three_pi_two = 3.0 * pi / 2.0 + real, parameter :: const = three_pi_two * sqrtpi1 + +! *** these are Qext/alfa and Qscat/alfv at the abscissas + real :: qalfip_e, qalfim_e ! extinction + real :: qalfip_s, qalfim_s ! scattering + real :: gsalfp, gsalfm ! scattering times asymmetry factor + integer :: IGH ! index for GH quadrature + integer :: i + real :: sum_e,sum_s, xi,wxi,xf, temp + real :: sum_sg + +! Gauss-Hermite abscissas and weights +! *** the following weights and abscissas are from Abramowitz +! Stegun, Table 25.10 page 924 +! FSB full precision from Table 25.10 + +! FSB ten-point - IGH = 5 + real, parameter :: ghxi_10(5) = (/ 0.342901327223705, & + 1.036610829789514, & + 1.756683649299882, & + 2.532731674232790, & + 3.436159118837738 /) + + real, parameter :: ghwi_10(5) = (/ 6.108626337353e-01, & + 2.401386110823e-01, & + 3.387439445548e-02, & + 1.343645746781e-03, & + 7.640432855233e-06 /) + +! FSB six-point - IGH = 3 + real, parameter :: ghxi_6(3) = (/ 0.436077411927617, & + 1.335849074013597, & + 2.350604973674492 /) + + real, parameter :: ghwi_6(3) = (/ 7.246295952244e-01, & + 1.570673203229e-01, & + 4.530009905509e-03 /) + +! FSB two-point - IGH = 1 + real, parameter :: ghxi_2(1) = (/ 0.707106781186548 /) + + real, parameter :: ghwi_2(1) = (/ 8.862269254528e-01 /) + + real GHXI(5), GHWI(5) ! weight and abscissas + integer NMAX ! number of weights and abscissa + +! start code +! FSB now choose IGH. These choices are designed to improve +! the computational efficiency without sacrificing accuracy. + + nr = real(RSHELL) + + IGH=3 ! default value; six_point is sufficient generally +! six point + NMAX = 3 + + if (nr .ge. 1.7) then +! 10 point + IGH = 5 ! more points needed here + NMAX = 5 + end if + + if ( XX .gt. 20.0 .or. XX .lt. 0.5 ) then + IGH = 1 ! in this range fewer points are needed + NMAX = 1 + end if + + if (IGH == 1) then +! two point + GHXI(1) = ghxi_2(1) + GHWI(1) = ghwi_2(1) + else if (IGH == 3) then + do i = 1, NMAX + GHXI(i) = ghxi_6(i) + GHWI(i) = ghwi_6(i) + end do + else + do i = 1,NMAX + GHXI(i) = ghxi_10(i) + GHWI(i) = ghwi_10(i) + end do + end if ! set up number of abscissas and weights + +! FSB now start the integration code + aa1 = sqrt2 * xlnsig ! This 1.0 / Sqrt( A ) in derivation of the integral + ! where A = 1.0 / ( 2.0 * xlnsg**2 ) + +! Then alpha = alfv * exp[ u / sqrt(A) ] +! For Gauss-Hermite Quadrature u = xi +! Therefore, xf = exp( xi / sqrt(A) ), +! or xf = exp( xi * aa1 ) + sum_e = 0.0 + sum_s = 0.0 + sum_sg = 0.0 +! FSB do NMAX calls to the MIE codes + do i = 1,NMAX + xi = GHXI(i) + wxi = GHWI(i) + xf = exp( xi * aa1 ) + temp = 1.0 / xf + XXP = XX * xf + XXM = XX * temp ! division cheaper than another exp() + YYP = YY * xf + YYM = YY * temp ! division cheaper than another exp() +! *** call subroutine to fetch the effficiencies + + call getqsgBHCS(XXP,YYP,RCORE,RSHELL,qalfip_e,qalfip_s,gsalfp, success) + call getqsgBHCS(XXM,YYM,RCORE,RSHELL,qalfim_e,qalfim_s,gsalfm, success) + + sum_e = sum_e + wxi * ( qalfip_e + qalfim_e ) + sum_s = sum_s + wxi * ( qalfip_s + qalfim_s ) + sum_sg = sum_sg + wxi * ( gsalfp + gsalfm ) + end do + + g_GH = sum_sg / sum_s ! this is + Qext_GH = const * sum_e ! + Qscat_GH = const * sum_s + + end subroutine ghintBH_CS_even + +! ------------------------------------------------------------------ + subroutine getqsgBHCS (XX,YY,RRFRL1,RRFRL2,qxtalf,qscalf,qsgalf, success) + implicit none + + real, intent(in) :: XX, YY + real, intent(out) :: qxtalf, qscalf, qsgalf + complex, intent(in) :: RRFRL1,RRFRL2 ! refractive indices Core , Shell + logical, intent(out) :: success ! flag for successful calculation + + real :: QEXT, QSCA, QBACK, G_MIE + real :: xx1 + character (len = 20) :: mystr1, mystr2, mystr3, mystr4 + + xx1 = 1.0 / YY + +! if ( (xx * real(RRFRL1) >= 30.0) & +! .or. (xx * aimag(RRFRL1) >= 30.0) & +! .or. (yy * aimag(RRFRL2) >= 30.0)) then +! print *, ' ==d== bhcoat error' +! end if + + call BHCOAT (XX,YY,RRFRL1,RRFRL2,QEXT,QSCA,QBACK,G_MIE, SUCCESS) + +! if ((trim(mystr1) == ' NaN') .or. & +! (trim(mystr2) == ' NaN') .or. & +! (trim(mystr3) == ' NaN') .or. & +! (trim(mystr4) == ' NaN')) then +! call BHCOAT (XX,YY,RRFRL1,RRFRL2,QEXT,QSCA,QBACK,G_MIE) +! end if + + qxtalf = QEXT * xx1 + qscalf = QSCA * xx1 + qsgalf = qscalf * G_MIE + + END subroutine getqsgBHCS + +! ------------------------------------------------------------------ + SUBROUTINE BHCOAT (XX, YY, RRFRL1, RRFRL2, QQEXT, QQSCA, QBACK, GGSCA, SUCCESS) + + use complex_number_module + + implicit none ! added by FSB + +! Arguments: + real, intent(in) :: XX,YY ! Defined below + complex, intent(in) :: RRFRL1,RRFRL2 ! Defined below + real, intent(out) :: QQEXT,QQSCA,QBACK ! Defined below + real, intent(out) :: GGSCA ! asymmetry factor added by FSB + logical,intent(out) :: success + +! Local variables: + + real*8, parameter :: DEL = 1.0D-08 + real*8, parameter :: ONE = 1.0D0, TWO = 2.0D0 +! complex*16, save :: II +! data II/(0.D0,1.D0)/ + type(complex_number) :: II + + integer :: IFLAG,N,NSTOP + + character (len = 128) :: mystr + +! ----------------------------------------------------------- +! del is the inner sphere convergence criterion +! ----------------------------------------------------------- + + real*8 :: CHI0Y,CHI1Y,CHIY,PSI0Y,PSI1Y,PSIY,QEXT,RN,QSCA,X,Y,YSTOP,GSCA + real*8 :: TWO_N_M_ONE, TWO_N_P_ONE + real*8 :: RY, RYY, RNRY, RN1, factor + +! complex*16 :: AMESS1,AMESS2,AMESS3,AMESS4,AN,ANCAP,AN1, BN,BNCAP,BN1, BRACK, & + type(complex_number) :: AMESS1,AMESS2,AMESS3,AMESS4,AN,ANCAP,AN1, BN,BNCAP,BN1, BRACK, & + CHI0X2,CHI0Y2,CHI1X2,CHI1Y2,CHIX2,CHIPX2,CHIPY2,CHIY2,CRACK, & + D0X1,D0X2,D0Y2,D1X1,D1X2,D1Y2,DNBAR,GNBAR, & + REFREL,RFREL1,RFREL2, XBACK,XI0Y,XI1Y,XIY, & + X1,X2,Y2,RCX1, RCX2,RCY2, FAC1, FAC2 + +!*********************************************************************** +! NOTES from Prof. Bruce T. Draine, Princeton University +! Subroutine BHCOAT calculates Q_ext, Q_sca, Q_back for coated sphere. +! All bessel functions computed by upward recurrence. +! Input: +! XX = 2*PI*RCORE*REFMED/WAVEL +! YY = 2*PI*RMANT*REFMED/WAVEL +! RFREL1 = REFCOR/REFMED +! RFREL2 = REFMAN/REFMED +! where REFCOR = complex refr.index of core) +! REFMAN = complex refr.index of mantle) +! REFMED = real refr.index of medium) +! RCORE = radius of core +! RMANT = radius of mantle +! WAVEL = wavelength of light in ambient medium +! +! Routine BHCOAT is taken from Bohren & Huffman (1983) +! Obtained from C.L.Joseph +! +! History: +! 92/11/24 (BTD) Explicit declaration of all variables +! April 30,2012 (FSB) added additional code to optimize +! run time by finding common terms and replacing multiple +! divisions by multiplication by a reciprocal. +! April 09, 2012 code transferred from BTD's BMHMIE to +! calculate the asymmetry factor by Prof. Francis S. Binkowski of +! The University of North Carolina at Chapel Hill. +! April 30,2012 (FSB) added additional code to optimize +! run time by finding common terms and replacing multiple +! divisions by multiplication by a reciprocal. +! July 16, 2010 more optimization by Dr. David Wong (DW) at US EPA + +! REFERENCE: +! Bohren, Craig F. and Donald R. Huffman, Absorption and +! Scattering of Light by Small Particles, Wiley-Interscience +! copyright 1983. Paperback Published 1998. +! This code was originally listed in Appendix B. pp 483-489. +! As noted above , the original code was subsequently +! modified by Prof. Bruce T. Draine of Princeton University. +! +! FSB The background for this code is discussed in Borhen & Huffman (1983) +! on pages 181-183 ( Equations 8.2 ) and on pages 483-484. +!*********************************************************************** +! +! Start Code + + SUCCESS = .TRUE. + + II = c_set(0.0D0, 1.0D0) + +! this technique will make the second 4 byte in the 8 byte variable be 0 +! rather than arbitrary digits to increase accuracy + write (mystr, *) xx, yy, real(RRFRL1), aimag(RRFRL1), real(RRFRL2), aimag(RRFRL2) + read (mystr, *) x, y, RFREL1, RFREL2 + +! X = XX +! Y = YY + RY = ONE / Y + RYY = RY * RY +! RFREL1%real_part = real(RRFRL1) +! RFREL1%imag_part = aimag(RRFRL1) +! RFREL2%real_part = real(RRFRL2) +! RFREL2%imag_part = aimag(RRFRL2) + x1 = c_mul(x, rfrel1) + x2 = c_mul(x, rfrel2) + y2 = c_mul(y, rfrel2) + RCX1 = c_div(ONE, X1) + RCX2 = c_div(ONE, X2) + RCY2 = c_div(ONE, Y2) + refrel = c_div(rfrel2, rfrel1) + ystop = y + 4.0 * y**0.3333 + 2.0 + nstop = INT( ystop ) + +! ----------------------------------------------------------- +! series terminated after nstop terms +! ----------------------------------------------------------- + +! initialize variables + d0x1 = c_div(c_cos(x1), c_sin(x1)) + d0x2 = c_div(c_cos(x2), c_sin(x2)) + d0y2 = c_div(c_cos(y2), c_sin(y2)) + + psi0y = cos(y) + psi1y = sin(y) + chi0y = -sin(y) + chi1y = cos(y) + + xi0y = c_sub(psi0y, c_mul(chi0y, II)) + xi1y = c_sub(psi1y, c_mul(chi1y, II)) + + chi0y2 = c_mul(-1.0d0, c_SIN(y2)) + chi1y2 = c_COS(y2) + chi0x2 = c_mul(-1.0d0, c_SIN(x2)) + chi1x2 = c_COS(x2) + qsca = 0.0d0 + qext = 0.0d0 + GSCA = 0.0d0 + xback = c_set(0.0d0, 0.0d0) + iflag = 0 + factor = 1.0d0 + +! FSB Start main loop + DO n = 1, nstop + rn = REAL( n, 8 ) + RN1 = ONE / RN + TWO_N_M_ONE = TWO * RN - ONE + TWO_N_P_ONE = TWO * RN + ONE + psiy = (TWO_N_M_ONE)*psi1y*RY - psi0y + chiy = (TWO_N_M_ONE)*chi1y*RY - chi0y + xiy = c_sub(psiy, c_mul(chiy, II)) + d1y2 = c_sub(c_div(ONE, c_sub(c_mul(rn, RCY2), d0y2)), c_mul(rn, RCY2)) + + IF (iflag .eq. 0) THEN +! *** Calculate inner sphere ancap, bncap +! and brack and crack + d1x1 = c_sub(c_div(ONE, c_sub(c_mul(rn, RCX1), d0x1)), c_mul(rn, RCX1)) + d1x2 = c_sub(c_div(ONE, c_sub(c_mul(rn, RCX2), d0x2)), c_mul(rn, RCX2)) + + chix2 = c_sub(c_mul(c_mul(TWO*rn - ONE, chi1x2), RCX2), chi0x2) + chiy2 = c_sub(c_mul(c_mul(TWO*rn - ONE, chi1y2), RCY2), chi0y2) + + chipx2 = c_sub(chi1x2, c_mul(c_mul(rn, chix2), RCX2)) + chipy2 = c_sub(chi1y2, c_mul(c_mul(rn, chiy2), RCY2)) + + ANCAP = c_sub(c_mul(c_mul(REFREL, D1X1), CHIX2), CHIPX2) + ANCAP = c_mul(ANCAP, c_sub(c_mul(CHIX2, D1X2), CHIPX2)) + ANCAP = c_div(c_sub(c_mul(REFREL, D1X1), D1X2), ANCAP) + + brack = c_mul(ancap, c_sub(c_mul(chiy2, d1y2), chipy2)) + + bncap = c_sub(c_mul(refrel, d1x2), d1x1) + bncap = c_div(bncap, c_sub(c_mul(refrel, chipx2), c_mul(d1x1, chix2))) + bncap = c_div(bncap, c_sub(c_mul(chix2, d1x2), chipx2)) + + crack = c_mul(bncap, c_sub(c_mul(chiy2, d1y2), chipy2)) +! *** calculate convergence test expressions +! for inner sphere. +! *** see pages 483-485 of Bohren & Huffman for +! definitions. + amess1 = c_mul(brack, chipy2) + amess2 = c_mul(brack, chiy2) + amess3 = c_mul(crack, chipy2) + amess4 = c_mul(crack, chiy2) + +! Now test for convergence for inner sphere +! All four criteria must be satisfied. See page 484 of B & H + IF (c_ABS(amess1) .LE. del*c_ABS(d1y2) .AND. & + (c_ABS(amess2) .LE. del) .AND. & + (c_ABS(amess3) .LE. del*c_ABS(d1y2)) .AND. & + (c_ABS(amess4) .LE. del) ) THEN +! convergence for inner sphere + brack = c_set(0.0D0,0.0D0) + crack = c_set(0.0D0,0.0D0) + iflag = 1 +! ELSE +! no convergence yet +! iflag = 0 + END IF + END IF ! test on iflag .eq. 0 + +! *** note usage of brack and crack See equations on +! Page 485 and discussion on pages 486 -487 of B & H + dnbar = c_sub(d1y2, c_mul(brack, chipy2)) + dnbar = c_div(dnbar, c_sub(ONE, c_mul(brack, chiy2))) + gnbar = c_sub(d1y2, c_mul(crack, chipy2)) + gnbar = c_div(gnbar, c_sub(ONE, c_mul(crack, chiy2))) +!*** Store previous values of an and bn for use +! in computation of g= + IF (N .GT. 1) THEN + AN1 = an + BN1 = bn + END IF +! *** update an and bn + RNRY = rn * RY + FAC1 = c_add(c_div(dnbar, rfrel2), RNRY) + + an = c_sub(c_mul(psiy, FAC1), psi1y) + an = c_div(an, c_sub(c_mul(FAC1, xiy), xi1y)) + FAC2 = c_add(c_mul(rfrel2, gnbar), RNRY) + bn = c_sub(c_mul(psiy, FAC2), psi1y) + bn = c_div(bn, c_sub(c_mul(FAC2, xiy), xi1y)) + +! *** Calculate sums for qsca, qext, xback + qsca = qsca + (TWO_N_P_ONE) * (c_ABS(an)**2 + c_ABS(bn)**2) + + qext = qext + TWO_N_P_ONE * (an%real_part + bn%real_part) + + FACTOR = FACTOR * (-1.0D0) + XBACK = c_add(XBACK, c_mul(TWO_N_P_ONE * FACTOR, c_sub(AN, BN))) + +! FSB calculate the sum for the asymmetry factor + + GSCA = GSCA + ((TWO_N_P_ONE)/(RN* (RN + ONE)))* & + (an%real_part*bn%real_part + an%imag_part*bn%imag_part) + + IF (n .GT. 1) THEN + + GSCA = GSCA + (RN - RN1) * & + (AN1%real_part*AN%real_part + AN1%imag_part*AN%imag_part + & + BN1%real_part*BN%real_part + BN1%imag_part*BN%imag_part) + + END IF +! continue update for next interation + psi0y = psi1y + psi1y = psiy + chi0y = chi1y + chi1y = chiy + xi1y = c_sub(psi1y, c_mul(chi1y, II)) + chi0x2 = chi1x2 + chi1x2 = chix2 + chi0y2 = chi1y2 + chi1y2 = chiy2 + d0x1 = d1x1 + d0x2 = d1x2 + d0y2 = d1y2 + END DO ! end of main loop + +!*** Have summed sufficient terms. +! Now compute QQSCA,QQEXT,QBACK,and GSCA + GGSCA = REAL( TWO * GSCA / qsca ) + QQSCA = REAL( TWO * qsca * RYY ) + QQEXT = REAL( TWO * qext * RYY ) + + QBACK = 0.5 * real((xback%real_part**2 + xback%imag_part**2) * RYY) + + end subroutine BHCOAT + +! ------------------------------------------------------------------ + subroutine ghintBH_Odd (INIT, crefin,alfv,xlnsig,Qext_GH,Qscat_GH,g_gh, success ) + +! *************** REVISED VERSION < NOTE +! FSB *********** This is the newest (04_14_2012) version of GhintBH +! this version does the Mie method and calculates the optimum set of +! set of Gauss-Hermite abscissas and weights. +! Dr. Francis S. Binkowski, The University of North Carolina +! at Chapel Hill +! FSB this code file now contains all of the necessary subroutines that +! are called to perform an integral of the Bohren and Huffman +! Mie codes ( as updated by Prof. Bruce C. Drain of Princeton) +! calculates the extinction and scattering coefficients +! normalized by wavelength and total particle volume +! concentration for a log normal particle distribution +! with the logarithm of the geometric standard deviation +! given by xlnsig. The integral of the +! asymmetry factor g is also calculated. +! FSB Change 12/20/2011 This code now has a choice of IGH based +! upon alfv and nr. +! FBB Changes Simplified code. Eliminated Penndorf code +! *** Does Gauss-Hermite quadrature of Qext / alfa & Qscat / alfa +! and asymmetry factor over log normal distribution using +! symmetric points. +! + implicit none + + logical, intent(INOUT) :: INIT ! initialize number of qudraure points + complex, intent(in) :: crefin ! complex index of refraction + real, intent(in) :: alfv ! Mie parameter for dgv + real, intent(in) :: xlnsig ! log of geometric standard deviation + real, intent(out) :: Qext_GH ! normalized extinction efficiency + real, intent(out) :: Qscat_GH ! normalized scattering efficiency + real, intent(out) :: g_GH ! asymmetry factor + logical, intent(out) :: success ! flag for successful calculation + + real :: nr ! real part of refractive index + real :: aa1 ! see below for definition + real :: alfaip, alfaim ! Mie parameters at abscissas + +! *** these are Qext/alfa and Qscat/alfv at the abscissas + real :: qalfip_e, qalfim_e ! extinction + real :: qalfip_s, qalfim_s ! scattering + real :: gsalfp, gsalfm ! scattering times asymmetry factor + +! FSB define parameters + real, parameter :: pi = 3.14159265 + real, parameter :: sqrtpi = 1.772454 + real, parameter :: sqrtpi1 = 1.0 / sqrtpi + real, parameter :: sqrt2 = 1.414214 + real, parameter :: three_pi_two = 3.0 * pi / 2.0 + real, parameter :: const = three_pi_two * sqrtpi1 + + integer :: i + real :: sum_e,sum_s, xi,wxi,xf + real :: sum_sg + + real, allocatable, save :: GHXI(:), GHWI(:) ! weight and abscissas + integer, save :: IGH ! number of weights and abscissa + integer, save :: NMAX ! optimumized number of weights and abscissa + + +! start code +! FSB now choose IGH. These choices are designed to improve +! the computational efficiency without sacrificing accuracy. + + If( INIT )Then + + Select Case( Quadrature_Points ) + Case( 1,3,9 ) + IGH = Quadrature_Points + Case Default + IGH = 3 + End Select + + NMAX = Max( Int( IGH / 2 ), 0) + + If( Allocated( GHXI ) .Or. Allocated( GHWI ) )Then + Success = .False. + Return + End If + + Allocate( GHXI( NMAX + 1 ), GHWI( NMAX + 1 ) ) + + Select Case ( IGH ) + Case ( 1 ) + GHXI(1) = ghxi_1(1) + GHWI(1) = ghwi_1(1) + Case ( 3 ) + do i = 1, NMAX + 1 + GHXI(i) = ghxi_3(i) + GHWI(i) = ghwi_3(i) + end do + Case ( 9 ) + do i = 1, NMAX + 1 + GHXI(i) = ghxi_9(i) + GHWI(i) = ghwi_9(i) + end do + end select + + If( AERO_UTIL_LOG .GT. 0 )Then + write(AERO_UTIL_LOG,*)'BHMIE: IGH,(NMAX + 1) = ',IGH,(NMAX + 1) + do i = 1, NMAX + 1 + write(AERO_UTIL_LOG,*)'BHMIE: i, GHXI(i), GHWI(i) = ',i, GHXI(i), GHWI(i) + end do + End If + + INIT = .False. + Else + If( .Not. Allocated( GHXI ) .Or. .Not. Allocated( GHWI ) )Then + Success = .False. + Return + End If + End If ! set up number of abscissas and weights + + nr = real(crefin) + +! FSB now start the integration code + aa1 = sqrt2 * xlnsig ! This 1.0 / Sqrt( A ) in derivation of the integral + ! where A = 1.0 / ( 2.0 * xlnsg**2 ) + +! Then alpha = alfv * exp[ u / sqrt(A) ] +! For Gauss-Hermite Quadrature u = xi +! Therefore, xf = exp( xi / sqrt(A) ), +! or xf = exp( xi * aa1 ) + +!start integration at zero point + xi = 0.0 + wxi = GHWI(NMAX+1) + xf = 1.0 + alfaip = alfv +! fetch the effficiencies at zero point + + call getqext_BH(alfaip,crefin,qalfip_e,qalfip_s, gsalfp, success) + + sum_e = wxi * qalfip_e + sum_s = wxi * qalfip_s + sum_sg = wxi * gsalfp + +! FSB do NMAX calls to the MIE codes + do i = 1, NMAX + xi = GHXI(i) + wxi = GHWI(i) + xf = exp( xi * aa1 ) + alfaip = alfv * xf + alfaim = alfv / xf ! division cheaper than another exp() +! *** call subroutine to fetch the effficiencies + + call getqext_BH(alfaip,crefin,qalfip_e,qalfip_s, gsalfp, success) + call getqext_BH(alfaim,crefin,qalfim_e,qalfim_s, gsalfm, success) + + sum_e = sum_e + wxi * ( qalfip_e + qalfim_e ) + sum_s = sum_s + wxi * ( qalfip_s + qalfim_s ) + sum_sg = sum_sg + wxi * ( gsalfp + gsalfm ) + + end do + + g_GH = sum_sg / sum_s ! this is + Qext_GH = const * sum_e ! + Qscat_GH = const * sum_s + + end subroutine ghintBH_Odd + +! ------------------------------------------------------------------ + subroutine ghintBH_CS_Odd (INIT, RCORE, RSHELL , XX, YY, xlnsig, & + Qext_GH,Qscat_GH, g_gh, success) + +! FSB code for coated-sphere (core-shell) version + +! *************** REVISED VERSION < NOTE +! FSB *********** This is the newest (04_14_2012) version of ghintBH_CS +! for the coated-sphere (core-shell) method using BHCOAT +! this version does the Mie method and calculates the optimum set of +! set of Gauss-Hermite abscissas and weights. +! Dr. Francis S. Binkowski, The University of North Carolina +! at Chapel Hill + +! FSB this code file now contains all of the necessary subroutines that +! are called to perform an integral of the Bohren and Huffman +! Mie codes ( as updated by Prof. Bruce C. Drain of Princeton) +! calculates the extinction and scattering coefficients +! normalized by wavelength and total particle volume +! concentration for a log normal particle distribution +! with the logarithm of the geometric standard deviation +! given by xlnsig. The integral of the +! asymmetry factor g is also calculated. +! FSB Change 12/20/2011 This code now has a choice of IGH based +! upon alfv and nr. +! FBB Changes Simplified code. Eliminated Penndorf code +! *** Does Gauss-Hermite quadrature of Qext / alfa & Qscat / alfa +! and asymmetry factor over log normal distribution using +! symmetric points. +! + implicit none + + logical, intent(inout) :: INIT ! initialize number of qudraure points + complex, intent(in) :: RCORE ! refractive index of core + complex, intent(in) :: RSHELL ! refractive index of shell + real, intent(in) :: XX ! Mie parameter for core + real, intent(in) :: YY ! Mie parameter for shell + real, intent(in) :: xlnsig ! log of geometric standard deviation + real, intent(out) :: Qext_GH ! normalized extinction efficiency + real, intent(out) :: Qscat_GH ! normalized scattering efficiency + real, intent(out) :: g_GH ! asymmetry factor + logical, intent(out) :: success ! flag for successful calculation + + real :: nr ! real part of refractive index + real :: aa1 ! see below for definition + real :: XXP, XXM ! Mie parameters at abscissas - CORE + real :: YYP, YYM ! Mie parameters at abscissas - SHELL + +! FSB define parameters + real, parameter :: pi = 3.14159265 + real, parameter :: sqrtpi = 1.772454 + real, parameter :: sqrtpi1 = 1.0 / sqrtpi + real, parameter :: sqrt2 = 1.414214 + real, parameter :: three_pi_two = 3.0 * pi / 2.0 + real, parameter :: const = three_pi_two * sqrtpi1 + +! *** these are Qext/alfa and Qscat/alfv at the abscissas + real :: qalfip_e, qalfim_e ! extinction + real :: qalfip_s, qalfim_s ! scattering + real :: gsalfp, gsalfm ! scattering times asymmetry factor + integer :: i + real :: sum_e,sum_s, xi,wxi,xf, temp + real :: sum_sg + + real, allocatable, save :: GHXI(:), GHWI(:) ! weight and abscissas + integer, save :: IGH ! number of weights and abscissa + integer, save :: NMAX ! optimized number of weights and abscissa + +! start code +! FSB now choose IGH. These choices are designed to improve +! the computational efficiency without sacrificing accuracy. + + If( INIT )Then + + Select Case( Quadrature_Points ) + Case( 1,3,9 ) + IGH = Quadrature_Points + Case Default + IGH = 3 + End Select + + If( Allocated( GHXI ) .Or. Allocated( GHWI ) )Then + Success = .False. + Return + End If + + NMAX = Max( Int( IGH / 2 ), 0) + + Allocate( GHXI( NMAX + 1 ), GHWI( NMAX + 1 ) ) + + Select Case ( IGH ) + Case ( 1 ) + GHXI(1) = ghxi_1(1) + GHWI(1) = ghwi_1(1) + Case ( 3 ) + do i = 1, NMAX + 1 + GHXI(i) = ghxi_3(i) + GHWI(i) = ghwi_3(i) + end do + Case ( 9 ) + do i = 1, NMAX + 1 + GHXI(i) = ghxi_9(i) + GHWI(i) = ghwi_9(i) + end do + end select + + If( AERO_UTIL_LOG .GT. 0 )Then + write(AERO_UTIL_LOG,*)'BHCoat: IGH,(NMAX + 1) = ',IGH,(NMAX + 1) + do i = 1, NMAX + 1 + write(AERO_UTIL_LOG,*)'BHCoat: i, GHXI(i), GHWI(i) = ',i, GHXI(i), GHWI(i) + end do + End If + + INIT = .False. + + Else + If( .Not. Allocated( GHXI ) .Or. .Not. Allocated( GHWI ) )Then + Success = .False. + Return + End If + End If ! set up number of abscissas and weights + + nr = real(RSHELL) + +! FSB now start the integration code + aa1 = sqrt2 * xlnsig ! This 1.0 / Sqrt( A ) in derivation of the integral + ! where A = 1.0 / ( 2.0 * xlnsg**2 ) + +! Then alpha = alfv * exp[ u / sqrt(A) ] +! For Gauss-Hermite Quadrature u = xi +! Therefore, xf = exp( xi / sqrt(A) ), +! or xf = exp( xi * aa1 ) + +!start integration at zero point + + xi = 0.0 + wxi = GHWI(NMAX+1) + xf = 1.0 + XXP = XX + YYP = YY + +! fetch the effficiencies at zero point + + call getqsgBHCS(XXP,YYP,RCORE,RSHELL,qalfip_e,qalfip_s,gsalfp, success) + + sum_e = wxi * qalfip_e + sum_s = wxi * qalfip_s + sum_sg = wxi * gsalfp + +! FSB do NMAX calls to the MIE codes + do i = 1, NMAX + xi = GHXI(i) + wxi = GHWI(i) + xf = exp( xi * aa1 ) + temp = 1.0 / xf + XXP = XX * xf + XXM = XX * temp ! division cheaper than another exp() + YYP = YY * xf + YYM = YY * temp ! division cheaper than another exp() +! *** call subroutine to fetch the effficiencies + + call getqsgBHCS(XXP,YYP,RCORE,RSHELL,qalfip_e,qalfip_s,gsalfp, success) + call getqsgBHCS(XXM,YYM,RCORE,RSHELL,qalfim_e,qalfim_s,gsalfm, success) + + sum_e = sum_e + wxi * ( qalfip_e + qalfim_e ) + sum_s = sum_s + wxi * ( qalfip_s + qalfim_s ) + sum_sg = sum_sg + wxi * ( gsalfp + gsalfm ) + end do + + g_GH = sum_sg / sum_s ! this is + Qext_GH = const * sum_e ! + Qscat_GH = const * sum_s + + end subroutine ghintBH_CS_Odd + +! ------------------------------------------------------------------ + SUBROUTINE BHMIE_FLEXI (X, NMX, NSTOP, REFREL, QQEXT, QQSCA, QBACK, GSCA, SUCCESS) + +! FSB Changed the call vector to return only QEXT, QSCAT QBACK GSCA +! and ignore NANG, S1 and S2 and all calculations for them + + implicit none + +! Arguments: + real, intent(in) :: X ! X = pi*particle_diameter / Wavelength + integer, intent(in) :: NMX ! maximum number of terms in Mie series + integer, intent(in) :: NSTOP ! minumum number of terms in Mie series + complex, intent(in) :: REFREL ! refractive index + +! REFREL = (complex refr. index of sphere)/(real index of medium) +! in the current use the index of refraction of the the medium +! i taken at 1.0 real. +! +! Output + + real, intent(out) :: QQEXT, QQSCA, QBACK, GSCA + logical, intent(out) :: SUCCESS + +! QQEXT Efficiency factor for extinction +! QQSCA Efficiency factor for scattering +! QQBACK Efficiency factor for back scatter +! GSCA asymmetry factor +! SUCCESS flag for successful calculation +! REFERENCE: +! Bohren, Craig F. and Donald R. Huffman, Absorption and +! Scattering of Light by Small Particles, Wiley-Interscience +! copyright 1983. Paperback Published 1998. +! FSB +! This code was originally listed in Appendix A. pp 477-482. +! As noted below, the original code was subsequently +! modified by Prof. Bruce T. Drain of Princetion University. +! The code was further modified for a specific application +! in a large three-dimensional code requiring as much +! computational efficiency as possible. +! Prof. Francis S. Binkowski of The University of North +! Carolina at Chapel Hill. + +! Declare parameters: +! Note: important that MXNANG be consistent with dimension of S1 and S2 +! in calling routine! + + integer, parameter :: MXNANG=10, NMXX=150000 ! FSB new limits + integer, parameter :: NANG = 2 + real*8, parameter :: PII = 3.1415916536D0 + real*8, parameter :: ONE = 1.0D0, TWO = 2.0D0 + complex*16, parameter :: COMPLEX_DZERO = (0.0D0,0.0D0) + complex, parameter :: COMPLEX_ZERO = (0.0,0.0) + +! Local variables: + integer :: N, NN + real*8 :: QSCA, QEXT, DX1, DXX1 + real*8 :: CHI,CHI0,CHI1,DX,EN,P,PSI,PSI0,PSI1,XSTOP,YMOD + real*8 :: TWO_N_M_ONE, TWO_N_P_ONE, EN1, FACTOR + complex*16 :: AN,AN1,BN,BN1,DREFRL,XI,XI1,Y, Y1, DREFRL1 + complex*16 :: D(NMX) + complex*16 :: FAC1, FAC2 + complex*16 :: XBACK + +!*********************************************************************** +! Subroutine BHMIE is the Bohren-Huffman Mie scattering subroutine +! to calculate scattering and absorption by a homogenous isotropic +! sphere. +! Given: +! X = 2*pi*a/lambda +! REFREL = (complex refr. index of sphere)/(real index of medium) +! real refractive index of medium taken as 1.0 +! Returns: +! QEXT = efficiency factor for extinction +! QSCA = efficiency factor for scattering +! QBACK = efficiency factor for backscatter +! see Bohren & Huffman 1983 p. 122 +! GSCA = asymmetry for scattering +! +! Original program taken from Bohren and Huffman (1983), Appendix A +! Modified by Prof. Bruce T.Draine, Princeton Univ. Obs., 90/10/26 +! in order to compute +! 91/05/07 (BTD): Modified to allow NANG=1 +! 91/08/15 (BTD): Corrected error (failure to initialize P) +! 91/08/15 (BTD): Modified to enhance vectorizability. +! 91/08/15 (BTD): Modified to make NANG=2 if called with NANG=1 +! 91/08/15 (BTD): Changed definition of QBACK. +! 92/01/08 (BTD): Converted to full double precision and double complex +! eliminated 2 unneed lines of code +! eliminated redundant variables (e.g. APSI,APSI0) +! renamed RN -> EN = double precision N +! Note that DOUBLE COMPLEX and DCMPLX are not part +! of f77 standard, so this version may not be fully +! portable. In event that portable version is +! needed, use src/bhmie_f77.f +! 93/06/01 (BTD): Changed AMAX1 to generic function MAX +! FSB April 09,2012 This code was modified by: +! Prof. Francis S. Binkowski University of North Carolina at +! Chapel Hill, Institue for the Environment. +! +! The modifications were made to enhance computation speed +! for use in a three-dimensional code. This was done by +! removing code that calculated angular scattering. The method +! of calculating QEXT, QBACK was also changed. + +!*********************************************************************** +!*** Safety checks + + SUCCESS = .TRUE. +! NANG = 2 ! FSB only this value +! IF(NANG.GT.MXNANG)STOP'***Error: NANG > MXNANG in bhmie' +! IF (NANG .LT. 2) NANG = 2 + + DX = REAL( X, 8 ) +! FSB Define reciprocals so that divisions can be replaced by multiplications. + DX1 = ONE / DX + DXX1 = DX1 * DX1 + DREFRL = DCMPLX( REAL( REFREL ), IMAG( REFREL ) ) + DREFRL1 = ONE / DREFRL + Y = DX * DREFRL + Y1 = ONE / Y +! YMOD = ABS(Y) + +!*** Series expansion terminated after NSTOP terms +! Logarithmic derivatives calculated from NMX on down +! XSTOP = X + 4.0 * X**0.3333 + 2.0 +! NMX = MAX(XSTOP,YMOD) + 15 + +! BTD experiment 91/1/15: add one more term to series and compare results +! NMX=AMAX1(XSTOP,YMOD)+16 +! test: compute 7001 wavelengths between .0001 and 1000 micron +! for a=1.0micron SiC grain. When NMX increased by 1, only a single +! computed number changed (out of 4*7001) and it only changed by 1/8387 +! conclusion: we are indeed retaining enough terms in series! + + FACTOR = 1.0D0 + +! IF (NMX .GT. NMXX) THEN +! WRITE(6,*)'Error: NMX > NMXX=',NMXX,' for |m|x=',YMOD +! SUCCESS = .FALSE. +! RETURN +! END IF + +! FSB all code relating to scattering angles is removed out for +! reasons of efficiency when running in a three-dimensional +! code. We only need QQSCA, QQEXT, GSCA AND QBACK + +!*** Logarithmic derivative D(J) calculated by downward recurrence +! beginning with initial value (0.,0.) + + D(NMX) = COMPLEX_DZERO + NN = NMX - 1 + DO N = 1,NN + EN = REAL( NMX - N + 1, 8 ) +! FSB In the following division by Y has been replaced by +! multiplication by Y1, the reciprocal of Y. + D(NMX-N) = ( EN * Y1 ) - (ONE / ( D(NMX-N+1) + EN * Y1)) + END DO + +!*** Riccati-Bessel functions with real argument X +! calculated by upward recurrence + + PSI0 = COS(DX) + PSI1 = SIN(DX) + CHI0 = -SIN(DX) + CHI1 = PSI0 + XI1 = DCMPLX(PSI1,-CHI1) + QSCA = 0.0D0 + GSCA = 0.0D0 + QEXT = 0.0D0 + P = -ONE + XBACK = COMPLEX_DZERO + +! FSB Start main loop + DO N = 1,NSTOP + EN = REAL( N, 8 ) + EN1 = ONE / EN + TWO_N_M_ONE = TWO * EN - ONE +! for given N, PSI = psi_n CHI = chi_n +! PSI1 = psi_{n-1} CHI1 = chi_{n-1} +! PSI0 = psi_{n-2} CHI0 = chi_{n-2} +! Calculate psi_n and chi_n + PSI = TWO_N_M_ONE * PSI1 * DX1 - PSI0 + CHI = TWO_N_M_ONE * CHI1 * DX1 - CHI0 + XI = DCMPLX(PSI,-CHI) + +!*** Compute AN and BN: +! FSB Rearrange to get common terms + FAC1 = D(N) * DREFRL1 + EN * DX1 + AN = (FAC1) * PSI - PSI1 + AN = AN / ( (FAC1 )* XI - XI1 ) + FAC2 = ( DREFRL * D(N) + EN * DX1) + BN = ( FAC2) * PSI -PSI1 + BN = BN / ((FAC2) * XI - XI1 ) + +! FSB calculate sum for QEXT as done by Wiscombe +! get common factor + TWO_N_P_ONE = (TWO * EN + ONE) + QEXT = QEXT + (TWO_N_P_ONE) * (REAL(AN) + REAL(BN) ) + QSCA = QSCA + (TWO_N_P_ONE) * ( ABS(AN)**2 + ABS(BN)**2 ) + +! FSB calculate XBACK from B & H Page 122 + FACTOR = -1.0d0 * FACTOR ! calculate (-1.0 ** N) + XBACK = XBACK + (TWO_N_P_ONE) * factor * (AN - BN) + +! FSB calculate asymmetry factor + + GSCA = GSCA + REAL((TWO_N_P_ONE)/(EN * (EN + ONE)) * & + (REAL(AN)*REAL(BN)+IMAG(AN)*IMAG(BN))) + + IF (N .GT. 1)THEN + GSCA = GSCA + REAL((EN - EN1) * & + (REAL(AN1)*REAL(AN) + IMAG(AN1)*IMAG(AN) + & + REAL(BN1)*REAL(BN) + IMAG(BN1)*IMAG(BN))) + ENDIF + +!*** Store previous values of AN and BN for use in computation of g= + AN1 = AN + BN1 = BN + +! FSB set up for next iteration + PSI0 = PSI1 + PSI1 = PSI + CHI0 = CHI1 + CHI1 = CHI + XI1 = DCMPLX(PSI1,-CHI1) + + END DO ! main loop on n + +!*** Have summed sufficient terms. + +! Now compute QQSCA,QQEXT,QBACK,and GSCA + GSCA = REAL( TWO / QSCA ) * GSCA + +! FSB in the following, divisions by DX * DX has been replaced by +! multiplication by DXX1 the reciprocal of 1.0 / (DX *DX) + QQSCA = REAL( TWO * QSCA * DXX1 ) + QQEXT = REAL( TWO * QEXT * DXX1 ) + QBACK = REAL( REAL( 0.5D0 * XBACK * CONJG(XBACK), 8 ) * DXX1 ) ! B&H Page 122 + + END subroutine BHMIE_FLEXI + +END MODULE rrtmg_aero_optical_util_module + +! ------------------------------------------------------------------ +! FSB REvised Mie calculations 02/09/2011 + +MODULE module_twoway_ra_rrtmg_sw + +contains + +!------------------------------------------------------------------ + Subroutine get_aerosol_Optics_RRTMG_SW ( ns, nmode,delta_z, INMASS_ws, & + INMASS_in, INMASS_ec, INMASS_ss, & + INMASS_h2o, INDGN, INSIG, & + tauaer, waer, gaer ) + +!FSB This version switches between BHCOAT to BHMIE depending upon whether +! EC is present or not. 04/15/2012. + +!FSB this version does a core-shell calculation with BHCOAT 04/11/2012 +! This version is set up to be used with RRTMG_SW <<<<<<<< +! wavelenght is calculated internally +! FSB This routine calculates the aerosol information ( tauaer, waer, +! gaer, needed to calculate the solar radiation) The calling +! program specifies the location ( row, column, layer, +! layer thicknes, and wave length for the calculation. +! FSB 02/09/2011 Modifications made to subroutine ghintBH. +! FSB 04/14/2012 REmoved MODULUS, made changes to ghintBH. +! Put in option for core-shell (coated-sphere). 2 + +! FSB Input variables: + + use rrtmg_aero_optical_util_module + + implicit none + + integer,intent(in) :: ns ! index for wavelength should be + ! between 1 and 14. <<< RRTMG_SW + integer,intent(in) :: nmode ! should be 3 for WRF/CMAQ calculation + real,intent(in) :: delta_z ! layer thickness [m] +! FSB mode types for WRF/CMAQ +! nmode = 1 Aitken +! nmode = 2 accumulation +! nmode = 3 coarse +! FSB modal mass concentration by species [ ug / m**3] NOTE: MKS + real, intent(in) :: INMASS_ws(nmode) ! water soluble + real, intent(in) :: INMASS_in(nmode) ! insolugle + real, intent(in) :: INMASS_ec(nmode) ! elemental carbon or soot like + real, intent(in) :: INMASS_ss(nmode) ! sea salt + real, intent(in) :: INMASS_h2o(nmode) ! water +! FSB particle size-distribution information + real, intent(in) :: INDGN( nmode) ! geometric mean diameter [ m ] NOTE: MKS + real, intent(in) :: INSIG( nmode) ! geometric standard deviation + +!FSB output aerosol radiative properties [dimensionless] + real, intent(out) :: tauaer ! aerosol extinction optical depth + real, intent(out) :: waer ! aerosol single scattering albedo + real, intent(out) :: gaer ! aerosol assymetry parameter + +! FSB Internal variables + + real :: NR(nmode), NI(nmode) ! refractive indices + complex :: refcor(nmode), refshell(nmode) ! complex refracive indices + complex :: crefin(nmode) ! complex refractive index + +! FSB special values for EC CORE-shell calculation + real :: DGNSHELL(nmode) ! modal geometric mean diameter [m] + real :: DGNCORE (nmode) ! modal geometric mean diameter [m] + +! FSB Modal volumes [ m**3 / m**3 ] + real :: MVOL_ws(nmode) ! water soluble + real :: MVOL_in(nmode) ! insolugle + real :: MVOL_ec(nmode) ! soot like + real :: MVOL_ss(nmode) !sea salt + real :: MVOL_h2o(nmode) ! water +! real :: VOL(nmode) ! total modal volume [m** 3 / m**3] +! FSB special values for EC CORE-shell calculation + real :: VOLCOR(nmode) ! volume of EC core [m** 3 / m**3] + real :: VOLSHELL(nmode) ! volume of shell [m** 3 / m**3] + + integer :: m ! loop index + real :: bext ! extinction coefficient [1 / m] + real :: bscat ! scattering coefficient [1 / m] + real :: gfac ! asymmetry factor + + real :: bextsum, bscatsum, bsgsum + +! FSB History variables by wavelength and mode +! real :: bext_wm(ns,nmode) +! real :: bscat_wm(ns,nmode) +! real :: gfac_wm(ns,nmode) + + real, parameter :: one3rd = 1.0 / 3.0 + real :: dfac ! ratio of (volcor/vol) ** one3rd + ! used for calculating the diameter + ! of the EC core + + logical :: succesS + +!...component densities [ g/ cm**3 ] <<<<< cgs + + real, parameter :: rhows = 1.8 ! bulk density of water soluble aerosol + + real, parameter :: rhoin = 2.2 ! bulk density forinsoluble aerosol + +! real, parameter :: rhoec = 1.7 ! bulk density for soot aerosol + real, parameter :: rhoec = 1.8 ! new value + + real, parameter :: rhoh2o = 1.0 ! bulk density of aerosol water + + real, parameter :: rhoss = 2.2 ! bulk density of seasalt + +! FSB scale factor for volume calculation +! 1.0d-12 * [ cm**3 / g] -> [ m** 3 / ug ] + real, parameter :: scalefactor = 1.0e-12 + +! FSB scale factor for [1/g] to [1/ug] + real, parameter :: cug2g = 1.0e-06 + +! FSB reciprocal component densities[ m ** 3 / ug ] + + real, parameter :: rhows1 = scalefactor / rhows ! water soluble aerosol + + real, parameter :: rhoin1 = scalefactor / rhoin ! insoluble aerosol + + real, parameter :: rhoec1 = scalefactor / rhoec ! soot aerosol + + real, parameter :: rhoh2o1 = scalefactor / rhoh2o ! aerosol water + + real, parameter :: rhoss1 = scalefactor / rhoss ! seasalt + + integer,parameter :: nspint_sw = 14 ! number of spectral intervals for RRTMG_SW + +! FSB Band numbers and wavelengths for RRTMG_SW + integer, parameter :: Band(nspint_sw) = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 /) + + real, parameter :: LAMDA_SW(nspint_sw) = (/ 3.4615, 2.7885, 2.325, 2.046, 1.784, & + 1.4625, 1.2705, 1.0101, 0.7016, 0.53325, & + 0.38815, 0.299, 0.2316, 8.24 /) ! wavelength [ um ] + +! *** refractive indices + +! *** Except as otherwise noted reference values of refractive +! indices for aerosol particles are from the OPAC Data base. +! Hess, Koepke, and Schult, Optical properties of +! aerosols and clouds: The software package OPAC, Bulletan of +! the American Meteorological Society, Vol 79, No 5, +! pp 831 - 844, May 1998. +! OPAC is a downloadable data set of optical properties of +! 10 aerosol components, 6 water clouds and 3 cirrus clouds +! at UV, visible and IR wavelengths +! www.lrz-muenchen.de/~uh234an/www/radaer/opac.htm + + +! FSB water soluble + real, parameter :: xnreal_ws(nspint_sw) = (/ 1.443, 1.420, 1.420, 1.420, 1.463, 1.510, 1.510, & + 1.520, 1.530, 1.530, 1.530, 1.530, 1.530, 1.710 /) + real, parameter :: xnimag_ws(nspint_sw) = (/ 5.718E-3, 1.777E-2, 1.060E-2, 8.368E-3, 1.621E-2, & + 2.198E-2, 1.929E-2, 1.564E-2, 7.000E-3, 5.666E-3, & + 5.000E-3, 8.440E-3, 3.000E-2, 1.100E-1 /) + +! FSB sea salt + real, parameter :: xnreal_ss(nspint_sw) = (/ 1.480, 1.534, 1.437, 1.448, 1.450, 1.462, 1.469, & + 1.470, 1.490, 1.500, 1.502, 1.510, 1.510, 1.510 /) + + real, parameter :: xnimag_ss(nspint_sw) = (/ 1.758E-3, 7.462E-3, 2.950E-3, 1.276E-3, 7.944E-4, & + 5.382E-4, 3.754E-4, 1.498E-4, 2.050E-7, 1.184E-8, & + 9.938E-8, 2.060E-6, 5.000E-6, 1.000E-2 /) + +! FSB insoluble + real, parameter :: xnreal_in(nspint_sw) = (/ 1.272, 1.168, 1.208, 1.253, 1.329, 1.418, 1.456, & + 1.518, 1.530, 1.530, 1.530, 1.530, 1.530, 1.470 /) + real, parameter :: xnimag_in(nspint_sw) = (/ 1.165E-2, 1.073E-2, 8.650E-3, 8.092E-3, 8.000E-3, & + 8.000E-3, 8.000E-3, 8.000E-3, 8.000E-3, 8.000E-3, & + 8.000E-3, 8.440E-3, 3.000E-2,9.000E-2 /) + +! FSB 02/11/2012 These values are replaced. +! data xnreal_ec /1.877, 1.832, 1.813, 1.802, 1.791, 1.769, 1.761, & +! 1.760, 1.750, 1.740, 1.750, 1.738, 1.620, 2.120/ +! data xnimag_ec/ 5.563E-1, 5.273E-1, 5.030E-1, 4.918E-1, 4.814E-1, & +! 4.585E-1, 4.508E-1, 4.404E-1, 4.300E-1, 4.400E-1, & +! 4.600E-1, 4.696E-1, 4.500E-1, 5.700E-1/ + +! New Refractive indices for EC at RRTMG Wavelengths +! Source lamda xnreal_ec xnimag_ec +! C&C Values +! 3.4615 2.089 1.070 +! 2.7885 2.014 0.939 +! 2.325 1.962 0.843 +! 2.046 1.950 0.784 +! Bond values +! 1.784 1.940 0.760 +! 1.4625 1.930 0.749 +! 1.2705 1.905 0.737 +! 1.0101 1.870 0.726 +! B&B Values +! 0.7016 1.85 0.71 +! 0.53325 1.85 0.71 +! 0.38815 1.85 0.71 +! 0.299 1.85 0.71 +! 0.2316 1.85 0.71 +! C & C values +! 8.24 2.589 1.771 +!References: +! Bond Personal Communication from Tami Bond +! B&B Bond, T.C. & R.W. Bergstrom (2006) Light absorption by +! Carbonaceous Particles: An investigative review, +! Aerosol Science and Technology. Vol. 40. pp 27-67 +! +! C&C Chang,H and T.T. Charalmpopoulos (1990) Determination of the +! wavelength dependence of refractive indices of flame soot, +! Proceeding of the Royal Society of London A, Vol. 430, pp 577-591. +! FSB new values + +! FSB elemental carbon - soot like + + real, parameter :: xnreal_ec(nspint_sw) = (/ 2.089, 2.014, 1.962, 1.950, 1.940, 1.930, 1.905, & + 1.870, 1.85, 1.85, 1.85, 1.85, 1.85, 2.589 /) + real, parameter :: xnimag_ec(nspint_sw) = (/ 1.070, 0.939, 0.843, 0.784, 0.760, 0.749, 0.737, & + 0.726, 0.71, 0.71, 0.71, 0.71, 0.71, 1.771 /) + +! FSB water + real, save :: xnreal_h2o(nspint_sw) = (/ 1.408, 1.324, 1.277, 1.302, 1.312, 1.321, 1.323, & + 1.327, 1.331, 1.334, 1.340, 1.349, 1.362, 1.260 /) + real, save :: xnimag_h2o(nspint_sw) = (/ 1.420E-2, 1.577E-1, 1.516E-3, 1.159E-3, 2.360E-4, & + 1.713E-4, 2.425E-5, 3.125E-6, 3.405E-8, 1.639E-9, & + 2.955E-9, 1.635E-8, 3.350E-8, 6.220E-2 /) + +! FSB Begin code ====================================================== + + bextsum = 0.0 + bscatsum = 0.0 + bsgsum = 0.0 + do m = 1, nmode +! FSB calculate volumes [ m**3 / m**3 ] +! FSB the reciprocal densities have been scaled to [ m**3 / ug ] + + MVOL_ws(m) = rhows1 * INMASS_ws(m) + MVOL_in(m) = rhoin1 * INMASS_in(m) + MVOL_ec(m) = rhoec1 * INMASS_ec(m) + MVOL_ss(m) = rhoss1 * INMASS_ss(m) + MVOL_h2o(m) = rhoh2o1 * INMASS_h2o(m) + + VOLSHELL(m) = MVOL_ws(m) + MVOL_in(m) + MVOL_ss(m) + MVOL_h2o(m) + VOLCOR(m) = MVOL_ec(m) +! VOL(m) = VOLSHELL(m) + VOLCOR(m) ! VOL is total volume + + if ( VOLCOR(m) .gt. 0.0 ) then +! FSB EC is present +! calculate the ratio of core to shell volume +! take cube root for scaling the diameter of +! the core to that of the shell. + +! dfac = ( VOLCOR(m) / VOL(m) ) ** one3rd + dfac = ( VOLCOR(m) / ( VOLSHELL(m) + VOLCOR(m) ) ) ** one3rd +! dfac = ( VOLCOR(m) / ( VOLSHELL(m) + VOLCOR(m) ) ) +! FSB Set shell and core diameters + DGNSHELL(m) = INDGN(m) + DGNCORE(M) = dfac * INDGN(m) +! FSB note that VOLSHELL(m) is the total volume when EC is not present + end if + +! internal mixture of non-EC species. + +! modal real refractive index No EC + nr(m) = (MVOL_ws(m) * xnreal_ws(ns) + & + MVOL_in(m) * xnreal_in(ns) + & + MVOL_ss(m) * xnreal_ss(ns) + & +! MVOL_h2o(m) * xnreal_h2o(ns)) / VOL(m) + MVOL_h2o(m) * xnreal_h2o(ns)) / VOLSHELL(m) + +! modal imaginary refractive index no EC + ni(m) = (MVOL_ws(m) * xnimag_ws(ns) + & + MVOL_in(m) * xnimag_in(ns) + & + MVOL_ss(m) * xnimag_ss(ns) + & +! MVOL_h2o(m) * xnimag_h2o(ns)) / VOL(m) + MVOL_h2o(m) * xnimag_h2o(ns)) / VOLSHELL(m) + + if ( VOLCOR(m) .gt. 0.0) then + +! FSB calculate the complex refractive indices for the CORE and +! the SHELL for case when and EC core is assumed to exist. + + refcor(m) = cmplx( xnreal_ec(ns), xnimag_ec(ns) ) + refshell(m) = cmplx(nr(m), ni(m) ) +! FSB do BHCOAT case + CALL aero_optical_CS( LAMDA_SW(ns), refcor(m), refshell(m), & + VOLCOR(m),VOLSHELL(m), DGNCORE(m), & + DGNSHELL(m), INSIG(m), & + bext, bscat, gfac, succesS ) +! else if ( VOL(m) .gt. 0.0) then + else if ( VOLSHELL(m) .gt. 0.0) then +! FSB do BHMIE case for the case when EC is not present. + crefin(m) = cmplx(nr(m), ni(m) ) +! CALL aero_optical2( LAMDA_SW(ns), crefin(m), VOL(m), & + CALL aero_optical2( LAMDA_SW(ns), crefin(m), VOLSHELL(m), & + INDGN(m), INSIG(m), & + bext, bscat, gfac, success ) + else + bext = 0.0 + bscat = 0.0 + gfac = 0.0 + end if + +! FSB sum for total values + bextsum = bextsum + bext + bscatsum = bscatsum +bscat + bsgsum = bsgsum + bscat * gfac +! FSB get history +! bext_wm(ns,m) = bext +! bscat_wm(ns,m) = bscat +! gfac_wm(ns,m) = gfac + end do ! loop over modes + +! FSB construct output variables + tauaer = bextsum * delta_z + waer = bscatsum / bextsum + gaer = bsgsum / bscatsum + + end subroutine get_aerosol_Optics_RRTMG_SW + +END MODULE module_twoway_ra_rrtmg_sw From 275c59df20dbd383031e4d538c967d098f85000a Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Sun, 19 Dec 2021 21:21:10 +0000 Subject: [PATCH 02/15] fix compilation error with solve rho --- dyn_em/solve_em.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dyn_em/solve_em.F b/dyn_em/solve_em.F index 47a38f9af1..a84668c3cb 100644 --- a/dyn_em/solve_em.F +++ b/dyn_em/solve_em.F @@ -4710,7 +4710,7 @@ END SUBROUTINE CMAQ_DRIVER CALL after_all_rk_steps ( grid, config_flags, & moist, chem, tracer, scalar, & - th_phy, pi_phy, p_phy, rho_phy, & + th_phy, pi_phy, p_phy, & p8w, t8w, dz8w, & REAL(curr_secs,8), curr_secs2, & diag_flag, & From 3cd5dcb22b0546982aaa6956ff27b07f7f232344 Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Sun, 19 Dec 2021 22:18:50 +0000 Subject: [PATCH 03/15] Fix inconsistent double precision complex functions modified: phys/module_twoway_rrtmg_aero_optical_util.F --- phys/module_twoway_rrtmg_aero_optical_util.F | 21 ++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/phys/module_twoway_rrtmg_aero_optical_util.F b/phys/module_twoway_rrtmg_aero_optical_util.F index e074d34fa5..dbe48cefd8 100644 --- a/phys/module_twoway_rrtmg_aero_optical_util.F +++ b/phys/module_twoway_rrtmg_aero_optical_util.F @@ -1687,7 +1687,11 @@ SUBROUTINE BHCOAT (XX, YY, RRFRL1, RRFRL2, QQEXT, QQSCA, QBACK, GGSCA, SUCCESS) SUCCESS = .TRUE. +#if ( RWORDSIZE == 8 ) + II = c_set(0.0, 1.0) +#else II = c_set(0.0D0, 1.0D0) +#endif ! this technique will make the second 4 byte in the 8 byte variable be 0 ! rather than arbitrary digits to increase accuracy @@ -1729,14 +1733,26 @@ SUBROUTINE BHCOAT (XX, YY, RRFRL1, RRFRL2, QQEXT, QQSCA, QBACK, GGSCA, SUCCESS) xi0y = c_sub(psi0y, c_mul(chi0y, II)) xi1y = c_sub(psi1y, c_mul(chi1y, II)) +#if ( RWORDSIZE == 8 ) + chi0y2 = c_mul(-1.0, c_SIN(y2)) +#else chi0y2 = c_mul(-1.0d0, c_SIN(y2)) +#endif chi1y2 = c_COS(y2) +#if ( RWORDSIZE == 8 ) + chi0x2 = c_mul(-1.0, c_SIN(x2)) +#else chi0x2 = c_mul(-1.0d0, c_SIN(x2)) +#endif chi1x2 = c_COS(x2) qsca = 0.0d0 qext = 0.0d0 GSCA = 0.0d0 +#if ( RWORDSIZE == 8 ) + xback = c_set(0.0, 0.0) +#else xback = c_set(0.0d0, 0.0d0) +#endif iflag = 0 factor = 1.0d0 @@ -1790,8 +1806,13 @@ SUBROUTINE BHCOAT (XX, YY, RRFRL1, RRFRL2, QQEXT, QQSCA, QBACK, GGSCA, SUCCESS) (c_ABS(amess3) .LE. del*c_ABS(d1y2)) .AND. & (c_ABS(amess4) .LE. del) ) THEN ! convergence for inner sphere +#if ( RWORDSIZE == 8 ) + brack = c_set(0.0,0.0) + crack = c_set(0.0,0.0) +#else brack = c_set(0.0D0,0.0D0) crack = c_set(0.0D0,0.0D0) +#endif iflag = 1 ! ELSE ! no convergence yet From 1572dfa61621f25ad4d4f9e36efe13baea2e5898 Mon Sep 17 00:00:00 2001 From: dwongepa Date: Mon, 3 Jan 2022 10:20:59 -0500 Subject: [PATCH 04/15] small tweaks from developer (switched back to use lioapi_temp variable otherwise gcc compiler will fail) --- arch/Config.pl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/arch/Config.pl b/arch/Config.pl index 8bb169ebe8..5b51789782 100644 --- a/arch/Config.pl +++ b/arch/Config.pl @@ -448,13 +448,13 @@ while ( ) { if ( $_ =~ /ifort compiler/ ) - { $lioapi = 'Linux2_x86_64ifort'; + { $lioapi_temp = 'Linux2_x86_64ifort'; } elsif ( $_ =~ /PGI compiler/ ) - { $lioapi = 'Linux2_x86_64pg'; + { $lioapi_temp = 'Linux2_x86_64pg'; } elsif ( $_ =~ /gfortran compiler/ ) - { $lioapi = 'Linux2_x86_64gfort'; + { $lioapi_temp = 'Linux2_x86_64gfort'; } if ( substr( $_, 0, 5 ) eq "#ARCH" && $latchon == 1 ) @@ -842,6 +842,7 @@ printf "Compile for nesting? (1=basic, 2=preset moves, 3=vortex following) [default 1]: " ; } $response = ; + $lioapi = $lioapi_temp; } printf "\n" ; lc $response ; From 6ab6adb2742ddfdc90410b3196d88793ee02f30e Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Wed, 19 Jan 2022 10:09:55 -0700 Subject: [PATCH 05/15] Last minute cleanups 1. Moved a number of cmaq options from the Registry.EM to the CMAQ registry (where they rightly belong). 2. Removed the hardcoded -DWRFIO_NCD_LARGE_FILE_SUPPORT in the standard netcdf makefile. If this is required by CMAQ, then an env variable is required to use the correct part of the if test (that already supports this capability). 3. Added a couple of tests in check_a_mundo: a) If CMAQ coupling and direct_sw_feedback is activated, then RRTMG SW is mandatory b) If NOT CMAQ coupling, then if ( direct_sw_feedback==T) or ( wrf_cmaq_option==1), stop modified: Registry/Registry.EM modified: Registry/registry.WRF-CMAQ-twoway modified: external/io_netcdf/makefile modified: share/module_check_a_mundo.F --- Registry/Registry.EM | 8 -------- Registry/registry.WRF-CMAQ-twoway | 8 ++++++++ external/io_netcdf/makefile | 2 +- share/module_check_a_mundo.F | 29 +++++++++++++++++++++++++++++ 4 files changed, 38 insertions(+), 9 deletions(-) diff --git a/Registry/Registry.EM b/Registry/Registry.EM index ee7a57d806..abfc5030bd 100644 --- a/Registry/Registry.EM +++ b/Registry/Registry.EM @@ -4,14 +4,6 @@ include registry.dimspec include registry.em_shared_collection -# WRF-CMAQ coupled model -rconfig integer wrf_cmaq_option namelist,wrf_cmaq 1 0 -rconfig integer wrf_cmaq_freq namelist,wrf_cmaq 1 1 -rconfig integer met_file_tstep namelist,wrf_cmaq 1 10000 - -rconfig logical direct_sw_feedback namelist,wrf_cmaq 1 .false. -rconfig logical feedback_restart namelist,wrf_cmaq 1 .false. - # added to output 5 for ESMF state real landmask ij misc 1 - i0125rh056d=(interp_fcnm_imask)u=(copy_fcnm) "LANDMASK" "LAND MASK (1 FOR LAND, 0 FOR WATER)" "" state real lakemask ij misc 1 - i012rhd=(interp_fcnm_imask)u=(copy_fcnm) "LAKEMASK" "LAKE MASK (1 FOR LAKE, 0 FOR NON-LAKE)" "" diff --git a/Registry/registry.WRF-CMAQ-twoway b/Registry/registry.WRF-CMAQ-twoway index fcbd101e72..b3aaba1d70 100644 --- a/Registry/registry.WRF-CMAQ-twoway +++ b/Registry/registry.WRF-CMAQ-twoway @@ -1,3 +1,11 @@ +# WRF-CMAQ coupled model +rconfig integer wrf_cmaq_option namelist,wrf_cmaq 1 0 +rconfig integer wrf_cmaq_freq namelist,wrf_cmaq 1 1 +rconfig integer met_file_tstep namelist,wrf_cmaq 1 10000 + +rconfig logical direct_sw_feedback namelist,wrf_cmaq 1 .false. +rconfig logical feedback_restart namelist,wrf_cmaq 1 .false. + #state real mass_ws_i ikj twoway_feedback_data 1 - r "mass_ws_i" "Water soluble i mode" "ug/m**3" #state real mass_ws_j ikj twoway_feedback_data 1 - r "mass_ws_j" "Water soluble j mode" "ug/m**3" #state real mass_ws_k ikj twoway_feedback_data 1 - r "mass_ws_k" "Water soluble k mode" "ug/m**3" diff --git a/external/io_netcdf/makefile b/external/io_netcdf/makefile index 9ac0d94eab..fa638d455f 100644 --- a/external/io_netcdf/makefile +++ b/external/io_netcdf/makefile @@ -29,7 +29,7 @@ wrf_io.o: wrf_io.F90 $(CODE) if [ $$a -a "$$WRFIO_NCD_LARGE_FILE_SUPPORT" = "1" ] ; then \ $(CPP1) -DWRFIO_NCD_LARGE_FILE_SUPPORT -I../ioapi_share wrf_io.F90 | $(M4) - > wrf_io.f ; \ else \ - $(CPP1) -DWRFIO_NCD_LARGE_FILE_SUPPORT -I../ioapi_share wrf_io.F90 | $(M4) - > wrf_io.f ; \ + $(CPP1) -I../ioapi_share wrf_io.F90 | $(M4) - > wrf_io.f ; \ fi $(FC) -o $@ $(FFLAGS) -c wrf_io.f diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index a9ee67bd91..7265350a41 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -625,6 +625,35 @@ END FUNCTION bep_bem_ngr_u # endif +!----------------------------------------------------------------------- +! With CMAQ coupling, if the option "direct_sw_feedback" is activated, +! then the only SW radiation scheme set up to support this is RRTMG. +!----------------------------------------------------------------------- +# if ( WRF_CMAQ == 1 ) + cmaq : DO i = 1, model_config_rec % max_dom + IF ( .NOT. model_config_rec % grid_allowed(i) ) CYCLE + IF ( ( model_config_rec % direct_sw_feedback ) .AND. & + ( model_config_rec % wrf_cmaq_option .EQ. 1 ) .AND. & + ( model_config_rec % ra_sw_physics(i) .NE. rrtmg_swscheme ) ) THEN + wrf_err_message = '--- ERROR: With CMAQ coupling, "direct_sw_feedback=T" requires RRTMG SW' + CALL wrf_message ( wrf_err_message ) + count_fatal_error = count_fatal_error + 1 + EXIT cmaq + END IF + ENDDO cmaq +# else + cmaq : DO i = 1, model_config_rec % max_dom + IF ( .NOT. model_config_rec % grid_allowed(i) ) CYCLE + IF ( ( model_config_rec % direct_sw_feedback ) .OR. & + ( model_config_rec % wrf_cmaq_option .EQ. 1 ) ) THEN + wrf_err_message = '--- ERROR: The option "direct_sw_feedback=T" and "wrf_cmaq_option==1" require CMAQ coupling' + CALL wrf_message ( wrf_err_message ) + count_fatal_error = count_fatal_error + 1 + EXIT cmaq + END IF + ENDDO cmaq +# endif + !----------------------------------------------------------------------- ! Print a warning message for not using a combination of radiation and microphysics from Goddard !----------------------------------------------------------------------- From bf8373cff6f61123740422d15bfa91904e02118b Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Wed, 19 Jan 2022 23:44:03 -0700 Subject: [PATCH 06/15] Get rid of the string "twoway" in CMAQ file names renamed: Registry/registry.WRF-CMAQ-twoway -> Registry/registry.CMAQ modified: Registry/registry.em_shared_collection modified: main/depend.common modified: phys/Makefile renamed: phys/module_twoway_rrtmg_aero_optical_util.F -> phys/module_ra_rrtmg_aero_optical_util_cmaq.F modified: phys/module_ra_rrtmg_sw.F renamed: phys/module_twoway_ra_rrtmg_sw.F -> phys/module_ra_rrtmg_sw_cmaq.F modified: phys/module_radiation_driver.F --- ...{registry.WRF-CMAQ-twoway => registry.CMAQ} | 0 Registry/registry.em_shared_collection | 2 +- main/depend.common | 6 +++--- phys/Makefile | 3 ++- ...> module_ra_rrtmg_aero_optical_util_cmaq.F} | 4 ++-- phys/module_ra_rrtmg_sw.F | 18 +++++++++--------- ...ra_rrtmg_sw.F => module_ra_rrtmg_sw_cmaq.F} | 4 ++-- phys/module_radiation_driver.F | 12 ++++++------ 8 files changed, 25 insertions(+), 24 deletions(-) rename Registry/{registry.WRF-CMAQ-twoway => registry.CMAQ} (100%) rename phys/{module_twoway_rrtmg_aero_optical_util.F => module_ra_rrtmg_aero_optical_util_cmaq.F} (99%) rename phys/{module_twoway_ra_rrtmg_sw.F => module_ra_rrtmg_sw_cmaq.F} (99%) diff --git a/Registry/registry.WRF-CMAQ-twoway b/Registry/registry.CMAQ similarity index 100% rename from Registry/registry.WRF-CMAQ-twoway rename to Registry/registry.CMAQ diff --git a/Registry/registry.em_shared_collection b/Registry/registry.em_shared_collection index c58105e805..484514c482 100644 --- a/Registry/registry.em_shared_collection +++ b/Registry/registry.em_shared_collection @@ -28,4 +28,4 @@ include registry.new3d_wif include registry.trad_fields include registry.solar_fields include registry.diags -include registry.WRF-CMAQ-twoway +include registry.CMAQ diff --git a/main/depend.common b/main/depend.common index e7eec48c15..ddd5943722 100644 --- a/main/depend.common +++ b/main/depend.common @@ -487,9 +487,9 @@ module_sf_ruclsm.o: ../frame/module_wrf_error.o module_data_gocart_dust.o module_sf_pxlsm.o: ../share/module_model_constants.o module_sf_pxlsm_data.o -module_twoway_ra_rrtmg_sw.o: module_twoway_rrtmg_aero_optical_util.o +module_ra_rrtmg_sw_cmaq.o: module_ra_rrtmg_aero_optical_util_cmaq.o -module_ra_rrtmg_sw.o: module_ra_rrtmg_lw.o +module_ra_rrtmg_sw.o: module_ra_rrtmg_sw_cmaq.o module_ra_rrtmg_lw.o module_ra_rrtmg_swf.o: module_ra_rrtmg_lwf.o module_ra_rrtmg_swk.o: module_ra_rrtmg_lwk.o module_ra_effective_radius.o @@ -705,7 +705,7 @@ module_radiation_driver.o: \ module_ra_rrtm.o \ module_ra_rrtmg_lw.o \ module_ra_rrtmg_sw.o \ - module_twoway_ra_rrtmg_sw.o \ + module_ra_rrtmg_sw_cmaq.o \ module_ra_rrtmg_lwf.o \ module_ra_rrtmg_swf.o \ module_ra_rrtmg_lwk.o \ diff --git a/phys/Makefile b/phys/Makefile index 3bdf54d7ff..1334011c88 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -9,7 +9,6 @@ MODULES = \ module_bep_bem_helper.o \ complex_number_module.o \ module_cam_shr_kind_mod.o \ - module_twoway_rrtmg_aero_optical_util.o \ module_cam_support.o \ module_cam_shr_const_mod.o \ module_cam_physconst.o \ @@ -118,6 +117,8 @@ MODULES = \ module_ra_rrtm.o \ module_ra_rrtmg_lw.o \ module_ra_rrtmg_sw.o \ + module_ra_rrtmg_sw_cmaq.o \ + module_ra_rrtmg_aero_optical_util_cmaq.o \ module_ra_rrtmg_lwf.o \ module_ra_rrtmg_swf.o \ module_ra_rrtmg_lwk.o \ diff --git a/phys/module_twoway_rrtmg_aero_optical_util.F b/phys/module_ra_rrtmg_aero_optical_util_cmaq.F similarity index 99% rename from phys/module_twoway_rrtmg_aero_optical_util.F rename to phys/module_ra_rrtmg_aero_optical_util_cmaq.F index dbe48cefd8..7ec0864a4e 100644 --- a/phys/module_twoway_rrtmg_aero_optical_util.F +++ b/phys/module_ra_rrtmg_aero_optical_util_cmaq.F @@ -2471,7 +2471,7 @@ END MODULE rrtmg_aero_optical_util_module ! ------------------------------------------------------------------ ! FSB REvised Mie calculations 02/09/2011 -MODULE module_twoway_ra_rrtmg_sw +MODULE module_ra_rrtmg_sw_cmaq contains @@ -2786,4 +2786,4 @@ Subroutine get_aerosol_Optics_RRTMG_SW ( ns, nmode,delta_z, INMASS_ws, & end subroutine get_aerosol_Optics_RRTMG_SW -END MODULE module_twoway_ra_rrtmg_sw +END MODULE module_ra_rrtmg_sw_cmaq diff --git a/phys/module_ra_rrtmg_sw.F b/phys/module_ra_rrtmg_sw.F index 0bf899c077..d475a3f84d 100644 --- a/phys/module_ra_rrtmg_sw.F +++ b/phys/module_ra_rrtmg_sw.F @@ -8855,9 +8855,9 @@ subroutine rrtmg_sw & sibvisdir, sibvisdif, sibnirdir, sibnirdif, & ! ---------------------- End, Zhenxin 2011-06-20 --------------------------------! swdkdir,swdkdif, & ! jararias, 2013/08/10 - swdkdirc & ! PAJ - ,calc_clean_atm_diag & - ,sw_zbbcddir, sw_dirdflux, sw_difdflux & ! WRF-CMAQ twoway coupled model + swdkdirc, & ! PAJ + calc_clean_atm_diag, & + sw_zbbcddir, sw_dirdflux, sw_difdflux & ! WRF-CMAQ twoway coupled model ) @@ -10012,7 +10012,7 @@ MODULE module_ra_rrtmg_sw use module_ra_rrtmg_lw, only : inirad, o3data, relcalc, reicalc, retab ! mcica_random_numbers, randomNumberSequence, & ! new_RandomNumberSequence, getRandomReal -use module_twoway_ra_rrtmg_sw +use module_ra_rrtmg_sw_cmaq CONTAINS @@ -10065,7 +10065,7 @@ SUBROUTINE RRTMG_SWRAD( & swdownc, swddnic, swddirc, & ! PAJ xcoszen,yr,julian, & ! jararias 2013/08 obscur, & ! amontornes-bcodina 2015/09 solar eclipses - proceed_twoway_sw, & ! WRF-CMAQ twoway coupled model + proceed_cmaq_sw, & ! WRF-CMAQ twoway coupled model nmode, & ! WRF-CMAQ twoway coupled model mass_ws_i, mass_ws_j, mass_ws_k, & ! WRF-CMAQ twoway coupled model mass_in_i, mass_in_j, mass_in_k, & ! WRF-CMAQ twoway coupled model @@ -10254,7 +10254,7 @@ SUBROUTINE RRTMG_SWRAD( & REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: obscur ! begin WRF-CMAQ twoway coupled model block - LOGICAL, INTENT(IN) :: proceed_twoway_sw + LOGICAL, INTENT(IN) :: proceed_cmaq_sw ! ** FSB items needed for new aerosol code from CMAQ integer, optional, intent(in) :: nmode ! number of log-normal modes @@ -11239,7 +11239,7 @@ SUBROUTINE RRTMG_SWRAD( & ssaaer(ncol,k,nb) = 1. asmaer(ncol,k,nb) = 0. - if (proceed_twoway_sw) then + if (proceed_cmaq_sw) then INMASS_ws(1) = mass_ws_i(i,k,j) INMASS_ws(2) = mass_ws_j(i,k,j) INMASS_ws(3) = mass_ws_k(i,k,j) @@ -11309,7 +11309,7 @@ SUBROUTINE RRTMG_SWRAD( & asmaer(ncol,k,nb) = gaer end if enddo ! loop over layers - if (proceed_twoway_sw) then + if (proceed_cmaq_sw) then ! No aerosols in top layer above model top (kte+1). tauaer(ncol, kte+1 ,nb) = 0. ssaaer(ncol, kte+1 ,nb) = 1. @@ -11534,7 +11534,7 @@ SUBROUTINE RRTMG_SWRAD( & enddo else - if (proceed_twoway_sw) then ! this is for WRF-CMAQ twoway coupled model + if (proceed_cmaq_sw) then ! this is for WRF-CMAQ twoway coupled model gtauxar_01 (i,:,j) = 0.0 gtauxar_02 (i,:,j) = 0.0 gtauxar_03 (i,:,j) = 0.0 diff --git a/phys/module_twoway_ra_rrtmg_sw.F b/phys/module_ra_rrtmg_sw_cmaq.F similarity index 99% rename from phys/module_twoway_ra_rrtmg_sw.F rename to phys/module_ra_rrtmg_sw_cmaq.F index 3a397e49a8..5e167bfaf1 100644 --- a/phys/module_twoway_ra_rrtmg_sw.F +++ b/phys/module_ra_rrtmg_sw_cmaq.F @@ -1,4 +1,4 @@ -MODULE module_twoway_ra_rrtmg_sw +MODULE module_ra_rrtmg_sw_cmaq contains @@ -324,4 +324,4 @@ Subroutine get_aerosol_Optics_RRTMG_SW ( ns, nmode,delta_z, INMASS_ws, & end subroutine get_aerosol_Optics_RRTMG_SW -END MODULE module_twoway_ra_rrtmg_sw +END MODULE module_ra_rrtmg_sw_cmaq diff --git a/phys/module_radiation_driver.F b/phys/module_radiation_driver.F index 3187b9b09e..23b4b3116d 100644 --- a/phys/module_radiation_driver.F +++ b/phys/module_radiation_driver.F @@ -269,7 +269,7 @@ SUBROUTINE radiation_driver ( & ! *** add new modules of schemes here #if ( WRF_CMAQ == 1) - USE module_twoway_ra_rrtmg_sw + USE module_ra_rrtmg_sw_cmaq #endif USE module_ra_sw , ONLY : swrad @@ -994,7 +994,7 @@ SUBROUTINE radiation_driver ( & !---------- local test vars ! real, dimension(ims:ime, kms:kme, jms:jme) :: hlw, hsw - LOGICAL :: proceed_twoway_sw + LOGICAL :: proceed_cmaq_sw logical, save :: firstime = .true. logical, save :: feedback_restart, direct_sw_feedback @@ -2499,7 +2499,7 @@ SUBROUTINE radiation_driver ( & CALL wrf_debug(100, 'CALL rrtmg_sw') if ( direct_sw_feedback .and. feedback_is_ready ) then - proceed_twoway_sw = .true. + proceed_cmaq_sw = .true. if (present(mass_ws_i)) then @@ -2532,7 +2532,7 @@ SUBROUTINE radiation_driver ( & sig_k(:,kte+1,:) = sig_k(:,kte,:) end if else - proceed_twoway_sw = .false. + proceed_cmaq_sw = .false. end if CALL RRTMG_SWRAD( & @@ -2598,7 +2598,7 @@ SUBROUTINE radiation_driver ( & swdownc=swdownc, swddnic=swddnic, swddirc=swddirc, & ! PAJ obscur=obscur_loc, & xcoszen=coszen,yr=yr,julian=julian,mp_physics=mp_physics, & ! jararias 2013/08/14 - proceed_twoway_sw=proceed_twoway_sw, & ! WRF-CMAQ coupled model + proceed_cmaq_sw=proceed_cmaq_sw, & ! WRF-CMAQ coupled model nmode=3, & ! WRF-CMAQ coupled model mass_ws_i=mass_ws_i, & ! WRF-CMAQ coupled model mass_ws_j=mass_ws_j, & ! WRF-CMAQ coupled model @@ -2642,7 +2642,7 @@ SUBROUTINE radiation_driver ( & ) ! = WRF-CMAQ twoway coupled model - if (proceed_twoway_sw) then + if (proceed_cmaq_sw) then DO j=jts,jte DO i=its,ite sw_ttauxar_01(i,j) = sum(sw_gtauxar_01(i,:,j)) From 077bac6a5bbfd94ee141d534ced9ad2a7e5e01b8 Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Mon, 24 Jan 2022 15:46:21 -0700 Subject: [PATCH 07/15] Zap dupe file module_ra_rrtmg_sw_cmaq.F The contents of this file are completely contained within module_ra_rrtmg_aero_optical_util_cmaq.F (also in the phys dir). THis change requires modifying some physics Makefiles and dependency files. M main/depend.common M phys/Makefile D phys/module_ra_rrtmg_sw_cmaq.F --- main/depend.common | 6 +- phys/Makefile | 1 - phys/module_ra_rrtmg_sw_cmaq.F | 327 --------------------------------- 3 files changed, 2 insertions(+), 332 deletions(-) delete mode 100644 phys/module_ra_rrtmg_sw_cmaq.F diff --git a/main/depend.common b/main/depend.common index d9e9d726ff..60feaa754b 100644 --- a/main/depend.common +++ b/main/depend.common @@ -487,9 +487,7 @@ module_sf_ruclsm.o: ../frame/module_wrf_error.o module_data_gocart_dust.o module_sf_pxlsm.o: ../share/module_model_constants.o module_sf_pxlsm_data.o -module_ra_rrtmg_sw_cmaq.o: module_ra_rrtmg_aero_optical_util_cmaq.o - -module_ra_rrtmg_sw.o: module_ra_rrtmg_sw_cmaq.o module_ra_rrtmg_lw.o +module_ra_rrtmg_sw.o: module_ra_rrtmg_aero_optical_util_cmaq.o module_ra_rrtmg_lw.o module_ra_rrtmg_swf.o: module_ra_rrtmg_lwf.o module_ra_rrtmg_swk.o: module_ra_rrtmg_lwk.o module_ra_effective_radius.o @@ -706,7 +704,7 @@ module_radiation_driver.o: \ module_ra_rrtm.o \ module_ra_rrtmg_lw.o \ module_ra_rrtmg_sw.o \ - module_ra_rrtmg_sw_cmaq.o \ + module_ra_rrtmg_aero_optical_util_cmaq.o \ module_ra_rrtmg_lwf.o \ module_ra_rrtmg_swf.o \ module_ra_rrtmg_lwk.o \ diff --git a/phys/Makefile b/phys/Makefile index 86a9ae6f10..b91e1bcbab 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -117,7 +117,6 @@ MODULES = \ module_ra_rrtm.o \ module_ra_rrtmg_lw.o \ module_ra_rrtmg_sw.o \ - module_ra_rrtmg_sw_cmaq.o \ module_ra_rrtmg_aero_optical_util_cmaq.o \ module_ra_rrtmg_lwf.o \ module_ra_rrtmg_swf.o \ diff --git a/phys/module_ra_rrtmg_sw_cmaq.F b/phys/module_ra_rrtmg_sw_cmaq.F deleted file mode 100644 index 5e167bfaf1..0000000000 --- a/phys/module_ra_rrtmg_sw_cmaq.F +++ /dev/null @@ -1,327 +0,0 @@ -MODULE module_ra_rrtmg_sw_cmaq - -contains - -!------------------------------------------------------------------ - Subroutine get_aerosol_Optics_RRTMG_SW ( ns, nmode,delta_z, INMASS_ws, & - INMASS_in, INMASS_ec, INMASS_ss, & - INMASS_h2o, INDGN, INSIG, & - tauaer, waer, gaer ) - -!FSB This version switches between BHCOAT to BHMIE depending upon whether -! EC is present or not. 04/15/2012. - -!FSB this version does a core-shell calculation with BHCOAT 04/11/2012 -! This version is set up to be used with RRTMG_SW <<<<<<<< -! wavelenght is calculated internally -! FSB This routine calculates the aerosol information ( tauaer, waer, -! gaer, needed to calculate the solar radiation) The calling -! program specifies the location ( row, column, layer, -! layer thicknes, and wave length for the calculation. -! FSB 02/09/2011 Modifications made to subroutine ghintBH. -! FSB 04/14/2012 REmoved MODULUS, made changes to ghintBH. -! Put in option for core-shell (coated-sphere). 2 - -! FSB Input variables: - - use rrtmg_aero_optical_util_module - - implicit none - - integer,intent(in) :: ns ! index for wavelength should be - ! between 1 and 14. <<< RRTMG_SW - integer,intent(in) :: nmode ! should be 3 for WRF/CMAQ calculation - real,intent(in) :: delta_z ! layer thickness [m] -! FSB mode types for WRF/CMAQ -! nmode = 1 Aitken -! nmode = 2 accumulation -! nmode = 3 coarse -! FSB modal mass concentration by species [ ug / m**3] NOTE: MKS - real, intent(in) :: INMASS_ws(nmode) ! water soluble - real, intent(in) :: INMASS_in(nmode) ! insolugle - real, intent(in) :: INMASS_ec(nmode) ! elemental carbon or soot like - real, intent(in) :: INMASS_ss(nmode) ! sea salt - real, intent(in) :: INMASS_h2o(nmode) ! water -! FSB particle size-distribution information - real, intent(in) :: INDGN( nmode) ! geometric mean diameter [ m ] NOTE: MKS - real, intent(in) :: INSIG( nmode) ! geometric standard deviation - -!FSB output aerosol radiative properties [dimensionless] - real, intent(out) :: tauaer ! aerosol extinction optical depth - real, intent(out) :: waer ! aerosol single scattering albedo - real, intent(out) :: gaer ! aerosol assymetry parameter - -! FSB Internal variables - - real :: NR(nmode), NI(nmode) ! refractive indices - complex :: refcor(nmode), refshell(nmode) ! complex refracive indices - complex :: crefin(nmode) ! complex refractive index - -! FSB special values for EC CORE-shell calculation - real :: DGNSHELL(nmode) ! modal geometric mean diameter [m] - real :: DGNCORE (nmode) ! modal geometric mean diameter [m] - -! FSB Modal volumes [ m**3 / m**3 ] - real :: MVOL_ws(nmode) ! water soluble - real :: MVOL_in(nmode) ! insolugle - real :: MVOL_ec(nmode) ! soot like - real :: MVOL_ss(nmode) !sea salt - real :: MVOL_h2o(nmode) ! water -! real :: VOL(nmode) ! total modal volume [m** 3 / m**3] -! FSB special values for EC CORE-shell calculation - real :: VOLCOR(nmode) ! volume of EC core [m** 3 / m**3] - real :: VOLSHELL(nmode) ! volume of shell [m** 3 / m**3] - - integer :: m ! loop index - real :: bext ! extinction coefficient [1 / m] - real :: bscat ! scattering coefficient [1 / m] - real :: gfac ! asymmetry factor - - real :: bextsum, bscatsum, bsgsum - -! FSB History variables by wavelength and mode -! real :: bext_wm(ns,nmode) -! real :: bscat_wm(ns,nmode) -! real :: gfac_wm(ns,nmode) - - real, parameter :: one3rd = 1.0 / 3.0 - real :: dfac ! ratio of (volcor/vol) ** one3rd - ! used for calculating the diameter - ! of the EC core - - logical :: succesS - -!...component densities [ g/ cm**3 ] <<<<< cgs - - real, parameter :: rhows = 1.8 ! bulk density of water soluble aerosol - - real, parameter :: rhoin = 2.2 ! bulk density forinsoluble aerosol - -! real, parameter :: rhoec = 1.7 ! bulk density for soot aerosol - real, parameter :: rhoec = 1.8 ! new value - - real, parameter :: rhoh2o = 1.0 ! bulk density of aerosol water - - real, parameter :: rhoss = 2.2 ! bulk density of seasalt - -! FSB scale factor for volume calculation -! 1.0d-12 * [ cm**3 / g] -> [ m** 3 / ug ] - real, parameter :: scalefactor = 1.0e-12 - -! FSB scale factor for [1/g] to [1/ug] - real, parameter :: cug2g = 1.0e-06 - -! FSB reciprocal component densities[ m ** 3 / ug ] - - real, parameter :: rhows1 = scalefactor / rhows ! water soluble aerosol - - real, parameter :: rhoin1 = scalefactor / rhoin ! insoluble aerosol - - real, parameter :: rhoec1 = scalefactor / rhoec ! soot aerosol - - real, parameter :: rhoh2o1 = scalefactor / rhoh2o ! aerosol water - - real, parameter :: rhoss1 = scalefactor / rhoss ! seasalt - - integer,parameter :: nspint_sw = 14 ! number of spectral intervals for RRTMG_SW - -! FSB Band numbers and wavelengths for RRTMG_SW - integer, parameter :: Band(nspint_sw) = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 /) - - real, parameter :: LAMDA_SW(nspint_sw) = (/ 3.4615, 2.7885, 2.325, 2.046, 1.784, & - 1.4625, 1.2705, 1.0101, 0.7016, 0.53325, & - 0.38815, 0.299, 0.2316, 8.24 /) ! wavelength [ um ] - -! *** refractive indices - -! *** Except as otherwise noted reference values of refractive -! indices for aerosol particles are from the OPAC Data base. -! Hess, Koepke, and Schult, Optical properties of -! aerosols and clouds: The software package OPAC, Bulletan of -! the American Meteorological Society, Vol 79, No 5, -! pp 831 - 844, May 1998. -! OPAC is a downloadable data set of optical properties of -! 10 aerosol components, 6 water clouds and 3 cirrus clouds -! at UV, visible and IR wavelengths -! www.lrz-muenchen.de/~uh234an/www/radaer/opac.htm - - -! FSB water soluble - real, parameter :: xnreal_ws(nspint_sw) = (/ 1.443, 1.420, 1.420, 1.420, 1.463, 1.510, 1.510, & - 1.520, 1.530, 1.530, 1.530, 1.530, 1.530, 1.710 /) - real, parameter :: xnimag_ws(nspint_sw) = (/ 5.718E-3, 1.777E-2, 1.060E-2, 8.368E-3, 1.621E-2, & - 2.198E-2, 1.929E-2, 1.564E-2, 7.000E-3, 5.666E-3, & - 5.000E-3, 8.440E-3, 3.000E-2, 1.100E-1 /) - -! FSB sea salt - real, parameter :: xnreal_ss(nspint_sw) = (/ 1.480, 1.534, 1.437, 1.448, 1.450, 1.462, 1.469, & - 1.470, 1.490, 1.500, 1.502, 1.510, 1.510, 1.510 /) - - real, parameter :: xnimag_ss(nspint_sw) = (/ 1.758E-3, 7.462E-3, 2.950E-3, 1.276E-3, 7.944E-4, & - 5.382E-4, 3.754E-4, 1.498E-4, 2.050E-7, 1.184E-8, & - 9.938E-8, 2.060E-6, 5.000E-6, 1.000E-2 /) - -! FSB insoluble - real, parameter :: xnreal_in(nspint_sw) = (/ 1.272, 1.168, 1.208, 1.253, 1.329, 1.418, 1.456, & - 1.518, 1.530, 1.530, 1.530, 1.530, 1.530, 1.470 /) - real, parameter :: xnimag_in(nspint_sw) = (/ 1.165E-2, 1.073E-2, 8.650E-3, 8.092E-3, 8.000E-3, & - 8.000E-3, 8.000E-3, 8.000E-3, 8.000E-3, 8.000E-3, & - 8.000E-3, 8.440E-3, 3.000E-2,9.000E-2 /) - -! FSB 02/11/2012 These values are replaced. -! data xnreal_ec /1.877, 1.832, 1.813, 1.802, 1.791, 1.769, 1.761, & -! 1.760, 1.750, 1.740, 1.750, 1.738, 1.620, 2.120/ -! data xnimag_ec/ 5.563E-1, 5.273E-1, 5.030E-1, 4.918E-1, 4.814E-1, & -! 4.585E-1, 4.508E-1, 4.404E-1, 4.300E-1, 4.400E-1, & -! 4.600E-1, 4.696E-1, 4.500E-1, 5.700E-1/ - -! New Refractive indices for EC at RRTMG Wavelengths -! Source lamda xnreal_ec xnimag_ec -! C&C Values -! 3.4615 2.089 1.070 -! 2.7885 2.014 0.939 -! 2.325 1.962 0.843 -! 2.046 1.950 0.784 -! Bond values -! 1.784 1.940 0.760 -! 1.4625 1.930 0.749 -! 1.2705 1.905 0.737 -! 1.0101 1.870 0.726 -! B&B Values -! 0.7016 1.85 0.71 -! 0.53325 1.85 0.71 -! 0.38815 1.85 0.71 -! 0.299 1.85 0.71 -! 0.2316 1.85 0.71 -! C & C values -! 8.24 2.589 1.771 -!References: -! Bond Personal Communication from Tami Bond -! B&B Bond, T.C. & R.W. Bergstrom (2006) Light absorption by -! Carbonaceous Particles: An investigative review, -! Aerosol Science and Technology. Vol. 40. pp 27-67 -! -! C&C Chang,H and T.T. Charalmpopoulos (1990) Determination of the -! wavelength dependence of refractive indices of flame soot, -! Proceeding of the Royal Society of London A, Vol. 430, pp 577-591. -! FSB new values - -! FSB elemental carbon - soot like - - real, parameter :: xnreal_ec(nspint_sw) = (/ 2.089, 2.014, 1.962, 1.950, 1.940, 1.930, 1.905, & - 1.870, 1.85, 1.85, 1.85, 1.85, 1.85, 2.589 /) - real, parameter :: xnimag_ec(nspint_sw) = (/ 1.070, 0.939, 0.843, 0.784, 0.760, 0.749, 0.737, & - 0.726, 0.71, 0.71, 0.71, 0.71, 0.71, 1.771 /) - -! FSB water - real, save :: xnreal_h2o(nspint_sw) = (/ 1.408, 1.324, 1.277, 1.302, 1.312, 1.321, 1.323, & - 1.327, 1.331, 1.334, 1.340, 1.349, 1.362, 1.260 /) - real, save :: xnimag_h2o(nspint_sw) = (/ 1.420E-2, 1.577E-1, 1.516E-3, 1.159E-3, 2.360E-4, & - 1.713E-4, 2.425E-5, 3.125E-6, 3.405E-8, 1.639E-9, & - 2.955E-9, 1.635E-8, 3.350E-8, 6.220E-2 /) - - -! FSB Begin code ====================================================== - - bextsum = 0.0 - bscatsum = 0.0 - bsgsum = 0.0 - do m = 1, nmode -! FSB calculate volumes [ m**3 / m**3 ] -! FSB the reciprocal densities have been scaled to [ m**3 / ug ] - - MVOL_ws(m) = rhows1 * INMASS_ws(m) - MVOL_in(m) = rhoin1 * INMASS_in(m) - MVOL_ec(m) = rhoec1 * INMASS_ec(m) - MVOL_ss(m) = rhoss1 * INMASS_ss(m) - MVOL_h2o(m) = rhoh2o1 * INMASS_h2o(m) - - VOLSHELL(m) = MVOL_ws(m) + MVOL_in(m) + MVOL_ss(m) + MVOL_h2o(m) - VOLCOR(m) = MVOL_ec(m) -! VOL(m) = VOLSHELL(m) + VOLCOR(m) ! VOL is total volume - - if ( VOLCOR(m) .gt. 0.0 ) then -! FSB EC is present -! calculate the ratio of core to shell volume -! take cube root for scaling the diameter of -! the core to that of the shell. - -! dfac = ( VOLCOR(m) / VOL(m) ) ** one3rd - dfac = ( VOLCOR(m) / ( VOLSHELL(m) + VOLCOR(m) ) ) ** one3rd -! dfac = ( VOLCOR(m) / ( VOLSHELL(m) + VOLCOR(m) ) ) -! FSB Set shell and core diameters - DGNSHELL(m) = INDGN(m) - DGNCORE(M) = dfac * INDGN(m) -! FSB note that VOLSHELL(m) is the total volume when EC is not present - end if - -! internal mixture of non-EC species. - -! modal real refractive index No EC - nr(m) = (MVOL_ws(m) * xnreal_ws(ns) + & - MVOL_in(m) * xnreal_in(ns) + & - MVOL_ss(m) * xnreal_ss(ns) + & -! MVOL_h2o(m) * xnreal_h2o(ns)) / VOL(m) - MVOL_h2o(m) * xnreal_h2o(ns)) / VOLSHELL(m) - -! modal imaginary refractive index no EC - ni(m) = (MVOL_ws(m) * xnimag_ws(ns) + & - MVOL_in(m) * xnimag_in(ns) + & - MVOL_ss(m) * xnimag_ss(ns) + & -! MVOL_h2o(m) * xnimag_h2o(ns)) / VOL(m) - MVOL_h2o(m) * xnimag_h2o(ns)) / VOLSHELL(m) - - if ( VOLCOR(m) .gt. 0.0) then - -! FSB calculate the complex refractive indices for the CORE and -! the SHELL for case when and EC core is assumed to exist. - - refcor(m) = cmplx( xnreal_ec(ns), xnimag_ec(ns) ) - refshell(m) = cmplx(nr(m), ni(m) ) -! FSB do BHCOAT case - CALL aero_optical_CS( LAMDA_SW(ns), refcor(m), refshell(m), & - VOLCOR(m),VOLSHELL(m), DGNCORE(m), & - DGNSHELL(m), INSIG(m), & - bext, bscat, gfac, succesS ) -! else if ( VOL(m) .gt. 0.0) then - else if ( VOLSHELL(m) .gt. 0.0) then -! FSB do BHMIE case for the case when EC is not present. - crefin(m) = cmplx(nr(m), ni(m) ) -! CALL aero_optical2( LAMDA_SW(ns), crefin(m), VOL(m), & - CALL aero_optical2( LAMDA_SW(ns), crefin(m), VOLSHELL(m), & - INDGN(m), INSIG(m), & - bext, bscat, gfac, success ) - else - bext = 0.0 - bscat = 0.0 - gfac = 0.0 - end if - -! FSB sum for total values - bextsum = bextsum + bext - bscatsum = bscatsum +bscat - bsgsum = bsgsum + bscat * gfac -! FSB get history -! bext_wm(ns,m) = bext -! bscat_wm(ns,m) = bscat -! gfac_wm(ns,m) = gfac - end do ! loop over modes - -! FSB construct output variables - tauaer = bextsum * delta_z - waer = bscatsum / bextsum - gaer = bsgsum / bscatsum - -! Write out modal values. - -! write(100,*) ' lamda mode bext bscat gaer ' -! m=1 -! write(100,*) lamda(ns), m, bext_wm(ns,m),bscat_wm(ns,m), gfac_wm(ns,m) -! m=2 -! write(100,*) lamda(ns), m, bext_wm(ns,m),bscat_wm(ns,m), gfac_wm(ns,m) -! m=3 -! write(100,*) lamda(ns), m, bext_wm(ns,m),bscat_wm(ns,m), gfac_wm(ns,m) - - end subroutine get_aerosol_Optics_RRTMG_SW - -END MODULE module_ra_rrtmg_sw_cmaq From ed69f2af8ad87f77aa5315cabb950222ac51172d Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Mon, 24 Jan 2022 17:15:07 -0700 Subject: [PATCH 08/15] CMAQ fix registry name 17:15 --- arch/Config.pl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/arch/Config.pl b/arch/Config.pl index e08e8e0e06..f1803e8a72 100644 --- a/arch/Config.pl +++ b/arch/Config.pl @@ -479,8 +479,8 @@ } close (FH); - # determine whether declarations in Registry/registry.WRF-CMAQ-twoway is commented out or not - my $file = "Registry/registry.WRF-CMAQ-twoway"; + # determine whether declarations in Registry/registry.CMAQ is commented out or not + my $file = "Registry/registry.CMAQ"; open (FH, $file) or die("File $file not found"); $registry_wo_cmaq = 0; From f9d165cee403348f813ea7dbf6c4fc937db18e8c Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Mon, 24 Jan 2022 18:49:06 -0700 Subject: [PATCH 09/15] TLADJ mods for CMAQ - arg from solve to rk part1 --- wrftladj/module_first_rk_step_part1_ad.F | 2 ++ wrftladj/module_first_rk_step_part1_tl.F | 2 ++ wrftladj/solve_em_ad.F | 4 ++++ wrftladj/solve_em_tl.F | 4 ++++ 4 files changed, 12 insertions(+) diff --git a/wrftladj/module_first_rk_step_part1_ad.F b/wrftladj/module_first_rk_step_part1_ad.F index 197bbf078b..2051a41d9c 100644 --- a/wrftladj/module_first_rk_step_part1_ad.F +++ b/wrftladj/module_first_rk_step_part1_ad.F @@ -32,6 +32,7 @@ SUBROUTINE a_first_rk_step_part1 ( grid , config_flags & , ipsy,ipey,jpsy,jpey,kpsy,kpey & , k_start , k_end & , f_flux & + , feedback_is_ready & ) USE module_state_description USE module_model_constants @@ -133,6 +134,7 @@ SUBROUTINE a_first_rk_step_part1 ( grid , config_flags & INTEGER, INTENT(IN) :: k_start, k_end LOGICAL, INTENT(IN), OPTIONAL :: f_flux + LOGICAL, INTENT(IN) :: feedback_is_ready ! For WRF-CMAQ coupled model, indicates feedback information is available ! Local REAL, DIMENSION( ims:ime, jms:jme ) :: exch_temf ! 1/7/09 WA diff --git a/wrftladj/module_first_rk_step_part1_tl.F b/wrftladj/module_first_rk_step_part1_tl.F index 6590fe24ca..3c290c889b 100644 --- a/wrftladj/module_first_rk_step_part1_tl.F +++ b/wrftladj/module_first_rk_step_part1_tl.F @@ -32,6 +32,7 @@ SUBROUTINE g_first_rk_step_part1 ( grid , config_flags & , ipsy,ipey,jpsy,jpey,kpsy,kpey & , k_start , k_end & , f_flux & + , feedback_is_ready & ) USE module_state_description USE module_model_constants @@ -134,6 +135,7 @@ SUBROUTINE g_first_rk_step_part1 ( grid , config_flags & INTEGER, INTENT(IN) :: k_start, k_end LOGICAL, INTENT(IN), OPTIONAL :: f_flux + LOGICAL, INTENT(IN) :: feedback_is_ready ! For WRF-CMAQ coupled model, indicates feedback information is available ! Local REAL, DIMENSION( ims:ime, jms:jme ) :: exch_temf ! 1/7/09 WA diff --git a/wrftladj/solve_em_ad.F b/wrftladj/solve_em_ad.F index 5b88f217d2..4ae309e574 100644 --- a/wrftladj/solve_em_ad.F +++ b/wrftladj/solve_em_ad.F @@ -138,6 +138,8 @@ SUBROUTINE solve_em_ad ( grid , config_flags & REAL :: t_new + LOGICAL :: feedback_is_ready ! CMAQ + ! Changes in tendency at this timestep real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,grid%sm33:grid%em33) :: h_tendency, & z_tendency @@ -251,6 +253,7 @@ SUBROUTINE solve_em_ad ( grid , config_flags & ! Initialize linkedlist !CALL linkedlist_initialize + feedback_is_ready = .false. ! set runge-kutta solver (2nd or 3rd order) dynamics_option = config_flags%rk_ord @@ -686,6 +689,7 @@ SUBROUTINE solve_em_ad ( grid , config_flags & , f_flux & , aerocu & , restart_flag & + , feedback_is_ready & ) #ifdef DM_PARALLEL diff --git a/wrftladj/solve_em_tl.F b/wrftladj/solve_em_tl.F index 75790e024a..957b46789c 100644 --- a/wrftladj/solve_em_tl.F +++ b/wrftladj/solve_em_tl.F @@ -115,6 +115,8 @@ SUBROUTINE solve_em_tl ( grid , config_flags & REAL :: t_new + LOGICAL :: feedback_is_ready ! CMAQ + ! Changes in tendency at this timestep real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,grid%sm33:grid%em33) :: h_tendency, & z_tendency @@ -219,6 +221,7 @@ SUBROUTINE solve_em_tl ( grid , config_flags & ! Initialize timers if compiled with -DBENCH #include "bench_solve_em_init.h" + feedback_is_ready = .false. ! set runge-kutta solver (2nd or 3rd order) dynamics_option = config_flags%rk_ord @@ -662,6 +665,7 @@ SUBROUTINE solve_em_tl ( grid , config_flags & , ipsy,ipey,jpsy,jpey,kpsy,kpey & , k_start , k_end & , f_flux & + , feedback_is_ready & ) #ifdef DM_PARALLEL From 92f657d7b4cbc0eaf8dc00c5f49657c13b65af1e Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Mon, 24 Jan 2022 19:38:11 -0700 Subject: [PATCH 10/15] Play with optional args in TLADJ solve -> first rk part1 modified: dyn_em/module_first_rk_step_part1.F modified: dyn_em/solve_em.F modified: wrftladj/module_first_rk_step_part1_ad.F modified: wrftladj/module_first_rk_step_part1_tl.F modified: wrftladj/solve_em_ad.F modified: wrftladj/solve_em_tl.F --- dyn_em/module_first_rk_step_part1.F | 6 +++--- dyn_em/solve_em.F | 6 +++--- wrftladj/module_first_rk_step_part1_ad.F | 2 +- wrftladj/module_first_rk_step_part1_tl.F | 2 +- wrftladj/solve_em_ad.F | 6 ++---- wrftladj/solve_em_tl.F | 4 ++-- 6 files changed, 12 insertions(+), 14 deletions(-) diff --git a/dyn_em/module_first_rk_step_part1.F b/dyn_em/module_first_rk_step_part1.F index 3f970ce63a..e0eafeac6b 100644 --- a/dyn_em/module_first_rk_step_part1.F +++ b/dyn_em/module_first_rk_step_part1.F @@ -90,7 +90,7 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & REAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_fdda3d),INTENT(INOUT) :: fdda3d REAL ,DIMENSION(ims:ime,1:1,jms:jme,num_fdda2d),INTENT(INOUT) :: fdda2d REAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_aerod),INTENT(INOUT) :: aerod - REAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_aerocu),INTENT(INOUT) ::aerocu + REAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_aerocu),INTENT(INOUT), OPTIONAL ::aerocu REAL ,DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: psim REAL ,DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: psih REAL ,DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: gz1oz0 @@ -118,9 +118,9 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & INTEGER, INTENT(IN) :: k_start, k_end LOGICAL, INTENT(IN), OPTIONAL :: f_flux - LOGICAL, INTENT(IN) :: restart_flag + LOGICAL, INTENT(IN), OPTIONAL :: restart_flag - LOGICAL, INTENT(IN) :: feedback_is_ready ! For WRF-CMAQ coupled model, indicates feedback information is available + LOGICAL, INTENT(IN), OPTIONAL :: feedback_is_ready ! For WRF-CMAQ coupled model, indicates feedback information is available ! Local real :: HYDRO_dt diff --git a/dyn_em/solve_em.F b/dyn_em/solve_em.F index a1c7690c22..0d1c743c5e 100644 --- a/dyn_em/solve_em.F +++ b/dyn_em/solve_em.F @@ -825,9 +825,9 @@ END SUBROUTINE CMAQ_DRIVER , imsy, imey, jmsy, jmey, kmsy, kmey & , ipsy, ipey, jpsy, jpey, kpsy, kpey & , k_start , k_end & - , f_flux & - , aerocu & - , restart_flag & + , f_flux=f_flux & + , aerocu=aerocu & + , restart_flag=restart_flag & , feedback_is_ready=feedback_is_ready & ) diff --git a/wrftladj/module_first_rk_step_part1_ad.F b/wrftladj/module_first_rk_step_part1_ad.F index 2051a41d9c..927c5cdabe 100644 --- a/wrftladj/module_first_rk_step_part1_ad.F +++ b/wrftladj/module_first_rk_step_part1_ad.F @@ -134,7 +134,7 @@ SUBROUTINE a_first_rk_step_part1 ( grid , config_flags & INTEGER, INTENT(IN) :: k_start, k_end LOGICAL, INTENT(IN), OPTIONAL :: f_flux - LOGICAL, INTENT(IN) :: feedback_is_ready ! For WRF-CMAQ coupled model, indicates feedback information is available + LOGICAL, INTENT(IN), OPTIONAL :: feedback_is_ready ! For WRF-CMAQ coupled model, indicates feedback information is available ! Local REAL, DIMENSION( ims:ime, jms:jme ) :: exch_temf ! 1/7/09 WA diff --git a/wrftladj/module_first_rk_step_part1_tl.F b/wrftladj/module_first_rk_step_part1_tl.F index 3c290c889b..cb61d747d3 100644 --- a/wrftladj/module_first_rk_step_part1_tl.F +++ b/wrftladj/module_first_rk_step_part1_tl.F @@ -135,7 +135,7 @@ SUBROUTINE g_first_rk_step_part1 ( grid , config_flags & INTEGER, INTENT(IN) :: k_start, k_end LOGICAL, INTENT(IN), OPTIONAL :: f_flux - LOGICAL, INTENT(IN) :: feedback_is_ready ! For WRF-CMAQ coupled model, indicates feedback information is available + LOGICAL, INTENT(IN), OPTIONAL :: feedback_is_ready ! For WRF-CMAQ coupled model, indicates feedback information is available ! Local REAL, DIMENSION( ims:ime, jms:jme ) :: exch_temf ! 1/7/09 WA diff --git a/wrftladj/solve_em_ad.F b/wrftladj/solve_em_ad.F index 4ae309e574..06cfec4ec6 100644 --- a/wrftladj/solve_em_ad.F +++ b/wrftladj/solve_em_ad.F @@ -686,10 +686,8 @@ SUBROUTINE solve_em_ad ( grid , config_flags & , imsy, imey, jmsy, jmey, kmsy, kmey & , ipsy, ipey, jpsy, jpey, kpsy, kpey & , k_start , k_end & - , f_flux & - , aerocu & - , restart_flag & - , feedback_is_ready & + , f_flux=f_flux & + , feedback_is_ready=feedback_is_ready & ) #ifdef DM_PARALLEL diff --git a/wrftladj/solve_em_tl.F b/wrftladj/solve_em_tl.F index 957b46789c..789c62d591 100644 --- a/wrftladj/solve_em_tl.F +++ b/wrftladj/solve_em_tl.F @@ -664,8 +664,8 @@ SUBROUTINE solve_em_tl ( grid , config_flags & , imsy,imey,jmsy,jmey,kmsy,kmey & , ipsy,ipey,jpsy,jpey,kpsy,kpey & , k_start , k_end & - , f_flux & - , feedback_is_ready & + , f_flux=f_flux & + , feedback_is_ready=feedback_is_ready & ) #ifdef DM_PARALLEL From 63a6111c6525cb0aba1da0cd9e58e3b1f9098986 Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Fri, 21 Jan 2022 22:41:49 -0700 Subject: [PATCH 11/15] BF: NETCDFPAR "Error" in build log when option not used TYPE: bug fix KEYWORDS: netcdfpar, Error SOURCE: internal DESCRIPTION OF CHANGES: Problem: With PR #1552 "Parallel netcdf-4 IO option" (SHA1 3cd4713), when then code was built without the new parallel NetCDF4 compression, the build log had an `Error`. The problem was related to constructing the object files in the io_netcdfpar directory. When the option is not selected at compile time, then we do not care about errors in the directory that will never be used. Solution: If the NETCDFPAR option is not selected at compile time, then do SKIP going into the io_netcdfpar directory all together. LIST OF MODIFIED FILES: m Makefile m arch/Config.pl m arch/configure.defaults TESTS CONDUCTED: 1. Without the NETCDFPAR parameter set, the build for the io_netcdfpar directory is bypassed: ``` cd ../io_netcdfpar ; \ echo SKIPPING make -i -r NETCDFPARPATH="/glade/u/apps/ch/opt/netcdf/4.7.3/gnu/9.1.0" \ cd ../io_netcdfpar ; \ echo SKIPPING make -i -r NETCDFPARPATH="/glade/u/apps/ch/opt/netcdf/4.7.3/gnu/9.1.0" \ ``` 2. When the NETCDFPAR env variable is set, the build includes the io_netcdfpar directory: cd ../io_netcdfpar ; \ make -i -r NETCDFPARPATH="/glade/u/apps/ch/opt/netcdf/4.8.0/gnu/9.1.0" \ cd ../io_netcdfpar ; \ make -i -r NETCDFPARPATH="/glade/u/apps/ch/opt/netcdf/4.8.0/gnu/9.1.0" \ ``` 3. Jenkins tests are all PASS. RELEASE NOTE: Include a stand-alone message suitable for the inclusion in the minor and annual releases. A publication citation is appropriate. --- Makefile | 4 ++-- arch/Config.pl | 2 ++ arch/configure.defaults | 44 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 48 insertions(+), 2 deletions(-) diff --git a/Makefile b/Makefile index 499884a078..2a2004570b 100644 --- a/Makefile +++ b/Makefile @@ -921,13 +921,13 @@ framework : LIB_LOCAL="$(LIB_LOCAL)" \ ESMF_MOD_DEPENDENCE="$(ESMF_MOD_DEPENDENCE)" AR="INTERNAL_BUILD_ERROR_SHOULD_NOT_NEED_AR"; \ cd ../io_netcdfpar ; \ - $(MAKE) NETCDFPARPATH="$(NETCDFPATH)" \ + $(NETCDFPAR_BUILD) $(MAKE) NETCDFPARPATH="$(NETCDFPATH)" \ FC="$(FC) $(FCBASEOPTS) $(PROMOTION) $(FCDEBUG) $(OMP)" RANLIB="$(RANLIB)" \ CPP="$(CPP)" LDFLAGS="$(LDFLAGS)" TRADFLAG="$(TRADFLAG)" ESMF_IO_LIB_EXT="$(ESMF_IO_LIB_EXT)" \ LIB_LOCAL="$(LIB_LOCAL)" \ ESMF_MOD_DEPENDENCE="$(ESMF_MOD_DEPENDENCE)" AR="INTERNAL_BUILD_ERROR_SHOULD_NOT_NEED_AR" diffwrf; \ cd ../io_netcdfpar ; \ - $(MAKE) NETCDFPARPATH="$(NETCDFPATH)" \ + $(NETCDFPAR_BUILD) $(MAKE) NETCDFPARPATH="$(NETCDFPATH)" \ FC="$(SFC) $(FCBASEOPTS) $(PROMOTION) $(FCDEBUG) $(OMP)" RANLIB="$(RANLIB)" \ CPP="$(CPP)" LDFLAGS="$(LDFLAGS)" TRADFLAG="$(TRADFLAG)" ESMF_IO_LIB_EXT="$(ESMF_IO_LIB_EXT)" \ LIB_LOCAL="$(LIB_LOCAL)" \ diff --git a/arch/Config.pl b/arch/Config.pl index f1803e8a72..fb4a0b4ca4 100644 --- a/arch/Config.pl +++ b/arch/Config.pl @@ -656,6 +656,7 @@ if ( $sw_netcdfpar_path ) { $_ =~ s/CONFIGURE_WRFIO_NFPAR/wrfio_nfpar/g ; $_ =~ s:CONFIGURE_NETCDFPAR_FLAG:-DNETCDFPAR: ; + $_ =~ s:CONFIGURE_NETCDFPAR_BUILD:: ; if ( $ENV{NETCDFPAR_LDFLAGS} ) { $_ =~ s:CONFIGURE_NETCDFPAR_LIB_PATH:\$\(WRF_SRC_ROOT_DIR\)/external/io_netcdfpar/libwrfio_nfpar.a $ENV{NETCDFPAR_LDFLAGS} : ; } elsif ( $sw_os eq "Interix" ) { @@ -667,6 +668,7 @@ else { $_ =~ s/CONFIGURE_WRFIO_NFPAR//g ; $_ =~ s:CONFIGURE_NETCDFPAR_FLAG::g ; + $_ =~ s:CONFIGURE_NETCDFPAR_BUILD:echo SKIPPING: ; $_ =~ s:CONFIGURE_NETCDFPAR_LIB_PATH::g ; } diff --git a/arch/configure.defaults b/arch/configure.defaults index cf307a502f..7f5740a2a3 100644 --- a/arch/configure.defaults +++ b/arch/configure.defaults @@ -42,6 +42,7 @@ RANLIB = ls RLFLAGS = #ranlib CC_TOOLS = cc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux i486 i586 i686 armv7l aarch64, gfortran compiler with gcc #serial smpar dmpar dm+sm @@ -86,6 +87,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux i486 i586 i686, g95 compiler with gcc #serial dmpar @@ -129,6 +131,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le i486 i586 i686, PGI compiler with gcc # serial smpar dmpar dm+sm @@ -172,6 +175,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le, PGI compiler with pgcc, SGI MPT # serial smpar dmpar dm+sm @@ -215,6 +219,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le, PGI accelerator compiler with gcc # serial smpar dmpar dm+sm @@ -257,6 +262,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le i486 i586 i686, ifort compiler with icc #serial smpar dmpar dm+sm @@ -334,6 +340,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le i486 i586 i686, Xeon Phi (MIC architecture) ifort compiler with icc # dm+sm @@ -381,6 +388,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = gcc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le i486 i586 i686, Xeon (SNB with AVX mods) ifort compiler with icc # serial smpar dmpar dm+sm @@ -428,6 +436,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = gcc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le i486 i586 i686, ifort compiler with icc, SGI MPT #serial smpar dmpar dm+sm @@ -499,6 +508,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le i486 i586 i686, ifort compiler with icc, IBM POE #serial smpar dmpar dm+sm @@ -548,6 +558,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH ia64 Linux ifort compiler with icc 9.x,10.x #serial smpar dmpar dm+sm @@ -629,6 +640,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux SGI Altix, ifort compiler with icc 9.x,10.x #serial smpar dmpar dm+sm @@ -712,6 +724,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux i486 i586 i686 x86_64 ppc64le, PathScale compiler with pathcc #serial dmpar @@ -755,6 +768,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le, gfortran compiler with gcc #serial smpar dmpar dm+sm @@ -799,6 +813,7 @@ M4 = m4 -G RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Darwin (MACOS) PGI compiler with pgcc #serial smpar dmpar dm+sm @@ -842,6 +857,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = -c CC_TOOLS = cc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Darwin (MACOS) intel compiler with icc #serial smpar dmpar dm+sm @@ -888,6 +904,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = -c CC_TOOLS = cc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Darwin (MACOS) intel compiler with clang EDIT FOR OPENMPI #serial smpar dmpar dm+sm @@ -933,6 +950,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = -c CC_TOOLS = cc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Darwin (MACOS) g95 with gcc #serial dmpar @@ -977,6 +995,7 @@ M4 = m4 -B 14000 RANLIB = ranlib -c RLFLAGS = -c CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Darwin (MACOS) gfortran with gcc #serial smpar dmpar dm+sm @@ -1021,6 +1040,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = -c CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Darwin (MACOS) gfortran with clang #serial smpar dmpar dm+sm @@ -1065,6 +1085,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = -c CC_TOOLS = clang +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Darwin (MACOS) xlf #serial dmpar @@ -1110,6 +1131,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = -c CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH AIX xlf compiler with xlc #serial smpar dmpar dm+sm @@ -1158,6 +1180,7 @@ M4 = m4 -B 20000 RANLIB = ranlib RLFLAGS = CC_TOOLS = cc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le i486 i586 i686, xlf compiler with xlc # serial smpar dmpar dm+sm @@ -1209,6 +1232,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = cc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Cray XC CLE/Linux x86_64, PGI compiler with gcc # serial dmpar smpar dm+sm @@ -1271,6 +1295,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Cray XE and XC CLE/Linux x86_64, Cray CCE compiler # serial dmpar smpar dm+sm @@ -1321,6 +1346,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = gcc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Cray XC CLE/Linux x86_64, Xeon ifort compiler # serial dmpar smpar dm+sm @@ -1370,6 +1396,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = gcc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### @@ -1422,6 +1449,8 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = cc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD + ########################################################### #ARCH Linux ppc64 BG /P xlf compiler with xlc # smpar dmpar dm+sm # developed on surveyor.alcf.anl.gov (thanks to ANL/ALCF) @@ -1469,6 +1498,8 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = cc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD + ########################################################### #ARCH Linux ppc64 IBM Blade Server xlf compiler with xlc # dmpar # provided by Luis C. Cana Cascallar for IBM JS21 blade server, May 2009 @@ -1513,6 +1544,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = xlc -q64 +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le i486 i586 i686, PGI compiler with pgcc # serial smpar dmpar dm+sm @@ -1556,6 +1588,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH CYGWIN_NT i686, PGI compiler on Windows # serial smpar dmpar dm+sm @@ -1599,6 +1632,7 @@ M4 = NA RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD LIB_EXTERNAL = \ ../external/io_netcdf/libwrfio_nf.a CONFIGURE_NETCDF_PATH/lib/libnetcdf.lib \ @@ -1656,6 +1690,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Darwin (MACOS) PGI compiler with pgcc -f90= #serial smpar dmpar dm+sm @@ -1699,6 +1734,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = -c CC_TOOLS = cc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Darwin (MACOS) intel compiler with icc #serial smpar dmpar dm+sm @@ -1745,6 +1781,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = -c CC_TOOLS = cc +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Darwin (MACOS) gfortran with gcc openmpi #serial smpar dmpar dm+sm @@ -1789,6 +1826,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = -c CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux x86_64 ppc64le i486 i586 i686, PGI compiler with pgcc -f90= # serial smpar dmpar dm+sm @@ -1832,6 +1870,7 @@ M4 = m4 -B 14000 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux HSW/BDW x86_64 ppc64le i486 i586 i686, ifort compiler with icc #serial smpar dmpar dm+sm @@ -1876,6 +1915,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux KNL x86_64 ppc64le i486 i586 i686, ifort compiler with icc #serial smpar dmpar dm+sm @@ -1920,6 +1960,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH CYGWIN_NT i686 x86_64 Cygwin, gfortran compiler with gcc #serial smpar dmpar dm+sm @@ -1964,6 +2005,7 @@ M4 = m4 -G RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD LIB_EXTERNAL = \ $(WRF_SRC_ROOT_DIR)/external/io_netcdf/libwrfio_nf.a CONFIGURE_NETCDF_PATH/lib/libnetcdf.dll.a \ @@ -2022,6 +2064,7 @@ M4 = m4 -G RANLIB = ranlib RLFLAGS = CC_TOOLS = $(SCC) +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD #insert new stanza here @@ -2067,6 +2110,7 @@ M4 = m4 RANLIB = ranlib RLFLAGS = CC_TOOLS = /usr/bin/gcc -Wall +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD #insert new stanza before the Fujitsu block, keep Fujitsu at the end of the list ########################################################### From b55d9906433ec28add0e582f70f3e817c5efbd97 Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Fri, 21 Jan 2022 23:09:25 -0700 Subject: [PATCH 12/15] Bug in the configure file, does not print out NC version When the NETCDFPAR option was not selected, there is no need to know the version number for NetCDF. The version was never assigned, and a later echo of the version was blank. This is now rectified. ``` NetCDF version: 4.7.3 Enabled NetCDF-4/HDF-5: yes NetCDF built with PnetCDF: no ``` M configure --- configure | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/configure b/configure index 5f0e780465..1c7b6d9e9f 100755 --- a/configure +++ b/configure @@ -159,6 +159,12 @@ if test -n "$PERL" ; then fi +if [ -e $NETCDF/bin/nc-config ] ; then + ncversion=`nc-config --version | awk '{print $2}'` +else + ncversion=OLD +fi + PROBS=FALSE if [ -n "$NETCDFPAR" ] ; then NETCDF="$NETCDFPAR" From 213893d3be21c60d393d32140a5b92c709207a77 Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Sat, 22 Jan 2022 23:42:19 -0700 Subject: [PATCH 13/15] NCPAR FIX nix DA 23:42 --- arch/configure.defaults | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/arch/configure.defaults b/arch/configure.defaults index 7f5740a2a3..081f64cfba 100644 --- a/arch/configure.defaults +++ b/arch/configure.defaults @@ -42,7 +42,7 @@ RANLIB = ls RLFLAGS = #ranlib CC_TOOLS = cc -NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD +NETCDFPAR_BUILD = CONFIGURE_NETCDFPAR_BUILD ########################################################### #ARCH Linux i486 i586 i686 armv7l aarch64, gfortran compiler with gcc #serial smpar dmpar dm+sm From 8827d3e347eb806542a5c5cbc33ca6b1caae7796 Mon Sep 17 00:00:00 2001 From: dwongepa Date: Tue, 25 Jan 2022 11:54:16 -0500 Subject: [PATCH 14/15] modified arch/Config.pl to accommodate file name shortening in the WRF-CMAQ coupled model --- arch/Config.pl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/arch/Config.pl b/arch/Config.pl index fb4a0b4ca4..f31b33827e 100644 --- a/arch/Config.pl +++ b/arch/Config.pl @@ -504,9 +504,7 @@ close (FH); if ( $wrf_cmaq_option eq 1 ) # build WRF-CMAQ coupled model - { system("cp Registry/registry.WRF-CMAQ-twoway_full Registry/registry.WRF-CMAQ-twoway"); - - if ( $lib_path_wo_cmaq == 1 ) + { if ( $lib_path_wo_cmaq == 1 ) { open (FILE, "; close (FILE); @@ -521,7 +519,7 @@ } if ( $registry_wo_cmaq == 1 ) - { open (FILE, "; close (FILE); @@ -529,7 +527,7 @@ { $_ =~ s/#state/state/g; } - open (FILE, ">Registry/registry.WRF-CMAQ-twoway") || die "File not found"; + open (FILE, ">Registry/registry.CMAQ") || die "File not found"; print FILE @lines; close (FILE); } @@ -574,7 +572,7 @@ } if ( $registry_wo_cmaq == 0 ) - { open (FILE, "; close (FILE); @@ -582,7 +580,7 @@ { $_ =~ s/state/#state/g; } - open (FILE, ">Registry/registry.WRF-CMAQ-twoway") || die "File not found"; + open (FILE, ">Registry/registry.CMAQ") || die "File not found"; print FILE @lines; close (FILE); } From 7ee699a9a90831324efa8a9d2970f402bce81d95 Mon Sep 17 00:00:00 2001 From: Dave Gill Date: Tue, 25 Jan 2022 13:38:00 -0700 Subject: [PATCH 15/15] Replace nonstandard ISNAN() with IEEE_IS_NAN() This fixes the PGI compiler, which frowned on non standard math functions. modified: phys/module_firebrand_spotting.F --- phys/module_firebrand_spotting.F | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/phys/module_firebrand_spotting.F b/phys/module_firebrand_spotting.F index 28e6c13022..5e0d64b4f7 100644 --- a/phys/module_firebrand_spotting.F +++ b/phys/module_firebrand_spotting.F @@ -9,6 +9,7 @@ MODULE module_firebrand_spotting + USE, INTRINSIC :: IEEE_ARITHMETIC USE module_domain, ONLY : get_ijk_from_grid, domain ! grid USE module_configure, ONLY : grid_config_rec_type ! config_flags USE module_symbols_util, ONLY : WRFU_TimeInterval, WRFU_TimeIntervalGet, WRFU_TimeIntervalSet @@ -2332,9 +2333,9 @@ SUBROUTINE firebrand_spotting_em_driver( & !------------------------------------------------------------------------------- sparse_mask = .FALSE. - WHERE (ISNAN(fs_p_mass) .OR. ISNAN(fs_p_diam) .OR. & - ISNAN(fs_p_effd) .OR. ISNAN(fs_p_temp) .OR. & - ISNAN(fs_p_tvel)) fs_p_effd = ZERO_dp ! ediam=zero will be removed in remove_br + WHERE (IEEE_IS_NAN(fs_p_mass) .OR. IEEE_IS_NAN(fs_p_diam) .OR. & + IEEE_IS_NAN(fs_p_effd) .OR. IEEE_IS_NAN(fs_p_temp) .OR. & + IEEE_IS_NAN(fs_p_tvel)) fs_p_effd = ZERO_dp ! ediam=zero will be removed in remove_br !=============================================================================== @@ -3141,7 +3142,7 @@ SUBROUTINE advect_xyz_m(grid, xp, yp, hgt, dtp, idp, & CYCLE ! Skip physics ENDIF - IF (ISNAN(zout(pp)) .OR. ABS(zout(pp)) > HUGE(1.0)) THEN + IF (IEEE_IS_NAN(zout(pp)) .OR. ABS(zout(pp)) > HUGE(1.0)) THEN fs_p_prop(pp) = p_properties(ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP) zout(pp) = ZERO_dp xout(pp) = ZERO_dp @@ -3192,7 +3193,7 @@ SUBROUTINE advect_xyz_m(grid, xp, yp, hgt, dtp, idp, & loc_d = loc_d, & fbprop = fs_p_prop(pp)) - IF (ISNAN(fs_p_prop(pp)%p_tvel)) THEN + IF (IEEE_IS_NAN(fs_p_prop(pp)%p_tvel)) THEN fs_p_prop(pp) = p_properties(ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP) zout(pp) = ZERO_dp