Changeset 4489 for LMDZ6/trunk/libf/phylmd/ecrad/radiation_config.F90
- Timestamp:
- Mar 31, 2023, 8:42:57 PM (17 months ago)
- Location:
- LMDZ6/trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk
- Property svn:mergeinfo changed
/LMDZ6/branches/LMDZ_ECRad (added) merged: 4175,4177-4183,4188,4192,4200-4203,4355,4366,4387-4388,4390,4444,4482,4486,4488
- Property svn:mergeinfo changed
-
LMDZ6/trunk/libf/phylmd/ecrad/radiation_config.F90
r4115 r4489 26 26 ! 2019-02-03 R. Hogan Added ability to fix out-of-physical-bounds inputs 27 27 ! 2019-02-10 R. Hogan Renamed "encroachment" to "entrapment" 28 ! 2020-05-18 R. Hogan Moved out_of_bounds_* to radiation_check.F90 29 ! 2021-07-04 R. Hogan Numerous changes for ecCKD and general cloud/aerosol optics 28 30 ! 29 31 ! Note: The aim is for ecRad in the IFS to be as similar as possible … … 37 39 38 40 use radiation_cloud_optics_data, only : cloud_optics_type 41 use radiation_general_cloud_optics_data, only : general_cloud_optics_type 39 42 use radiation_aerosol_optics_data, only : aerosol_optics_type 40 43 use radiation_pdf_sampler, only : pdf_sampler_type 41 44 use radiation_cloud_cover, only : OverlapName, & 42 45 & IOverlapMaximumRandom, IOverlapExponentialRandom, IOverlapExponential 46 use radiation_ecckd, only : ckd_model_type 43 47 44 48 implicit none … … 70 74 & IEntrapmentExplicitNonFractal, & ! As above but ignore fractal nature of clouds 71 75 & IEntrapmentMaximum ! Complete horizontal homogenization within regions (old SPARTACUS assumption) 72 73 76 end enum 74 77 … … 94 97 ! Gas models 95 98 enum, bind(c) 96 enumerator IGasModelMonochromatic, IGasModelIFSRRTMG 99 enumerator IGasModelMonochromatic, IGasModelIFSRRTMG, IGasModelECCKD 97 100 end enum 98 character(len=*), parameter :: GasModelName(0:1) = (/ 'Monochromatic', & 99 & 'RRTMG-IFS ' /) 101 character(len=*), parameter :: GasModelName(0:2) = (/ 'Monochromatic', & 102 & 'RRTMG-IFS ', & 103 & 'ECCKD '/) 100 104 101 105 ! Hydrometeor scattering models … … 130 134 integer, parameter :: NMaxAerosolTypes = 256 131 135 136 ! Maximum number of different cloud types that can be provided 137 integer, parameter :: NMaxCloudTypes = 12 138 132 139 ! Maximum number of shortwave albedo and longwave emissivity 133 140 ! intervals … … 155 162 character(len=511) :: directory_name = '.' 156 163 164 ! If this is true then support arbitrary hydrometeor types (not 165 ! just ice and liquid) and arbitrary spectral discretization (not 166 ! just RRTMG). It is required that this is true if the ecCKD gas 167 ! optics model is selected. General cloud optics has only been 168 ! available from ecRad version 1.5. 169 logical :: use_general_cloud_optics = .true. 170 171 ! If this is true then support aerosol properties at an arbitrary 172 ! spectral discretization (not just RRTMG). It is required that 173 ! this is true if the ecCKD gas optics model is selected. 174 logical :: use_general_aerosol_optics = .true. 175 157 176 ! Cloud is deemed to be present in a layer if cloud fraction 158 177 ! exceeds this value … … 168 187 ! (2000)? 169 188 logical :: use_beta_overlap = .false. 189 190 ! Use a more vectorizable McICA cloud generator, at the expense of 191 ! more random numbers being generated? This is the default on NEC 192 ! SX. 193 #ifdef __SX__ 194 logical :: use_vectorizable_generator = .true. 195 #else 196 logical :: use_vectorizable_generator = .false. 197 #endif 170 198 171 199 ! Shape of sub-grid cloud water PDF … … 236 264 logical :: do_sw_delta_scaling_with_gases = .false. 237 265 238 ! Codes describing the gas and cloud scattering models to use, the 239 ! latter of which is currently not used 266 ! Codes describing the gas model 240 267 integer :: i_gas_model = IGasModelIFSRRTMG 241 ! integer :: i_cloud_model242 268 243 269 ! Optics if i_gas_model==IGasModelMonochromatic. … … 270 296 ! according to the spectral overlap of each interval with each 271 297 ! band 272 logical :: do_nearest_spectral_sw_albedo = . true.273 logical :: do_nearest_spectral_lw_emiss = . true.298 logical :: do_nearest_spectral_sw_albedo = .false. 299 logical :: do_nearest_spectral_lw_emiss = .false. 274 300 275 301 ! User-defined monotonically increasing wavelength bounds (m) 276 302 ! between input surface albedo/emissivity intervals. Implicitly 277 ! the first interval starts at zero and the last ends at infinity. 303 ! the first interval starts at zero and the last ends at 304 ! infinity. These must be set with define_sw_albedo_intervals and 305 ! define_lw_emiss_intervals. 278 306 real(jprb) :: sw_albedo_wavelength_bound(NMaxAlbedoIntervals-1) = -1.0_jprb 279 real(jprb) :: lw_emiss_wavelength_bound( NMaxAlbedoIntervals-1) 307 real(jprb) :: lw_emiss_wavelength_bound( NMaxAlbedoIntervals-1) = -1.0_jprb 280 308 281 309 ! The index to the surface albedo/emissivity intervals for each of … … 296 324 logical :: do_3d_effects = .true. 297 325 326 character(len=511) :: cloud_type_name(NMaxCloudTypes) = ["","","","","","","","","","","",""] 327 ! & 328 ! & = ["mie_droplet ", & 329 ! & "baum-general-habit-mixture_ice"] 330 331 ! Spectral averaging method to use with generalized cloud optics; 332 ! see Edwards & Slingo (1996) for definition. Experimentation 333 ! with ecRad suggests that "thick" averaging is more accurate for 334 ! both liquid and ice clouds. 335 logical :: use_thick_cloud_spectral_averaging(NMaxCloudTypes) & 336 & = [.true.,.true.,.true.,.true.,.true.,.true., & 337 & .true.,.true.,.true.,.true.,.true.,.true.] 338 298 339 ! To what extent do we include "entrapment" effects in the 299 340 ! SPARTACUS solver? This essentially means that in a situation … … 420 461 ! doesn't start with a '/' character then it will be prepended by 421 462 ! the contents of directory_name. 422 character(len=511) :: ice_optics_override_file_name = ''423 character(len=511) :: liq_optics_override_file_name = ''463 character(len=511) :: ice_optics_override_file_name = '' 464 character(len=511) :: liq_optics_override_file_name = '' 424 465 character(len=511) :: aerosol_optics_override_file_name = '' 466 character(len=511) :: gas_optics_sw_override_file_name = '' 467 character(len=511) :: gas_optics_lw_override_file_name = '' 425 468 426 469 ! Optionally override the look-up table file for the cloud-water … … 428 471 character(len=511) :: cloud_pdf_override_file_name = '' 429 472 473 ! Do we compute cloud, aerosol and surface optical properties per 474 ! g point? Not available with RRTMG gas optics model. 475 logical :: do_cloud_aerosol_per_sw_g_point = .true. 476 logical :: do_cloud_aerosol_per_lw_g_point = .true. 477 478 ! Do we weight the mapping from surface emissivity/albedo to 479 ! g-point/band weighting by a reference Planck function (more 480 ! accurate) or weighting each wavenumber equally (less accurate 481 ! but consistent with IFS Cycle 48r1 and earlier)? 482 logical :: do_weighted_surface_mapping = .true. 483 484 ! COMPUTED PARAMETERS 485 486 ! Users of this library should not edit these parameters directly; 487 ! they are set by the "consolidate" routine 488 430 489 ! Has "consolidate" been called? 431 490 logical :: is_consolidated = .false. 432 491 433 ! COMPUTED PARAMETERS 434 ! Users of this library should not edit these parameters directly; 435 ! they are set by the "consolidate" routine 436 437 ! Wavenumber range for each band, in cm-1, which will be allocated 438 ! to be of length n_bands_sw or n_bands_lw 439 real(jprb), allocatable, dimension(:) :: wavenumber1_sw 440 real(jprb), allocatable, dimension(:) :: wavenumber2_sw 441 real(jprb), allocatable, dimension(:) :: wavenumber1_lw 442 real(jprb), allocatable, dimension(:) :: wavenumber2_lw 492 ! Fraction of each g point in each wavenumber interval, 493 ! dimensioned (n_wav_frac_[l|s]w, n_g_[l|s]w) 494 real(jprb), allocatable, dimension(:,:) :: g_frac_sw, g_frac_lw 443 495 444 496 ! If the nearest surface albedo/emissivity interval is to be used … … 490 542 integer :: n_canopy_bands_lw = 1 491 543 544 ! Data structures containing gas optics description in the case of 545 ! ecCKD 546 type(ckd_model_type) :: gas_optics_sw, gas_optics_lw 547 492 548 ! Data structure containing cloud scattering data 493 549 type(cloud_optics_type) :: cloud_optics 550 551 ! Number of general cloud types, default liquid and ice 552 integer :: n_cloud_types = 2 553 554 ! List of data structures (one per cloud type) containing cloud 555 ! scattering data 556 type(general_cloud_optics_type), allocatable :: cloud_optics_sw(:) 557 type(general_cloud_optics_type), allocatable :: cloud_optics_lw(:) 494 558 495 559 ! Data structure containing aerosol scattering data … … 502 566 character(len=511) :: ice_optics_file_name, & 503 567 & liq_optics_file_name, & 504 & aerosol_optics_file_name 568 & aerosol_optics_file_name, & 569 & gas_optics_sw_file_name, & 570 & gas_optics_lw_file_name 505 571 506 572 ! McICA PDF look-up table file name … … 515 581 ! g points or the number of bands 516 582 integer :: n_spec_sw = 0, n_spec_lw = 0 583 584 ! Number of wavenumber intervals used to describe the mapping from 585 ! g-points to wavenumber space 586 integer :: n_wav_frac_sw = 0, n_wav_frac_lw = 0 517 587 518 588 ! Dimensions to store variables that are only needed if longwave … … 539 609 procedure :: define_sw_albedo_intervals 540 610 procedure :: define_lw_emiss_intervals 541 procedure :: consolidate_intervals 611 procedure :: set_aerosol_wavelength_mono 612 procedure :: consolidate_sw_albedo_intervals 613 procedure :: consolidate_lw_emiss_intervals 542 614 543 615 end type config_type … … 574 646 logical :: do_sw, do_lw, do_clear, do_sw_direct 575 647 logical :: do_3d_effects, use_expm_everywhere, use_aerosols 648 logical :: use_general_cloud_optics, use_general_aerosol_optics 576 649 logical :: do_lw_side_emissivity 577 650 logical :: do_3d_lw_multilayer_effects, do_fu_lw_ice_optics_bug … … 579 652 logical :: do_save_radiative_properties, do_save_spectral_flux 580 653 logical :: do_save_gpoint_flux, do_surface_sw_spectral_flux 581 logical :: use_beta_overlap, do_lw_derivatives 654 logical :: use_beta_overlap, do_lw_derivatives, use_vectorizable_generator 582 655 logical :: do_sw_delta_scaling_with_gases 583 656 logical :: do_canopy_fluxes_sw, do_canopy_fluxes_lw 584 657 logical :: use_canopy_full_spectrum_sw, use_canopy_full_spectrum_lw 585 658 logical :: do_canopy_gases_sw, do_canopy_gases_lw 659 logical :: do_cloud_aerosol_per_sw_g_point, do_cloud_aerosol_per_lw_g_point 660 logical :: do_weighted_surface_mapping 586 661 integer :: n_regions, iverbose, iverbosesetup, n_aerosol_types 587 662 real(jprb):: mono_lw_wavelength, mono_lw_total_od, mono_sw_total_od … … 596 671 character(511) :: liq_optics_override_file_name, ice_optics_override_file_name 597 672 character(511) :: cloud_pdf_override_file_name 673 character(511) :: gas_optics_sw_override_file_name, gas_optics_lw_override_file_name 598 674 character(63) :: liquid_model_name, ice_model_name, gas_model_name 599 675 character(63) :: sw_solver_name, lw_solver_name, overlap_scheme_name 600 676 character(63) :: sw_entrapment_name, sw_encroachment_name, cloud_pdf_shape_name 677 character(len=511) :: cloud_type_name(NMaxCloudTypes) = ["","","","","","","","","","","",""] 678 logical :: use_thick_cloud_spectral_averaging(NMaxCloudTypes) & 679 & = [.false.,.false.,.false.,.false.,.false.,.false., & 680 & .false.,.false.,.false.,.false.,.false.,.false.] 601 681 integer :: i_aerosol_type_map(NMaxAerosolTypes) ! More than 256 is an error 602 682 603 logical :: do_nearest_spectral_sw_albedo = .true.604 logical :: do_nearest_spectral_lw_emiss = .true.683 logical :: do_nearest_spectral_sw_albedo 684 logical :: do_nearest_spectral_lw_emiss 605 685 real(jprb) :: sw_albedo_wavelength_bound(NMaxAlbedoIntervals-1) 606 686 real(jprb) :: lw_emiss_wavelength_bound( NMaxAlbedoIntervals-1) … … 609 689 610 690 integer :: iunit ! Unit number of namelist file 611 612 logical :: lldeb_conf = .false.613 691 614 692 namelist /radiation/ do_sw, do_lw, do_sw_direct, & … … 622 700 & ice_optics_override_file_name, liq_optics_override_file_name, & 623 701 & aerosol_optics_override_file_name, cloud_pdf_override_file_name, & 702 & gas_optics_sw_override_file_name, gas_optics_lw_override_file_name, & 624 703 & liquid_model_name, ice_model_name, max_3d_transfer_rate, & 625 704 & min_cloud_effective_size, overhang_factor, encroachment_scaling, & … … 627 706 & do_canopy_fluxes_sw, do_canopy_fluxes_lw, & 628 707 & do_canopy_gases_sw, do_canopy_gases_lw, & 708 & use_general_cloud_optics, use_general_aerosol_optics, & 629 709 & do_sw_delta_scaling_with_gases, overlap_scheme_name, & 630 & sw_solver_name, lw_solver_name, use_beta_overlap, &710 & sw_solver_name, lw_solver_name, use_beta_overlap, use_vectorizable_generator, & 631 711 & use_expm_everywhere, iverbose, iverbosesetup, & 632 712 & cloud_inhom_decorr_scaling, cloud_fraction_threshold, & … … 637 717 & mono_lw_single_scattering_albedo, mono_sw_single_scattering_albedo, & 638 718 & mono_lw_asymmetry_factor, mono_sw_asymmetry_factor, & 639 & cloud_pdf_shape_name, &719 & cloud_pdf_shape_name, cloud_type_name, use_thick_cloud_spectral_averaging, & 640 720 & do_nearest_spectral_sw_albedo, do_nearest_spectral_lw_emiss, & 641 721 & sw_albedo_wavelength_bound, lw_emiss_wavelength_bound, & 642 & i_sw_albedo_index, i_lw_emiss_index 643 722 & i_sw_albedo_index, i_lw_emiss_index, & 723 & do_cloud_aerosol_per_lw_g_point, & 724 & do_cloud_aerosol_per_sw_g_point, do_weighted_surface_mapping 725 644 726 real(jprb) :: hook_handle 645 727 … … 670 752 ice_optics_override_file_name = this%ice_optics_override_file_name 671 753 aerosol_optics_override_file_name = this%aerosol_optics_override_file_name 754 gas_optics_sw_override_file_name = this%gas_optics_sw_override_file_name 755 gas_optics_lw_override_file_name = this%gas_optics_lw_override_file_name 672 756 use_expm_everywhere = this%use_expm_everywhere 673 757 use_aerosols = this%use_aerosols … … 679 763 iverbose = this%iverbose 680 764 iverbosesetup = this%iverbosesetup 765 use_general_cloud_optics = this%use_general_cloud_optics 766 use_general_aerosol_optics = this%use_general_aerosol_optics 681 767 cloud_fraction_threshold = this%cloud_fraction_threshold 682 768 cloud_mixing_ratio_threshold = this%cloud_mixing_ratio_threshold 683 769 use_beta_overlap = this%use_beta_overlap 770 use_vectorizable_generator = this%use_vectorizable_generator 684 771 cloud_inhom_decorr_scaling = this%cloud_inhom_decorr_scaling 685 772 clear_to_thick_fraction = this%clear_to_thick_fraction … … 689 776 max_3d_transfer_rate = this%max_3d_transfer_rate 690 777 min_cloud_effective_size = this%min_cloud_effective_size 778 cloud_type_name = this%cloud_type_name 779 use_thick_cloud_spectral_averaging = this%use_thick_cloud_spectral_averaging 780 691 781 overhang_factor = this%overhang_factor 692 782 encroachment_scaling = -1.0_jprb … … 715 805 i_sw_albedo_index = this%i_sw_albedo_index 716 806 i_lw_emiss_index = this%i_lw_emiss_index 807 do_cloud_aerosol_per_lw_g_point = this%do_cloud_aerosol_per_lw_g_point 808 do_cloud_aerosol_per_sw_g_point = this%do_cloud_aerosol_per_sw_g_point 809 do_weighted_surface_mapping = this%do_weighted_surface_mapping 717 810 718 811 if (present(file_name) .and. present(unit)) then … … 746 839 end if 747 840 else 841 842 ! This version exits correctly, but provides less information 843 ! about how the namelist was incorrect 748 844 read(unit=iunit, iostat=iosread, nml=radiation) 845 846 ! Depending on compiler this version provides more information 847 ! about the error in the namelist 848 !read(unit=iunit, nml=radiation) 849 749 850 if (iosread /= 0) then 750 851 ! An error occurred reading the file … … 812 913 this%mono_sw_asymmetry_factor = mono_sw_asymmetry_factor 813 914 this%use_beta_overlap = use_beta_overlap 915 this%use_vectorizable_generator = use_vectorizable_generator 814 916 this%cloud_inhom_decorr_scaling = cloud_inhom_decorr_scaling 815 917 this%clear_to_thick_fraction = clear_to_thick_fraction … … 819 921 this%max_3d_transfer_rate = max_3d_transfer_rate 820 922 this%min_cloud_effective_size = max(1.0e-6_jprb, min_cloud_effective_size) 923 this%cloud_type_name = cloud_type_name 924 this%use_thick_cloud_spectral_averaging = use_thick_cloud_spectral_averaging 821 925 if (encroachment_scaling >= 0.0_jprb) then 822 926 this%overhang_factor = encroachment_scaling … … 832 936 this%ice_optics_override_file_name = ice_optics_override_file_name 833 937 this%aerosol_optics_override_file_name = aerosol_optics_override_file_name 938 this%gas_optics_sw_override_file_name = gas_optics_sw_override_file_name 939 this%gas_optics_lw_override_file_name = gas_optics_lw_override_file_name 940 this%use_general_cloud_optics = use_general_cloud_optics 941 this%use_general_aerosol_optics = use_general_aerosol_optics 834 942 this%cloud_fraction_threshold = cloud_fraction_threshold 835 943 this%cloud_mixing_ratio_threshold = cloud_mixing_ratio_threshold … … 845 953 this%i_sw_albedo_index = i_sw_albedo_index 846 954 this%i_lw_emiss_index = i_lw_emiss_index 847 848 ! AI mars 2022 849 if (lldeb_conf) then 850 print*,'**************PARAMETRES DE CONFIGURATION OFFLINE*******************' 851 print*,'config%iverbosesetup = ', iverbosesetup 852 print*,'config%do_lw = ', do_lw 853 print*,'config%do_sw = ', do_sw 854 print*,'config%do_clear = ', do_clear 855 print*,'config%do_sw_direct = ', do_sw_direct 856 print*,'config%do_3d_effects = ', do_3d_effects 857 print*,'config%do_3d_lw_multilayer_effects = ', do_3d_lw_multilayer_effects 858 print*,'config%do_lw_side_emissivity = ', do_lw_side_emissivity 859 print*,'config%use_expm_everywhere = ', use_expm_everywhere 860 print*,'config%use_aerosols = ', use_aerosols 861 print*,'config%do_lw_cloud_scattering = ', do_lw_cloud_scattering 862 print*,'config%do_lw_aerosol_scattering = ', do_lw_aerosol_scattering 863 print*,'config%nregions = ', n_regions 864 print*,'config%do_surface_sw_spectral_flux = ', do_surface_sw_spectral_flux 865 print*,'config%do_sw_delta_scaling_with_gases = ', & 866 do_sw_delta_scaling_with_gases 867 print*,'config%do_fu_lw_ice_optics_bug = ', do_fu_lw_ice_optics_bug 868 print*,'config%do_canopy_fluxes_sw = ', do_canopy_fluxes_sw 869 print*,'config%do_canopy_fluxes_lw = ', do_canopy_fluxes_lw 870 print*,'config%use_canopy_full_spectrum_sw = ', use_canopy_full_spectrum_sw 871 print*,'config%use_canopy_full_spectrum_lw = ', use_canopy_full_spectrum_lw 872 print*,'config%do_canopy_gases_sw = ', do_canopy_gases_sw 873 print*,'config%do_canopy_gases_lw = ', do_canopy_gases_lw 874 print*,'config%mono_lw_wavelength = ', mono_lw_wavelength 875 print*,'config%mono_lw_total_od = ', mono_lw_total_od 876 print*,'config%mono_sw_total_od = ', mono_sw_total_od 877 print*,'config%mono_lw_single_scattering_albedo = ', & 878 mono_lw_single_scattering_albedo 879 print*,'config%mono_sw_single_scattering_albedo = ', & 880 mono_sw_single_scattering_albedo 881 print*,'config%mono_lw_asymmetry_factor = ', mono_lw_asymmetry_factor 882 print*,'config%mono_sw_asymmetry_factor = ', mono_sw_asymmetry_factor 883 print*,'config%use_beta_overlap = ', use_beta_overlap 884 print*,'config%cloud_inhom_decorr_scaling = ', cloud_inhom_decorr_scaling 885 print*,'config%clear_to_thick_fraction = ', clear_to_thick_fraction 886 print*,'config%overhead_sun_factor = ', overhead_sun_factor 887 print*,'config%max_gas_od_3d = ', max_gas_od_3d 888 print*,'config%max_cloud_od = ', max_cloud_od 889 print*,'config%max_3d_transfer_rate = ', max_3d_transfer_rate 890 print*,'config%min_cloud_effective_size = ', & 891 max(1.0e-6_jprb,min_cloud_effective_size) 892 print*,'config%overhang_factor = ', encroachment_scaling 893 894 print*,'config%directory_name = ',directory_name 895 print*,'config%cloud_pdf_override_file_name = ',cloud_pdf_override_file_name 896 print*,'config%liq_optics_override_file_name = ',liq_optics_override_file_name 897 print*,'config%ice_optics_override_file_name = ',ice_optics_override_file_name 898 print*,'config%aerosol_optics_override_file_name = ', & 899 aerosol_optics_override_file_name 900 print*,'config%cloud_fraction_threshold = ',cloud_fraction_threshold 901 print*,'config%cloud_mixing_ratio_threshold = ',cloud_mixing_ratio_threshold 902 print*,'config%n_aerosol_types = ',n_aerosol_types 903 print*,'config%do_save_radiative_properties = ',do_save_radiative_properties 904 print*,'config%do_lw_derivatives = ',do_lw_derivatives 905 print*,'config%do_save_spectral_flux = ',do_save_spectral_flux 906 print*,'config%do_save_gpoint_flux = ',do_save_gpoint_flux 907 print*,'config%do_nearest_spectral_sw_albedo = ',do_nearest_spectral_sw_albedo 908 print*,'config%do_nearest_spectral_lw_emiss = ',do_nearest_spectral_lw_emiss 909 print*,'config%sw_albedo_wavelength_bound = ',sw_albedo_wavelength_bound 910 print*,'config%lw_emiss_wavelength_bound = ',lw_emiss_wavelength_bound 911 print*,'config%i_sw_albedo_index = ',i_sw_albedo_index 912 print*,'config%i_lw_emiss_index = ',i_lw_emiss_index 913 print*,'************************************************************************' 914 endif 955 this%do_cloud_aerosol_per_lw_g_point = do_cloud_aerosol_per_lw_g_point 956 this%do_cloud_aerosol_per_sw_g_point = do_cloud_aerosol_per_sw_g_point 957 this%do_weighted_surface_mapping = do_weighted_surface_mapping 958 915 959 if (do_save_gpoint_flux) then 916 960 ! Saving the fluxes every g-point overrides saving as averaged … … 919 963 ! save anything 920 964 this%do_save_spectral_flux = .true. 921 print*,'config%do_save_spectral_flux = .true.'922 965 end if 923 966 … … 925 968 call get_enum_code(liquid_model_name, LiquidModelName, & 926 969 & 'liquid_model_name', this%i_liq_model) 927 print*,'config%i_liq_model =', this%i_liq_model928 970 929 971 ! Determine ice optics model 930 972 call get_enum_code(ice_model_name, IceModelName, & 931 973 & 'ice_model_name', this%i_ice_model) 932 print*,'config%i_ice_model =', this%i_ice_model 974 933 975 ! Determine gas optics model 934 976 call get_enum_code(gas_model_name, GasModelName, & 935 977 & 'gas_model_name', this%i_gas_model) 936 print*,'config%%i_gas_model = ', this%i_gas_model937 978 938 979 ! Determine solvers 939 980 call get_enum_code(sw_solver_name, SolverName, & 940 981 & 'sw_solver_name', this%i_solver_sw) 941 print*,'config%i_solver_sw = ', this%i_solver_sw942 982 call get_enum_code(lw_solver_name, SolverName, & 943 983 & 'lw_solver_name', this%i_solver_lw) 944 print*,'config%i_solver_lw = ', this%i_solver_lw 984 945 985 if (len_trim(sw_encroachment_name) > 1) then 946 986 call get_enum_code(sw_encroachment_name, EncroachmentName, & … … 950 990 call get_enum_code(sw_entrapment_name, EntrapmentName, & 951 991 & 'sw_entrapment_name', this%i_3d_sw_entrapment) 952 print*,'config%i_3d_sw_entrapment = ', this%i_3d_sw_entrapment953 992 end if 954 993 … … 956 995 call get_enum_code(overlap_scheme_name, OverlapName, & 957 996 & 'overlap_scheme_name', this%i_overlap_scheme) 958 print*,'config%i_overlap_scheme = ', this%i_overlap_scheme997 959 998 ! Determine cloud PDF shape 960 999 call get_enum_code(cloud_pdf_shape_name, PdfShapeName, & 961 1000 & 'cloud_pdf_shape_name', this%i_cloud_pdf_shape) 962 print*,'config%i_cloud_pdf_shape = ', this%i_cloud_pdf_shape 1001 963 1002 this%i_aerosol_type_map = 0 964 1003 if (this%use_aerosols) then 965 1004 this%i_aerosol_type_map(1:n_aerosol_types) & 966 1005 & = i_aerosol_type_map(1:n_aerosol_types) 967 print*,'config%i_aerosol_type_map = ', this%i_aerosol_type_map968 1006 end if 969 1007 … … 975 1013 this%do_clouds = .false. 976 1014 end if 977 print*,'config%do_clouds = ', this%do_clouds 1015 1016 if (this%i_gas_model == IGasModelIFSRRTMG & 1017 & .and. (this%use_general_cloud_optics & 1018 & .or. this%use_general_aerosol_optics)) then 1019 if (this%do_sw .and. this%do_cloud_aerosol_per_sw_g_point) then 1020 write(nulout,'(a)') 'Warning: RRTMG SW only supports cloud/aerosol/surface optical properties per band, not per g-point' 1021 this%do_cloud_aerosol_per_sw_g_point = .false. 1022 end if 1023 if (this%do_lw .and. this%do_cloud_aerosol_per_lw_g_point) then 1024 write(nulout,'(a)') 'Warning: RRTMG LW only supports cloud/aerosol/surface optical properties per band, not per g-point' 1025 this%do_cloud_aerosol_per_lw_g_point = .false. 1026 end if 1027 end if 1028 978 1029 979 1030 ! Normal subroutine exit … … 993 1044 subroutine consolidate_config(this) 994 1045 1046 use parkind1, only : jprd 995 1047 use yomhook, only : lhook, dr_hook 996 1048 use radiation_io, only : nulout, nulerr, radiation_abort … … 1026 1078 write(nulerr,'(a)') '*** Error: SPARTACUS/Tripleclouds solvers can only do Exponential-Random overlap' 1027 1079 call radiation_abort('Radiation configuration error') 1080 end if 1081 1082 if (jprb < jprd .and. this%iverbosesetup >= 1 & 1083 & .and. (this%i_solver_sw == ISolverSPARTACUS & 1084 & .or. this%i_solver_lw == ISolverSPARTACUS)) then 1085 write(nulout,'(a)') 'Warning: the SPARTACUS solver may be unstable in single precision' 1086 end if 1087 1088 ! If ecCKD gas optics model is being used set relevant file names 1089 if (this%i_gas_model == IGasModelECCKD) then 1090 1091 ! This gas optics model requires the general cloud and 1092 ! aerosol optics settings 1093 if (.not. this%use_general_cloud_optics) then 1094 write(nulerr,'(a)') '*** Error: ecCKD gas optics model requires general cloud optics' 1095 call radiation_abort('Radiation configuration error') 1096 end if 1097 if (.not. this%use_general_aerosol_optics) then 1098 write(nulerr,'(a)') '*** Error: ecCKD gas optics model requires general aerosol optics' 1099 call radiation_abort('Radiation configuration error') 1100 end if 1101 1102 if (len_trim(this%gas_optics_sw_override_file_name) > 0) then 1103 if (this%gas_optics_sw_override_file_name(1:1) == '/') then 1104 this%gas_optics_sw_file_name = trim(this%gas_optics_sw_override_file_name) 1105 else 1106 this%gas_optics_sw_file_name = trim(this%directory_name) & 1107 & // '/' // trim(this%gas_optics_sw_override_file_name) 1108 end if 1109 else 1110 ! In the IFS, the gas optics files should be specified in 1111 ! ifs/module/radiation_setup.F90, not here 1112 this%gas_optics_sw_file_name = trim(this%directory_name) & 1113 & // "/ecckd-1.0_sw_climate_rgb-32b_ckd-definition.nc" 1114 end if 1115 1116 if (len_trim(this%gas_optics_lw_override_file_name) > 0) then 1117 if (this%gas_optics_lw_override_file_name(1:1) == '/') then 1118 this%gas_optics_lw_file_name = trim(this%gas_optics_lw_override_file_name) 1119 else 1120 this%gas_optics_lw_file_name = trim(this%directory_name) & 1121 & // '/' // trim(this%gas_optics_lw_override_file_name) 1122 end if 1123 else 1124 ! In the IFS, the gas optics files should be specified in 1125 ! ifs/module/radiation_setup.F90, not here 1126 this%gas_optics_lw_file_name = trim(this%directory_name) & 1127 & // "/ecckd-1.0_lw_climate_fsck-32b_ckd-definition.nc" 1128 end if 1028 1129 1029 1130 end if … … 1040 1141 ! In the IFS, the aerosol optics file should be specified in 1041 1142 ! ifs/module/radiation_setup.F90, not here 1042 this%aerosol_optics_file_name & 1043 & = trim(this%directory_name) // "/aerosol_ifs_rrtm_45R2.nc" 1143 if (this%use_general_aerosol_optics) then 1144 this%aerosol_optics_file_name & 1145 & = trim(this%directory_name) // "/aerosol_ifs_48R1.nc" 1146 else 1147 this%aerosol_optics_file_name & 1148 & = trim(this%directory_name) // "/aerosol_ifs_rrtm_46R1_with_NI_AM.nc" 1149 end if 1044 1150 end if 1045 1151 … … 1229 1335 & this%i_gas_model) 1230 1336 call print_logical(' Aerosols are', 'use_aerosols', this%use_aerosols) 1231 call print_logical(' Clouds are', 'do_clouds', this%do_clouds) 1337 if (this%use_aerosols) then 1338 call print_logical(' General aerosol optics', & 1339 & 'use_general_aerosol_optics', this%use_general_aerosol_optics) 1340 end if 1341 if (this%do_clouds) then 1342 write(nulout,'(a)') ' Clouds are ON' 1343 else 1344 write(nulout,'(a)') ' Clouds are OFF' 1345 end if 1346 if (this%do_sw) then 1347 call print_logical(' Do cloud/aerosol/surface SW properties per g-point', & 1348 & 'do_cloud_aerosol_per_sw_g_point', this%do_cloud_aerosol_per_sw_g_point) 1349 end if 1350 if (this%do_lw) then 1351 call print_logical(' Do cloud/aerosol/surface LW properties per g-point', & 1352 & 'do_cloud_aerosol_per_lw_g_point', this%do_cloud_aerosol_per_lw_g_point) 1353 end if 1232 1354 1233 1355 !--------------------------------------------------------------------- … … 1253 1375 & 'do_nearest_spectral_lw_emiss', this%do_nearest_spectral_lw_emiss) 1254 1376 end if 1377 call print_logical(' Planck-weighted surface albedo/emiss mapping', & 1378 & 'do_weighted_surface_mapping', this%do_weighted_surface_mapping) 1379 1255 1380 !--------------------------------------------------------------------- 1256 1381 if (this%do_clouds) then … … 1260 1385 call print_real(' Cloud mixing-ratio threshold', & 1261 1386 & 'cloud_mixing_ratio_threshold', this%cloud_mixing_ratio_threshold) 1262 call print_enum(' Liquid optics scheme is', LiquidModelName, & 1263 & 'i_liq_model',this%i_liq_model) 1264 call print_enum(' Ice optics scheme is', IceModelName, & 1265 & 'i_ice_model',this%i_ice_model) 1266 if (this%i_ice_model == IIceModelFu) then 1267 call print_logical(' Longwave ice optics bug in Fu scheme is', & 1268 & 'do_fu_lw_ice_optics_bug',this%do_fu_lw_ice_optics_bug) 1387 call print_logical(' General cloud optics', & 1388 & 'use_general_cloud_optics', this%use_general_cloud_optics) 1389 if (.not. this%use_general_cloud_optics) then 1390 call print_enum(' Liquid optics scheme is', LiquidModelName, & 1391 & 'i_liq_model',this%i_liq_model) 1392 call print_enum(' Ice optics scheme is', IceModelName, & 1393 & 'i_ice_model',this%i_ice_model) 1394 if (this%i_ice_model == IIceModelFu) then 1395 call print_logical(' Longwave ice optics bug in Fu scheme is', & 1396 & 'do_fu_lw_ice_optics_bug',this%do_fu_lw_ice_optics_bug) 1397 end if 1269 1398 end if 1270 1399 call print_enum(' Cloud overlap scheme is', OverlapName, & … … 1360 1489 & 'overhang_factor', this%overhang_factor) 1361 1490 end if 1491 1492 else if (this%i_solver_sw == ISolverMcICA & 1493 & .or. this%i_solver_lw == ISolverMcICA) then 1494 call print_logical(' Use vectorizable McICA cloud generator', & 1495 & 'use_vectorizable_generator', this%use_vectorizable_generator) 1362 1496 end if 1363 1497 … … 1383 1517 use parkind1, only : jprb 1384 1518 use radiation_io, only : nulout, nulerr, radiation_abort 1519 use radiation_spectral_definition, only : SolarReferenceTemperature 1385 1520 1386 1521 class(config_type), intent(in) :: this … … 1396 1531 character(len=*), optional, intent(in) :: weighting_name 1397 1532 1533 real(jprb), allocatable :: mapping(:,:) 1534 1398 1535 ! Internally we deal with wavenumber 1399 1536 real(jprb) :: wavenumber1, wavenumber2 ! cm-1 1400 1537 1538 real(jprb) :: wavenumber1_band, wavenumber2_band ! cm-1 1539 1401 1540 integer :: jband ! Loop index for spectral band 1402 1541 1403 1542 if (this%n_bands_sw <= 0) then 1404 1543 write(nulerr,'(a)') '*** Error: get_sw_weights called before number of shortwave bands set' 1405 call radiation_abort( )1544 call radiation_abort('Radiation configuration error') 1406 1545 end if 1407 1546 … … 1410 1549 wavenumber2 = 0.01_jprb / wavelength1 1411 1550 1551 call this%gas_optics_sw%spectral_def%calc_mapping_from_bands(SolarReferenceTemperature, & 1552 & [wavelength1, wavelength2], [1, 2, 3], mapping, & 1553 & use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point), use_fluxes=.true.) 1554 1555 ! "mapping" now contains a 3*nband matrix, where mapping(2,:) 1556 ! contains the weights of interest. We now find the non-zero weights 1412 1557 nweights = 0 1413 1414 do jband = 1,this%n_bands_sw 1415 if (wavenumber1 < this%wavenumber2_sw(jband) & 1416 & .and. wavenumber2 > this%wavenumber1_sw(jband)) then 1558 do jband = 1,size(mapping,2) 1559 if (mapping(2,jband) > 0.0_jprb) then 1417 1560 nweights = nweights+1 1418 iband(nweights) = jband 1419 weight(nweights) = (min(wavenumber2,this%wavenumber2_sw(jband)) & 1420 & - max(wavenumber1,this%wavenumber1_sw(jband))) & 1421 & / (this%wavenumber2_sw(jband) - this%wavenumber1_sw(jband)) 1561 iband(nweights) = jband; 1562 weight(nweights) = mapping(2,jband) 1422 1563 end if 1423 1564 end do … … 1426 1567 write(nulerr,'(a,e8.4,a,e8.4,a)') '*** Error: wavelength range ', & 1427 1568 & wavelength1, ' to ', wavelength2, ' m is outside shortwave band' 1428 call radiation_abort( )1569 call radiation_abort('Radiation configuration error') 1429 1570 else if (this%iverbosesetup >= 2 .and. present(weighting_name)) then 1430 1571 write(nulout,'(a,a,a,f6.0,a,f6.0,a)') 'Spectral weights for ', & 1431 1572 & weighting_name, ' (', wavenumber1, ' to ', & 1432 1573 & wavenumber2, ' cm-1):' 1433 do jband = 1, nweights 1434 write(nulout, '(a,i0,a,f6.0,a,f6.0,a,f8.4)') ' Shortwave band ', & 1435 & iband(jband), ' (', this%wavenumber1_sw(iband(jband)), ' to ', & 1436 & this%wavenumber2_sw(iband(jband)), ' cm-1): ', weight(jband) 1437 end do 1574 if (this%do_cloud_aerosol_per_sw_g_point) then 1575 do jband = 1, nweights 1576 write(nulout, '(a,i0,a,f8.4)') ' Shortwave g point ', iband(jband), ': ', weight(jband) 1577 end do 1578 else 1579 do jband = 1, nweights 1580 wavenumber1_band = this%gas_optics_sw%spectral_def%wavenumber1_band(iband(jband)) 1581 wavenumber2_band = this%gas_optics_sw%spectral_def%wavenumber2_band(iband(jband)) 1582 write(nulout, '(a,i0,a,f6.0,a,f6.0,a,f8.4)') ' Shortwave band ', & 1583 & iband(jband), ' (', wavenumber1_band, ' to ', & 1584 & wavenumber2_band, ' cm-1): ', weight(jband) 1585 end do 1586 end if 1438 1587 end if 1439 1588 … … 1452 1601 1453 1602 use radiation_io, only : nulerr, radiation_abort 1603 use radiation_spectral_definition, only : SolarReferenceTemperature 1454 1604 1455 1605 class(config_type), intent(inout) :: this … … 1467 1617 write(nulerr,'(a,i0,a,i0)') '*** Error: ', ninterval, & 1468 1618 & ' albedo intervals exceeds maximum of ', NMaxAlbedoIntervals 1469 call radiation_abort( );1619 call radiation_abort('Radiation configuration error') 1470 1620 end if 1471 1621 … … 1480 1630 this%i_sw_albedo_index(ninterval+1:) = 0 1481 1631 1632 ! If this routine is called before setup_radiation then the 1633 ! spectral intervals are not yet known 1634 ! consolidate_sw_albedo_intervals is called later. Otherwise it 1635 ! is called immediately and overwrites any existing mapping. 1482 1636 if (this%is_consolidated) then 1483 call this%consolidate_intervals(.true., & 1484 & this%do_nearest_spectral_sw_albedo, & 1485 & this%sw_albedo_wavelength_bound, this%i_sw_albedo_index, & 1486 & this%wavenumber1_sw, this%wavenumber2_sw, & 1487 & this%i_albedo_from_band_sw, this%sw_albedo_weights) 1637 call this%consolidate_sw_albedo_intervals 1488 1638 end if 1489 1639 … … 1497 1647 1498 1648 use radiation_io, only : nulerr, radiation_abort 1649 use radiation_spectral_definition, only : TerrestrialReferenceTemperature 1499 1650 1500 1651 class(config_type), intent(inout) :: this … … 1512 1663 write(nulerr,'(a,i0,a,i0)') '*** Error: ', ninterval, & 1513 1664 & ' emissivity intervals exceeds maximum of ', NMaxAlbedoIntervals 1514 call radiation_abort( );1665 call radiation_abort('Radiation configuration error') 1515 1666 end if 1516 1667 … … 1526 1677 1527 1678 if (this%is_consolidated) then 1528 call this%consolidate_intervals(.false., & 1529 & this%do_nearest_spectral_lw_emiss, & 1530 & this%lw_emiss_wavelength_bound, this%i_lw_emiss_index, & 1531 & this%wavenumber1_lw, this%wavenumber2_lw, & 1532 & this%i_emiss_from_band_lw, this%lw_emiss_weights) 1679 call this%consolidate_lw_emiss_intervals 1533 1680 end if 1534 1681 … … 1537 1684 1538 1685 !--------------------------------------------------------------------- 1539 ! This routine consolidates either the input shortwave albedo 1540 ! intervals with the shortwave bands, or the input longwave 1541 ! emissivity intervals with the longwave bands, depending on the 1542 ! arguments provided. 1543 subroutine consolidate_intervals(this, is_sw, do_nearest, & 1544 & wavelength_bound, i_intervals, wavenumber1, wavenumber2, & 1545 & i_mapping, weights) 1546 1547 use radiation_io, only : nulout, nulerr, radiation_abort 1686 ! Set the wavelengths (m) at which monochromatic aerosol properties 1687 ! are required. This routine must be called before consolidation of 1688 ! settings. 1689 subroutine set_aerosol_wavelength_mono(this, wavelength_mono) 1690 1691 use radiation_io, only : nulerr, radiation_abort 1692 1693 class(config_type), intent(inout) :: this 1694 real(jprb), intent(in) :: wavelength_mono(:) 1695 1696 if (this%is_consolidated) then 1697 write(nulerr,'(a)') '*** Errror: set_aerosol_wavelength_mono must be called before setup_radiation' 1698 call radiation_abort('Radiation configuration error') 1699 end if 1700 1701 if (allocated(this%aerosol_optics%wavelength_mono)) then 1702 deallocate(this%aerosol_optics%wavelength_mono) 1703 end if 1704 allocate(this%aerosol_optics%wavelength_mono(size(wavelength_mono))) 1705 this%aerosol_optics%wavelength_mono = wavelength_mono 1706 1707 end subroutine set_aerosol_wavelength_mono 1708 1709 1710 !--------------------------------------------------------------------- 1711 ! Consolidate the surface shortwave albedo intervals with the 1712 ! band/g-point intervals 1713 subroutine consolidate_sw_albedo_intervals(this) 1714 1715 use radiation_io, only : nulout 1716 use radiation_spectral_definition, only : SolarReferenceTemperature 1548 1717 1549 1718 class(config_type), intent(inout) :: this 1550 ! Is this the shortwave? Otherwise longwave 1551 logical, intent(in) :: is_sw 1552 ! Do we find the nearest albedo interval to the centre of each 1553 ! band, or properly weight the contributions? This can be modified 1554 ! if there is only one albedo intervals. 1555 logical, intent(inout) :: do_nearest 1556 ! Monotonically increasing wavelength bounds between intervals, 1557 ! not including the outer bounds (which are assumed to be zero and 1558 ! infinity) 1559 real(jprb), intent(in) :: wavelength_bound(NMaxAlbedoIntervals-1) 1560 ! The albedo band indices corresponding to each interval 1561 integer, intent(in) :: i_intervals(NMaxAlbedoIntervals) 1562 ! Start and end wavenumber bounds for the ecRad bands (cm-1) 1563 real(jprb), intent(in) :: wavenumber1(:), wavenumber2(:) 1564 1565 ! if do_nearest is TRUE then the result is expressed in i_mapping, 1566 ! which will be allocated to have the same length as wavenumber1, 1567 ! and contain the index of the albedo interval corresponding to 1568 ! that band 1569 integer, allocatable, intent(inout) :: i_mapping(:) 1570 ! ...otherwise the result is expressed in "weights", of 1571 ! size(n_intervals, n_bands) containing how much of each interval 1572 ! contributes to each band. 1573 real(jprb), allocatable, intent(inout) :: weights(:,:) 1574 1575 ! Number and loop index of ecRad bands 1576 integer :: nband, jband 1577 ! Number and index of albedo/emissivity intervals 1578 integer :: ninterval, iinterval 1579 ! Sometimes an albedo or emissivity value will be used in more 1580 ! than one interval, so nvalue indicates how many values will 1581 ! actually be provided 1582 integer :: nvalue 1583 ! Wavenumber bounds of the albedo/emissivity interval 1584 real(jprb) :: wavenumber1_albedo, wavenumber2_albedo 1585 ! Reciprocal of the wavenumber range of the ecRad band 1586 real(jprb) :: recip_dwavenumber ! cm 1587 ! Midpoint/bound of wavenumber band 1588 real(jprb) :: wavenumber_mid, wavenumber_bound ! cm-1 1589 1590 nband = size(wavenumber1) 1719 1720 integer :: ninterval, jint, jband 1591 1721 1592 1722 ! Count the number of albedo/emissivity intervals 1593 1723 ninterval = 0 1594 do iinterval= 1,NMaxAlbedoIntervals1595 if ( i_intervals(iinterval) > 0) then1596 ninterval = iinterval1724 do jint = 1,NMaxAlbedoIntervals 1725 if (this%i_sw_albedo_index(jint) > 0) then 1726 ninterval = jint 1597 1727 else 1598 1728 exit … … 1600 1730 end do 1601 1731 1602 if (ninterval < 2) then 1603 ! Zero or one albedo/emissivity intervals found, so we index all 1604 ! bands to one interval 1605 if (allocated(i_mapping)) then 1606 deallocate(i_mapping) 1607 end if 1608 allocate(i_mapping(nband)) 1609 i_mapping(:) = 1 1610 do_nearest = .true. 1732 if (ninterval < 1) then 1733 ! The user has not specified shortwave albedo bands - assume 1734 ! only one 1611 1735 ninterval = 1 1612 nvalue = 1 1613 else 1614 ! Check wavelength is monotonically increasing 1615 do jband = 2,ninterval-1 1616 if (wavelength_bound(jband) <= wavelength_bound(jband-1)) then 1617 if (is_sw) then 1618 write(nulerr, '(a,a)') '*** Error: wavelength bounds for shortwave albedo intervals ', & 1619 & 'must be monotonically increasing' 1620 else 1621 write(nulerr, '(a,a)') '*** Error: wavelength bounds for longwave emissivity intervals ', & 1622 & 'must be monotonically increasing' 1623 end if 1624 call radiation_abort() 1625 end if 1626 end do 1627 1628 ! What is the maximum index, indicating the number of 1629 ! albedo/emissivity values to expect? 1630 nvalue = maxval(i_intervals(1:ninterval)) 1631 1632 if (do_nearest) then 1633 ! Simpler nearest-neighbour mapping from band to 1634 ! albedo/emissivity interval 1635 if (allocated(i_mapping)) then 1636 deallocate(i_mapping) 1637 end if 1638 allocate(i_mapping(nband)) 1639 1640 ! Loop over bands 1641 do jband = 1,nband 1642 ! Compute mid-point of band in wavenumber space (cm-1) 1643 wavenumber_mid = 0.5_jprb * (wavenumber1(jband) & 1644 & + wavenumber2(jband)) 1645 iinterval = 1 1646 ! Convert wavelength (m) into wavenumber (cm-1) at the lower 1647 ! bound of the albedo interval 1648 wavenumber_bound = 0.01_jprb / wavelength_bound(iinterval) 1649 ! Find the albedo interval that has the largest overlap with 1650 ! the band; this approach assumes that the albedo intervals 1651 ! are larger than the spectral bands 1652 do while (wavenumber_bound >= wavenumber_mid & 1653 & .and. iinterval < ninterval) 1654 iinterval = iinterval + 1 1655 if (iinterval < ninterval) then 1656 wavenumber_bound = 0.01_jprb / wavelength_bound(iinterval) 1657 else 1658 ! For the last interval there is no lower bound 1659 wavenumber_bound = 0.0_jprb 1660 end if 1661 end do 1662 ! Save the index of the band corresponding to the albedo 1663 ! interval and move onto the next band 1664 i_mapping(jband) = i_intervals(iinterval) 1665 end do 1666 else 1667 ! More accurate weighting 1668 if (allocated(weights)) then 1669 deallocate(weights) 1670 end if 1671 allocate(weights(nvalue,nband)) 1672 weights(:,:) = 0.0_jprb 1673 1674 ! Loop over bands 1675 do jband = 1,nband 1676 recip_dwavenumber = 1.0_jprb / (wavenumber2(jband) & 1677 & - wavenumber1(jband)) 1678 ! Find the first overlapping albedo band 1679 iinterval = 1 1680 ! Convert wavelength (m) into wavenumber (cm-1) at the lower 1681 ! bound (in wavenumber space) of the albedo/emissivty interval 1682 wavenumber1_albedo = 0.01_jprb / wavelength_bound(iinterval) 1683 do while (wavenumber1_albedo >= wavenumber2(jband) & 1684 & .and. iinterval < ninterval) 1685 iinterval = iinterval + 1 1686 wavenumber1_albedo = 0.01_jprb / wavelength_bound(iinterval) 1687 end do 1688 1689 wavenumber2_albedo = wavenumber2(jband) 1690 1691 ! Add all overlapping bands 1692 do while (wavenumber2_albedo > wavenumber1(jband) & 1693 & .and. iinterval <= ninterval) 1694 weights(i_intervals(iinterval),jband) & 1695 & = weights(i_intervals(iinterval),jband) & 1696 & + recip_dwavenumber & 1697 & * (min(wavenumber2_albedo,wavenumber2(jband)) & 1698 & - max(wavenumber1_albedo,wavenumber1(jband))) 1699 wavenumber2_albedo = wavenumber1_albedo 1700 iinterval = iinterval + 1 1701 if (iinterval < ninterval) then 1702 wavenumber1_albedo = 0.01_jprb / wavelength_bound(iinterval) 1703 else 1704 wavenumber1_albedo = 0.0_jprb 1705 end if 1706 end do 1707 end do 1708 end if 1709 end if 1710 1711 ! Define how many bands to use for reporting surface downwelling 1712 ! fluxes for canopy radiation scheme 1713 if (is_sw) then 1736 this%i_sw_albedo_index(1) = 1 1737 this%i_sw_albedo_index(2:) = 0 1714 1738 if (this%use_canopy_full_spectrum_sw) then 1715 1739 this%n_canopy_bands_sw = this%n_g_sw 1716 1740 else 1717 this%n_canopy_bands_sw = nvalue 1741 this%n_canopy_bands_sw = 1 1742 end if 1743 else 1744 if (this%use_canopy_full_spectrum_sw) then 1745 this%n_canopy_bands_sw = this%n_g_sw 1746 else 1747 this%n_canopy_bands_sw = maxval(this%i_sw_albedo_index(1:ninterval)) 1748 end if 1749 end if 1750 1751 if (this%do_weighted_surface_mapping) then 1752 call this%gas_optics_sw%spectral_def%calc_mapping_from_bands(SolarReferenceTemperature, & 1753 & this%sw_albedo_wavelength_bound(1:ninterval-1), this%i_sw_albedo_index(1:ninterval), & 1754 & this%sw_albedo_weights, use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point)) 1755 else 1756 ! Weight each wavenumber equally as in IFS Cycles 48 and earlier 1757 call this%gas_optics_sw%spectral_def%calc_mapping_from_bands(-1.0_jprb, & 1758 & this%sw_albedo_wavelength_bound(1:ninterval-1), this%i_sw_albedo_index(1:ninterval), & 1759 & this%sw_albedo_weights, use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point)) 1760 end if 1761 1762 ! Legacy method uses input band with largest weight 1763 if (this%do_nearest_spectral_sw_albedo) then 1764 allocate(this%i_albedo_from_band_sw(this%n_bands_sw)) 1765 this%i_albedo_from_band_sw = maxloc(this%sw_albedo_weights, dim=1) 1766 end if 1767 1768 if (this%iverbosesetup >= 2) then 1769 write(nulout, '(a)') 'Surface shortwave albedo' 1770 if (.not. this%do_nearest_spectral_sw_albedo) then 1771 call this%gas_optics_sw%spectral_def%print_mapping_from_bands(this%sw_albedo_weights, & 1772 & use_bands=(.not. this%do_cloud_aerosol_per_sw_g_point)) 1773 else if (ninterval <= 1) then 1774 write(nulout, '(a)') 'All shortwave bands will use the same albedo' 1775 else 1776 write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_albedo_from_band_sw), & 1777 & ' shortwave intervals to albedo intervals:' 1778 do jband = 1,size(this%i_albedo_from_band_sw) 1779 write(nulout,'(a,i0)',advance='no') ' ', this%i_albedo_from_band_sw(jband) 1780 end do 1781 write(nulout, '()') 1782 end if 1783 end if 1784 1785 end subroutine consolidate_sw_albedo_intervals 1786 1787 1788 !--------------------------------------------------------------------- 1789 ! Consolidate the surface longwave emissivity intervals with the 1790 ! band/g-point intervals 1791 subroutine consolidate_lw_emiss_intervals(this) 1792 1793 use radiation_io, only : nulout 1794 use radiation_spectral_definition, only : TerrestrialReferenceTemperature 1795 1796 class(config_type), intent(inout) :: this 1797 1798 integer :: ninterval, jint, jband 1799 1800 ! Count the number of albedo/emissivity intervals 1801 ninterval = 0 1802 do jint = 1,NMaxAlbedoIntervals 1803 if (this%i_lw_emiss_index(jint) > 0) then 1804 ninterval = jint 1805 else 1806 exit 1807 end if 1808 end do 1809 1810 if (ninterval < 1) then 1811 ! The user has not specified longwave emissivity bands - assume 1812 ! only one 1813 ninterval = 1 1814 this%i_lw_emiss_index(1) = 1 1815 this%i_lw_emiss_index(2:) = 0 1816 if (this%use_canopy_full_spectrum_sw) then 1817 this%n_canopy_bands_lw = this%n_g_lw 1818 else 1819 this%n_canopy_bands_lw = 1 1718 1820 end if 1719 1821 else … … 1721 1823 this%n_canopy_bands_lw = this%n_g_lw 1722 1824 else 1723 this%n_canopy_bands_lw = nvalue 1724 end if 1825 this%n_canopy_bands_lw = maxval(this%i_lw_emiss_index(1:ninterval)) 1826 end if 1827 end if 1828 1829 if (this%do_weighted_surface_mapping) then 1830 call this%gas_optics_lw%spectral_def%calc_mapping_from_bands(TerrestrialReferenceTemperature, & 1831 & this%lw_emiss_wavelength_bound(1:ninterval-1), this%i_lw_emiss_index(1:ninterval), & 1832 & this%lw_emiss_weights, use_bands=(.not. this%do_cloud_aerosol_per_lw_g_point)) 1833 else 1834 ! Weight each wavenumber equally as in IFS Cycles 48 and earlier 1835 call this%gas_optics_lw%spectral_def%calc_mapping_from_bands(-1.0_jprb, & 1836 & this%lw_emiss_wavelength_bound(1:ninterval-1), this%i_lw_emiss_index(1:ninterval), & 1837 & this%lw_emiss_weights, use_bands=(.not. this%do_cloud_aerosol_per_lw_g_point)) 1838 end if 1839 1840 ! Legacy method uses input band with largest weight 1841 if (this%do_nearest_spectral_lw_emiss) then 1842 allocate(this%i_emiss_from_band_lw(this%n_bands_lw)) 1843 this%i_emiss_from_band_lw = maxloc(this%lw_emiss_weights, dim=1) 1725 1844 end if 1726 1845 1727 1846 if (this%iverbosesetup >= 2) then 1728 if (.not. do_nearest) then 1729 if (is_sw) then 1730 write(nulout, '(a,i0,a,i0,a)') 'Weighting of ', nvalue, ' albedo values in ', & 1731 & nband, ' shortwave bands (wavenumber ranges in cm-1):' 1732 else 1733 write(nulout, '(a,i0,a,i0,a)') 'Weighting of ', nvalue, ' emissivity values in ', & 1734 & nband, ' longwave bands (wavenumber ranges in cm-1):' 1735 end if 1736 do jband = 1,nband 1737 write(nulout,'(i6,a,i6,a)',advance='no') nint(wavenumber1(jband)), ' to', & 1738 & nint(wavenumber2(jband)), ':' 1739 do iinterval = 1,nvalue 1740 write(nulout,'(f5.2)',advance='no') weights(iinterval,jband) 1741 end do 1742 write(nulout, '()') 1743 end do 1847 write(nulout, '(a)') 'Surface longwave emissivity' 1848 if (.not. this%do_nearest_spectral_lw_emiss) then 1849 call this%gas_optics_lw%spectral_def%print_mapping_from_bands(this%lw_emiss_weights, & 1850 & use_bands=(.not. this%do_cloud_aerosol_per_lw_g_point)) 1744 1851 else if (ninterval <= 1) then 1745 if (is_sw) then 1746 write(nulout, '(a)') 'All shortwave bands will use the same albedo' 1747 else 1748 write(nulout, '(a)') 'All longwave bands will use the same emissivty' 1749 end if 1852 write(nulout, '(a)') 'All longwave bands will use the same emissivty' 1750 1853 else 1751 if (is_sw) then 1752 write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', nband, & 1753 & ' shortwave bands to albedo intervals:' 1754 else 1755 write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', nband, & 1756 & ' longwave bands to emissivity intervals:' 1757 end if 1758 do jband = 1,nband 1759 write(nulout,'(a,i0)',advance='no') ' ', i_mapping(jband) 1854 write(nulout, '(a,i0,a)',advance='no') 'Mapping from ', size(this%i_emiss_from_band_lw), & 1855 & ' longwave intervals to emissivity intervals:' 1856 do jband = 1,size(this%i_emiss_from_band_lw) 1857 write(nulout,'(a,i0)',advance='no') ' ', this%i_emiss_from_band_lw(jband) 1760 1858 end do 1761 1859 write(nulout, '()') … … 1763 1861 end if 1764 1862 1765 end subroutine consolidate_ intervals1863 end subroutine consolidate_lw_emiss_intervals 1766 1864 1767 1865 … … 1866 1964 end subroutine print_enum 1867 1965 1868 1869 !---------------------------------------------------------------------1870 ! Return .true. if 1D allocatable array "var" is out of physical1871 ! range specified by boundmin and boundmax, and issue a warning.1872 ! "do_fix" determines whether erroneous values are fixed to lie1873 ! within the physical range. To check only a subset of the array,1874 ! specify i1 and i2 for the range.1875 function out_of_bounds_1d(var, var_name, boundmin, boundmax, do_fix, i1, i2) result (is_bad)1876 1877 use radiation_io, only : nulout1878 1879 real(jprb), allocatable, intent(inout) :: var(:)1880 character(len=*), intent(in) :: var_name1881 real(jprb), intent(in) :: boundmin, boundmax1882 logical, intent(in) :: do_fix1883 integer, optional, intent(in) :: i1, i21884 1885 logical :: is_bad1886 1887 real(jprb) :: varmin, varmax1888 1889 is_bad = .false.1890 1891 if (allocated(var)) then1892 1893 if (present(i1) .and. present(i2)) then1894 varmin = minval(var(i1:i2))1895 varmax = maxval(var(i1:i2))1896 else1897 varmin = minval(var)1898 varmax = maxval(var)1899 end if1900 1901 if (varmin < boundmin .or. varmax > boundmax) then1902 write(nulout,'(a,a,a,g12.4,a,g12.4,a,g12.4,a,g12.4)',advance='no') &1903 & '*** Warning: ', var_name, ' range', varmin, ' to', varmax, &1904 & ' is out of physical range', boundmin, 'to', boundmax1905 is_bad = .true.1906 if (do_fix) then1907 if (present(i1) .and. present(i2)) then1908 var(i1:i2) = max(boundmin, min(boundmax, var(i1:i2)))1909 else1910 var = max(boundmin, min(boundmax, var))1911 end if1912 write(nulout,'(a)') ': corrected'1913 else1914 write(nulout,'(1x)')1915 end if1916 end if1917 1918 end if1919 1920 end function out_of_bounds_1d1921 1922 1923 !---------------------------------------------------------------------1924 ! Return .true. if 2D allocatable array "var" is out of physical1925 ! range specified by boundmin and boundmax, and issue a warning. To1926 ! check only a subset of the array, specify i1 and i2 for the range1927 ! of the first dimension and j1 and j2 for the range of the second.1928 function out_of_bounds_2d(var, var_name, boundmin, boundmax, do_fix, &1929 & i1, i2, j1, j2) result (is_bad)1930 1931 use radiation_io, only : nulout1932 1933 real(jprb), allocatable, intent(inout) :: var(:,:)1934 character(len=*), intent(in) :: var_name1935 real(jprb), intent(in) :: boundmin, boundmax1936 logical, intent(in) :: do_fix1937 integer, optional, intent(in) :: i1, i2, j1, j21938 1939 ! Local copies of indices1940 integer :: ii1, ii2, jj1, jj21941 1942 logical :: is_bad1943 1944 real(jprb) :: varmin, varmax1945 1946 is_bad = .false.1947 1948 if (allocated(var)) then1949 1950 if (present(i1) .and. present(i2)) then1951 ii1 = i11952 ii2 = i21953 else1954 ii1 = lbound(var,1)1955 ii2 = ubound(var,1)1956 end if1957 if (present(j1) .and. present(j2)) then1958 jj1 = j11959 jj2 = j21960 else1961 jj1 = lbound(var,2)1962 jj2 = ubound(var,2)1963 end if1964 varmin = minval(var(ii1:ii2,jj1:jj2))1965 varmax = maxval(var(ii1:ii2,jj1:jj2))1966 1967 if (varmin < boundmin .or. varmax > boundmax) then1968 write(nulout,'(a,a,a,g12.4,a,g12.4,a,g12.4,a,g12.4)',advance='no') &1969 & '*** Warning: ', var_name, ' range', varmin, ' to', varmax,&1970 & ' is out of physical range', boundmin, 'to', boundmax1971 is_bad = .true.1972 if (do_fix) then1973 var(ii1:ii2,jj1:jj2) = max(boundmin, min(boundmax, var(ii1:ii2,jj1:jj2)))1974 write(nulout,'(a)') ': corrected'1975 else1976 write(nulout,'(1x)')1977 end if1978 end if1979 1980 end if1981 1982 end function out_of_bounds_2d1983 1984 1985 !---------------------------------------------------------------------1986 ! Return .true. if 3D allocatable array "var" is out of physical1987 ! range specified by boundmin and boundmax, and issue a warning. To1988 ! check only a subset of the array, specify i1 and i2 for the range1989 ! of the first dimension, j1 and j2 for the second and k1 and k2 for1990 ! the third.1991 function out_of_bounds_3d(var, var_name, boundmin, boundmax, do_fix, &1992 & i1, i2, j1, j2, k1, k2) result (is_bad)1993 1994 use radiation_io, only : nulout1995 1996 real(jprb), allocatable, intent(inout) :: var(:,:,:)1997 character(len=*), intent(in) :: var_name1998 real(jprb), intent(in) :: boundmin, boundmax1999 logical, intent(in) :: do_fix2000 integer, optional, intent(in) :: i1, i2, j1, j2, k1, k22001 2002 ! Local copies of indices2003 integer :: ii1, ii2, jj1, jj2, kk1, kk22004 2005 logical :: is_bad2006 2007 real(jprb) :: varmin, varmax2008 2009 is_bad = .false.2010 2011 if (allocated(var)) then2012 2013 if (present(i1) .and. present(i2)) then2014 ii1 = i12015 ii2 = i22016 else2017 ii1 = lbound(var,1)2018 ii2 = ubound(var,1)2019 end if2020 if (present(j1) .and. present(j2)) then2021 jj1 = j12022 jj2 = j22023 else2024 jj1 = lbound(var,2)2025 jj2 = ubound(var,2)2026 end if2027 if (present(k1) .and. present(k2)) then2028 kk1 = k12029 kk2 = k22030 else2031 kk1 = lbound(var,3)2032 kk2 = ubound(var,3)2033 end if2034 varmin = minval(var(ii1:ii2,jj1:jj2,kk1:kk2))2035 varmax = maxval(var(ii1:ii2,jj1:jj2,kk1:kk2))2036 2037 if (varmin < boundmin .or. varmax > boundmax) then2038 write(nulout,'(a,a,a,g12.4,a,g12.4,a,g12.4,a,g12.4)',advance='no') &2039 & '*** Warning: ', var_name, ' range', varmin, ' to', varmax,&2040 & ' is out of physical range', boundmin, 'to', boundmax2041 is_bad = .true.2042 if (do_fix) then2043 var(ii1:ii2,jj1:jj2,kk1:kk2) = max(boundmin, min(boundmax, &2044 & var(ii1:ii2,jj1:jj2,kk1:kk2)))2045 write(nulout,'(a)') ': corrected'2046 else2047 write(nulout,'(1x)')2048 end if2049 end if2050 2051 end if2052 2053 end function out_of_bounds_3d2054 2055 2056 1966 end module radiation_config
Note: See TracChangeset
for help on using the changeset viewer.