 Timestamp:
 Mar 10, 2020 8:15:32 PM (11 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/radiation_model_mod.f90
r4442 r4452 28 28 !  29 29 ! $Id$ 30 !  Change order of dimension in surface arrays %frac, %emissivity and %albedo 30 ! Bugfix in calc_albedo 31 ! 32 ! 4442 20200304 19:21:13Z suehring 33 !  Change order of dimension in surface arrays %frac, %emissivity and %albedo 31 34 ! to allow for better vectorization in the radiation interactions. 32 35 !  Minor formatting issues 33 ! 36 ! 34 37 ! 4441 20200304 19:20:35Z suehring 35 38 ! bugfixes: cppdirectives for serial mode moved, small changes to get serial mode compiled 36 ! 39 ! 37 40 ! 4400 20200210 20:32:41Z suehring 38 41 ! Initialize radiation arrays with zero 39 ! 42 ! 40 43 ! 4392 20200131 16:14:57Z pavelkrc 41 44 !  Add debug tracing of large radiative fluxes (option trace_fluxes_above) 42 45 !  Print exact counts of SVF and CSF if debut_output is enabled 43 46 !  Update descriptions of RTM 3.0 and related comments 44 ! 47 ! 45 48 ! 4360 20200107 11:25:50Z suehring 46 49 ! Renamed pc_heating_rate, pc_transpiration_rate, pc_transpiration_rate to 47 50 ! pcm_heating_rate, pcm_latent_rate, pcm_transpiration_rate 48 ! 51 ! 49 52 ! 4340 20191216 08:17:03Z Giersch 50 ! Albedo indices for building_surface_pars are now declared as parameters to 53 ! Albedo indices for building_surface_pars are now declared as parameters to 51 54 ! prevent an error if the gfortran compiler with Werror=unusedvalue is used 52 ! 55 ! 53 56 ! 4291 20191111 12:36:54Z moh.hefny 54 ! Enabled RTM in case of biometeorology even if there is no vertical 57 ! Enabled RTM in case of biometeorology even if there is no vertical 55 58 ! surfaces or 3D vegetation in the domain 56 ! 59 ! 57 60 ! 4286 20191030 16:01:14Z resler 58 61 !  Fix wrong treating of time_rad during interpolation in radiation_model_mod … … 62 65 ! 4271 20191023 10:46:41Z maronga 63 66 ! Bugfix: missing parentheses in calculation of snow albedo 64 ! 67 ! 65 68 ! 4245 20190930 08:40:37Z pavelkrc 66 69 ! Initialize explicit persurface albedos from building_surface_pars 67 ! 70 ! 68 71 ! 4238 20190925 16:06:01Z suehring 69 72 ! Modify check in order to avoid equality comparisons of floating points 70 ! 73 ! 71 74 ! 4227 20190910 18:04:34Z gronemeier 72 75 ! implement new palm_date_time_mod 73 ! 76 ! 74 77 ! 4226 20190910 17:03:24Z suehring 75 78 !  Netcdf input routine for dimension length renamed 76 !  Define time variable for external radiation input relative to time_utc_init 77 ! 79 !  Define time variable for external radiation input relative to time_utc_init 80 ! 78 81 ! 4210 20190902 13:07:09Z suehring 79 82 !  Revise steering of splitting diffuse and direct radiation … … 81 84 !  Optimize mapping of radiation components onto 2D arrays, avoid unnecessary 82 85 ! operations 83 ! 86 ! 84 87 ! 4208 20190902 09:01:07Z suehring 85 88 ! Bugfix in accessing albedo_pars in the clearsky branch 86 ! (merge from branch resler) 87 ! 89 ! (merge from branch resler) 90 ! 88 91 ! 4198 20190829 15:17:48Z gronemeier 89 92 ! Prohibit execution of radiation model if rotation_angle is not zero 90 ! 93 ! 91 94 ! 4197 20190829 14:33:32Z suehring 92 95 ! Revise steering of surface albedo initialization when albedo_pars is provided 93 ! 96 ! 94 97 ! 4190 20190827 15:42:37Z suehring 95 ! Implement external radiation forcing also for levelofdetail = 2 98 ! Implement external radiation forcing also for levelofdetail = 2 96 99 ! (horizontally 2D radiation) 97 ! 100 ! 98 101 ! 4188 20190826 14:15:47Z suehring 99 102 ! Minor adjustment in error message 100 ! 103 ! 101 104 ! 4187 20190826 12:43:15Z suehring 102 !  Take external radiation from root domain dynamic input if not provided for 105 !  Take external radiation from root domain dynamic input if not provided for 103 106 ! each nested domain 104 107 !  Combine MPI_ALLREDUCE calls to reduce mpi overhead 105 ! 108 ! 106 109 ! 4182 20190822 15:20:23Z scharf 107 110 ! Corrected "Former revisions" section 108 ! 111 ! 109 112 ! 4179 20190821 11:16:12Z suehring 110 113 ! Remove debug prints 111 ! 114 ! 112 115 ! 4178 20190821 11:13:06Z suehring 113 116 ! External radiation forcing implemented. 114 ! 117 ! 115 118 ! 4168 20190816 13:50:17Z suehring 116 119 ! Replace function get_topography_top_index by topo_top_ind 117 ! 120 ! 118 121 ! 4157 20190814 09:19:12Z suehring 119 122 ! Give informative message on raytracing distance only by core zero 120 ! 123 ! 121 124 ! 4148 20190808 11:26:00Z suehring 122 125 ! Comments added 123 ! 126 ! 124 127 ! 4134 20190802 18:39:57Z suehring 125 128 ! Bugfix in formatted write statement 126 ! 129 ! 127 130 ! 4127 20190730 14:47:10Z suehring 128 131 ! Remove unused pch_index (merge from branch resler) 129 ! 132 ! 130 133 ! 4089 20190711 14:30:27Z suehring 131 134 !  Correct level 2 initialization of spectral albedos in rrtmg branch, long and 132 ! shortwave albedos were mixedup. 135 ! shortwave albedos were mixedup. 133 136 !  Change order of albedo_pars so that it is now consistent with the defined 134 ! order of albedo_pars in PIDS 135 ! 137 ! order of albedo_pars in PIDS 138 ! 136 139 ! 4069 20190701 14:05:51Z Giersch 137 ! Masked output running index mid has been introduced as a local variable to 140 ! Masked output running index mid has been introduced as a local variable to 138 141 ! avoid runtime error (Loop variable has been modified) in time_integration 139 ! 142 ! 140 143 ! 4067 20190701 13:29:25Z suehring 141 144 ! Bugfix, pass dummy string to MPI_INFO_SET (J. Resler) 142 ! 145 ! 143 146 ! 4039 20190618 10:32:41Z suehring 144 147 ! Bugfix for masked data output 145 ! 148 ! 146 149 ! 4008 20190530 09:50:11Z moh.hefny 147 150 ! Bugfix in check variable when a variable's string is less than 3 … … 150 153 ! 151 154 ! 3992 20190522 16:49:38Z suehring 152 ! Bugfix in rrtmg radiation branch in a nested run when the lowest prognistic 153 ! grid points in a child domain are all inside topography 154 ! 155 ! Bugfix in rrtmg radiation branch in a nested run when the lowest prognistic 156 ! grid points in a child domain are all inside topography 157 ! 155 158 ! 3987 20190522 09:52:13Z kanani 156 159 ! Introduce alternative switch for debug output during timestepping 157 ! 160 ! 158 161 ! 3943 20190502 09:50:41Z maronga 159 162 ! Missing blank characteer added. 160 ! 163 ! 161 164 ! 3900 20190416 15:17:43Z suehring 162 165 ! Fixed initialization problem 163 ! 166 ! 164 167 ! 3885 20190411 11:29:34Z kanani 165 ! Changes related to global restructuring of location messages and introduction 168 ! Changes related to global restructuring of location messages and introduction 166 169 ! of additional debug messages 167 ! 170 ! 168 171 ! 3881 20190410 09:31:22Z suehring 169 172 ! Output of albedo and emissivity moved from USM, bugfixes in initialization 170 173 ! of albedo 171 ! 174 ! 172 175 ! 3861 20190404 06:27:41Z maronga 173 176 ! Bugfix: factor of 4.0 required instead of 3.0 in calculation of rad_lw_out_change_0 174 ! 177 ! 175 178 ! 3859 20190403 20:30:31Z maronga 176 179 ! Added some descriptions 177 ! 180 ! 178 181 ! 3847 20190401 14:51:44Z suehring 179 182 ! Implement check for dt_radiation (must be > 0) 180 ! 183 ! 181 184 ! 3846 20190401 13:55:30Z suehring 182 185 ! unused variable removed 183 ! 186 ! 184 187 ! 3814 20190326 08:40:31Z pavelkrc 185 188 ! Change zenith(0:0) and others to scalar. 186 189 ! Code review. 187 190 ! Rename exported nzu, nzp and related variables due to name conflict 188 ! 191 ! 189 192 ! 3771 20190228 12:19:33Z raasch 190 193 ! rrtmg preprocessor for directives moved/added, save attribute added to temporary 191 194 ! pointers to avoid compiler warnings about outlived pointer targets, 192 195 ! statement added to avoid compiler warning about unused variable 193 ! 196 ! 194 197 ! 3769 20190228 10:16:49Z moh.hefny 195 198 ! removed unused variables and subroutine radiation_radflux_gridbox … … 197 200 ! 3767 20190227 08:18:02Z raasch 198 201 ! unused variable for file index removed from rrdsubroutines parameter list 199 ! 202 ! 200 203 ! 3760 20190221 18:47:35Z moh.hefny 201 204 ! Bugfix: initialized simulated_time before calculating solar position … … 205 208 ! (resler, pavelkrc) 206 209 ! Bugfixes: add further required MRT factors to read/write_svf, 207 ! fix for aggregating view factors to eliminate local noise in reflected 208 ! irradiance at mutually close surfaces (corners, presence of trees) in the 210 ! fix for aggregating view factors to eliminate local noise in reflected 211 ! irradiance at mutually close surfaces (corners, presence of trees) in the 209 212 ! angular discretization scheme. 210 213 ! … … 213 216 ! 214 217 ! 3705 20190129 19:56:39Z suehring 215 ! Make variables that are sampled in virtual measurement module public 216 ! 218 ! Make variables that are sampled in virtual measurement module public 219 ! 217 220 ! 3704 20190129 19:51:41Z suehring 218 221 ! Some interface calls moved to module_interface + cleanup 219 ! 222 ! 220 223 ! 3667 20190110 14:26:24Z schwenkel 221 224 ! Modified check for rrtmg input files 222 ! 225 ! 223 226 ! 3655 20190107 16:51:22Z knoop 224 227 ! nopointer option removed 225 ! 228 ! 226 229 ! 1496 20141202 17:25:50Z maronga 227 230 ! Initial revision 228 ! 231 ! 229 232 ! 230 233 ! Description: … … 247 250 !! 248 251 MODULE radiation_model_mod 249 252 250 253 USE arrays_3d, & 251 254 ONLY: dzw, hyp, nc, pt, p, q, ql, u, v, w, zu, zw, exner, d_exner … … 276 279 277 280 USE grid_variables, & 278 ONLY: ddx, ddy, dx, dy 281 ONLY: ddx, ddy, dx, dy 279 282 280 283 USE indices, & … … 366 369 ! Predefined Land surface classes (albedo_type) after Briegleb (1992) 367 370 CHARACTER(37), DIMENSION(0:33), PARAMETER :: albedo_type_name = (/ & 368 'user defined ', & ! 0 371 'user defined ', & ! 0 369 372 'ocean ', & ! 1 370 373 'mixed farming, tall grassland ', & ! 2 371 'tall/medium grassland ', & ! 3 374 'tall/medium grassland ', & ! 3 372 375 'evergreen shrubland ', & ! 4 373 376 'short grassland/meadow/shrubland ', & ! 5 … … 377 380 'tropical evergreen broadleaved forest', & ! 9 378 381 'medium/tall grassland/woodland ', & ! 10 379 'desert, sandy ', & ! 11 380 'desert, rocky ', & ! 12 382 'desert, sandy ', & ! 11 383 'desert, rocky ', & ! 12 381 384 'tundra ', & ! 13 382 'land ice ', & ! 14 383 'sea ice ', & ! 15 385 'land ice ', & ! 14 386 'sea ice ', & ! 15 384 387 'snow ', & ! 16 385 388 'bare soil ', & ! 17 … … 415 418 416 419 INTEGER(iwp) :: albedo_type = 9999999_iwp, & !< Albedo surface type 417 dots_rad = 0_iwp !< starting index for timeseries output 420 dots_rad = 0_iwp !< starting index for timeseries output 418 421 419 422 LOGICAL :: unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed … … 467 470 REAL(wp), PARAMETER :: emissivity_atm_clsky = 0.8_wp !< emissivity of the clearsky atmosphere 468 471 ! 469 ! Land surface albedos for solar zenith angle of 60degree after Briegleb (1992) 472 ! Land surface albedos for solar zenith angle of 60degree after Briegleb (1992) 470 473 ! (broadband, longwave, shortwave ): bb, lw, sw, 471 REAL(wp), DIMENSION(0:2,1:33), PARAMETER :: albedo_pars = RESHAPE( (/& 474 REAL(wp), DIMENSION(0:2,1:33), PARAMETER :: albedo_pars = RESHAPE( (/& 472 475 0.06_wp, 0.06_wp, 0.06_wp, & ! 1 473 476 0.19_wp, 0.28_wp, 0.09_wp, & ! 2 … … 541 544 542 545 ! 543 ! The following variables should be only changed with care, as this will 546 ! The following variables should be only changed with care, as this will 544 547 ! require further setting of some variables, which is currently not 545 548 ! implemented (aerosols, ice phase). … … 570 573 rrtm_cliqwp, & !< incloud liquid water path (g/m2) 571 574 rrtm_co2vmr, & !< CO2 volume mixing ratio (g/mol) 572 rrtm_emis, & !< surface emissivity (01) 575 rrtm_emis, & !< surface emissivity (01) 573 576 rrtm_h2ovmr, & !< H2O volume mixing ratio 574 577 rrtm_n2ovmr, & !< N2O volume mixing ratio … … 582 585 rrtm_tlev, & !< actual temperature (K, zwgrid) 583 586 rrtm_lwdflx, & !< RRTM output of incoming longwave radiation flux (W/m2) 584 rrtm_lwdflxc, & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2) 587 rrtm_lwdflxc, & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2) 585 588 rrtm_lwuflx, & !< RRTM output of outgoing longwave radiation flux (W/m2) 586 589 rrtm_lwuflxc, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2) … … 590 593 rrtm_lwhrc, & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d) 591 594 rrtm_swdflx, & !< RRTM output of incoming shortwave radiation flux (W/m2) 592 rrtm_swdflxc, & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2) 595 rrtm_swdflxc, & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2) 593 596 rrtm_swuflx, & !< RRTM output of outgoing shortwave radiation flux (W/m2) 594 597 rrtm_swuflxc, & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2) 595 598 rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d) 596 rrtm_swhrc, & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 597 rrtm_dirdflux, & !< RRTM output of incoming direct shortwave (W/m2) 598 rrtm_difdflux !< RRTM output of incoming diffuse shortwave (W/m2) 599 rrtm_swhrc, & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 600 rrtm_dirdflux, & !< RRTM output of incoming direct shortwave (W/m2) 601 rrtm_difdflux !< RRTM output of incoming diffuse shortwave (W/m2) 599 602 600 603 REAL(wp), DIMENSION(1) :: rrtm_aldif, & !< surface albedo for longwave diffuse radiation … … 734 737 ! type for calculation of svf 735 738 TYPE t_svf 736 INTEGER(iwp) :: isurflt !< 737 INTEGER(iwp) :: isurfs !< 738 REAL(wp) :: rsvf !< 739 REAL(wp) :: rtransp !< 739 INTEGER(iwp) :: isurflt !< 740 INTEGER(iwp) :: isurfs !< 741 REAL(wp) :: rsvf !< 742 REAL(wp) :: rtransp !< 740 743 END TYPE 741 744 742 745 ! type for calculation of csf 743 746 TYPE t_csf 744 INTEGER(iwp) :: ip !< 745 INTEGER(iwp) :: itx !< 746 INTEGER(iwp) :: ity !< 747 INTEGER(iwp) :: itz !< 747 INTEGER(iwp) :: ip !< 748 INTEGER(iwp) :: itx !< 749 INTEGER(iwp) :: ity !< 750 INTEGER(iwp) :: itz !< 748 751 INTEGER(iwp) :: isurfs !< Idx of source face / 1 for sky 749 752 REAL(wp) :: rcvf !< Canopy view factor for faces / … … 788 791 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif !< array of diffuse sw radiation from sky and model boundary falling to local surface 789 792 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif !< array of diffuse lw radiation from sky and model boundary falling to local surface 790 793 791 794 !< Outward radiation is only valid for nonvirtual surfaces 792 795 REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsl !< array of reflected sw radiation for local surface in ith reflection … … 835 838 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: boxes !< coordinates of gridboxes being crossed by ray 836 839 REAL(wp), DIMENSION(:), ALLOCATABLE :: crlens !< array of crossing lengths of ray for particular grid boxes 837 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: lad_ip !< array of numbers of process where lad is stored 840 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: lad_ip !< array of numbers of process where lad is stored 838 841 #if defined( __parallel ) 839 842 INTEGER(kind=MPI_ADDRESS_KIND), & … … 876 879 REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf !< emissivity of the wall surface 877 880 ! 878 ! External radiation. Depending on the given level of detail either a 1D or 881 ! External radiation. Depending on the given level of detail either a 1D or 879 882 ! a 3D array will be allocated. 880 TYPE( real_1d_3d ) :: rad_lw_in_f !< external incoming longwave radiation, from observation or model 881 TYPE( real_1d_3d ) :: rad_sw_in_f !< external incoming shortwave radiation, from observation or model 883 TYPE( real_1d_3d ) :: rad_lw_in_f !< external incoming longwave radiation, from observation or model 884 TYPE( real_1d_3d ) :: rad_sw_in_f !< external incoming shortwave radiation, from observation or model 882 885 TYPE( real_1d_3d ) :: rad_sw_in_dif_f !< external incoming shortwave radiation, diffuse part, from observation or model 883 TYPE( real_1d_3d ) :: time_rad_f !< time dimension for external radiation, from observation or model 886 TYPE( real_1d_3d ) :: time_rad_f !< time dimension for external radiation, from observation or model 884 887 885 888 INTERFACE radiation_check_data_output … … 894 897 MODULE PROCEDURE radiation_check_data_output_pr 895 898 END INTERFACE radiation_check_data_output_pr 896 899 897 900 INTERFACE radiation_check_parameters 898 901 MODULE PROCEDURE radiation_check_parameters 899 902 END INTERFACE radiation_check_parameters 900 903 901 904 INTERFACE radiation_clearsky 902 905 MODULE PROCEDURE radiation_clearsky 903 906 END INTERFACE radiation_clearsky 904 907 905 908 INTERFACE radiation_constant 906 909 MODULE PROCEDURE radiation_constant 907 910 END INTERFACE radiation_constant 908 911 909 912 INTERFACE radiation_control 910 913 MODULE PROCEDURE radiation_control … … 933 936 INTERFACE radiation_header 934 937 MODULE PROCEDURE radiation_header 935 END INTERFACE radiation_header 936 938 END INTERFACE radiation_header 939 937 940 INTERFACE radiation_init 938 941 MODULE PROCEDURE radiation_init … … 942 945 MODULE PROCEDURE radiation_parin 943 946 END INTERFACE radiation_parin 944 947 945 948 INTERFACE radiation_rrtmg 946 949 MODULE PROCEDURE radiation_rrtmg … … 969 972 MODULE PROCEDURE radiation_interaction_init 970 973 END INTERFACE radiation_interaction_init 971 974 972 975 INTERFACE radiation_presimulate_solar_pos 973 976 MODULE PROCEDURE radiation_presimulate_solar_pos … … 1005 1008 radiation_read_svf, radiation_presimulate_solar_pos 1006 1009 1007 1010 1008 1011 ! 1009 1012 ! Public variables and constants / NEEDS SORTING … … 1042 1045 !! 1043 1046 SUBROUTINE radiation_control 1044 1045 1047 1048 1046 1049 IMPLICIT NONE 1047 1050 … … 1054 1057 CASE ( 'constant' ) 1055 1058 CALL radiation_constant 1056 1057 CASE ( 'clearsky' ) 1059 1060 CASE ( 'clearsky' ) 1058 1061 CALL radiation_clearsky 1059 1062 1060 1063 CASE ( 'rrtmg' ) 1061 1064 CALL radiation_rrtmg 1062 1065 1063 1066 CASE ( 'external' ) 1064 1067 ! … … 1084 1087 !! 1085 1088 SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k ) 1086 1087 1089 1090 1088 1091 USE control_parameters, & 1089 1092 ONLY: data_output, message_string … … 1262 1265 !  1263 1266 !> Check data output of profiles for radiation model 1264 !! 1267 !! 1265 1268 SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, & 1266 1269 dopr_unit ) 1267 1270 1268 1271 USE arrays_3d, & 1269 1272 ONLY: zu … … 1279 1282 1280 1283 IMPLICIT NONE 1281 1282 CHARACTER (LEN=*) :: unit !< 1283 CHARACTER (LEN=*) :: variable !< 1284 1285 CHARACTER (LEN=*) :: unit !< 1286 CHARACTER (LEN=*) :: variable !< 1284 1287 CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit 1285 1286 INTEGER(iwp) :: var_count !< 1288 1289 INTEGER(iwp) :: var_count !< 1287 1290 1288 1291 SELECT CASE ( TRIM( variable ) ) 1289 1292 1290 1293 CASE ( 'rad_net' ) 1291 1294 IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )& … … 1315 1318 dopr_unit = 'W/m2' 1316 1319 hom(:,2,100,:) = SPREAD( zw, 2, statistic_regions+1 ) 1317 unit = dopr_unit 1320 unit = dopr_unit 1318 1321 ENDIF 1319 1322 … … 1330 1333 dopr_unit = 'W/m2' 1331 1334 hom(:,2,101,:) = SPREAD( zw, 2, statistic_regions+1 ) 1332 unit = dopr_unit 1335 unit = dopr_unit 1333 1336 ENDIF 1334 1337 … … 1431 1434 1432 1435 END SUBROUTINE radiation_check_data_output_pr 1433 1434 1436 1437 1435 1438 !! 1436 1439 ! Description: … … 1444 1447 1445 1448 USE netcdf_data_input_mod, & 1446 ONLY: input_pids_static 1447 1449 ONLY: input_pids_static 1450 1448 1451 IMPLICIT NONE 1449 1450 ! 1451 ! In case no urbansurface or landsurface model is applied, usage of 1452 ! a radiation model make no sense. 1452 1453 ! 1454 ! In case no urbansurface or landsurface model is applied, usage of 1455 ! a radiation model make no sense. 1453 1456 IF ( .NOT. land_surface .AND. .NOT. urban_surface ) THEN 1454 1457 message_string = 'Usage of radiation module is only allowed if ' // & … … 1466 1469 ELSEIF ( radiation_scheme == 'rrtmg' ) THEN 1467 1470 #if ! defined ( __rrtmg ) 1468 message_string = 'radiation_scheme = "rrtmg" requires ' // & 1471 message_string = 'radiation_scheme = "rrtmg" requires ' // & 1469 1472 'compilation of PALM with preprocessor ' // & 1470 1473 'directive D__rrtmg' … … 1472 1475 #endif 1473 1476 #if defined ( __rrtmg ) && ! defined( __netcdf ) 1474 message_string = 'radiation_scheme = "rrtmg" requires ' // & 1477 message_string = 'radiation_scheme = "rrtmg" requires ' // & 1475 1478 'the use of NetCDF (preprocessor directive ' // & 1476 1479 'D__netcdf' … … 1480 1483 ENDIF 1481 1484 ! 1482 ! Checks performed only if data is given via namelist only. 1485 ! Checks performed only if data is given via namelist only. 1483 1486 IF ( .NOT. input_pids_static ) THEN 1484 1487 IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND. & 1485 1488 radiation_scheme == 'clearsky') THEN 1486 message_string = 'radiation_scheme = "clearsky" in combination'//& 1489 message_string = 'radiation_scheme = "clearsky" in combination'//& 1487 1490 'with albedo_type = 0 requires setting of'// & 1488 1491 'albedo /= 9999999.9' … … 1492 1495 IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND. & 1493 1496 ( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp& 1494 .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 1497 .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 1495 1498 ) ) THEN 1496 message_string = 'radiation_scheme = "rrtmg" in combination' // & 1499 message_string = 'radiation_scheme = "rrtmg" in combination' // & 1497 1500 'with albedo_type = 0 requires setting of ' // & 1498 1501 'albedo_lw_dif /= 9999999.9' // & … … 1506 1509 ! Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented 1507 1510 ! Serial mode does not allow mpi_rma 1508 #if defined( __parallel ) 1511 #if defined( __parallel ) 1509 1512 IF ( rad_angular_discretization .AND. .NOT. raytrace_mpi_rma ) THEN 1510 1513 message_string = 'rad_angular_discretization can only be used ' // & … … 1522 1525 IF ( cloud_droplets .AND. radiation_scheme == 'rrtmg' .AND. & 1523 1526 average_radiation ) THEN 1524 message_string = 'average_radiation = .T. with radiation_scheme'// & 1527 message_string = 'average_radiation = .T. with radiation_scheme'// & 1525 1528 '= "rrtmg" in combination cloud_droplets = .T.'// & 1526 1529 'is not implementd' … … 1539 1542 ! Check for dt_radiation 1540 1543 IF ( dt_radiation <= 0.0 ) THEN 1541 message_string = 'dt_radiation must be > 0.0' 1542 CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 ) 1544 message_string = 'dt_radiation must be > 0.0' 1545 CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 ) 1543 1546 ENDIF 1544 1547 ! … … 1549 1552 'model.&Using rotation_angle /= 0.0 is not allowed in combination ' // & 1550 1553 'with the radiation model at the moment!' 1551 CALL message( 'check_parameters', 'PA0675', 1, 2, 0, 6, 0 ) 1554 CALL message( 'check_parameters', 'PA0675', 1, 2, 0, 6, 0 ) 1552 1555 ENDIF 1553 1554 END SUBROUTINE radiation_check_parameters 1555 1556 1556 1557 END SUBROUTINE radiation_check_parameters 1558 1559 1557 1560 !! 1558 1561 ! Description: … … 1561 1564 !! 1562 1565 SUBROUTINE radiation_init 1563 1566 1564 1567 IMPLICIT NONE 1565 1568 1566 INTEGER(iwp) :: i !< running index xdirection 1569 INTEGER(iwp) :: i !< running index xdirection 1567 1570 INTEGER(iwp) :: is !< running index for input surface elements 1568 INTEGER(iwp) :: ioff !< offset in x between surface element reference grid point in atmosphere and actual surface 1569 INTEGER(iwp) :: j !< running index ydirection 1570 INTEGER(iwp) :: joff !< offset in y between surface element reference grid point in atmosphere and actual surface 1571 INTEGER(iwp) :: ioff !< offset in x between surface element reference grid point in atmosphere and actual surface 1572 INTEGER(iwp) :: j !< running index ydirection 1573 INTEGER(iwp) :: joff !< offset in y between surface element reference grid point in atmosphere and actual surface 1571 1574 INTEGER(iwp) :: k !< running index zdirection 1572 INTEGER(iwp) :: l !< running index for orientation of vertical surfaces 1575 INTEGER(iwp) :: l !< running index for orientation of vertical surfaces 1573 1576 INTEGER(iwp) :: m !< running index for surface elements 1574 1577 INTEGER(iwp) :: ntime = 0 !< number of available external radiation timesteps … … 1584 1587 ! or if biometeorology output is required for flat surfaces. 1585 1588 ! The namelist parameter radiation_interactions_on can override this behavior. 1586 ! (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in 1589 ! (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in 1587 1590 ! init_surface_arrays.) 1588 1591 IF ( radiation_interactions_on ) THEN … … 1591 1594 average_radiation = .TRUE. 1592 1595 ELSE 1593 radiation_interactions_on = .FALSE. !< reset namelist parameter: no interactions 1596 radiation_interactions_on = .FALSE. !< reset namelist parameter: no interactions 1594 1597 !< calculations necessary in case of flat surface 1595 1598 ENDIF 1596 1599 ELSEIF ( vertical_surfaces_exist .OR. plant_canopy .OR. biometeorology ) THEN 1597 1600 message_string = 'radiation_interactions_on is set to .FALSE. although ' // & 1598 'vertical surfaces and/or trees or biometeorology exist ' // & 1601 'vertical surfaces and/or trees or biometeorology exist ' // & 1599 1602 'is ON. The model will run without RTM (no shadows, no ' // & 1600 1603 'radiation reflections)' … … 1607 1610 1608 1611 ! 1609 ! If required, initialize radiation interactions between surfaces 1612 ! If required, initialize radiation interactions between surfaces 1610 1613 ! via skyview factors. This must be done before radiation is initialized. 1611 1614 IF ( radiation_interactions ) CALL radiation_interaction_init … … 1615 1618 surf_lsm_h%ns > 0 ) THEN 1616 1619 ALLOCATE( surf_lsm_h%rad_net(1:surf_lsm_h%ns) ) 1617 surf_lsm_h%rad_net = 0.0_wp 1620 surf_lsm_h%rad_net = 0.0_wp 1618 1621 ENDIF 1619 1622 IF ( .NOT. ALLOCATED ( surf_usm_h%rad_net ) .AND. & 1620 1623 surf_usm_h%ns > 0 ) THEN 1621 1624 ALLOCATE( surf_usm_h%rad_net(1:surf_usm_h%ns) ) 1622 surf_usm_h%rad_net = 0.0_wp 1625 surf_usm_h%rad_net = 0.0_wp 1623 1626 ENDIF 1624 1627 DO l = 0, 3 … … 1626 1629 surf_lsm_v(l)%ns > 0 ) THEN 1627 1630 ALLOCATE( surf_lsm_v(l)%rad_net(1:surf_lsm_v(l)%ns) ) 1628 surf_lsm_v(l)%rad_net = 0.0_wp 1631 surf_lsm_v(l)%rad_net = 0.0_wp 1629 1632 ENDIF 1630 1633 IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_net ) .AND. & 1631 1634 surf_usm_v(l)%ns > 0 ) THEN 1632 1635 ALLOCATE( surf_usm_v(l)%rad_net(1:surf_usm_v(l)%ns) ) 1633 surf_usm_v(l)%rad_net = 0.0_wp 1636 surf_usm_v(l)%rad_net = 0.0_wp 1634 1637 ENDIF 1635 1638 ENDDO … … 1641 1644 surf_lsm_h%ns > 0 ) THEN 1642 1645 ALLOCATE( surf_lsm_h%rad_lw_out_change_0(1:surf_lsm_h%ns) ) 1643 surf_lsm_h%rad_lw_out_change_0 = 0.0_wp 1646 surf_lsm_h%rad_lw_out_change_0 = 0.0_wp 1644 1647 ENDIF 1645 1648 IF ( .NOT. ALLOCATED ( surf_usm_h%rad_lw_out_change_0 ) .AND. & 1646 1649 surf_usm_h%ns > 0 ) THEN 1647 1650 ALLOCATE( surf_usm_h%rad_lw_out_change_0(1:surf_usm_h%ns) ) 1648 surf_usm_h%rad_lw_out_change_0 = 0.0_wp 1651 surf_usm_h%rad_lw_out_change_0 = 0.0_wp 1649 1652 ENDIF 1650 1653 DO l = 0, 3 … … 1652 1655 surf_lsm_v(l)%ns > 0 ) THEN 1653 1656 ALLOCATE( surf_lsm_v(l)%rad_lw_out_change_0(1:surf_lsm_v(l)%ns) ) 1654 surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp 1657 surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp 1655 1658 ENDIF 1656 1659 IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_lw_out_change_0 ) .AND. & 1657 1660 surf_usm_v(l)%ns > 0 ) THEN 1658 1661 ALLOCATE( surf_usm_v(l)%rad_lw_out_change_0(1:surf_usm_v(l)%ns) ) 1659 surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp 1662 surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp 1660 1663 ENDIF 1661 1664 ENDDO … … 1676 1679 ALLOCATE( surf_lsm_h%rad_lw_ref(1:surf_lsm_h%ns) ) 1677 1680 ALLOCATE( surf_lsm_h%rad_lw_res(1:surf_lsm_h%ns) ) 1678 surf_lsm_h%rad_sw_in = 0.0_wp 1679 surf_lsm_h%rad_sw_out = 0.0_wp 1680 surf_lsm_h%rad_sw_dir = 0.0_wp 1681 surf_lsm_h%rad_sw_dif = 0.0_wp 1682 surf_lsm_h%rad_sw_ref = 0.0_wp 1683 surf_lsm_h%rad_sw_res = 0.0_wp 1684 surf_lsm_h%rad_lw_in = 0.0_wp 1685 surf_lsm_h%rad_lw_out = 0.0_wp 1686 surf_lsm_h%rad_lw_dif = 0.0_wp 1687 surf_lsm_h%rad_lw_ref = 0.0_wp 1688 surf_lsm_h%rad_lw_res = 0.0_wp 1681 surf_lsm_h%rad_sw_in = 0.0_wp 1682 surf_lsm_h%rad_sw_out = 0.0_wp 1683 surf_lsm_h%rad_sw_dir = 0.0_wp 1684 surf_lsm_h%rad_sw_dif = 0.0_wp 1685 surf_lsm_h%rad_sw_ref = 0.0_wp 1686 surf_lsm_h%rad_sw_res = 0.0_wp 1687 surf_lsm_h%rad_lw_in = 0.0_wp 1688 surf_lsm_h%rad_lw_out = 0.0_wp 1689 surf_lsm_h%rad_lw_dif = 0.0_wp 1690 surf_lsm_h%rad_lw_ref = 0.0_wp 1691 surf_lsm_h%rad_lw_res = 0.0_wp 1689 1692 ENDIF 1690 1693 IF ( .NOT. ALLOCATED ( surf_usm_h%rad_sw_in ) .AND. & … … 1701 1704 ALLOCATE( surf_usm_h%rad_lw_ref(1:surf_usm_h%ns) ) 1702 1705 ALLOCATE( surf_usm_h%rad_lw_res(1:surf_usm_h%ns) ) 1703 surf_usm_h%rad_sw_in = 0.0_wp 1704 surf_usm_h%rad_sw_out = 0.0_wp 1705 surf_usm_h%rad_sw_dir = 0.0_wp 1706 surf_usm_h%rad_sw_dif = 0.0_wp 1707 surf_usm_h%rad_sw_ref = 0.0_wp 1708 surf_usm_h%rad_sw_res = 0.0_wp 1709 surf_usm_h%rad_lw_in = 0.0_wp 1710 surf_usm_h%rad_lw_out = 0.0_wp 1711 surf_usm_h%rad_lw_dif = 0.0_wp 1712 surf_usm_h%rad_lw_ref = 0.0_wp 1713 surf_usm_h%rad_lw_res = 0.0_wp 1706 surf_usm_h%rad_sw_in = 0.0_wp 1707 surf_usm_h%rad_sw_out = 0.0_wp 1708 surf_usm_h%rad_sw_dir = 0.0_wp 1709 surf_usm_h%rad_sw_dif = 0.0_wp 1710 surf_usm_h%rad_sw_ref = 0.0_wp 1711 surf_usm_h%rad_sw_res = 0.0_wp 1712 surf_usm_h%rad_lw_in = 0.0_wp 1713 surf_usm_h%rad_lw_out = 0.0_wp 1714 surf_usm_h%rad_lw_dif = 0.0_wp 1715 surf_usm_h%rad_lw_ref = 0.0_wp 1716 surf_usm_h%rad_lw_res = 0.0_wp 1714 1717 ENDIF 1715 1718 DO l = 0, 3 … … 1729 1732 ALLOCATE( surf_lsm_v(l)%rad_lw_res(1:surf_lsm_v(l)%ns) ) 1730 1733 1731 surf_lsm_v(l)%rad_sw_in = 0.0_wp 1734 surf_lsm_v(l)%rad_sw_in = 0.0_wp 1732 1735 surf_lsm_v(l)%rad_sw_out = 0.0_wp 1733 1736 surf_lsm_v(l)%rad_sw_dir = 0.0_wp … … 1736 1739 surf_lsm_v(l)%rad_sw_res = 0.0_wp 1737 1740 1738 surf_lsm_v(l)%rad_lw_in = 0.0_wp 1739 surf_lsm_v(l)%rad_lw_out = 0.0_wp 1740 surf_lsm_v(l)%rad_lw_dif = 0.0_wp 1741 surf_lsm_v(l)%rad_lw_ref = 0.0_wp 1742 surf_lsm_v(l)%rad_lw_res = 0.0_wp 1741 surf_lsm_v(l)%rad_lw_in = 0.0_wp 1742 surf_lsm_v(l)%rad_lw_out = 0.0_wp 1743 surf_lsm_v(l)%rad_lw_dif = 0.0_wp 1744 surf_lsm_v(l)%rad_lw_ref = 0.0_wp 1745 surf_lsm_v(l)%rad_lw_res = 0.0_wp 1743 1746 ENDIF 1744 1747 IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_sw_in ) .AND. & … … 1755 1758 ALLOCATE( surf_usm_v(l)%rad_lw_ref(1:surf_usm_v(l)%ns) ) 1756 1759 ALLOCATE( surf_usm_v(l)%rad_lw_res(1:surf_usm_v(l)%ns) ) 1757 surf_usm_v(l)%rad_sw_in = 0.0_wp 1760 surf_usm_v(l)%rad_sw_in = 0.0_wp 1758 1761 surf_usm_v(l)%rad_sw_out = 0.0_wp 1759 1762 surf_usm_v(l)%rad_sw_dir = 0.0_wp … … 1761 1764 surf_usm_v(l)%rad_sw_ref = 0.0_wp 1762 1765 surf_usm_v(l)%rad_sw_res = 0.0_wp 1763 surf_usm_v(l)%rad_lw_in = 0.0_wp 1764 surf_usm_v(l)%rad_lw_out = 0.0_wp 1765 surf_usm_v(l)%rad_lw_dif = 0.0_wp 1766 surf_usm_v(l)%rad_lw_ref = 0.0_wp 1767 surf_usm_v(l)%rad_lw_res = 0.0_wp 1766 surf_usm_v(l)%rad_lw_in = 0.0_wp 1767 surf_usm_v(l)%rad_lw_out = 0.0_wp 1768 surf_usm_v(l)%rad_lw_dif = 0.0_wp 1769 surf_usm_v(l)%rad_lw_ref = 0.0_wp 1770 surf_usm_v(l)%rad_lw_res = 0.0_wp 1768 1771 ENDIF 1769 1772 ENDDO … … 1855 1858 ENDDO 1856 1859 ! 1857 ! Level 2 initialization of broadband albedo via given albedo_type. 1858 ! Only if albedo_type is nonzero. In case of urban surface and 1860 ! Level 2 initialization of broadband albedo via given albedo_type. 1861 ! Only if albedo_type is nonzero. In case of urban surface and 1859 1862 ! input data is read from ASCII file, albedo_type will be zero, so that 1860 ! albedo won't be overwritten. 1863 ! albedo won't be overwritten. 1861 1864 DO m = 1, surf_lsm_h%ns 1862 1865 IF ( surf_lsm_h%albedo_type(m,ind_veg_wall) /= 0 ) & … … 1931 1934 surf_usm_h%albedo(m,ind_wat_win) = albedo_pars_f%pars_xy(0,j,i) 1932 1935 ENDIF 1933 ENDDO 1934 ! 1935 ! Vertical surfaces 1936 ENDDO 1937 ! 1938 ! Vertical surfaces 1936 1939 DO l = 0, 3 1937 1940 … … 1961 1964 ENDDO 1962 1965 1963 ENDIF 1966 ENDIF 1964 1967 ! 1965 1968 ! Read explicit albedo values from building surface pars. If present, … … 2051 2054 ! Allocate albedos for short/longwave radiation, horizontal surfaces 2052 2055 ! for wall/green/window (USM) or vegetation/pavement/water surfaces 2053 ! (LSM). 2056 ! (LSM). 2054 2057 ALLOCATE ( surf_lsm_h%aldif(1:surf_lsm_h%ns,0:2) ) 2055 2058 ALLOCATE ( surf_lsm_h%aldir(1:surf_lsm_h%ns,0:2) ) … … 2071 2074 2072 2075 ! 2073 ! Allocate broadband albedo (temporary for the current radiation 2076 ! Allocate broadband albedo (temporary for the current radiation 2074 2077 ! implementations) 2075 2078 IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) ) & … … 2079 2082 2080 2083 ! 2081 ! Allocate albedos for short/longwave radiation, vertical surfaces 2084 ! Allocate albedos for short/longwave radiation, vertical surfaces 2082 2085 DO l = 0, 3 2083 2086 … … 2111 2114 ENDDO 2112 2115 ! 2113 ! Level 1 initialization of spectral albedos via namelist 2114 ! paramters. Please note, this case all surface tiles are initialized 2115 ! the same. 2116 ! Level 1 initialization of spectral albedos via namelist 2117 ! paramters. Please note, this case all surface tiles are initialized 2118 ! the same. 2116 2119 IF ( surf_lsm_h%ns > 0 ) THEN 2117 2120 surf_lsm_h%aldif = albedo_lw_dif … … 2162 2165 2163 2166 ! 2164 ! Level 2 initialization of spectral albedos via albedo_type. 2165 ! Please note, for natural and urbantype surfaces, a tile approach 2166 ! is applied so that the resulting albedo is calculated via the weighted 2167 ! average of respective surface fractions. 2167 ! Level 2 initialization of spectral albedos via albedo_type. 2168 ! Please note, for natural and urbantype surfaces, a tile approach 2169 ! is applied so that the resulting albedo is calculated via the weighted 2170 ! average of respective surface fractions. 2168 2171 DO m = 1, surf_lsm_h%ns 2169 2172 ! … … 2186 2189 ENDDO 2187 2190 ! 2188 ! For urban surface only if albedo has not been already initialized 2191 ! For urban surface only if albedo has not been already initialized 2189 2192 ! in the urbansurface model via the ASCII file. 2190 2193 IF ( .NOT. surf_usm_h%albedo_from_ascii ) THEN … … 2231 2234 ENDDO 2232 2235 ! 2233 ! For urban surface only if albedo has not been already initialized 2236 ! For urban surface only if albedo has not been already initialized 2234 2237 ! in the urbansurface model via the ASCII file. 2235 2238 IF ( .NOT. surf_usm_v(l)%albedo_from_ascii ) THEN … … 2269 2272 IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill ) & 2270 2273 surf_lsm_h%albedo(m,ind_type) = & 2271 albedo_pars_f%pars_xy(0,j,i) 2274 albedo_pars_f%pars_xy(0,j,i) 2272 2275 IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill ) & 2273 2276 surf_lsm_h%aldir(m,ind_type) = & 2274 albedo_pars_f%pars_xy(1,j,i) 2277 albedo_pars_f%pars_xy(1,j,i) 2275 2278 IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill ) & 2276 2279 surf_lsm_h%aldif(m,ind_type) = & 2277 albedo_pars_f%pars_xy(1,j,i) 2280 albedo_pars_f%pars_xy(1,j,i) 2278 2281 IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill ) & 2279 2282 surf_lsm_h%asdir(m,ind_type) = & 2280 albedo_pars_f%pars_xy(2,j,i) 2283 albedo_pars_f%pars_xy(2,j,i) 2281 2284 IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill ) & 2282 2285 surf_lsm_h%asdif(m,ind_type) = & … … 2285 2288 ENDDO 2286 2289 ! 2287 ! For urban surface only if albedo has not been already initialized 2290 ! For urban surface only if albedo has not been already initialized 2288 2291 ! in the urbansurface model via the ASCII file. 2289 2292 IF ( .NOT. surf_usm_h%albedo_from_ascii ) THEN … … 2378 2381 ENDDO 2379 2382 ! 2380 ! For urban surface only if albedo has not been already initialized 2383 ! For urban surface only if albedo has not been already initialized 2381 2384 ! in the urbansurface model via the ASCII file. 2382 2385 IF ( .NOT. surf_usm_v(l)%albedo_from_ascii ) THEN … … 2411 2414 albedo_pars_f%pars_xy(2,j+joff,i+ioff) 2412 2415 ENDIF 2413 ! 2416 ! 2414 2417 ! Spectral albedos especially for building green surfaces 2415 2418 IF ( albedo_pars_f%pars_xy(3,j+joff,i+ioff) /= & … … 2427 2430 albedo_pars_f%pars_xy(4,j+joff,i+ioff) 2428 2431 ENDIF 2429 ! 2432 ! 2430 2433 ! Spectral albedos especially for building window surfaces 2431 2434 IF ( albedo_pars_f%pars_xy(5,j+joff,i+ioff) /= & … … 2640 2643 2641 2644 ! 2642 ! Calculate initial values of current (cosine of) the zenith angle and 2645 ! Calculate initial values of current (cosine of) the zenith angle and 2643 2646 ! whether the sun is up 2644 2647 CALL get_date_time( time_since_reference_point, & … … 2785 2788 2786 2789 ! 2787 ! Allocate 1element array for surface temperature 2790 ! Allocate 1element array for surface temperature 2788 2791 ! (RRTMG anticipates an array as passed argument). 2789 2792 ALLOCATE ( rrtm_tsfc(1) ) 2790 2793 ! 2791 ! Allocate surface emissivity. 2794 ! Allocate surface emissivity. 2792 2795 ! Values will be given directly before calling rrtm_lw. 2793 2796 ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) ) … … 2801 2804 '&Please provide <jobname>_lsw file in the INPUT directory.' 2802 2805 CALL message( 'radiation_init', 'PA0583', 1, 2, 0, 6, 0 ) 2803 ENDIF 2806 ENDIF 2804 2807 INQUIRE( FILE='rrtmg_sw.nc', EXIST=sw_exists ) 2805 2808 IF ( .NOT. sw_exists ) THEN … … 2808 2811 '&Please provide <jobname>_rsw file in the INPUT directory.' 2809 2812 CALL message( 'radiation_init', 'PA0584', 1, 2, 0, 6, 0 ) 2810 ENDIF 2811 2813 ENDIF 2814 2812 2815 IF ( lw_radiation ) CALL rrtmg_lw_ini ( c_p ) 2813 2816 IF ( sw_radiation ) CALL rrtmg_sw_ini ( c_p ) 2814 2817 2815 2818 ! 2816 2819 ! Set input files for RRTMG 2817 INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists) 2820 INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists) 2818 2821 IF ( .NOT. snd_exists ) THEN 2819 2822 rrtm_input_file = "rrtmg_lw.nc" … … 2823 2826 ! Read vertical layers for RRTMG from sounding data 2824 2827 ! The routine provides nzt_rad, hyp_snd(1:nzt_rad), 2825 ! t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1), 2828 ! t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1), 2826 2829 ! rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1) 2827 2830 CALL read_sounding_data … … 2834 2837 ENDIF 2835 2838 ! 2836 ! Initializaion actions exclusively required for external 2839 ! Initializaion actions exclusively required for external 2837 2840 ! radiation forcing 2838 2841 IF ( radiation_scheme == 'external' ) THEN 2839 2842 ! 2840 ! Open the radiation input file. Note, for child domain, a dynamic 2841 ! input file is often not provided. In order to do not need to 2843 ! Open the radiation input file. Note, for child domain, a dynamic 2844 ! input file is often not provided. In order to do not need to 2842 2845 ! duplicate the dynamic input file just for the radiation input, take 2843 ! it from the dynamic file for the parent if not available for the 2844 ! child domain(s). In this case this is possible because radiation 2846 ! it from the dynamic file for the parent if not available for the 2847 ! child domain(s). In this case this is possible because radiation 2845 2848 ! input should be the same for each model. 2846 2849 INQUIRE( FILE = TRIM( input_file_dynamic ), & 2847 2850 EXIST = radiation_input_root_domain ) 2848 2851 2849 2852 IF ( .NOT. input_pids_dynamic .AND. & 2850 2853 .NOT. radiation_input_root_domain ) THEN … … 2859 2862 ! 2860 2863 ! Open dynamic input file for child domain if available, else, open 2861 ! dynamic input file for the root domain. 2864 ! dynamic input file for the root domain. 2862 2865 IF ( input_pids_dynamic ) THEN 2863 2866 CALL open_read_file( TRIM( input_file_dynamic ) // & … … 2868 2871 pids_id ) 2869 2872 ENDIF 2870 2873 2871 2874 CALL inquire_num_variables( pids_id, num_var_pids ) 2872 ! 2875 ! 2873 2876 ! Allocate memory to store variable names and read them 2874 2877 ALLOCATE( vars_pids(1:num_var_pids) ) 2875 2878 CALL inquire_variable_names( pids_id, vars_pids ) 2876 ! 2879 ! 2877 2880 ! Input time dimension. 2878 2881 IF ( check_existence( vars_pids, 'time_rad' ) ) THEN 2879 2882 CALL get_dimension_length( pids_id, ntime, 'time_rad' ) 2880 2883 2881 2884 ALLOCATE( time_rad_f%var1d(0:ntime1) ) 2882 ! 2883 ! Read variable 2885 ! 2886 ! Read variable 2884 2887 CALL get_variable( pids_id, 'time_rad', time_rad_f%var1d ) 2885 2888 2886 2889 time_rad_f%from_file = .TRUE. 2887 ENDIF 2888 ! 2890 ENDIF 2891 ! 2889 2892 ! Input shortwave downwelling. 2890 2893 IF ( check_existence( vars_pids, 'rad_sw_in' ) ) THEN 2891 ! 2894 ! 2892 2895 ! Get _FillValue attribute 2893 2896 CALL get_attribute( pids_id, char_fill, rad_sw_in_f%fill, & 2894 2897 .FALSE., 'rad_sw_in' ) 2895 ! 2898 ! 2896 2899 ! Get levelofdetail 2897 2900 CALL get_attribute( pids_id, char_lod, rad_sw_in_f%lod, & … … 2905 2908 ! 2906 2909 ! Levelofdetail 2  radiation depends on time_rad, y, x 2907 ELSEIF ( rad_sw_in_f%lod == 2 ) THEN 2910 ELSEIF ( rad_sw_in_f%lod == 2 ) THEN 2908 2911 ALLOCATE( rad_sw_in_f%var3d(0:ntime1,nys:nyn,nxl:nxr) ) 2909 2912 2910 2913 CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var3d, & 2911 2914 nxl, nxr, nys, nyn, 0, ntime1 ) 2912 2915 2913 2916 rad_sw_in_f%from_file = .TRUE. 2914 ELSE 2917 ELSE 2915 2918 message_string = '"rad_sw_in" has no valid lod attribute' 2916 2919 CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 ) 2917 2920 ENDIF 2918 2921 ENDIF 2919 ! 2922 ! 2920 2923 ! Input longwave downwelling. 2921 2924 IF ( check_existence( vars_pids, 'rad_lw_in' ) ) THEN 2922 ! 2925 ! 2923 2926 ! Get _FillValue attribute 2924 2927 CALL get_attribute( pids_id, char_fill, rad_lw_in_f%fill, & 2925 2928 .FALSE., 'rad_lw_in' ) 2926 ! 2929 ! 2927 2930 ! Get levelofdetail 2928 2931 CALL get_attribute( pids_id, char_lod, rad_lw_in_f%lod, & … … 2936 2939 ! 2937 2940 ! Levelofdetail 2  radiation depends on time_rad, y, x 2938 ELSEIF ( rad_lw_in_f%lod == 2 ) THEN 2941 ELSEIF ( rad_lw_in_f%lod == 2 ) THEN 2939 2942 ALLOCATE( rad_lw_in_f%var3d(0:ntime1,nys:nyn,nxl:nxr) ) 2940 2943 2941 2944 CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var3d, & 2942 2945 nxl, nxr, nys, nyn, 0, ntime1 ) 2943 2946 2944 2947 rad_lw_in_f%from_file = .TRUE. 2945 ELSE 2948 ELSE 2946 2949 message_string = '"rad_lw_in" has no valid lod attribute' 2947 2950 CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 ) 2948 2951 ENDIF 2949 2952 ENDIF 2950 ! 2953 ! 2951 2954 ! Input shortwave downwelling, diffuse part. 2952 2955 IF ( check_existence( vars_pids, 'rad_sw_in_dif' ) ) THEN 2953 ! 2956 ! 2954 2957 ! Read _FillValue attribute 2955 2958 CALL get_attribute( pids_id, char_fill, rad_sw_in_dif_f%fill, & 2956 2959 .FALSE., 'rad_sw_in_dif' ) 2957 ! 2960 ! 2958 2961 ! Get levelofdetail 2959 2962 CALL get_attribute( pids_id, char_lod, rad_sw_in_dif_f%lod, & … … 2968 2971 ! 2969 2972 ! Levelofdetail 2  radiation depends on time_rad, y, x 2970 ELSEIF ( rad_sw_in_dif_f%lod == 2 ) THEN 2973 ELSEIF ( rad_sw_in_dif_f%lod == 2 ) THEN 2971 2974 ALLOCATE( rad_sw_in_dif_f%var3d(0:ntime1,nys:nyn,nxl:nxr) ) 2972 2975 2973 2976 CALL get_variable( pids_id, 'rad_sw_in_dif', & 2974 2977 rad_sw_in_dif_f%var3d, & 2975 2978 nxl, nxr, nys, nyn, 0, ntime1 ) 2976 2979 2977 2980 rad_sw_in_dif_f%from_file = .TRUE. 2978 ELSE 2981 ELSE 2979 2982 message_string = '"rad_sw_in_dif" has no valid lod attribute' 2980 2983 CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 ) 2981 2984 ENDIF 2982 2985 ENDIF 2983 ! 2986 ! 2984 2987 ! Finally, close the input file and deallocate temporary arrays 2985 2988 DEALLOCATE( vars_pids ) 2986 2989 2987 2990 CALL close_input_file( pids_id ) 2988 2991 #endif … … 2995 2998 CALL message( 'radiation_init', 'PA0195', 1, 2, 0, 6, 0 ) 2996 2999 ENDIF 2997 3000 2998 3001 IF ( .NOT. time_rad_f%from_file ) THEN 2999 3002 message_string = 'In case of external radiation forcing ' // & … … 3001 3004 CALL message( 'radiation_init', 'PA0196', 1, 2, 0, 6, 0 ) 3002 3005 ENDIF 3003 3006 3004 3007 CALL get_date_time( 0.0_wp, second_of_day=second_of_day ) 3005 3008 … … 3018 3021 ENDIF 3019 3022 ENDIF 3020 3023 3021 3024 IF ( ALLOCATED( rad_lw_in_f%var1d ) ) THEN 3022 3025 IF ( ANY( rad_lw_in_f%var1d == rad_lw_in_f%fill ) ) THEN … … 3026 3029 ENDIF 3027 3030 ENDIF 3028 3031 3029 3032 IF ( ALLOCATED( rad_sw_in_dif_f%var1d ) ) THEN 3030 3033 IF ( ANY( rad_sw_in_dif_f%var1d == rad_sw_in_dif_f%fill ) ) THEN … … 3034 3037 ENDIF 3035 3038 ENDIF 3036 3039 3037 3040 IF ( ALLOCATED( rad_sw_in_f%var3d ) ) THEN 3038 3041 IF ( ANY( rad_sw_in_f%var3d == rad_sw_in_f%fill ) ) THEN … … 3042 3045 ENDIF 3043 3046 ENDIF 3044 3047 3045 3048 IF ( ALLOCATED( rad_lw_in_f%var3d ) ) THEN 3046 3049 IF ( ANY( rad_lw_in_f%var3d == rad_lw_in_f%fill ) ) THEN … … 3050 3053 ENDIF 3051 3054 ENDIF 3052 3055 3053 3056 IF ( ALLOCATED( rad_sw_in_dif_f%var3d ) ) THEN 3054 3057 IF ( ANY( rad_sw_in_dif_f%var3d == rad_sw_in_dif_f%fill ) ) THEN … … 3068 3071 ENDIF 3069 3072 ! 3070 ! All radiation input should have the same level of detail. The sum 3073 ! All radiation input should have the same level of detail. The sum 3071 3074 ! of lods divided by the number of available radiation arrays must be 3072 ! 1 (if all are lod = 1) or 2 (if all are lod = 2). 3075 ! 1 (if all are lod = 1) or 2 (if all are lod = 2). 3073 3076 IF ( REAL( MERGE( rad_lw_in_f%lod, 0, rad_lw_in_f%from_file ) + & 3074 3077 MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) + & … … 3109 3112 CASE ( 'constant' ) 3110 3113 CALL radiation_constant 3111 3114 3112 3115 CASE ( 'external' ) 3113 3116 ! … … 3128 3131 3129 3132 ! 3130 ! If required, read or calculate and write out the SVF 3133 ! If required, read or calculate and write out the SVF 3131 3134 IF ( radiation_interactions .AND. read_svf) THEN 3132 3135 ! … … 3147 3150 3148 3151 ! 3149 ! Adjust radiative fluxes. In case of urban and land surfaces, also 3152 ! Adjust radiative fluxes. In case of urban and land surfaces, also 3150 3153 ! call an initial interaction. 3151 3154 IF ( radiation_interactions ) THEN … … 3171 3174 INTEGER(iwp) :: l !< running index for surface orientation 3172 3175 INTEGER(iwp) :: t !< index of current timestep 3173 INTEGER(iwp) :: tm !< index of previous timestep 3174 3176 INTEGER(iwp) :: tm !< index of previous timestep 3177 3175 3178 LOGICAL :: horizontal !< flag indicating treatment of horinzontal surfaces 3176 3177 REAL(wp) :: fac_dt !< interpolation factor 3179 3180 REAL(wp) :: fac_dt !< interpolation factor 3178 3181 REAL(wp) :: second_of_day_init !< second of the day at model start 3179 3182 3180 TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine 3183 TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine 3181 3184 3182 3185 ! … … 3187 3190 CALL calc_zenith( day_of_year, second_of_day ) 3188 3191 ! 3189 ! Interpolate external radiation on current timestep 3192 ! Interpolate external radiation on current timestep 3190 3193 IF ( time_since_reference_point <= 0.0_wp ) THEN 3191 3194 t = 0 … … 3198 3201 t = t + 1 3199 3202 ENDDO 3200 3203 3201 3204 tm = MAX( t1, 0 ) 3202 3205 3203 3206 fac_dt = ( time_since_reference_point & 3204 3207  time_rad_f%var1d(tm) + dt_3d ) & … … 3207 3210 ENDIF 3208 3211 ! 3209 ! Call clearsky calculation for each surface orientation. 3212 ! Call clearsky calculation for each surface orientation. 3210 3213 ! First, horizontal surfaces 3211 3214 horizontal = .TRUE. … … 3223 3226 CALL radiation_external_surf 3224 3227 ENDDO 3225 3228 3226 3229 CONTAINS 3227 3230 … … 3229 3232 3230 3233 USE control_parameters 3231 3234 3232 3235 IMPLICIT NONE 3233 3236 3234 INTEGER(iwp) :: i !< grid index along xdimension 3235 INTEGER(iwp) :: j !< grid index along ydimension 3236 INTEGER(iwp) :: k !< grid index along zdimension 3237 INTEGER(iwp) :: i !< grid index along xdimension 3238 INTEGER(iwp) :: j !< grid index along ydimension 3239 INTEGER(iwp) :: k !< grid index along zdimension 3237 3240 INTEGER(iwp) :: m !< running index for surface elements 3238 3239 REAL(wp) :: lw_in !< downwelling longwave radiation, interpolated value 3241 3242 REAL(wp) :: lw_in !< downwelling longwave radiation, interpolated value 3240 3243 REAL(wp) :: sw_in !< downwelling shortwave radiation, interpolated value 3241 REAL(wp) :: sw_in_dif !< downwelling diffuse shortwave radiation, interpolated value 3244 REAL(wp) :: sw_in_dif !< downwelling diffuse shortwave radiation, interpolated value 3242 3245 3243 3246 IF ( surf%ns < 1 ) RETURN … … 3245 3248 ! levelofdetail = 1. Note, here it must be distinguished between 3246 3249 ! averaged radiation and nonaveraged radiation for the upwelling 3247 ! fluxes. 3250 ! fluxes. 3248 3251 IF ( rad_sw_in_f%lod == 1 ) THEN 3249 3252 3250 3253 sw_in = ( 1.0_wp  fac_dt ) * rad_sw_in_f%var1d(tm) & 3251 3254 + fac_dt * rad_sw_in_f%var1d(t) 3252 3255 3253 3256 lw_in = ( 1.0_wp  fac_dt ) * rad_lw_in_f%var1d(tm) & 3254 3257 + fac_dt * rad_lw_in_f%var1d(t) 3255 3258 ! 3256 ! Limit shortwave incoming radiation to positive values, in order 3259 ! Limit shortwave incoming radiation to positive values, in order 3257 3260 ! to overcome possible observation errors. 3258 3261 sw_in = MAX( 0.0_wp, sw_in ) 3259 3262 sw_in = MERGE( sw_in, 0.0_wp, sun_up ) 3260 3261 surf%rad_sw_in = sw_in 3263 3264 surf%rad_sw_in = sw_in 3262 3265 surf%rad_lw_in = lw_in 3263 3266 3264 3267 IF ( average_radiation ) THEN 3265 3268 surf%rad_sw_out = albedo_urb * surf%rad_sw_in 3266 3269 3267 3270 surf%rad_lw_out = emissivity_urb * sigma_sb * t_rad_urb**4 & 3268 3271 + ( 1.0_wp  emissivity_urb ) * surf%rad_lw_in 3269 3272 3270 3273 surf%rad_net = surf%rad_sw_in  surf%rad_sw_out & 3271 3274 + surf%rad_lw_in  surf%rad_lw_out 3272 3275 3273 3276 surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb & 3274 3277 * sigma_sb & … … 3284 3287 surf%albedo(m,ind_wat_win) ) & 3285 3288 * surf%rad_sw_in(m) 3286 3289 3287 3290 surf%rad_lw_out(m) = ( surf%frac(m,ind_veg_wall) * & 3288 3291 surf%emissivity(m,ind_veg_wall) & … … 3294 3297 * sigma_sb & 3295 3298 * ( surf%pt_surface(m) * exner(k) )**4 3296 3299 3297 3300 surf%rad_lw_out_change_0(m) = & 3298 3301 ( surf%frac(m,ind_veg_wall) * & … … 3305 3308 * ( surf%pt_surface(m) * exner(k) )**3 3306 3309 ENDDO 3307 3310 3308 3311 ENDIF 3309 3312 ! 3310 ! If diffuse shortwave radiation is available, store it on 3313 ! If diffuse shortwave radiation is available, store it on 3311 3314 ! the respective files. 3312 IF ( rad_sw_in_dif_f%from_file ) THEN 3315 IF ( rad_sw_in_dif_f%from_file ) THEN 3313 3316 sw_in_dif= ( 1.0_wp  fac_dt ) * rad_sw_in_dif_f%var1d(tm) & 3314 3317 + fac_dt * rad_sw_in_dif_f%var1d(t) 3315 3318 3316 3319 IF ( ALLOCATED( rad_sw_in_diff ) ) rad_sw_in_diff = sw_in_dif 3317 3320 IF ( ALLOCATED( rad_sw_in_dir ) ) rad_sw_in_dir = sw_in & 3318 3321  sw_in_dif 3319 ! 3320 ! Diffuse longwave radiation equals the total downwelling 3322 ! 3323 ! Diffuse longwave radiation equals the total downwelling 3321 3324 ! longwave radiation 3322 3325 IF ( ALLOCATED( rad_lw_in_diff ) ) rad_lw_in_diff = lw_in … … 3330 3333 j = surf%j(m) 3331 3334 k = surf%k(m) 3332 3335 3333 3336 surf%rad_sw_in(m) = ( 1.0_wp  fac_dt ) & 3334 3337 * rad_sw_in_f%var3d(tm,j,i) & 3335 + fac_dt * rad_sw_in_f%var3d(t,j,i) 3336 ! 3337 ! Limit shortwave incoming radiation to positive values, in 3338 + fac_dt * rad_sw_in_f%var3d(t,j,i) 3339 ! 3340 ! Limit shortwave incoming radiation to positive values, in 3338 3341 ! order to overcome possible observation errors. 3339 3342 surf%rad_sw_in(m) = MAX( 0.0_wp, surf%rad_sw_in(m) ) 3340 3343 surf%rad_sw_in(m) = MERGE( surf%rad_sw_in(m), 0.0_wp, sun_up ) 3341 3344 3342 3345 surf%rad_lw_in(m) = ( 1.0_wp  fac_dt ) & 3343 3346 * rad_lw_in_f%var3d(tm,j,i) & 3344 + fac_dt * rad_lw_in_f%var3d(t,j,i) 3345 ! 3346 ! Weighted average according to surface fraction. 3347 + fac_dt * rad_lw_in_f%var3d(t,j,i) 3348 ! 3349 ! Weighted average according to surface fraction. 3347 3350 surf%rad_sw_out(m) = ( surf%frac(m,ind_veg_wall) * & 3348 3351 surf%albedo(m,ind_veg_wall) & … … 3376 3379 + surf%rad_lw_in(m)  surf%rad_lw_out(m) 3377 3380 ! 3378 ! If diffuse shortwave radiation is available, store it on 3379 ! the respective files. 3381 ! If diffuse shortwave radiation is available, store it on 3382 ! the respective files. 3380 3383 IF ( rad_sw_in_dif_f%from_file ) THEN 3381 3384 IF ( ALLOCATED( rad_sw_in_diff ) ) & … … 3384 3387 + fac_dt * rad_sw_in_dif_f%var3d(t,j,i) 3385 3388 ! 3386 ! dir = sw_in  sw_in_dif. 3389 ! dir = sw_in  sw_in_dif. 3387 3390 IF ( ALLOCATED( rad_sw_in_dir ) ) & 3388 3391 rad_sw_in_dir(j,i) = surf%rad_sw_in(m)  & 3389 3392 rad_sw_in_diff(j,i) 3390 ! 3391 ! Diffuse longwave radiation equals the total downwelling 3393 ! 3394 ! Diffuse longwave radiation equals the total downwelling 3392 3395 ! longwave radiation 3393 3396 IF ( ALLOCATED( rad_lw_in_diff ) ) & … … 3399 3402 ENDIF 3400 3403 ! 3401 ! Store radiation also on 2D arrays, which are still used for 3404 ! Store radiation also on 2D arrays, which are still used for 3402 3405 ! directdiffuse splitting. Note, this is only required 3403 3406 ! for horizontal surfaces, which covers all x,y position. … … 3406 3409 i = surf%i(m) 3407 3410 j = surf%j(m) 3408 3411 3409 3412 rad_sw_in(0,j,i) = surf%rad_sw_in(m) 3410 3413 rad_lw_in(0,j,i) = surf%rad_lw_in(m) … … 3413 3416 ENDDO 3414 3417 ENDIF 3415 3418 3416 3419 END SUBROUTINE radiation_external_surf 3417 3420 3418 END SUBROUTINE radiation_external 3419 3421 END SUBROUTINE radiation_external 3422 3420 3423 !! 3421 3424 ! Description: … … 3431 3434 LOGICAL :: horizontal !< flag indicating treatment of horinzontal surfaces 3432 3435 3433 REAL(wp) :: pt1 !< potential temperature at first grid level or mean value at urban layer top 3434 REAL(wp) :: pt1_l !< potential temperature at first grid level or mean value at urban layer top at local subdomain 3435 REAL(wp) :: ql1 !< liquid water mixing ratio at first grid level or mean value at urban layer top 3436 REAL(wp) :: ql1_l !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain 3437 3438 TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine 3436 REAL(wp) :: pt1 !< potential temperature at first grid level or mean value at urban layer top 3437 REAL(wp) :: pt1_l !< potential temperature at first grid level or mean value at urban layer top at local subdomain 3438 REAL(wp) :: ql1 !< liquid water mixing ratio at first grid level or mean value at urban layer top 3439 REAL(wp) :: ql1_l !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain 3440 3441 TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine 3439 3442 3440 3443 ! … … 3452 3455 ! Calculate value of the Exner function at model surface 3453 3456 ! 3454 ! In case averaged radiation is used, calculate mean temperature and 3457 ! In case averaged radiation is used, calculate mean temperature and 3455 3458 ! liquid water mixing ratio at the urbanlayer top. 3456 3459 IF ( average_radiation ) THEN … … 3461 3464 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) ) 3462 3465 3463 #if defined( __parallel ) 3466 #if defined( __parallel ) 3464 3467 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 3465 3468 CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) … … 3477 3480 ENDIF 3478 3481 #else 3479 pt1 = pt1_l 3482 pt1 = pt1_l 3480 3483 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = ql1_l 3481 3484 #endif … … 3487 3490 ENDIF 3488 3491 ! 3489 ! Call clearsky calculation for each surface orientation. 3492 ! Call clearsky calculation for each surface orientation. 3490 3493 ! First, horizontal surfaces 3491 3494 horizontal = .TRUE. … … 3518 3521 3519 3522 ! 3520 ! Calculate radiation fluxes and net radiation (rad_net) assuming 3521 ! homogeneous urban radiation conditions. 3522 IF ( average_radiation ) THEN 3523 ! Calculate radiation fluxes and net radiation (rad_net) assuming 3524 ! homogeneous urban radiation conditions. 3525 IF ( average_radiation ) THEN 3523 3526 3524 3527 k = nz_urban_t … … 3526 3529 surf%rad_sw_in = solar_constant * sky_trans * cos_zenith 3527 3530 surf%rad_sw_out = albedo_urb * surf%rad_sw_in 3528 3531 3529 3532 surf%rad_lw_in = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4 3530 3533 … … 3539 3542 3540 3543 ! 3541 ! Calculate radiation fluxes and net radiation (rad_net) for each surface 3544 ! Calculate radiation fluxes and net radiation (rad_net) for each surface 3542 3545 ! element. 3543 3546 ELSE … … 3551 3554 3552 3555 ! 3553 ! Weighted average according to surface fraction. 3554 ! ATTENTION: when radiation interactions are switched on the 3555 ! calculated fluxes below are not actually used as they are 3556 ! Weighted average according to surface fraction. 3557 ! ATTENTION: when radiation interactions are switched on the 3558 ! calculated fluxes below are not actually used as they are 3556 3559 ! overwritten in radiation_interaction. 3557 3560 surf%rad_sw_out(m) = ( surf%frac(m,ind_veg_wall) * & … … 3611 3614 ENDDO 3612 3615 ENDIF 3613 3616 3614 3617 END SUBROUTINE radiation_clearsky_surf 3615 3618 … … 3628 3631 3629 3632 INTEGER(iwp) :: l !< running index for surface orientation 3630 3633 3631 3634 LOGICAL :: horizontal !< flag indicating treatment of horinzontal surfaces 3632 3635 3633 REAL(wp) :: pt1 !< potential temperature at first grid level or mean value at urban layer top 3634 REAL(wp) :: pt1_l !< potential temperature at first grid level or mean value at urban layer top at local subdomain 3635 REAL(wp) :: ql1 !< liquid water mixing ratio at first grid level or mean value at urban layer top 3636 REAL(wp) :: ql1_l !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain 3636 REAL(wp) :: pt1 !< potential temperature at first grid level or mean value at urban layer top 3637 REAL(wp) :: pt1_l !< potential temperature at first grid level or mean value at urban layer top at local subdomain 3638 REAL(wp) :: ql1 !< liquid water mixing ratio at first grid level or mean value at urban layer top 3639 REAL(wp) :: ql1_l !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain 3637 3640 3638 3641 TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine 3639 3642 3640 3643 ! 3641 ! In case averaged radiation is used, calculate mean temperature and 3644 ! In case averaged radiation is used, calculate mean temperature and 3642 3645 ! liquid water mixing ratio at the urbanlayer top. 3643 IF ( average_radiation ) THEN 3646 IF ( average_radiation ) THEN 3644 3647 pt1 = 0.0_wp 3645 3648 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = 0.0_wp … … 3648 3651 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) ) 3649 3652 3650 #if defined( __parallel ) 3653 #if defined( __parallel ) 3651 3654 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 3652 3655 CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) … … 3729 3732 ELSE 3730 3733 ! 3731 ! Determine index offset between surface element and adjacent 3734 ! Determine index offset between surface element and adjacent 3732 3735 ! atmospheric grid point 3733 3736 ioff = surf%ioff … … 3753 3756 3754 3757 ! 3755 ! Weighted average according to surface fraction. 3758 ! Weighted average according to surface fraction. 3756 3759 surf%rad_lw_out(m) = ( surf%frac(m,ind_veg_wall) * & 3757 3760 surf%emissivity(m,ind_veg_wall) & … … 3802 3805 3803 3806 END SUBROUTINE radiation_constant_surf 3804 3807 3805 3808 3806 3809 END SUBROUTINE radiation_constant … … 3815 3818 3816 3819 IMPLICIT NONE 3817 3820 3818 3821 INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file 3819 3820 3821 3822 3823 3824 3822 3825 ! 3823 3826 ! Write radiation model header … … 3834 3837 ELSEIF ( radiation_scheme == "external" ) THEN 3835 3838 WRITE( io, 14 ) 3836 ENDIF 3839 ENDIF 3837 3840 3838 3841 IF ( albedo_type_f%from_file .OR. vegetation_type_f%from_file .OR. & … … 3840 3843 building_type_f%from_file ) THEN 3841 3844 WRITE( io, 13 ) 3842 ELSE 3845 ELSE 3843 3846 IF ( albedo_type == 0 ) THEN 3844 3847 WRITE( io, 7 ) albedo … … 3850 3853 WRITE( io, 9 ) 3851 3854 ENDIF 3852 3855 3853 3856 WRITE( io, 12 ) dt_radiation 3854 3857 3855 3858 3856 3859 3 FORMAT (//' Radiation model information:'/ & … … 3873 3876 3874 3877 END SUBROUTINE radiation_header 3875 3878 3876 3879 3877 3880 !! … … 3885 3888 IMPLICIT NONE 3886 3889 3887 CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file 3888 3890 CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file 3891 3889 3892 NAMELIST /radiation_par/ albedo, albedo_lw_dif, albedo_lw_dir, & 3890 3893 albedo_sw_dif, albedo_sw_dir, albedo_type, & … … 3903 3906 unscheduled_radiation_calls 3904 3907 3905 3908 3906 3909 NAMELIST /radiation_parameters/ albedo, albedo_lw_dif, albedo_lw_dir, & 3907 3910 albedo_sw_dif, albedo_sw_dir, albedo_type, & … … 3919 3922 svfnorm_report_thresh, sw_radiation, & 3920 3923 unscheduled_radiation_calls 3921 3924 3922 3925 line = ' ' 3923 3926 3924 3927 ! 3925 3928 ! Try to find radiation model namelist … … 3979 3982 3980 3983 14 CONTINUE 3981 3984 3982 3985 END SUBROUTINE radiation_parin 3983 3986 … … 4039 4042 zenith(0) = cos_zenith 4040 4043 ! 4041 ! Calculate surface albedo. In case average radiation is applied, 4044 ! Calculate surface albedo. In case average radiation is applied, 4042 4045 ! this is not required. 4043 4046 #if defined( __netcdf ) … … 4071 4074 IF ( average_radiation ) THEN 4072 4075 ! 4073 ! Determine minimum topography top index. 4076 ! Determine minimum topography top index. 4074 4077 k_topo_l = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ) 4075 4078 #if defined( __parallel ) … … 4079 4082 k_topo = k_topo_l 4080 4083 #endif 4081 4084 4082 4085 rrtm_asdir(1) = albedo_urb 4083 4086 rrtm_asdif(1) = albedo_urb … … 4087 4090 rrtm_emis = emissivity_urb 4088 4091 ! 4089 ! Calculate mean pt profile. 4092 ! Calculate mean pt profile. 4090 4093 CALL calc_mean_profile( pt, 4 ) 4091 4094 pt_av = hom(:, 1, 4, 0) 4092 4095 4093 4096 IF ( humidity ) THEN 4094 4097 CALL calc_mean_profile( q, 41 ) … … 4104 4107 ! average ql is now in hom(:, 1, 54, 0) 4105 4108 ql_av = hom(:, 1, 54, 0) 4106 4109 4107 4110 DO k = nzb+1, nzt+1 4108 4111 rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp & … … 4126 4129 4127 4130 ! 4128 ! Avoid temperature/humidity jumps at the top of the PALM domain by 4131 ! Avoid temperature/humidity jumps at the top of the PALM domain by 4129 4132 ! linear interpolation from nzt+2 to nzt+7. Jumps are induced by 4130 4133 ! discrepancies between the values in the domain and those above that … … 4163 4166 rrtm_cliqwp(0,k) = ql_av(k) * 1000._wp * & 4164 4167 (rrtm_plev(0,k)  rrtm_plev(0,k+1)) & 4165 * 100._wp / g 4168 * 100._wp / g 4166 4169 4167 4170 IF ( rrtm_cliqwp(0,k) > 0._wp ) THEN … … 4189 4192 ! Set surface temperature 4190 4193 rrtm_tsfc = t_rad_urb 4191 4192 IF ( lw_radiation ) THEN 4194 4195 IF ( lw_radiation ) THEN 4193 4196 ! 4194 4197 ! Due to technical reasons, copy optical depth to dummy arguments 4195 4198 ! which are allocated on the exact size as the rrtmg_lw is called. 4196 4199 ! As one dimesion is allocated with zero size, compiler complains 4197 ! that rank of the array does not match that of the 4198 ! assumedshaped arguments in the RRTMG library. In order to 4199 ! avoid this, write to dummy arguments and give pass the entire 4200 ! dummy array. Seems to be the only existing workaround. 4200 ! that rank of the array does not match that of the 4201 ! assumedshaped arguments in the RRTMG library. In order to 4202 ! avoid this, write to dummy arguments and give pass the entire 4203 ! dummy array. Seems to be the only existing workaround. 4201 4204 ALLOCATE( rrtm_lw_taucld_dum(1:nbndlw+1,0:0,k_topo+1:nzt_rad+1) ) 4202 4205 ALLOCATE( rrtm_lw_tauaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1) ) … … 4206 4209 rrtm_lw_tauaer_dum = & 4207 4210 rrtm_lw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1) 4208 4209 CALL rrtmg_lw( 1, & 4211 4212 CALL rrtmg_lw( 1, & 4210 4213 nzt_radk_topo, & 4211 4214 rrtm_icld, & … … 4234 4237 rrtm_cicewp(:,k_topo+1:), & 4235 4238 rrtm_cliqwp(:,k_topo+1:), & 4236 rrtm_reice(:,k_topo+1:), & 4239 rrtm_reice(:,k_topo+1:), & 4237 4240 rrtm_reliq(:,k_topo+1:), & 4238 4241 rrtm_lw_tauaer_dum, & … … 4245 4248 rrtm_lwuflx_dt(:,k_topo:), & 4246 4249 rrtm_lwuflxc_dt(:,k_topo:) ) 4247 4250 4248 4251 DEALLOCATE ( rrtm_lw_taucld_dum ) 4249 4252 DEALLOCATE ( rrtm_lw_tauaer_dum ) … … 4256 4259 rad_lw_in_diff(:,:) = rad_lw_in(k_topo,:,:) 4257 4260 ! 4258 ! Save heating rates (convert from K/d to K/h). 4261 ! Save heating rates (convert from K/d to K/h). 4259 4262 ! Further, even though an aggregated radiation is computed, map 4260 ! signlecolumn profiles on top of any topography, in order to 4263 ! signlecolumn profiles on top of any topography, in order to 4261 4264 ! obtain correct near surface radiation heating/cooling rates. 4262 4265 DO i = nxl, nxr … … 4274 4277 IF ( sw_radiation .AND. sun_up ) THEN 4275 4278 ! 4276 ! Due to technical reasons, copy optical depths and other 4277 ! to dummy arguments which are allocated on the exact size as the 4279 ! Due to technical reasons, copy optical depths and other 4280 ! to dummy arguments which are allocated on the exact size as the 4278 4281 ! rrtmg_sw is called. 4279 4282 ! As one dimesion is allocated with zero size, compiler complains 4280 ! that rank of the array does not match that of the 4281 ! assumedshaped arguments in the RRTMG library. In order to 4282 ! avoid this, write to dummy arguments and give pass the entire 4283 ! dummy array. Seems to be the only existing workaround. 4283 ! that rank of the array does not match that of the 4284 ! assumedshaped arguments in the RRTMG library. In order to 4285 ! avoid this, write to dummy arguments and give pass the entire 4286 ! dummy array. Seems to be the only existing workaround. 4284 4287 ALLOCATE( rrtm_sw_taucld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) ) 4285 4288 ALLOCATE( rrtm_sw_ssacld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) ) … … 4290 4293 ALLOCATE( rrtm_sw_asmaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) ) 4291 4294 ALLOCATE( rrtm_sw_ecaer_dum(0:0,k_topo+1:nzt_rad+1,1:naerec+1) ) 4292 4295 4293 4296 rrtm_sw_taucld_dum = rrtm_sw_taucld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) 4294 4297 rrtm_sw_ssacld_dum = rrtm_sw_ssacld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) … … 4309 4312 rrtm_tlev(:,k_topo+1:nzt_rad+2), & 4310 4313 rrtm_tsfc, & 4311 rrtm_h2ovmr(:,k_topo+1:nzt_rad+1), & 4312 rrtm_o3vmr(:,k_topo+1:nzt_rad+1), & 4314 rrtm_h2ovmr(:,k_topo+1:nzt_rad+1), & 4315 rrtm_o3vmr(:,k_topo+1:nzt_rad+1), & 4313 4316 rrtm_co2vmr(:,k_topo+1:nzt_rad+1), & 4314 4317 rrtm_ch4vmr(:,k_topo+1:nzt_rad+1), & 4315 4318 rrtm_n2ovmr(:,k_topo+1:nzt_rad+1), & 4316 4319 rrtm_o2vmr(:,k_topo+1:nzt_rad+1), & 4317 rrtm_asdir, & 4320 rrtm_asdir, & 4318 4321 rrtm_asdif, & 4319 4322 rrtm_aldir, & … … 4339 4342 rrtm_sw_asmaer_dum, & 4340 4343 rrtm_sw_ecaer_dum, & 4341 rrtm_swuflx(:,k_topo:nzt_rad+1), & 4342 rrtm_swdflx(:,k_topo:nzt_rad+1), & 4343 rrtm_swhr(:,k_topo+1:nzt_rad+1), & 4344 rrtm_swuflxc(:,k_topo:nzt_rad+1), & 4344 rrtm_swuflx(:,k_topo:nzt_rad+1), & 4345 rrtm_swdflx(:,k_topo:nzt_rad+1), & 4346 rrtm_swhr(:,k_topo+1:nzt_rad+1), & 4347 rrtm_swuflxc(:,k_topo:nzt_rad+1), & 4345 4348 rrtm_swdflxc(:,k_topo:nzt_rad+1), & 4346 4349 rrtm_swhrc(:,k_topo+1:nzt_rad+1), & 4347 4350 rrtm_dirdflux(:,k_topo:nzt_rad+1), & 4348 4351 rrtm_difdflux(:,k_topo:nzt_rad+1) ) 4349 4352 4350 4353 DEALLOCATE( rrtm_sw_taucld_dum ) 4351 4354 DEALLOCATE( rrtm_sw_ssacld_dum ) … … 4356 4359 DEALLOCATE( rrtm_sw_asmaer_dum ) 4357 4360 DEALLOCATE( rrtm_sw_ecaer_dum ) 4358 4361 4359 4362 ! 4360 4363 ! Save radiation fluxes for the entire depth of the model domain … … 4382 4385 ENDIF 4383 4386 ! 4384 ! RRTMG is called for each (j,i) grid point separately, starting at the 4387 ! RRTMG is called for each (j,i) grid point separately, starting at the 4385 4388 ! highest topography level. Here no RTM is used since average_radiation is false 4386 4389 ELSE … … 4410 4413 rrtm_tlay(0,k) = pt(k,j,i) * exner(k) & 4411 4414 + lv_d_cp * ql(k,j,i) 4412 rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i) 4415 rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i) 4413 4416 ENDDO 4414 4417 ELSE … … 4420 4423 DO k = nzb+1, nzt+1 4421 4424 rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i) 4422 ENDDO 4425 ENDDO 4423 4426 ELSE 4424 4427 rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp … … 4427 4430 4428 4431 ! 4429 ! Avoid temperature/humidity jumps at the top of the LES domain by 4432 ! Avoid temperature/humidity jumps at the top of the LES domain by 4430 4433 ! linear interpolation from nzt+2 to nzt+7 4431 4434 DO k = nzt+2, nzt+7 … … 4463 4466 rrtm_cliqwp(0,k) = ql(k,j,i) * 1000.0_wp * & 4464 4467 (rrtm_plev(0,k)  rrtm_plev(0,k+1)) & 4465 * 100.0_wp / g 4468 * 100.0_wp / g 4466 4469 4467 4470 IF ( rrtm_cliqwp(0,k) > 0.0_wp ) THEN … … 4473 4476 IF ( bulk_cloud_model ) THEN 4474 4477 ! 4475 ! Calculete effective droplet radius. In case of using 4476 ! cloud_scheme = 'morrison' and a non reasonable number 4477 ! of cloud droplets the inital aerosol number 4478 ! Calculete effective droplet radius. In case of using 4479 ! cloud_scheme = 'morrison' and a non reasonable number 4480 ! of cloud droplets the inital aerosol number 4478 4481 ! concentration is considered. 4479 4482 IF ( microphysics_morrison ) THEN … … 4482 4485 ELSE 4483 4486 nc_rad = na_init 4484 ENDIF 4487 ENDIF 4485 4488 ELSE 4486 4489 nc_rad = nc_const 4487 ENDIF 4490 ENDIF 4488 4491 4489 4492 rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i) & … … 4525 4528 4526 4529 ! 4527 ! Write surface emissivity and surface temperature at current 4530 ! Write surface emissivity and surface temperature at current 4528 4531 ! surface element on RRTMGshaped array. 4529 4532 ! Please note, as RRTMG is a single column model, surface attributes 4530 ! are only obtained from horizontally aligned surfaces (for 4531 ! simplicity). Taking surface attributes from horizontal and 4532 ! vertical walls would lead to multiple solutions. 4533 ! are only obtained from horizontally aligned surfaces (for 4534 ! simplicity). Taking surface attributes from horizontal and 4535 ! vertical walls would lead to multiple solutions. 4533 4536 ! Moreover, for natural and urbantype surfaces, several surface 4534 ! classes can exist at a surface element next to each other. 4535 ! To obtain bulk parameters, apply a weighted average for these 4536 ! surfaces. 4537 ! classes can exist at a surface element next to each other. 4538 ! To obtain bulk parameters, apply a weighted average for these 4539 ! surfaces. 4537 4540 DO m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i) 4538 4541 rrtm_emis = surf_lsm_h%frac(m,ind_veg_wall) * & 4539 4542 surf_lsm_h%emissivity(m,ind_veg_wall) + & 4540 4543 surf_lsm_h%frac(m,ind_pav_green) * & 4541 surf_lsm_h%emissivity(m,ind_pav_green) + & 4544 surf_lsm_h%emissivity(m,ind_pav_green) + & 4542 4545 surf_lsm_h%frac(m,ind_wat_win) * & 4543 4546 surf_lsm_h%emissivity(m,ind_wat_win) 4544 4547 rrtm_tsfc = surf_lsm_h%pt_surface(m) * exner(nzb) 4545 ENDDO 4548 ENDDO 4546 4549 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) 4547 4550 rrtm_emis = surf_usm_h%frac(m,ind_veg_wall) * & 4548 4551 surf_usm_h%emissivity(m,ind_veg_wall) + & 4549 4552 surf_usm_h%frac(m,ind_pav_green) * & 4550 surf_usm_h%emissivity(m,ind_pav_green) + & 4553 surf_usm_h%emissivity(m,ind_pav_green) + & 4551 4554 surf_usm_h%frac(m,ind_wat_win) * & 4552 4555 surf_usm_h%emissivity(m,ind_wat_win) … … 4562 4565 ! which are allocated on the exact size as the rrtmg_lw is called. 4563 4566 ! As one dimesion is allocated with zero size, compiler complains 4564 ! that rank of the array does not match that of the 4565 ! assumedshaped arguments in the RRTMG library. In order to 4566 ! avoid this, write to dummy arguments and give pass the entire 4567 ! dummy array. Seems to be the only existing workaround. 4567 ! that rank of the array does not match that of the 4568 ! assumedshaped arguments in the RRTMG library. In order to 4569 ! avoid this, write to dummy arguments and give pass the entire 4570 ! dummy array. Seems to be the only existing workaround. 4568 4571 ALLOCATE( rrtm_lw_taucld_dum(1:nbndlw+1,0:0,k_topo+1:nzt_rad+1) ) 4569 4572 ALLOCATE( rrtm_lw_tauaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1) ) … … 4574 4577 rrtm_lw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1) 4575 4578 4576 CALL rrtmg_lw( 1, & 4579 CALL rrtmg_lw( 1, & 4577 4580 nzt_radk_topo, & 4578 4581 rrtm_icld, & … … 4601 4604 rrtm_cicewp(:,k_topo+1:nzt_rad+1), & 4602 4605 rrtm_cliqwp(:,k_topo+1:nzt_rad+1), & 4603 rrtm_reice(:,k_topo+1:nzt_rad+1), & 4606 rrtm_reice(:,k_topo+1:nzt_rad+1), & 4604 4607 rrtm_reliq(:,k_topo+1:nzt_rad+1), & 4605 4608 rrtm_lw_tauaer_dum, & … … 4630 4633 4631 4634 ! 4632 ! Save surface radiative fluxes and change in LW heating rate 4635 ! Save surface radiative fluxes and change in LW heating rate 4633 4636 ! onto respective surface elements 4634 4637 ! Horizontal surfaces … … 4638 4641 surf_lsm_h%rad_lw_out(m) = rrtm_lwuflx(0,k_topo) 4639 4642 surf_lsm_h%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo) 4640 ENDDO 4643 ENDDO 4641 4644 DO m = surf_usm_h%start_index(j,i), & 4642 4645 surf_usm_h%end_index(j,i) … … 4644 4647 surf_usm_h%rad_lw_out(m) = rrtm_lwuflx(0,k_topo) 4645 4648 surf_usm_h%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo) 4646 ENDDO 4647 ! 4648 ! Vertical surfaces. Fluxes are obtain at vertical level of the 4649 ENDDO 4650 ! 4651 ! Vertical surfaces. Fluxes are obtain at vertical level of the 4649 4652 ! respective surface element 4650 4653 DO l = 0, 3 … … 4655 4658 surf_lsm_v(l)%rad_lw_out(m) = rrtm_lwuflx(0,k) 4656 4659 surf_lsm_v(l)%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k) 4657 ENDDO 4660 ENDDO 4658 4661 DO m = surf_usm_v(l)%start_index(j,i), & 4659 4662 surf_usm_v(l)%end_index(j,i) … … 4662 4665 surf_usm_v(l)%rad_lw_out(m) = rrtm_lwuflx(0,k) 4663 4666 surf_usm_v(l)%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k) 4664 ENDDO 4667 ENDDO 4665 4668 ENDDO 4666 4669 … … 4669 4672 IF ( sw_radiation .AND. sun_up ) THEN 4670 4673 ! 4671 ! Get albedo for direct/diffusive long/shortwave radiation at 4672 ! current (y,x)location from surface variables. 4673 ! Only obtain it from horizontal surfaces, as RRTMG is a single 4674 ! Get albedo for direct/diffusive long/shortwave radiation at 4675 ! current (y,x)location from surface variables. 4676 ! Only obtain it from horizontal surfaces, as RRTMG is a single 4674 4677 ! column model 4675 ! (Please note, only one loop will entered, controlled by 4678 ! (Please note, only one loop will entered, controlled by 4676 4679 ! startend index.) 4677 4680 DO m = surf_lsm_h%start_index(j,i), & 4678 4681 surf_lsm_h%end_index(j,i) 4679 4682 rrtm_asdir(1) = SUM( surf_lsm_h%frac(m,:) * & 4680 surf_lsm_h%rrtm_asdir( :,m) )4683 surf_lsm_h%rrtm_asdir(m,:) ) 4681 4684 rrtm_asdif(1) = SUM( surf_lsm_h%frac(m,:) * & 4682 surf_lsm_h%rrtm_asdif( :,m) )4685 surf_lsm_h%rrtm_asdif(m,:) ) 4683 4686 rrtm_aldir(1) = SUM( surf_lsm_h%frac(m,:) * & 4684 surf_lsm_h%rrtm_aldir( :,m) )4687 surf_lsm_h%rrtm_aldir(m,:) ) 4685 4688 rrtm_aldif(1) = SUM( surf_lsm_h%frac(m,:) * & 4686 surf_lsm_h%rrtm_aldif( :,m) )4687 ENDDO 4689 surf_lsm_h%rrtm_aldif(m,:) ) 4690 ENDDO 4688 4691 DO m = surf_usm_h%start_index(j,i), & 4689 4692 surf_usm_h%end_index(j,i) 4690 4693 rrtm_asdir(1) = SUM( surf_usm_h%frac(m,:) * & 4691 surf_usm_h%rrtm_asdir( :,m) )4694 surf_usm_h%rrtm_asdir(m,:) ) 4692 4695 rrtm_asdif(1) = SUM( surf_usm_h%frac(m,:) * & 4693 surf_usm_h%rrtm_asdif( :,m) )4696 surf_usm_h%rrtm_asdif(m,:) ) 4694 4697 rrtm_aldir(1) = SUM( surf_usm_h%frac(m,:) * & 4695 surf_usm_h%rrtm_aldir( :,m) )4698 surf_usm_h%rrtm_aldir(m,:) ) 4696 4699 rrtm_aldif(1) = SUM( surf_usm_h%frac(m,:) * & 4697 surf_usm_h%rrtm_aldif( :,m) )4700 surf_usm_h%rrtm_aldif(m,:) ) 4698 4701 ENDDO 4699 4702 ! 4700 ! Due to technical reasons, copy optical depths and other 4701 ! to dummy arguments which are allocated on the exact size as the 4703 ! Due to technical reasons, copy optical depths and other 4704 ! to dummy arguments which are allocated on the exact size as the 4702 4705 ! rrtmg_sw is called. 4703 4706 ! As one dimesion is allocated with zero size, compiler complains 4704 ! that rank of the array does not match that of the 4705 ! assumedshaped arguments in the RRTMG library. In order to 4706 ! avoid this, write to dummy arguments and give pass the entire 4707 ! dummy array. Seems to be the only existing workaround. 4707 ! that rank of the array does not match that of the 4708 ! assumedshaped arguments in the RRTMG library. In order to 4709 ! avoid this, write to dummy arguments and give pass the entire 4710 ! dummy array. Seems to be the only existing workaround. 4708 4711 ALLOCATE( rrtm_sw_taucld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) ) 4709 4712 ALLOCATE( rrtm_sw_ssacld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) ) … … 4714 4717 ALLOCATE( rrtm_sw_asmaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) ) 4715 4718 ALLOCATE( rrtm_sw_ecaer_dum(0:0,k_topo+1:nzt_rad+1,1:naerec+1) ) 4716 4719 4717 4720 rrtm_sw_taucld_dum = rrtm_sw_taucld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) 4718 4721 rrtm_sw_ssacld_dum = rrtm_sw_ssacld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) … … 4733 4736 rrtm_tlev(:,k_topo+1:nzt_rad+2), & 4734 4737 rrtm_tsfc, & 4735 rrtm_h2ovmr(:,k_topo+1:nzt_rad+1), & 4736 rrtm_o3vmr(:,k_topo+1:nzt_rad+1), & 4738 rrtm_h2ovmr(:,k_topo+1:nzt_rad+1), & 4739 rrtm_o3vmr(:,k_topo+1:nzt_rad+1), & 4737 4740 rrtm_co2vmr(:,k_topo+1:nzt_rad+1), & 4738 4741 rrtm_ch4vmr(:,k_topo+1:nzt_rad+1), & 4739 4742 rrtm_n2ovmr(:,k_topo+1:nzt_rad+1), & 4740 4743 rrtm_o2vmr(:,k_topo+1:nzt_rad+1), & 4741 rrtm_asdir, & 4744 rrtm_asdir, & 4742 4745 rrtm_asdif, & 4743 4746 rrtm_aldir, & … … 4763 4766 rrtm_sw_asmaer_dum, & 4764 4767 rrtm_sw_ecaer_dum, & 4765 rrtm_swuflx(:,k_topo:nzt_rad+1), & 4766 rrtm_swdflx(:,k_topo:nzt_rad+1), & 4767 rrtm_swhr(:,k_topo+1:nzt_rad+1), & 4768 rrtm_swuflxc(:,k_topo:nzt_rad+1), & 4768 rrtm_swuflx(:,k_topo:nzt_rad+1), & 4769 rrtm_swdflx(:,k_topo:nzt_rad+1), & 4770 rrtm_swhr(:,k_topo+1:nzt_rad+1), & 4771 rrtm_swuflxc(:,k_topo:nzt_rad+1), & 4769 4772 rrtm_swdflxc(:,k_topo:nzt_rad+1), & 4770 4773 rrtm_swhrc(:,k_topo+1:nzt_rad+1), & … … 4800 4803 surf_lsm_h%rad_sw_in(m) = rrtm_swdflx(0,k_topo) 4801 4804 surf_lsm_h%rad_sw_out(m) = rrtm_swuflx(0,k_topo) 4802 ENDDO 4805 ENDDO 4803 4806 DO m = surf_usm_h%start_index(j,i), & 4804 4807 surf_usm_h%end_index(j,i) 4805 4808 surf_usm_h%rad_sw_in(m) = rrtm_swdflx(0,k_topo) 4806 4809 surf_usm_h%rad_sw_out(m) = rrtm_swuflx(0,k_topo) 4807 ENDDO 4808 ! 4809 ! Vertical surfaces. Fluxes are obtain at respective vertical 4810 ENDDO 4811 ! 4812 ! Vertical surfaces. Fluxes are obtain at respective vertical 4810 4813 ! level of the surface element 4811 4814 DO l = 0, 3 … … 4815 4818 surf_lsm_v(l)%rad_sw_in(m) = rrtm_swdflx(0,k) 4816 4819 surf_lsm_v(l)%rad_sw_out(m) = rrtm_swuflx(0,k) 4817 ENDDO 4820 ENDDO 4818 4821 DO m = surf_usm_v(l)%start_index(j,i), & 4819 4822 surf_usm_v(l)%end_index(j,i) … … 4821 4824 surf_usm_v(l)%rad_sw_in(m) = rrtm_swdflx(0,k) 4822 4825 surf_usm_v(l)%rad_sw_out(m) = rrtm_swuflx(0,k) 4823 ENDDO 4826 ENDDO 4824 4827 ENDDO 4825 4828 ! … … 4829 4832 rad_sw_out = 0.0_wp 4830 4833 ! !!!!!!!! ATTENSION !!!!!!!!!!!!!!! 4831 ! Surface radiative fluxes should be also set to zero here 4834 ! Surface radiative fluxes should be also set to zero here 4832 4835 ! Save surface radiative fluxes onto respective surface elements 4833 4836 ! Horizontal surfaces … … 4836 4839 surf_lsm_h%rad_sw_in(m) = 0.0_wp 4837 4840 surf_lsm_h%rad_sw_out(m) = 0.0_wp 4838 ENDDO 4841 ENDDO 4839 4842 DO m = surf_usm_h%start_index(j,i), & 4840 4843 surf_usm_h%end_index(j,i) 4841 4844 surf_usm_h%rad_sw_in(m) = 0.0_wp 4842 4845 surf_usm_h%rad_sw_out(m) = 0.0_wp 4843 ENDDO 4844 ! 4845 ! Vertical surfaces. Fluxes are obtain at respective vertical 4846 ENDDO 4847 ! 4848 ! Vertical surfaces. Fluxes are obtain at respective vertical 4846 4849 ! level of the surface element 4847 4850 DO l = 0, 3 … … 4851 4854 surf_lsm_v(l)%rad_sw_in(m) = 0.0_wp 4852 4855 surf_lsm_v(l)%rad_sw_out(m) = 0.0_wp 4853 ENDDO 4856 ENDDO 4854 4857 DO m = surf_usm_v(l)%start_index(j,i), & 4855 4858 surf_usm_v(l)%end_index(j,i) … … 4857 4860 surf_usm_v(l)%rad_sw_in(m) = 0.0_wp 4858 4861 surf_usm_v(l)%rad_sw_out(m) = 0.0_wp 4859 ENDDO 4862 ENDDO 4860 4863 ENDDO 4861 4864 ENDIF … … 4868 4871 ! Finally, calculate surface net radiation for surface elements. 4869 4872 IF ( .NOT. radiation_interactions ) THEN 4870 ! First, for horizontal surfaces 4873 ! First, for horizontal surfaces 4871 4874 DO m = 1, surf_lsm_h%ns 4872 4875 surf_lsm_h%rad_net(m) = surf_lsm_h%rad_sw_in(m) & … … 4882 4885 ENDDO 4883 4886 ! 4884 ! Vertical surfaces. 4887 ! Vertical surfaces. 4885 4888 ! Todo: weight with azimuth and zenith angle according to their orientation! 4886 DO l = 0, 3 4889 DO l = 0, 3 4887 4890 DO m = 1, surf_lsm_v(l)%ns 4888 4891 surf_lsm_v(l)%rad_net(m) = surf_lsm_v(l)%rad_sw_in(m) & … … 4907 4910 4908 4911 CALL exchange_horiz( rad_sw_in, nbgp ) 4909 CALL exchange_horiz( rad_sw_out, nbgp ) 4912 CALL exchange_horiz( rad_sw_out, nbgp ) 4910 4913 CALL exchange_horiz( rad_sw_hr, nbgp ) 4911 4914 CALL exchange_horiz( rad_sw_cs_hr, nbgp ) … … 4921 4924 !> Calculate the cosine of the zenith angle (variable is called zenith) 4922 4925 !! 4923 SUBROUTINE calc_zenith( day_of_year, second_of_day ) 4926 SUBROUTINE calc_zenith( day_of_year, second_of_day ) 4924 4927 4925 4928 USE palm_date_time_mod, & … … 4973 4976 ! Description: 4974 4977 !  4975 !> Calculates surface albedo components based on Briegleb (1992) and 4978 !> Calculates surface albedo components based on Briegleb (1992) and 4976 4979 !> Briegleb et al. (1986) 4977 4980 !! … … 4990 4993 ! 4991 4994 ! Loop over surface elements 4992 DO ind_type = 0, SIZE( surf%albedo_type, 1)  14993 4995 DO ind_type = 0, SIZE( surf%albedo_type, 2 )  1 4996 4994 4997 ! 4995 4998 ! Ocean … … 5081 5084 surf%rrtm_asdif = albedo_urb 5082 5085 surf%rrtm_aldir = albedo_urb 5083 surf%rrtm_aldif = albedo_urb 5086 surf%rrtm_aldif = albedo_urb 5084 5087 ! 5085 5088 ! Darkness … … 5150 5153 DEALLOCATE ( rrtm_sw_tauaer ) 5151 5154 DEALLOCATE ( rrtm_sw_ssaaer ) 5152 DEALLOCATE ( rrtm_sw_asmaer ) 5153 DEALLOCATE ( rrtm_sw_ecaer ) 5154 5155 DEALLOCATE ( rrtm_sw_asmaer ) 5156 DEALLOCATE ( rrtm_sw_ecaer ) 5157 5155 5158 DEALLOCATE ( rrtm_swdflx ) 5156 5159 DEALLOCATE ( rrtm_swdflxc ) … … 5236 5239 5237 5240 ! 5238 ! Calculate pressure levels on zu and zw grid. Sounding data is added at 5239 ! top of the LES domain. This routine does not consider horizontal or 5241 ! Calculate pressure levels on zu and zw grid. Sounding data is added at 5242 ! top of the LES domain. This routine does not consider horizontal or 5240 5243 ! vertical variability of pressure and temperature 5241 5244 ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1) ) … … 5263 5266 5264 5267 ! 5265 ! Calculate temperature/humidity levels at top of the LES domain. 5266 ! Currently, the temperature is taken from sounding data (might lead to a 5267 ! temperature jump at interface. To do: Humidity is currently not 5268 ! Calculate temperature/humidity levels at top of the LES domain. 5269 ! Currently, the temperature is taken from sounding data (might lead to a 5270 ! temperature jump at interface. To do: Humidity is currently not 5268 5271 ! calculated above the LES domain. 5269 5272 ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1) ) … … 5299 5302 ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 5300 5303 ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 5301 ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 5302 ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) ) 5304 ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 5305 ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) ) 5303 5306 5304 5307 ! … … 5332 5335 rrtm_swdflx = 0.0_wp 5333 5336 rrtm_swuflx = 0.0_wp 5334 rrtm_swhr = 0.0_wp 5337 rrtm_swhr = 0.0_wp 5335 5338 rrtm_swuflxc = 0.0_wp 5336 5339 rrtm_swdflxc = 0.0_wp … … 5348 5351 rrtm_lwdflx = 0.0_wp 5349 5352 rrtm_lwuflx = 0.0_wp 5350 rrtm_lwhr = 0.0_wp 5353 rrtm_lwhr = 0.0_wp 5351 5354 rrtm_lwuflxc = 0.0_wp 5352 5355 rrtm_lwdflxc = 0.0_wp … … 5367 5370 !> Read trace gas data from file and convert into trace gas paths / volume 5368 5371 !> mixing ratios. If a userdefined input file is provided it needs to follow 5369 !> the convections used in RRTMG (see respective netCDF files shipped with 5372 !> the convections used in RRTMG (see respective netCDF files shipped with 5370 5373 !> RRTMG) 5371 5374 !! … … 5422 5425 DEALLOCATE ( rrtm_cfc22vmr ) 5423 5426 DEALLOCATE ( rrtm_ccl4vmr ) 5424 DEALLOCATE ( rrtm_h2ovmr ) 5427 DEALLOCATE ( rrtm_h2ovmr ) 5425 5428 ENDIF 5426 5429 … … 5446 5449 nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim ) 5447 5450 CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 ) 5448 nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np) 5451 nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np) 5449 5452 CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 ) 5450 5453 5451 5454 nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim ) 5452 5455 CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 ) 5453 nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs ) 5456 nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs ) 5454 5457 CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 ) 5455 5456 5457 ! 5458 ! Allocate pressure, and trace gas arrays 5458 5459 5460 ! 5461 ! Allocate pressure, and trace gas arrays 5459 5462 ALLOCATE( p_mls(1:np) ) 5460 ALLOCATE( trace_mls(1:num_trace_gases,1:np) ) 5461 ALLOCATE( trace_mls_tmp(1:nabs,1:np) ) 5463 ALLOCATE( trace_mls(1:num_trace_gases,1:np) ) 5464 ALLOCATE( trace_mls_tmp(1:nabs,1:np) ) 5462 5465 5463 5466 … … 5482 5485 ! 5483 5486 ! Replace missing values by zero 5484 WHERE ( trace_mls(n,:) > 2.0_wp ) 5487 WHERE ( trace_mls(n,:) > 2.0_wp ) 5485 5488 trace_mls(n,:) = 0.0_wp 5486 5489 END WHERE … … 5497 5500 ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) ) 5498 5501 5499 rrtm_play_tmp(1:nzt_rad) = rrtm_play(0,1:nzt_rad) 5502 rrtm_play_tmp(1:nzt_rad) = rrtm_play(0,1:nzt_rad) 5500 5503 rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1) 5501 5504 rrtm_play_tmp(nzt_rad+1) = rrtm_plev(0,nzt_rad+1) * 0.5_wp 5502 5505 rrtm_plev_tmp(nzt_rad+2) = MIN( 1.0E4_wp, 0.25_wp & 5503 5506 * rrtm_plev(0,nzt_rad+1) ) 5504 5507 5505 5508 ! 5506 5509 ! Calculate trace gas path (zero at surface) with interpolation to the … … 5509 5512 5510 5513 trace_mls_path(nzb+1,:) = 0.0_wp 5511 5514 5512 5515 DO k = nzb+2, nzt_rad+2 5513 5516 DO m = 1, num_trace_gases … … 5516 5519 ! 5517 5520 ! When the pressure level is higher than the trace gas pressure 5518 ! level, assume that 5519 IF ( rrtm_plev_tmp(k1) > p_mls(1) ) THEN 5520 5521 ! level, assume that 5522 IF ( rrtm_plev_tmp(k1) > p_mls(1) ) THEN 5523 5521 5524 trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1) & 5522 5525 * ( rrtm_plev_tmp(k1) & … … 5526 5529 5527 5530 ! 5528 ! Integrate for each sounding level from the contributing p_mls 5531 ! Integrate for each sounding level from the contributing p_mls 5529 5532 ! levels 5530 5533 DO n = 2, np … … 5553 5556 ENDDO 5554 5557 5555 IF ( rrtm_plev_tmp(k) < p_mls(np) ) THEN 5558 IF ( rrtm_plev_tmp(k) < p_mls(np) ) THEN 5556 5559 trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np) & 5557 5560 * ( MIN( rrtm_plev_tmp(k1), p_mls(np) ) & 5558 5561  rrtm_plev_tmp(k) & 5559 ) / g 5560 ENDIF 5562 ) / g 5563 ENDIF 5561 5564 ENDDO 5562 5565 ENDDO … … 5617 5620 5618 5621 rrtm_h2ovmr(0,:) = trace_path_tmp(:) 5619 5622 5620 5623 CASE DEFAULT 5621 5624 … … 5772 5775 INTEGER(iwp) :: pc_box_dimshift !< transform for best accuracy 5773 5776 INTEGER(iwp), DIMENSION(0:3) :: reorder = (/ 1, 0, 3, 2 /) 5774 5777 5775 5778 REAL(wp), DIMENSION(3,3) :: mrot !< grid rotation matrix (zyx) 5776 5779 REAL(wp), DIMENSION(3,0:nsurf_type):: vnorm !< face direction normal vectors (zyx) … … 5797 5800 REAL(wp) :: area_surf !< total area of surfaces in all processor 5798 5801 REAL(wp) :: area_hor !< total horizontal area of domain in all processor 5799 #if defined( __parallel ) 5802 #if defined( __parallel ) 5800 5803 REAL(wp), DIMENSION(1:7) :: combine_allreduce !< dummy array used to combine several MPI_ALLREDUCE calls 5801 5804 REAL(wp), DIMENSION(1:7) :: combine_allreduce_l !< dummy array used to combine several MPI_ALLREDUCE calls … … 5850 5853 ENDIF 5851 5854 ! 5852 ! Split downwelling shortwave radiation into a diffuse and a direct part. 5855 ! Split downwelling shortwave radiation into a diffuse and a direct part. 5853 5856 ! Note, if radiation scheme is RRTMG or diffuse radiation is externally 5854 5857 ! prescribed, this is not required. Please note, in case of external 5855 ! radiation, the clearsky model is applied during spinup, so that 5856 ! radiation need to be split also in this case. 5858 ! radiation, the clearsky model is applied during spinup, so that 5859 ! radiation need to be split also in this case. 5857 5860 IF ( radiation_scheme == 'constant' .OR. & 5858 5861 radiation_scheme == 'clearsky' .OR. & … … 6141 6144 IF ( .NOT. surface_reflections ) THEN 6142 6145 ! 6143 ! Set nrefsteps to 0 to disable reflections 6146 ! Set nrefsteps to 0 to disable reflections 6144 6147 nrefsteps = 0 6145 6148 surfoutsl = albedo_surf * surfins … … 6549 6552 combine_allreduce_l(6) = emiss_sum_surfl 6550 6553 combine_allreduce_l(7) = area_surfl 6551 6554 6552 6555 CALL MPI_ALLREDUCE( combine_allreduce_l, & 6553 6556 combine_allreduce, & … … 6557 6560 comm2d, & 6558 6561 ierr ) 6559 6562 6560 6563 pinsw = combine_allreduce(1) 6561 6564 pinlw = combine_allreduce(2) … … 6580 6583 IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf 6581 6584 ! 6582 ! Temporally comment out calculation of effective radiative temperature. 6585 ! Temporally comment out calculation of effective radiative temperature. 6583 6586 ! See below for more explanation. 6584 6587 ! (3) temperature 6585 6588 ! first we calculate an effective horizontal area to account for 6586 6589 ! the effect of vertical surfaces (which contributes to LW emission) 6587 ! We simply use the ratio of the total LW to the incoming LW flux 6590 ! We simply use the ratio of the total LW to the incoming LW flux 6588 6591 area_hor = pinlw / rad_lw_in_diff(nyn,nxl) 6589 6592 t_rad_urb = ( ( pemitlw  pabslw + emissivity_urb * pinlw ) / & 6590 (emissivity_urb * sigma_sb * area_hor) )**0.25_wp 6593 (emissivity_urb * sigma_sb * area_hor) )**0.25_wp 6591 6594 6592 6595 IF ( debug_output_timestep ) CALL debug_message( 'radiation_interaction', 'end' ) … … 6675 6678 !> It follows Boland, Ridley & Brown (2008) 6676 6679 !! 6677 SUBROUTINE calc_diffusion_radiation 6680 SUBROUTINE calc_diffusion_radiation 6678 6681 6679 6682 USE palm_date_time_mod, & … … 6706 6709 0.000719_wp * cos(2.0_wp * year_angle) + & 6707 6710 0.000077_wp * sin(2.0_wp * year_angle)) 6708 6709 ! 6711 6712 ! 6710 6713 ! Under a very low angle, we keep extraterestrial radiation at 6711 6714 ! the last small value, therefore the clearness index will be pushed 6712 ! towards 0 while keeping full continuity. 6715 ! towards 0 while keeping full continuity. 6713 6716 IF ( cos_zenith <= lowest_solarUp ) THEN 6714 6717 corrected_solarUp = lowest_solarUp … … 6716 6719 corrected_solarUp = cos_zenith 6717 6720 ENDIF 6718 6721 6719 6722 horizontalETR = etr * corrected_solarUp 6720 6723 6721 6724 DO i = nxl, nxr 6722 6725 DO j = nys, nyn … … 6728 6731 ENDDO 6729 6732 ENDDO 6730 6733 6731 6734 END SUBROUTINE calc_diffusion_radiation 6732 6735 … … 6826 6829 6827 6830 END SUBROUTINE radiation_interaction 6828 6831 6829 6832 !! 6830 6833 ! Description: … … 6934 6937 #endif 6935 6938 ! 6936 ! Stretching (nonuniform grid spacing) is not considered in the radiation 6937 ! model. Therefore, vertical stretching has to be applied above the area 6939 ! Stretching (nonuniform grid spacing) is not considered in the radiation 6940 ! model. Therefore, vertical stretching has to be applied above the area 6938 6941 ! where the parts of the radiation model which assume constant grid spacing 6939 ! are active. ABS (...) is required because the default value of 6940 ! dz_stretch_level_start is 9999999.9_wp (negative). 6942 ! are active. ABS (...) is required because the default value of 6943 ! dz_stretch_level_start is 9999999.9_wp (negative). 6941 6944 IF ( ABS( dz_stretch_level_start(1) ) <= zw(nz_urban_t) ) THEN 6942 6945 WRITE( message_string, * ) 'The lowest level where vertical ', & … … 6944 6947 'greater than ', zw(nz_urban_t) 6945 6948 CALL message( 'radiation_interaction_init', 'PA0496', 1, 2, 0, 6, 0 ) 6946 ENDIF 6949 ENDIF 6947 6950 ! 6948 6951 ! global number of urban and plant layers … … 7365 7368 !! 7366 7369 SUBROUTINE radiation_calc_svf 7367 7370 7368 7371 IMPLICIT NONE 7369 7372 7370 7373 INTEGER(iwp) :: i, j, k, d, ip, jp 7371 7374 INTEGER(iwp) :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrt, imrtf, ipcgb … … 7411 7414 INTEGER(kind=MPI_ADDRESS_KIND) :: size_lad_rma 7412 7415 #endif 7413 ! 7416 ! 7414 7417 INTEGER(iwp), DIMENSION(0:svfnorm_report_num) :: svfnorm_counts 7415 7418 … … 7440 7443 ray_skip_maxdist = 0 7441 7444 ray_skip_minval = 0 7442 7445 7443 7446 ! initialize temporary terrain and plant canopy height arrays (global 2D array!) 7444 7447 ALLOCATE( nzterr(0:(nx+1)*(ny+1)1) ) … … 7482 7485 IF ( raytrace_mpi_rma ) THEN 7483 7486 ALLOCATE( lad_s_ray(maxboxesg) ) 7484 7487 7485 7488 ! set conditions for RMA communication 7486 7489 CALL MPI_Info_create(minfo, ierr) … … 7556 7559 FLUSH(9) 7557 7560 ENDIF 7558 7561 7559 7562 ELSE 7560 7563 ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nz_plant1) ) … … 7589 7592 ! 7590 7593 dsitrans(:,:) = 0._wp 7591 7594 7592 7595 DO isurflt = 1, nsurfl 7593 7596 ! determine face centers … … 7694 7697 SUM(ztransp(itarg0:itarg0+lowest_free_ray1) & 7695 7698 * vffrac(itarg0:itarg0+lowest_free_ray1)) 7696 7699 7697 7700 ! Save direct solar transparency 7698 7701 j = MODULO(NINT(azmid/ & … … 7760 7763 CALL debug_message( debug_string, 'info' ) 7761 7764 ENDIF 7762 7765 7763 7766 nsvfla = k 7764 7767 ENDIF … … 7785 7788 CYCLE 7786 7789 ENDIF 7787 7790 7788 7791 sd = surf(id, isurfs) 7789 7792 sa = (/ REAL(surf(iz, isurfs), wp)  0.5_wp * kdir(sd), & … … 7801 7804 CYCLE 7802 7805 ENDIF 7803 7806 7804 7807 difvf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction 7805 7808 * dot_product((/ kdir(td), jdir(td), idir(td) /), uv) & ! cosine of target normal and reverse direction … … 7845 7848 CALL debug_message( debug_string, 'info' ) 7846 7849 ENDIF 7847 7850 7848 7851 nsvfla = k 7849 7852 ENDIF … … 8082 8085 ! deallocate temporary global arrays 8083 8086 DEALLOCATE(nzterr) 8084 8087 8085 8088 IF ( plant_canopy ) THEN 8086 8089 ! finalize mpi_rma communication and deallocate temporary arrays … … 8202 8205 ! sort and merge csf for the last time, keeping the array size to minimum 8203 8206 CALL merge_and_grow_csf(1) 8204 8207 8205 8208 ! aggregate csb among processors 8206 8209 ! allocate necessary arrays … … 8214 8217 ALLOCATE( ipcsflt(0:numprocs1) ) 8215 8218 ALLOCATE( dpcsflt(0:numprocs1) ) 8216 8217 ! fill out arrays of csf values and 8219 8220 ! fill out arrays of csf values and 8218 8221 ! arrays of number of elements and displacements 8219 8222 ! for particular precessors … … 8257 8260 dcsflt(jp) = d 8258 8261 ENDDO 8259 8262 8260 8263 ! deallocate temporary acsf array 8261 8264 ! DEALLOCATE(acsf)  ifort has a problem with deallocation of allocatable target … … 8267 8270 DEALLOCATE(acsf2) 8268 8271 ENDIF 8269 8272 8270 8273 #if defined( __parallel ) 8271 8274 ! scatter and gather the number of elements to and from all processor … … 8309 8312 FLUSH(9) 8310 8313 ENDIF 8311 8314 8312 8315 #else 8313 8316 npcsfl = ncsfl … … 8370 8373 ENDDO 8371 8374 ENDIF 8372 8375 8373 8376 ! deallocation of temporary arrays 8374 8377 IF ( npcbl > 0 ) DEALLOCATE( gridpcbl ) … … 8379 8382 CALL debug_message( debug_string, 'info' ) 8380 8383 ENDIF 8381 8384 8382 8385 ENDIF 8383 8386 … … 8388 8391 8389 8392 RETURN !todo: remove 8390 8393 8391 8394 ! WRITE( message_string, * ) & 8392 8395 ! 'I/O error when processing shape view factors / ', & 8393 8396 ! 'plant canopy sink factors / direct irradiance factors.' 8394 8397 ! CALL message( 'init_urban_surface', 'PA0502', 2, 2, 0, 6, 0 ) 8395 8398 8396 8399 END SUBROUTINE radiation_calc_svf 8397 8400 8398 8401 8399 8402 !! 8400 8403 ! Description: … … 8440 8443 INTEGER(iwp), DIMENSION(3) :: dimnext !< next dimension increments along path 8441 8444 INTEGER(iwp), DIMENSION(3) :: dimdelta !< dimension direction = + 1 8442 INTEGER(iwp) :: px, py !< number of processors in x and y dir before 8445 INTEGER(iwp) :: px, py !< number of processors in x and y dir before 8443 8446 !< the processor in the question 8444 8447 INTEGER(iwp) :: ip !< number of processor where gridbox reside 8445 8448 INTEGER(iwp) :: ig !< 1D index of gridbox in global 2D array 8446 8449 8447 8450 REAL(wp) :: eps = 1E10_wp !< epsilon for value comparison 8448 8451 REAL(wp) :: lad_s_target !< recieved lad_s of particular grid box … … 8461 8464 CALL merge_and_grow_csf(k) 8462 8465 ENDIF 8463 8466 8464 8467 transparency = 1._wp 8465 8468 ncsb = 0 … … 8534 8537 dimnextdist(seldim) = (dimnext(seldim)  .5_wp  src(seldim)) / uvect(seldim) 8535 8538 ENDDO 8536 8539 8537 8540 IF ( plant_canopy ) THEN 8538 8541 #if defined( __parallel ) … … 8549 8552 ENDIF 8550 8553 ENDDO 8551 8554 8552 8555 ! wait for all pending local requests complete 8553 8556 CALL MPI_Win_flush_local_all(win_lad, ierr) … … 8557 8560 ENDIF 8558 8561 CALL cpu_log( log_point_s(77), 'rad_rma_lad', 'stop' ) 8559 8562 8560 8563 ENDIF 8561 8564 #endif … … 8586 8589 8587 8590 transparency = transparency * (1._wp  cursink) 8588 8591 8589 8592 ENDDO 8590 8593 ENDIF 8591 8594 8592 8595 visible = .TRUE. 8593 8596 8594 8597 END SUBROUTINE raytrace 8595 8596 8598 8599 8597 8600 !! 8598 8601 ! Description: … … 8640 8643 INTEGER(iwp), DIMENSION(2) :: dimnext !< next dimension increments along path 8641 8644 INTEGER(iwp), DIMENSION(2) :: dimdelta !< dimension direction = + 1 8642 INTEGER(iwp) :: px, py !< number of processors in x and y dir before 8645 INTEGER(iwp) :: px, py !< number of processors in x and y dir before 8643 8646 !< the processor in the question 8644 8647 INTEGER(iwp) :: ip !< number of processor where gridbox reside 8645 8648 INTEGER(iwp) :: ig !< 1D index of gridbox in global 2D array 8646 8649 INTEGER(iwp) :: maxboxes !< max no of CSF created 8647 INTEGER(iwp) :: nly !< maximum plant canopy height 8650 INTEGER(iwp) :: nly !< maximum plant canopy height 8648 8651 INTEGER(iwp) :: ntrack 8649 8652 8650 8653 INTEGER(iwp) :: zb0 8651 8654 INTEGER(iwp) :: zb1 … … 8661 8664 INTEGER(MPI_ADDRESS_KIND) :: wdisp !< RMA window displacement 8662 8665 #endif 8663 8666 8664 8667 REAL(wp) :: eps = 1E10_wp !< epsilon for value comparison 8665 8668 REAL(wp) :: zbottom, ztop !< urban surface boundary in real numbers … … 8669 8672 REAL(wp) :: dxxyy !< square of real horizontal distance 8670 8673 REAL(wp) :: curtrans !< transparency of current PC box crossing 8671 8672 8673 8674 8675 8676 8674 8677 yxorigin(:) = origin(2:3) 8675 8678 transparency(:) = 1._wp ! Preset the all rays to transparent before reducing … … 8816 8819 IF ( plant_canopy ) THEN 8817 8820 ! Request LAD WHERE applicable 8818 ! 8821 ! 8819 8822 #if defined( __parallel ) 8820 8823 IF ( raytrace_mpi_rma ) THEN … … 8898 8901 ENDIF 8899 8902 ENDDO 8900 8903 8901 8904 DEALLOCATE( target_surfl ) 8902 8905 8903 8906 ELSE 8904 8907 itarget(:) = 1 … … 8907 8910 IF ( plant_canopy ) THEN 8908 8911 ! Skip the PCB around origin if requested (for MRT, the PCB might not be there) 8909 ! 8912 ! 8910 8913 IF ( skip_1st_pcb .AND. NINT(origin(1)) <= plantt_max ) THEN 8911 8914 rt2_track_lad(NINT(origin(1), iwp), 1) = 0._wp … … 8913 8916 8914 8917 ! Assert that we have space allocated for CSFs 8915 ! 8918 ! 8916 8919 maxboxes = (ntrack + MAX(CEILING(origin(1).5_wp)  nz_urban_b, & 8917 8920 nz_urban_t  CEILING(origin(1).5_wp))) * nrays … … 8926 8929 8927 8930 ! Calculate transparencies and store new CSFs 8928 ! 8931 ! 8929 8932 zbottom = REAL(nz_urban_b, wp)  .5_wp 8930 8933 ztop = REAL(plantt_max, wp) + .5_wp 8931 8934 8932 8935 ! Reverse direction of radiation (face>sky), only when calc_svf 8933 ! 8936 ! 8934 8937 IF ( calc_svf ) THEN 8935 8938 DO i = 1, ntrack ! for each column … … 9002 9005 9003 9006 ! Forward direction of radiation (sky>face), always 9004 ! 9007 ! 9005 9008 DO i = ntrack, 1, 1 ! for each column backwards 9006 9009 dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)rt2_track_dist(i1))**2 … … 9107 9110 9108 9111 END SUBROUTINE raytrace_2d 9109 9112 9110 9113 9111 9114 !! … … 9132 9135 ndsidir = 0 9133 9136 sun_direction = .TRUE. 9134 9137 9135 9138 ! 9136 9139 ! Process spinup time if configured … … 9176 9179 CALL calc_zenith( day_of_year, second_of_day ) 9177 9180 IF ( cos_zenith > 0 ) THEN 9178 ! 9181 ! 9179 9182 ! Identify solar direction vector (discretized number) 1) 9180 9183 i = MODULO(NINT(ATAN2(sun_dir_lon, sun_dir_lat) & … … 9210 9213 IMPLICIT NONE 9211 9214 INTEGER(iwp), INTENT(in) :: x, y, z, d, x2, y2, z2, d2 9212 9215 9213 9216 surface_facing = .FALSE. 9214 9217 … … 9255 9258 9256 9259 surface_facing = .TRUE. 9257 9260 9258 9261 END FUNCTION surface_facing 9259 9262 … … 9270 9273 9271 9274 IMPLICIT NONE 9272 9275 9273 9276 CHARACTER(rad_version_len) :: rad_version_field 9274 9277 9275 9278 INTEGER(iwp) :: i 9276 9279 INTEGER(iwp) :: ndsidir_from_file = 0 … … 9296 9299 'that has written the svf data ',& 9297 9300 'and the one that will read it ',& 9298 'is not allowed' 9301 'is not allowed' 9299 9302 CALL message( 'check_open', 'PA0491', 1, 2, 0, 6, 0 ) 9300 9303 ENDIF 9301 9304 9302 9305 ENDIF 9303 9304 ! 9305 ! Open binary file 9306 9307 ! 9308 ! Open binary file 9306 9309 CALL check_open( 88 ) 9307 9310 … … 9315 9318 CALL message( 'radiation_read_svf', 'PA0482', 1, 2, 0, 6, 0 ) 9316 9319 ENDIF 9317 9320 9318 9321 ! 9319 9322 ! read nsvfl, ncsfl, nsurfl, nmrtf 9320 9323 READ ( 88 ) nsvfl, ncsfl, nsurfl_from_file, npcbl_from_file, & 9321 9324 ndsidir_from_file, nmrtbl_from_file, nmrtf 9322 9325 9323 9326 IF ( nsvfl < 0 .OR. ncsfl < 0 ) THEN 9324 9327 WRITE( message_string, * ) 'Wrong number of SVF or CSF' … … 9330 9333 IF ( debug_output ) CALL debug_message( debug_string, 'info' ) 9331 9334 ENDIF 9332 9335 9333 9336 IF ( nsurfl_from_file /= nsurfl ) THEN 9334 9337 WRITE( message_string, * ) 'nsurfl from SVF file does not ', & … … 9337 9340 CALL message( 'radiation_read_svf', 'PA0490', 1, 2, 0, 6, 0 ) 9338 9341 ENDIF 9339 9342 9340 9343 IF ( npcbl_from_file /= npcbl ) THEN 9341 9344 WRITE( message_string, * ) 'npcbl from SVF file does not ', & … … 9344 9347 CALL message( 'radiation_read_svf', 'PA0493', 1, 2, 0, 6, 0 ) 9345 9348 ENDIF 9346 9349 9347 9350 IF ( ndsidir_from_file /= ndsidir ) THEN 9348 9351 WRITE( message_string, * ) 'ndsidir from SVF file does not ', & … … 9360 9363 IF ( debug_output ) CALL debug_message( debug_string, 'info' ) 9361 9364 ENDIF 9362 9363 ! 9364 ! Arrays skyvf, skyvft, dsitrans and dsitransc are allready 9365 ! allocated in radiation_interaction_init and 9365 9366 ! 9367 ! Arrays skyvf, skyvft, dsitrans and dsitransc are allready 9368 ! allocated in radiation_interaction_init and 9366 9369 ! radiation_presimulate_solar_pos 9367 9370 IF ( nsurfl > 0 ) THEN 9368 9371 READ(88) skyvf 9369 9372 READ(88) skyvft 9370 READ(88) dsitrans 9373 READ(88) dsitrans 9371 9374 ENDIF 9372 9375 9373 9376 IF ( plant_canopy .AND. npcbl > 0 ) THEN 9374 9377 READ ( 88 ) dsitransc 9375 9378 ENDIF 9376 9377 ! 9378 ! The allocation of svf, svfsurf, csf, csfsurf, mrtf, mrtft, and 9379 ! mrtfsurf happens in routine radiation_calc_svf which is not 9380 ! called if the program enters radiation_read_svf. Therefore 9379 9380 ! 9381 ! The allocation of svf, svfsurf, csf, csfsurf, mrtf, mrtft, and 9382 ! mrtfsurf happens in routine radiation_calc_svf which is not 9383 ! called if the program enters radiation_read_svf. Therefore 9381 9384 ! these arrays has to allocate in the following 9382 9385 IF ( nsvfl > 0 ) THEN … … 9408 9411 READ(88) mrtfsurf 9409 9412 ENDIF 9410 9411 ! 9412 ! Close binary file 9413 9414 ! 9415 ! Close binary file 9413 9416 CALL close_file( 88 ) 9414 9417 9415 9418 ENDIF 9416 9419 #if defined( __parallel ) … … 9436 9439 9437 9440 IMPLICIT NONE 9438 9441 9439 9442 INTEGER(iwp) :: i 9440 9443 … … 9445 9448 IF ( i == io_group ) THEN 9446 9449 ! 9447 ! Open binary file 9450 ! Open binary file 9448 9451 CALL check_open( 89 ) 9449 9452 … … 9473 9476 IF ( nmrtf > 0 ) THEN 9474 9477 WRITE ( 89 ) mrtf 9475 WRITE ( 89 ) mrtft 9478 WRITE ( 89 ) mrtft 9476 9479 WRITE ( 89 ) mrtfsurf 9477 9480 ENDIF 9478 9481 ! 9479 ! Close binary file 9482 ! Close binary file 9480 9483 CALL close_file( 89 ) 9481 9484 … … 9627 9630 END SUBROUTINE quicksort_csf 9628 9631 9629 9632 9630 9633 !! 9631 9634 ! … … 9712 9715 END SUBROUTINE merge_and_grow_csf 9713 9716 9714 9717 9715 9718 ! quicksort.f *f90* 9716 9719 ! Author: tnissie, adaptation J.Resler … … 9746 9749 IF ( j+1 < last ) CALL quicksort_csf2(kpcsflt, pcsflt, j+1, last) 9747 9750 END SUBROUTINE quicksort_csf2 9748 9751 9749 9752 9750 9753 PURE FUNCTION csf_lt2(item1, item2) result(res) … … 9781 9784 !! 9782 9785 SUBROUTINE radiation_3d_data_averaging( mode, variable ) 9783 9786 9784 9787 9785 9788 USE control_parameters … … 9791 9794 IMPLICIT NONE 9792 9795 9793 CHARACTER (LEN=*) :: mode !< 9794 CHARACTER (LEN=*) :: variable !< 9796 CHARACTER (LEN=*) :: mode !< 9797 CHARACTER (LEN=*) :: variable !< 9795 9798 9796 9799 LOGICAL :: match_lsm !< flag indicating naturaltype surface 9797 9800 LOGICAL :: match_usm !< flag indicating urbantype surface 9798 9799 INTEGER(iwp) :: i !< 9800 INTEGER(iwp) :: j !< 9801 INTEGER(iwp) :: k !< 9801 9802 INTEGER(iwp) :: i !< 9803 INTEGER(iwp) :: j !< 9804 INTEGER(iwp) :: k !< 9802 9805 INTEGER(iwp) :: l, m !< index of current surface element 9803 9806 … … 9835 9838 ENDIF 9836 9839 rad_net_av = 0.0_wp 9837 9840 9838 9841 CASE ( 'rad_lw_in*' ) 9839 9842 IF ( .NOT. ALLOCATED( rad_lw_in_xy_av ) ) THEN … … 9841 9844 ENDIF 9842 9845 rad_lw_in_xy_av = 0.0_wp 9843 9846 9844 9847 CASE ( 'rad_lw_out*' ) 9845 9848 IF ( .NOT. ALLOCATED( rad_lw_out_xy_av ) ) THEN … … 9847 9850 ENDIF 9848 9851 rad_lw_out_xy_av = 0.0_wp 9849 9852 9850 9853 CASE ( 'rad_sw_in*' ) 9851 9854 IF ( .NOT. ALLOCATED( rad_sw_in_xy_av ) ) THEN … … 9853 9856 ENDIF 9854 9857 rad_sw_in_xy_av = 0.0_wp 9855 9858 9856 9859 CASE ( 'rad_sw_out*' ) 9857 9860 IF ( .NOT. ALLOCATED( rad_sw_out_xy_av ) ) THEN 9858 9861 ALLOCATE( rad_sw_out_xy_av(nysg:nyng,nxlg:nxrg) ) 9859 9862 ENDIF 9860 rad_sw_out_xy_av = 0.0_wp 9863 rad_sw_out_xy_av = 0.0_wp 9861 9864 9862 9865 CASE ( 'rad_lw_in' ) … … 10097 10100 ENDDO 10098 10101 ENDIF 10099 10102 10100 10103 CASE ( 'rad_lw_out*' ) 10101 10104 IF ( ALLOCATED( rad_lw_out_xy_av ) ) THEN … … 10119 10122 ENDDO 10120 10123 ENDIF 10121 10124 10122 10125 CASE ( 'rad_sw_in*' ) 10123 10126 IF ( ALLOCATED( rad_sw_in_xy_av ) ) THEN … … 10141 10144 ENDDO 10142 10145 ENDIF 10143 10146 10144 10147 CASE ( 'rad_sw_out*' ) 10145 10148 IF ( ALLOCATED( rad_sw_out_xy_av ) ) THEN … … 10163 10166 ENDDO 10164 10167 ENDIF 10165 10168 10166 10169 CASE ( 'rad_lw_in' ) 10167 10170 IF ( ALLOCATED( rad_lw_in_av ) ) THEN … … 10418 10421 ENDDO 10419 10422 ENDIF 10420 10423 10421 10424 CASE ( 'rad_lw_in*' ) 10422 10425 IF ( ALLOCATED( rad_lw_in_xy_av ) ) THEN … … 10428 10431 ENDDO 10429 10432 ENDIF 10430 10433 10431 10434 CASE ( 'rad_lw_out*' ) 10432 10435 IF ( ALLOCATED( rad_lw_out_xy_av ) ) THEN … … 10438 10441 ENDDO 10439 10442 ENDIF 10440 10443 10441 10444 CASE ( 'rad_sw_in*' ) 10442 10445 IF ( ALLOCATED( rad_sw_in_xy_av ) ) THEN … … 10448 10451 ENDDO 10449 10452 ENDIF 10450 10453 10451 10454 CASE ( 'rad_sw_out*' ) 10452 10455 IF ( ALLOCATED( rad_sw_out_xy_av ) ) THEN … … 10702 10705 !! 10703 10706 SUBROUTINE radiation_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z ) 10704 10707 10705 10708 IMPLICIT NONE 10706 10709 10707 10710 CHARACTER (LEN=*), INTENT(IN) :: variable !< 10708 LOGICAL, INTENT(OUT) :: found !< 10709 CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< 10710 CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< 10711 CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< 10711 LOGICAL, INTENT(OUT) :: found !< 10712 CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< 10713 CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< 10714 CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< 10712 10715 10713 10716 CHARACTER (len=varnamelength) :: var … … 10779 10782 SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode, & 10780 10783 local_pf, two_d, nzb_do, nzt_do ) 10781 10784 10782 10785 USE indices 10783 10786 … … 10787 10790 IMPLICIT NONE 10788 10791 10789 CHARACTER (LEN=*) :: grid !< 10790 CHARACTER (LEN=*) :: mode !< 10791 CHARACTER (LEN=*) :: variable !< 10792 10793 INTEGER(iwp) :: av !< 10794 INTEGER(iwp) :: i !< 10795 INTEGER(iwp) :: j !< 10796 INTEGER(iwp) :: k !< 10792 CHARACTER (LEN=*) :: grid !< 10793 CHARACTER (LEN=*) :: mode !< 10794 CHARACTER (LEN=*) :: variable !< 10795 10796 INTEGER(iwp) :: av !< 10797 INTEGER(iwp) :: i !< 10798 INTEGER(iwp) :: j !< 10799 INTEGER(iwp) :: k !< 10797 10800 INTEGER(iwp) :: m !< index of surface element at grid point (j,i) 10798 INTEGER(iwp) :: nzb_do !< 10799 INTEGER(iwp) :: nzt_do !< 10800 10801 LOGICAL :: found !< 10801 INTEGER(iwp) :: nzb_do !< 10802 INTEGER(iwp) :: nzt_do !< 10803 10804 LOGICAL :: found !< 10802 10805 LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections) 10803 10806 10804 10807 REAL(wp) :: fill_value = 999.0_wp !< value for the _FillValue attribute 10805 10808 10806 REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< 10809 REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< 10807 10810 10808 10811 found = .TRUE. … … 10818 10821 ! Naturaltype surfaces 10819 10822 DO m = surf_lsm_h%start_index(j,i), & 10820 surf_lsm_h%end_index(j,i) 10823 surf_lsm_h%end_index(j,i) 10821 10824 local_pf(i,j,nzb+1) = surf_lsm_h%rad_net(m) 10822 10825 ENDDO … … 10824 10827 ! Urbantype surfaces 10825 10828 DO m = surf_usm_h%start_index(j,i), & 10826 surf_usm_h%end_index(j,i) 10829 surf_usm_h%end_index(j,i) 10827 10830 local_pf(i,j,nzb+1) = surf_usm_h%rad_net(m) 10828 10831 ENDDO … … 10835 10838 ENDIF 10836 10839 DO i = nxl, nxr 10837 DO j = nys, nyn 10840 DO j = nys, nyn 10838 10841 local_pf(i,j,nzb+1) = rad_net_av(j,i) 10839 10842 ENDDO … … 10842 10845 two_d = .TRUE. 10843 10846 grid = 'zu1' 10844 10847 10845 10848 CASE ( 'rad_lw_in*_xy' ) ! 2darray 10846 10849 IF ( av == 0 ) THEN … … 10851 10854 ! Naturaltype surfaces 10852 10855 DO m = surf_lsm_h%start_index(j,i), & 10853 surf_lsm_h%end_index(j,i) 10856 surf_lsm_h%end_index(j,i) 10854 10857 local_pf(i,j,nzb+1) = surf_lsm_h%rad_lw_in(m) 10855 10858 ENDDO … … 10857 10860 ! Urbantype surfaces 10858 10861 DO m = surf_usm_h%start_index(j,i), & 10859 surf_usm_h%end_index(j,i) 10862 surf_usm_h%end_index(j,i) 10860 10863 local_pf(i,j,nzb+1) = surf_usm_h%rad_lw_in(m) 10861 10864 ENDDO … … 10868 10871 ENDIF 10869 10872 DO i = nxl, nxr 10870 DO j = nys, nyn 10873 DO j = nys, nyn 10871 10874 local_pf(i,j,nzb+1) = rad_lw_in_xy_av(j,i) 10872 10875 ENDDO … … 10875 10878 two_d = .TRUE. 10876 10879 grid = 'zu1' 10877 10880 10878 10881 CASE ( 'rad_lw_out*_xy' ) ! 2darray 10879 10882 IF ( av == 0 ) THEN … … 10884 10887 ! Naturaltype surfaces 10885 10888 DO m = surf_lsm_h%start_index(j,i), & 10886 surf_lsm_h%end_index(j,i) 10889 surf_lsm_h%end_index(j,i) 10887 10890 local_pf(i,j,nzb+1) = surf_lsm_h%rad_lw_out(m) 10888 10891 ENDDO … … 10890 10893 ! Urbantype surfaces 10891 10894 DO m = surf_usm_h%start_index(j,i), & 10892 surf_usm_h%end_index(j,i) 10895 surf_usm_h%end_index(j,i) 10893 10896 local_pf(i,j,nzb+1) = surf_usm_h%rad_lw_out(m) 10894 10897 ENDDO … … 10901 10904 ENDIF 10902 10905 DO i = nxl, nxr 10903 DO j = nys, nyn 10906 DO j = nys, nyn 10904 10907 local_pf(i,j,nzb+1) = rad_lw_out_xy_av(j,i) 10905 10908 ENDDO … … 10908 10911 two_d = .TRUE. 10909 10912 grid = 'zu1' 10910 10913 10911 10914 CASE ( 'rad_sw_in*_xy' ) ! 2darray 10912 10915 IF ( av == 0 ) THEN … … 10917 10920 ! Naturaltype surfaces 10918 10921 DO m = surf_lsm_h%start_index(j,i), & 10919 surf_lsm_h%end_index(j,i) 10922 surf_lsm_h%end_index(j,i) 10920 10923 local_pf(i,j,nzb+1) = surf_lsm_h%rad_sw_in(m) 10921 10924 ENDDO … … 10923 10926 ! Urbantype surfaces 10924 10927 DO m = surf_usm_h%start_index(j,i), & 10925 surf_usm_h%end_index(j,i) 10928 surf_usm_h%end_index(j,i) 10926 10929 local_pf(i,j,nzb+1) = surf_usm_h%rad_sw_in(m) 10927 10930 ENDDO … … 10934 10937 ENDIF 10935 10938 DO i = nxl, nxr 10936 DO j = nys, nyn 10939 DO j = nys, nyn 10937 10940 local_pf(i,j,nzb+1) = rad_sw_in_xy_av(j,i) 10938 10941 ENDDO … … 10941 10944 two_d = .TRUE. 10942 10945 grid = 'zu1' 10943 10946 10944 10947 CASE ( 'rad_sw_out*_xy' ) ! 2darray 10945 10948 IF ( av == 0 ) THEN … … 10950 10953 ! Naturaltype surfaces 10951 10954 DO m = surf_lsm_h%start_index(j,i), & 10952 surf_lsm_h%end_index(j,i) 10955 surf_lsm_h%end_index(j,i) 10953 10956 local_pf(i,j,nzb+1) = surf_lsm_h%rad_sw_out(m) 10954 10957 ENDDO … … 10956 10959 ! Urbantype surfaces 10957 10960 DO m = surf_usm_h%start_index(j,i), & 10958 surf_usm_h%end_index(j,i) 10961 surf_usm_h%end_index(j,i) 10959 10962 local_pf(i,j,nzb+1) = surf_usm_h%rad_sw_out(m) 10960 10963 ENDDO … … 10967 10970 ENDIF 10968 10971 DO i = nxl, nxr 10969 DO j = nys, nyn 10972 DO j = nys, nyn 10970 10973 local_pf(i,j,nzb+1) = rad_sw_out_xy_av(j,i) 10971 10974 ENDDO … … 10973 10976 ENDIF 10974 10977 two_d = .TRUE. 10975 grid = 'zu1' 10976 10978 grid = 'zu1' 10979 10977 10980 CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' ) 10978 10981 IF ( av == 0 ) THEN … … 10990 10993 ENDIF 10991 10994 DO i = nxl, nxr 10992 DO j = nys, nyn 10995 DO j = nys, nyn 10993 10996 DO k = nzb_do, nzt_do 10994 10997 local_pf(i,j,k) = rad_lw_in_av(k,j,i) … … 11014 11017 ENDIF 11015 11018 DO i = nxl, nxr 11016 DO j = nys, nyn 11019 DO j = nys, nyn 11017 11020 DO k = nzb_do, nzt_do 11018 11021 local_pf(i,j,k) = rad_lw_out_av(k,j,i) … … 11020 11023 ENDDO 11021 11024 ENDDO 11022 ENDIF 11025 ENDIF 11023 11026 IF ( mode == 'xy' ) grid = 'zu' 11024 11027 … … 11038 11041 ENDIF 11039 11042 DO i = nxl, nxr 11040 DO j = nys, nyn 11043 DO j = nys, nyn 11041 11044 DO k = nzb_do, nzt_do 11042 11045 local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i) … … 11062 11065 ENDIF 11063 11066 DO i = nxl, nxr 11064 DO j = nys, nyn 11067 DO j = nys, nyn 11065 11068 DO k = nzb_do, nzt_do 11066 11069 local_pf(i,j,k) = rad_lw_hr_av(k,j,i) … … 11086 11089 ENDIF 11087 11090 DO i = nxl, nxr 11088 DO j = nys, nyn 11091 DO j = nys, nyn 11089 11092 DO k = nzb_do, nzt_do 11090 11093 local_pf(i,j,k) = rad_sw_in_av(k,j,i) … … 11110 11113 ENDIF 11111 11114 DO i = nxl, nxr 11112 DO j = nys, nyn 11115 DO j = nys, nyn 11113 11116 DO k = nzb, nzt+1 11114 11117 local_pf(i,j,k) = rad_sw_out_av(k,j,i) … … 11134 11137 ENDIF 11135 11138 DO i = nxl, nxr 11136 DO j = nys, nyn 11139 DO j = nys, nyn 11137 11140 DO k = nzb_do, nzt_do 11138 11141 local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i) … … 11158 11161 ENDIF 11159 11162 DO i = nxl, nxr 11160 DO j = nys, nyn 11163 DO j = nys, nyn 11161 11164 DO k = nzb_do, nzt_do 11162 11165 local_pf(i,j,k) = rad_sw_hr_av(k,j,i) … … 11172 11175 11173 11176 END SELECT 11174 11177 11175 11178 END SUBROUTINE radiation_data_output_2d 11176 11179 … … 11183 11186 !! 11184 11187 SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) 11185 11188 11186 11189 11187 11190 USE indices … … 11192 11195 IMPLICIT NONE 11193 11196 11194 CHARACTER (LEN=*) :: variable !< 11197 CHARACTER (LEN=*) :: variable !< 11195 11198 11196 11199 INTEGER(iwp) :: av !< … … 11217 11220 RETURN 11218 11221 ENDIF 11219 11222 11220 11223 IF ( var(1:3) /= 'rad' .AND. var(1:3) /= 'rtm' ) THEN 11221 11224 found = .FALSE. … … 11684 11687 ENDIF 11685 11688 ENDIF 11686 ! 11689 ! 11687 11690 ! block of RTM output variables 11688 11691 ! variables are intended mainly for debugging and detailed analyse purposes 11689 11692 CASE ( 'rtm_skyvf' ) 11690 ! 11693 ! 11691 11694 ! sky view factor 11692 11695 DO isurf = dirstart(ids), dirend(ids) … … 11758 11761 !! 11759 11762 SUBROUTINE radiation_data_output_mask( av, variable, found, local_pf, mid ) 11760 11763 11761 11764 USE control_parameters 11762 11765 11763 11766 USE indices 11764 11767 11765 11768 USE kinds 11766 11769 11767 11770 11768 11771 IMPLICIT NONE 11769 11772 11770 CHARACTER (LEN=*) :: variable !< 11773 CHARACTER (LEN=*) :: variable !< 11771 11774 11772 11775 CHARACTER(LEN=5) :: grid !< flag to distinquish between staggered grids 11773 11776 11774 INTEGER(iwp) :: av !< 11775 INTEGER(iwp) :: i !< 11776 INTEGER(iwp) :: j !< 11777 INTEGER(iwp) :: k !< 11777 INTEGER(iwp) :: av !< 11778 INTEGER(iwp) :: i !< 11779 INTEGER(iwp) :: j !< 11780 INTEGER(iwp) :: k !< 11778 11781 INTEGER(iwp) :: mid !< masked output running index 11779 11782 INTEGER(iwp) :: topo_top_index !< k index of highest horizontal surface … … 11785 11788 REAL(wp), & 11786 11789 DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: & 11787 local_pf !< 11790 local_pf !< 11788 11791 11789 11792 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output … … 11918 11921 WRITE ( 14 ) rad_net_av 11919 11922 ENDIF 11920 11923 11921 11924 IF ( ALLOCATED( rad_lw_in_xy_av ) ) THEN 11922 11925 CALL wrd_write_string( 'rad_lw_in_xy_av' ) 11923 11926 WRITE ( 14 ) rad_lw_in_xy_av 11924 11927 ENDIF 11925 11928 11926 11929 IF ( ALLOCATED( rad_lw_out_xy_av ) ) THEN 11927 11930 CALL wrd_write_string( 'rad_lw_out_xy_av' ) 11928 11931 WRITE ( 14 ) rad_lw_out_xy_av 11929 11932 ENDIF 11930 11933 11931 11934 IF ( ALLOCATED( rad_sw_in_xy_av ) ) THEN 11932 11935 CALL wrd_write_string( 'rad_sw_in_xy_av' ) 11933 11936 WRITE ( 14 ) rad_sw_in_xy_av 11934 11937 ENDIF 11935 11938 11936 11939 IF ( ALLOCATED( rad_sw_out_xy_av ) ) THEN 11937 11940 CALL wrd_write_string( 'rad_sw_out_xy_av' ) … … 12030 12033 nxr_on_file, nynf, nync, nyn_on_file, nysf, & 12031 12034 nysc, nys_on_file, tmp_2d, tmp_3d, found ) 12032 12035 12033 12036 12034 12037 USE control_parameters 12035 12038 12036 12039 USE indices 12037 12040 12038 12041 USE kinds 12039 12042 12040 12043 USE pegrid 12041 12044 … … 12043 12046 IMPLICIT NONE 12044 12047 12045 INTEGER(iwp) :: k !< 12046 INTEGER(iwp) :: nxlc !< 12047 INTEGER(iwp) :: nxlf !< 12048 INTEGER(iwp) :: nxl_on_file !< 12049 INTEGER(iwp) :: nxrc !< 12050 INTEGER(iwp) :: nxrf !< 12051 INTEGER(iwp) :: nxr_on_file !< 12052 INTEGER(iwp) :: nync !< 12053 INTEGER(iwp) :: nynf !< 12054 INTEGER(iwp) :: nyn_on_file !< 12055 INTEGER(iwp) :: nysc !< 12056 INTEGER(iwp) :: nysf !< 12057 INTEGER(iwp) :: nys_on_file !< 12048 INTEGER(iwp) :: k !< 12049 INTEGER(iwp) :: nxlc !< 12050 INTEGER(iwp) :: nxlf !< 12051 INTEGER(iwp) :: nxl_on_file !< 12052 INTEGER(iwp) :: nxrc !< 12053 INTEGER(iwp) :: nxrf !< 12054 INTEGER(iwp) :: nxr_on_file !< 12055 INTEGER(iwp) :: nync !< 12056 INTEGER(iwp) :: nynf !< 12057 INTEGER(iwp) :: nyn_on_file !< 12058 INTEGER(iwp) :: nysc !< 12059 INTEGER(iwp) :: nysf !< 12060 INTEGER(iwp) :: nys_on_file !< 12058 12061 12059 12062 LOGICAL, INTENT(OUT) :: found 12060 12063 12061 REAL(wp), DIMENSION(nys_on_filenbgp:nyn_on_file+nbgp,nxl_on_filenbgp:nxr_on_file+nbgp) :: tmp_2d !< 12062 12063 REAL(wp), DIMENSION(nzb:nzt+1,nys_on_filenbgp:nyn_on_file+nbgp,nxl_on_filenbgp:nxr_on_file+nbgp) :: tmp_3d !< 12064 12065 REAL(wp), DIMENSION(0:0,nys_on_filenbgp:nyn_on_file+nbgp,nxl_on_filenbgp:nxr_on_file+nbgp) :: tmp_3d2 !< 12064 REAL(wp), DIMENSION(nys_on_filenbgp:nyn_on_file+nbgp,nxl_on_filenbgp:nxr_on_file+nbgp) :: tmp_2d !< 12065 12066 REAL(wp), DIMENSION(nzb:nzt+1,nys_on_filenbgp:nyn_on_file+nbgp,nxl_on_filenbgp:nxr_on_file+nbgp) :: tmp_3d !< 12067 12068 REAL(wp), DIMENSION(0:0,nys_on_filenbgp:nyn_on_file+nbgp,nxl_on_filenbgp:nxr_on_file+nbgp) :: tmp_3d2 !< 12066 12069 12067 12070 … … 12074 12077 IF ( .NOT. ALLOCATED( rad_net_av ) ) THEN 12075 12078 ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) ) 12076 ENDIF 12079 ENDIF 12077 12080 IF ( k == 1 ) READ ( 13 ) tmp_2d 12078 12081 rad_net_av(nyscnbgp:nync+nbgp,nxlcnbgp:nxrc+nbgp) = & 12079 12082 tmp_2d(nysfnbgp:nynf+nbgp,nxlfnbgp:nxrf+nbgp) 12080 12083 12081 12084 CASE ( 'rad_lw_in_xy_av' ) 12082 12085 IF ( .NOT. ALLOCATED( rad_lw_in_xy_av ) ) THEN 12083 12086 ALLOCATE( rad_lw_in_xy_av(nysg:nyng,nxlg:nxrg) ) 12084 ENDIF 12087 ENDIF 12085 12088 IF ( k == 1 ) READ ( 13 ) tmp_2d 12086 12089 rad_lw_in_xy_av(nyscnbgp:nync+nbgp,nxlcnbgp:nxrc+nbgp) = & 12087 12090 tmp_2d(nysfnbgp:nynf+nbgp,nxlfnbgp:nxrf+nbgp) 12088 12091 12089 12092 CASE ( 'rad_lw_out_xy_av' ) 12090 12093 IF ( .NOT. ALLOCATED( rad_lw_out_xy_av ) ) THEN 12091 12094 ALLOCATE( rad_lw_out_xy_av(nysg:nyng,nxlg:nxrg) ) 12092 ENDIF 12095 ENDIF 12093 12096 IF ( k == 1 ) READ ( 13 ) tmp_2d 12094 12097 rad_lw_out_xy_av(nyscnbgp:nync+nbgp,nxlcnbgp:nxrc+nbgp) = & 12095 12098 tmp_2d(nysfnbgp:nynf+nbgp,nxlfnbgp:nxrf+nbgp) 12096 12099 12097 12100 CASE ( 'rad_sw_in_xy_av' ) 12098 12101 IF ( .NOT. ALLOCATED( rad_sw_in_xy_av ) ) THEN 12099 12102 ALLOCATE( rad_sw_in_xy_av(nysg:nyng,nxlg:nxrg) ) 12100 ENDIF 12103 ENDIF 12101 12104 IF ( k == 1 ) READ ( 13 ) tmp_2d 12102 12105 rad_sw_in_xy_av(nyscnbgp:nync+nbgp,nxlcnbgp:nxrc+nbgp) = & 12103 12106 tmp_2d(nysfnbgp:nynf+nbgp,nxlfnbgp:nxrf+nbgp) 12104 12107 12105 12108 CASE ( 'rad_sw_out_xy_av' ) 12106 12109 IF ( .NOT. ALLOCATED( rad_sw_out_xy_av ) ) THEN 12107 12110 ALLOCATE( rad_sw_out_xy_av(nysg:nyng,nxlg:nxrg) ) 12108 ENDIF 12111 ENDIF 12109 12112 IF ( k == 1 ) READ ( 13 ) tmp_2d 12110 12113 rad_sw_out_xy_av(nyscnbgp:nync+nbgp,nxlcnbgp:nxrc+nbgp) = & 12111 12114 tmp_2d(nysfnbgp:nynf+nbgp,nxlfnbgp:nxrf+nbgp) 12112 12115 12113 12116 CASE ( 'rad_lw_in' ) 12114 12117 IF ( .NOT. ALLOCATED( rad_lw_in ) ) THEN … … 12120 12123 ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 12121 12124 ENDIF 12122 ENDIF 12125 ENDIF 12123 12126 IF ( k == 1 ) THEN 12124 12127 IF ( radiation_scheme == 'clearsky' .OR. & … … 12144 12147 ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 12145 12148 ENDIF 12146 ENDIF 12149 ENDIF 12147 12150 IF ( k == 1 ) THEN 12148 12151 IF ( radiation_scheme == 'clearsky' .OR. & … … 12168 12171 ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 12169 12172 ENDIF 12170 ENDIF 12173 ENDIF 12171 12174 IF ( k == 1 ) THEN 12172 12175 IF ( radiation_scheme == 'clearsky' .OR. & … … 12192 12195 ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 12193 12196 ENDIF 12194 ENDIF 12197 ENDIF 12195 12198 IF ( k == 1 ) THEN 12196 12199 IF ( radiation_scheme == 'clearsky' .OR. & … … 12248 12251 ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 12249 12252 ENDIF 12250 ENDIF 12253 ENDIF 12251 12254 IF ( k == 1 ) THEN 12252 12255 IF ( radiation_scheme == 'clearsky' .OR. & … … 12272 12275 ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 12273 12276 ENDIF 12274 ENDIF 12277 ENDIF 12275 12278 IF ( k == 1 ) THEN 12276 12279 IF ( radiation_scheme == 'clearsky' .OR. & … … 12296 12299 ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 12297 12300 ENDIF 12298 ENDIF 12301 ENDIF 12299 12302 IF ( k == 1 ) THEN 12300 12303 IF ( radiation_scheme == 'clearsky' .OR. & … … 12320 12323 ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 12321 12324 ENDIF 12322 ENDIF 12325 ENDIF 12323 12326 IF ( k == 1 ) THEN 12324 12327 IF ( radiation_scheme == 'clearsky' .OR. &
Note: See TracChangeset
for help on using the changeset viewer.