- Timestamp:
- Aug 21, 2019 11:13:06 AM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
netcdf_data_input_mod.f90 (modified) (5 diffs)
-
radiation_model_mod.f90 (modified) (34 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/netcdf_data_input_mod.f90
r4150 r4178 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Implement input of external radiation forcing. Therefore, provide public 28 ! subroutines and variables. 29 ! 30 ! 4150 2019-08-08 20:00:47Z suehring 27 31 ! Some variables are given the public attribute, in order to call netcdf input 28 32 ! from single routines … … 608 612 LOGICAL :: from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used 609 613 END TYPE int_2d_32bit 610 614 ! 615 !-- Define data type to read 1D real variables 616 TYPE real_1d 617 LOGICAL :: from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used 618 619 REAL(wp) :: fill = -9999.9_wp !< fill value 620 621 REAL(wp), DIMENSION(:), ALLOCATABLE :: var !< respective variable 622 END TYPE real_1d 611 623 ! 612 624 !-- Define data type to read 2D real variables … … 834 846 CHARACTER(LEN=100) :: input_file_vm = 'PIDS_VM' !< Name of file which comprises virtual measurement data 835 847 836 CHARACTER(LEN=25), ALLOCATABLE, DIMENSION(:) :: string_values !< output of string variables read from netcdf input files 848 CHARACTER(LEN=25), ALLOCATABLE, DIMENSION(:) :: string_values !< output of string variables read from netcdf input files 849 CHARACTER(LEN=50), DIMENSION(:), ALLOCATABLE :: vars_pids !< variable in input file 837 850 838 851 INTEGER(iwp) :: id_emis !< NetCDF id of input file for chemistry emissions: TBD: It has to be removed 839 852 840 853 INTEGER(iwp) :: nc_stat !< return value of nf90 function call 854 INTEGER(iwp) :: num_var_pids !< number of variables in file 855 INTEGER(iwp) :: pids_id !< file id 841 856 842 857 LOGICAL :: input_pids_static = .FALSE. !< Flag indicating whether Palm-input-data-standard file containing static information exists … … 949 964 ! 950 965 !-- Public data structures 951 PUBLIC real_2d 966 PUBLIC real_1d, & 967 real_2d 952 968 ! 953 969 !-- Public variables … … 956 972 chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type, & 957 973 coord_ref_sys, & 958 init_3d, init_model, input_file_atts, input_file_static, & 974 init_3d, init_model, input_file_atts, & 975 input_file_dynamic, & 976 input_file_static, & 959 977 input_pids_static, & 960 978 input_pids_dynamic, input_pids_vm, input_file_vm, & 961 979 leaf_area_density_f, nest_offl, & 980 num_var_pids, & 962 981 pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f, & 982 pids_id, & 963 983 root_area_density_lad_f, root_area_density_lsm_f, soil_pars_f, & 964 984 soil_type_f, street_crossing_f, street_type_f, surface_fraction_f, & 965 985 terrain_height_f, vegetation_pars_f, vegetation_type_f, & 986 vars_pids, & 966 987 water_pars_f, water_type_f 967 988 ! -
palm/trunk/SOURCE/radiation_model_mod.f90
r4168 r4178 28 28 ! ----------------- 29 29 ! $Id$ 30 ! External radiation forcing implemented. 31 ! 32 ! 4168 2019-08-16 13:50:17Z suehring 30 33 ! Replace function get_topography_top_index by topo_top_ind 31 34 ! … … 667 670 ONLY: cloud_droplets, coupling_char, & 668 671 debug_output, debug_output_timestep, debug_string, & 672 dt_3d, & 669 673 dz, dt_spinup, end_time, & 670 674 humidity, & … … 704 708 705 709 USE netcdf_data_input_mod, & 706 ONLY: albedo_type_f, albedo_pars_f, building_type_f, pavement_type_f, & 707 vegetation_type_f, water_type_f 710 ONLY: albedo_type_f, & 711 albedo_pars_f, & 712 building_type_f, & 713 pavement_type_f, & 714 vegetation_type_f, & 715 water_type_f, & 716 char_fill, & 717 check_existence, & 718 close_input_file, & 719 get_attribute, & 720 get_variable, & 721 inquire_num_variables, & 722 inquire_variable_names, & 723 input_file_dynamic, & 724 netcdf_data_input_get_dimension_length, & 725 num_var_pids, & 726 pids_id, & 727 open_read_file, & 728 real_1d, & 729 vars_pids 708 730 709 731 USE plant_canopy_model_mod, & … … 1236 1258 REAL(wp), DIMENSION(:), ALLOCATABLE :: albedo_surf !< albedo of the surface 1237 1259 REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf !< emissivity of the wall surface 1260 ! 1261 !-- External radiation 1262 TYPE( real_1d ) :: rad_lw_in_f !< external incoming longwave radiation, from observation or model 1263 TYPE( real_1d ) :: rad_sw_in_f !< external incoming shortwave radiation, from observation or model 1264 TYPE( real_1d ) :: rad_sw_in_dif_f !< external incoming shortwave radiation, diffuse part, from observation or model 1265 TYPE( real_1d ) :: time_rad_f !< time dimension for external radiation, from observation or model 1238 1266 1239 1267 … … 1414 1442 CASE ( 'rrtmg' ) 1415 1443 CALL radiation_rrtmg 1444 1445 CASE ( 'external' ) 1446 ! 1447 !-- During spinup apply clear-sky model 1448 IF ( time_since_reference_point < 0.0_wp ) THEN 1449 CALL radiation_clearsky 1450 ELSE 1451 CALL radiation_external 1452 ENDIF 1416 1453 1417 1454 CASE DEFAULT … … 1804 1841 IF ( radiation_scheme /= 'constant' .AND. & 1805 1842 radiation_scheme /= 'clear-sky' .AND. & 1806 radiation_scheme /= 'rrtmg' ) THEN 1843 radiation_scheme /= 'rrtmg' .AND. & 1844 radiation_scheme /= 'external' ) THEN 1807 1845 message_string = 'unknown radiation_scheme = '// & 1808 1846 TRIM( radiation_scheme ) … … 1899 1937 INTEGER(iwp) :: l !< running index for orientation of vertical surfaces 1900 1938 INTEGER(iwp) :: m !< running index for surface elements 1939 INTEGER(iwp) :: ntime = 0 !< number of available external radiation timesteps 1901 1940 #if defined( __rrtmg ) 1902 1941 INTEGER(iwp) :: ind_type !< running index for subgrid-surface tiles … … 2113 2152 2114 2153 IF ( radiation_scheme == 'clear-sky' .OR. & 2115 radiation_scheme == 'constant') THEN 2116 2117 2154 radiation_scheme == 'constant' .OR. & 2155 radiation_scheme == 'external' ) THEN 2118 2156 ! 2119 2157 !-- Allocate arrays for incoming/outgoing short/longwave radiation … … 2896 2934 #endif 2897 2935 ENDIF 2898 2936 ! 2937 !-- Initializaion actions exclusively required for external 2938 !-- radiation forcing 2939 IF ( radiation_scheme == 'external' ) THEN 2940 ! 2941 !-- Open the radiation input file 2942 #if defined( __netcdf ) 2943 CALL open_read_file( TRIM( input_file_dynamic ) // & 2944 TRIM( coupling_char ), & 2945 pids_id ) 2946 2947 CALL inquire_num_variables( pids_id, num_var_pids ) 2948 ! 2949 !-- Allocate memory to store variable names and read them 2950 ALLOCATE( vars_pids(1:num_var_pids) ) 2951 CALL inquire_variable_names( pids_id, vars_pids ) 2952 2953 ! 2954 !-- Input time dimension. 2955 IF ( check_existence( vars_pids, 'time_rad' ) ) THEN 2956 CALL netcdf_data_input_get_dimension_length( pids_id, & 2957 ntime, 'time_rad' ) 2958 2959 ALLOCATE( time_rad_f%var(0:ntime-1) ) 2960 ! 2961 !-- Read variable 2962 CALL get_variable( pids_id, 'time_rad', time_rad_f%var ) 2963 2964 time_rad_f%from_file = .TRUE. 2965 ENDIF 2966 ! 2967 !-- Input shortwave downwelling. 2968 IF ( check_existence( vars_pids, 'rad_sw_in' ) ) THEN 2969 ALLOCATE( rad_sw_in_f%var(0:ntime-1) ) 2970 ! 2971 !-- Read _FillValue attribute 2972 CALL get_attribute( pids_id, char_fill, rad_sw_in_f%fill, & 2973 .FALSE., 'rad_sw_in' ) 2974 ! 2975 !-- Read variable 2976 CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var ) 2977 2978 rad_sw_in_f%from_file = .TRUE. 2979 ENDIF 2980 ! 2981 !-- Input longwave downwelling. 2982 IF ( check_existence( vars_pids, 'rad_lw_in' ) ) THEN 2983 ALLOCATE( rad_lw_in_f%var(0:ntime-1) ) 2984 ! 2985 !-- Read _FillValue attribute 2986 CALL get_attribute( pids_id, char_fill, rad_lw_in_f%fill, & 2987 .FALSE., 'rad_lw_in' ) 2988 ! 2989 !-- Read variable 2990 CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var ) 2991 2992 rad_lw_in_f%from_file = .TRUE. 2993 ENDIF 2994 ! 2995 !-- Input shortwave downwelling, diffuse part. 2996 IF ( check_existence( vars_pids, 'rad_sw_in_dif' ) ) THEN 2997 ALLOCATE( rad_sw_in_dif_f%var(0:ntime-1) ) 2998 ! 2999 !-- Read _FillValue attribute 3000 CALL get_attribute( pids_id, char_fill, rad_sw_in_dif_f%fill, & 3001 .FALSE., 'rad_sw_in_dif' ) 3002 ! 3003 !-- Read variable 3004 CALL get_variable( pids_id, 'rad_sw_in_dif', rad_sw_in_dif_f%var ) 3005 3006 rad_sw_in_dif_f%from_file = .TRUE. 3007 ENDIF 3008 ! 3009 !-- Finally, close the input file. 3010 CALL close_input_file( pids_id ) 3011 #endif 3012 ! 3013 !-- Make some consistency checks. 3014 IF ( .NOT. rad_sw_in_f%from_file .OR. & 3015 .NOT. rad_lw_in_f%from_file ) THEN 3016 message_string = 'In case of external radiation forcing ' // & 3017 'both, rad_sw_in and rad_lw_in are required.' 3018 CALL message( 'radiation_init', 'PA0195', 1, 2, 0, 6, 0 ) 3019 ENDIF 3020 3021 IF ( .NOT. time_rad_f%from_file ) THEN 3022 message_string = 'In case of external radiation forcing ' // & 3023 'dimension time_rad is required.' 3024 CALL message( 'radiation_init', 'PA0196', 1, 2, 0, 6, 0 ) 3025 ENDIF 3026 3027 IF ( ANY( rad_sw_in_f%var == rad_sw_in_f%fill ) ) THEN 3028 message_string = 'External radiation forcing: rad_sw_in must ' // & 3029 'not contain any fill values' 3030 CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 ) 3031 ENDIF 3032 3033 IF ( ANY( rad_lw_in_f%var == rad_lw_in_f%fill ) ) THEN 3034 message_string = 'External radiation forcing: rad_lw_in must ' // & 3035 'not contain any fill values' 3036 CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 ) 3037 ENDIF 3038 3039 IF ( rad_sw_in_dif_f%from_file ) THEN 3040 IF ( ANY( rad_sw_in_dif_f%var == rad_sw_in_dif_f%fill ) ) THEN 3041 message_string = 'External radiation forcing: ' // & 3042 'rad_lw_in_dif must not contain any fill ' // & 3043 'values' 3044 CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 ) 3045 ENDIF 3046 ENDIF 3047 3048 IF ( time_rad_f%var(0) /= 0.0_wp ) THEN 3049 message_string = 'External radiation forcing: first point in ' // & 3050 'time is /= 0.0.' 3051 CALL message( 'radiation_init', 'PA0313', 1, 2, 0, 6, 0 ) 3052 ENDIF 3053 3054 IF ( end_time - spinup_time > time_rad_f%var(ntime-1) ) THEN 3055 message_string = 'External radiation forcing does not cover ' // & 3056 'the entire simulation time.' 3057 CALL message( 'radiation_init', 'PA0314', 1, 2, 0, 6, 0 ) 3058 ENDIF 3059 3060 ENDIF 2899 3061 ! 2900 3062 !-- Perform user actions if required … … 2913 3075 CASE ( 'constant' ) 2914 3076 CALL radiation_constant 3077 3078 CASE ( 'external' ) 3079 ! 3080 !-- During spinup apply clear-sky model 3081 IF ( time_since_reference_point < 0.0_wp ) THEN 3082 CALL radiation_clearsky 3083 ELSE 3084 CALL radiation_external 3085 ENDIF 2915 3086 2916 3087 CASE DEFAULT 2917 3088 2918 3089 END SELECT 2919 2920 ! readjust date and time to its initial value3090 ! 3091 !-- Readjust date and time to its initial value 2921 3092 CALL init_date_and_time 2922 3093 … … 2958 3129 2959 3130 3131 !------------------------------------------------------------------------------! 3132 ! Description: 3133 ! ------------ 3134 !> A simple clear sky radiation model 3135 !------------------------------------------------------------------------------! 3136 SUBROUTINE radiation_external 3137 3138 IMPLICIT NONE 3139 3140 INTEGER(iwp) :: l !< running index for surface orientation 3141 INTEGER(iwp) :: t !< index of current timestep 3142 INTEGER(iwp) :: tm !< index of previous timestep 3143 3144 REAL(wp) :: fac_dt !< interpolation factor 3145 REAL(wp) :: lw_in !< downwelling longwave radiation, interpolated value 3146 REAL(wp) :: sw_in !< downwelling shortwave radiation, interpolated value 3147 REAL(wp) :: sw_in_dif !< downwelling diffuse shortwave radiation, interpolated value 3148 3149 TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine 3150 3151 ! 3152 !-- Calculate current zenith angle 3153 CALL calc_zenith 3154 ! 3155 !-- Interpolate external radiation on current timestep 3156 IF ( time_since_reference_point <= 0.0_wp ) THEN 3157 t = 0 3158 tm = 0 3159 fac_dt = 0 3160 ELSE 3161 t = 0 3162 DO WHILE ( time_rad_f%var(t) <= time_since_reference_point ) 3163 t = t + 1 3164 ENDDO 3165 3166 tm = MAX( t-1, 0 ) 3167 3168 fac_dt = ( time_since_reference_point - time_rad_f%var(tm) + dt_3d ) & 3169 / ( time_rad_f%var(t) - time_rad_f%var(tm) ) 3170 fac_dt = MIN( 1.0_wp, fac_dt ) 3171 ENDIF 3172 3173 sw_in = ( 1.0_wp - fac_dt ) * rad_sw_in_f%var(tm) & 3174 + fac_dt * rad_sw_in_f%var(t) 3175 3176 lw_in = ( 1.0_wp - fac_dt ) * rad_lw_in_f%var(tm) & 3177 + fac_dt * rad_lw_in_f%var(t) 3178 ! 3179 !-- Limit shortwave incoming radiation to positive values, in order 3180 !-- to overcome possible observation errors. 3181 sw_in = MAX( 0.0_wp, sw_in ) 3182 sw_in = MERGE( sw_in, 0.0_wp, sun_up ) 3183 ! 3184 !-- Call clear-sky calculation for each surface orientation. 3185 !-- First, horizontal surfaces 3186 surf => surf_lsm_h 3187 CALL radiation_external_surf 3188 surf => surf_usm_h 3189 CALL radiation_external_surf 3190 ! 3191 !-- Vertical surfaces 3192 DO l = 0, 3 3193 surf => surf_lsm_v(l) 3194 CALL radiation_external_surf 3195 surf => surf_usm_v(l) 3196 CALL radiation_external_surf 3197 ENDDO 3198 ! 3199 !-- In case of average radiation check if external diffuse shortwave 3200 !-- radiation is available and, if required, interpolate it onto the 3201 !-- respective arrays. 3202 IF ( average_radiation ) THEN 3203 IF ( rad_sw_in_dif_f%from_file ) THEN 3204 sw_in_dif= ( 1.0_wp - fac_dt ) * rad_sw_in_dif_f%var(tm) & 3205 + fac_dt * rad_sw_in_dif_f%var(t) 3206 3207 IF ( ALLOCATED( rad_sw_in_diff ) ) rad_sw_in_diff = sw_in_dif 3208 IF ( ALLOCATED( rad_sw_in_dir ) ) rad_sw_in_dir = sw_in - sw_in_dif 3209 ! 3210 !-- Diffuse longwave radiation equals the total downwelling longwaver 3211 !-- radiation 3212 IF ( ALLOCATED( rad_lw_in_diff ) ) rad_lw_in_diff = lw_in 3213 ENDIF 3214 ENDIF 3215 3216 CONTAINS 3217 3218 SUBROUTINE radiation_external_surf 3219 3220 USE control_parameters 3221 3222 IMPLICIT NONE 3223 3224 INTEGER(iwp) :: i !< grid index along x-dimension 3225 INTEGER(iwp) :: j !< grid index along y-dimension 3226 INTEGER(iwp) :: k !< grid index along z-dimension 3227 INTEGER(iwp) :: m !< running index for surface elements 3228 3229 IF ( surf%ns < 1 ) RETURN 3230 3231 surf%rad_sw_in = sw_in 3232 surf%rad_lw_in = lw_in 3233 3234 IF ( average_radiation ) THEN 3235 3236 surf%rad_sw_out = albedo_urb * surf%rad_sw_in 3237 3238 surf%rad_lw_out = emissivity_urb * sigma_sb * t_rad_urb**4 & 3239 + ( 1.0_wp - emissivity_urb ) * surf%rad_lw_in 3240 3241 surf%rad_net = surf%rad_sw_in - surf%rad_sw_out & 3242 + surf%rad_lw_in - surf%rad_lw_out 3243 3244 surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb * sigma_sb & 3245 * t_rad_urb**3 3246 ! 3247 !-- Calculate upwelling radiation fluxes and net radiation for each 3248 !-- surface element depending on its properties. 3249 ELSE 3250 3251 DO m = 1, surf%ns 3252 i = surf%i(m) 3253 j = surf%j(m) 3254 k = surf%k(m) 3255 ! 3256 !-- Weighted average according to surface fraction. 3257 surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m) * & 3258 surf%albedo(ind_veg_wall,m) & 3259 + surf%frac(ind_pav_green,m) * & 3260 surf%albedo(ind_pav_green,m) & 3261 + surf%frac(ind_wat_win,m) * & 3262 surf%albedo(ind_wat_win,m) ) & 3263 * surf%rad_sw_in(m) 3264 3265 surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m) * & 3266 surf%emissivity(ind_veg_wall,m) & 3267 + surf%frac(ind_pav_green,m) * & 3268 surf%emissivity(ind_pav_green,m) & 3269 + surf%frac(ind_wat_win,m) * & 3270 surf%emissivity(ind_wat_win,m) & 3271 ) & 3272 * sigma_sb & 3273 * ( surf%pt_surface(m) * exner(k) )**4 3274 3275 surf%rad_lw_out_change_0(m) = & 3276 ( surf%frac(ind_veg_wall,m) * & 3277 surf%emissivity(ind_veg_wall,m) & 3278 + surf%frac(ind_pav_green,m) * & 3279 surf%emissivity(ind_pav_green,m) & 3280 + surf%frac(ind_wat_win,m) * & 3281 surf%emissivity(ind_wat_win,m) & 3282 ) * 4.0_wp * sigma_sb & 3283 * ( surf%pt_surface(m) * exner(k) )**3 3284 3285 surf%rad_net(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m) & 3286 + surf%rad_lw_in(m) - surf%rad_lw_out(m) 3287 3288 ENDDO 3289 3290 ENDIF 3291 3292 ! 3293 !-- Fill out values in radiation arrays 3294 DO m = 1, surf%ns 3295 i = surf%i(m) 3296 j = surf%j(m) 3297 rad_sw_out(0,j,i) = surf%rad_sw_out(m) 3298 rad_lw_out(0,j,i) = surf%rad_lw_out(m) 3299 ENDDO 3300 3301 rad_sw_in(0,:,:) = sw_in 3302 rad_lw_in(0,:,:) = lw_in 3303 3304 END SUBROUTINE radiation_external_surf 3305 3306 END SUBROUTINE radiation_external 3307 2960 3308 !------------------------------------------------------------------------------! 2961 3309 ! Description: … … 3061 3409 3062 3410 surf%rad_lw_in = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4 3411 3412 print*, "clear sky", emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4 3063 3413 3064 3414 surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4 & … … 3136 3486 i = surf%i(m) 3137 3487 j = surf%j(m) 3138 rad_sw_in(0,j,i) = surf%rad_sw_in(m)3488 rad_sw_in(0,j,i) = surf%rad_sw_in(m) 3139 3489 rad_sw_out(0,j,i) = surf%rad_sw_out(m) 3140 rad_lw_in(0,j,i) = surf%rad_lw_in(m)3490 rad_lw_in(0,j,i) = surf%rad_lw_in(m) 3141 3491 rad_lw_out(0,j,i) = surf%rad_lw_out(m) 3142 3492 ENDDO … … 3355 3705 IF ( .NOT. lw_radiation ) WRITE( io, 10 ) 3356 3706 IF ( .NOT. sw_radiation ) WRITE( io, 11 ) 3707 ELSEIF ( radiation_scheme == "external" ) THEN 3708 WRITE( io, 14 ) 3357 3709 ENDIF 3358 3710 … … 3390 3742 13 FORMAT (/' Albedo is set individually for each xy-location, according ', & 3391 3743 'to given surface type.') 3744 14 FORMAT (' --> External radiation forcing is used') 3392 3745 3393 3746 … … 5328 5681 ENDIF 5329 5682 ENDIF 5330 5331 !-- if radiation scheme is not RRTMG, split diffusion and direct part of the solar downward radiation 5332 !-- comming from radiation model and store it in 2D arrays 5333 IF ( radiation_scheme /= 'rrtmg' ) CALL calc_diffusion_radiation 5683 ! 5684 !-- Split downwelling shortwave radiation into a diffuse and a direct part. 5685 !-- Note, if radiation scheme is RRTMG or diffuse radiation is externally 5686 !-- prescribed, this is not required. 5687 IF ( radiation_scheme /= 'rrtmg' .AND. & 5688 .NOT. rad_sw_in_dif_f%from_file ) CALL calc_diffusion_radiation 5334 5689 5335 5690 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 6118 6473 SUBROUTINE calc_diffusion_radiation 6119 6474 6120 REAL(wp), PARAMETER :: lowest_solarUp = 0.1_wp !< limit the sun elevation to protect stability of the calculation 6121 INTEGER(iwp) :: i, j 6122 REAL(wp) :: year_angle !< angle 6123 REAL(wp) :: etr !< extraterestrial radiation 6124 REAL(wp) :: corrected_solarUp !< corrected solar up radiation 6125 REAL(wp) :: horizontalETR !< horizontal extraterestrial radiation 6126 REAL(wp) :: clearnessIndex !< clearness index 6127 REAL(wp) :: diff_frac !< diffusion fraction of the radiation 6128 6475 INTEGER(iwp) :: i !< grid index x-direction 6476 INTEGER(iwp) :: j !< grid index y-direction 6129 6477 6478 REAL(wp) :: year_angle !< angle 6479 REAL(wp) :: etr !< extraterestrial radiation 6480 REAL(wp) :: corrected_solarUp !< corrected solar up radiation 6481 REAL(wp) :: horizontalETR !< horizontal extraterestrial radiation 6482 REAL(wp) :: clearnessIndex !< clearness index 6483 REAL(wp) :: diff_frac !< diffusion fraction of the radiation 6484 6485 REAL(wp), PARAMETER :: lowest_solarUp = 0.1_wp !< limit the sun elevation to protect stability of the calculation 6486 print*, "in diffusion_radiation" 6487 ! 6130 6488 !-- Calculate current day and time based on the initial values and simulation time 6131 6489 year_angle = ( (day_of_year_init * 86400) + time_utc_init & … … 6138 6496 0.000719_wp * cos(2.0_wp * year_angle) + & 6139 6497 0.000077_wp * sin(2.0_wp * year_angle)) 6140 6498 ! 6141 6499 !-- 6142 6500 !-- Under a very low angle, we keep extraterestrial radiation at 6143 6501 !-- the last small value, therefore the clearness index will be pushed 6144 !-- towards 0 while keeping full continuity. 6145 !-- 6502 !-- towards 0 while keeping full continuity. 6146 6503 IF ( cos_zenith <= lowest_solarUp ) THEN 6147 6504 corrected_solarUp = lowest_solarUp … … 11451 11808 IF ( .NOT. ALLOCATED( rad_lw_in ) ) THEN 11452 11809 IF ( radiation_scheme == 'clear-sky' .OR. & 11453 radiation_scheme == 'constant') THEN 11810 radiation_scheme == 'constant' .OR. & 11811 radiation_scheme == 'external' ) THEN 11454 11812 ALLOCATE( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) ) 11455 11813 ELSE … … 11459 11817 IF ( k == 1 ) THEN 11460 11818 IF ( radiation_scheme == 'clear-sky' .OR. & 11461 radiation_scheme == 'constant') THEN 11819 radiation_scheme == 'constant' .OR. & 11820 radiation_scheme == 'external' ) THEN 11462 11821 READ ( 13 ) tmp_3d2 11463 11822 rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & … … 11473 11832 IF ( .NOT. ALLOCATED( rad_lw_in_av ) ) THEN 11474 11833 IF ( radiation_scheme == 'clear-sky' .OR. & 11475 radiation_scheme == 'constant') THEN 11834 radiation_scheme == 'constant' .OR. & 11835 radiation_scheme == 'external' ) THEN 11476 11836 ALLOCATE( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) ) 11477 11837 ELSE … … 11481 11841 IF ( k == 1 ) THEN 11482 11842 IF ( radiation_scheme == 'clear-sky' .OR. & 11483 radiation_scheme == 'constant') THEN 11843 radiation_scheme == 'constant' .OR. & 11844 radiation_scheme == 'external' ) THEN 11484 11845 READ ( 13 ) tmp_3d2 11485 11846 rad_lw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =& … … 11495 11856 IF ( .NOT. ALLOCATED( rad_lw_out ) ) THEN 11496 11857 IF ( radiation_scheme == 'clear-sky' .OR. & 11497 radiation_scheme == 'constant') THEN 11858 radiation_scheme == 'constant' .OR. & 11859 radiation_scheme == 'external' ) THEN 11498 11860 ALLOCATE( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) ) 11499 11861 ELSE … … 11503 11865 IF ( k == 1 ) THEN 11504 11866 IF ( radiation_scheme == 'clear-sky' .OR. & 11505 radiation_scheme == 'constant') THEN 11867 radiation_scheme == 'constant' .OR. & 11868 radiation_scheme == 'external' ) THEN 11506 11869 READ ( 13 ) tmp_3d2 11507 11870 rad_lw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & … … 11517 11880 IF ( .NOT. ALLOCATED( rad_lw_out_av ) ) THEN 11518 11881 IF ( radiation_scheme == 'clear-sky' .OR. & 11519 radiation_scheme == 'constant') THEN 11882 radiation_scheme == 'constant' .OR. & 11883 radiation_scheme == 'external' ) THEN 11520 11884 ALLOCATE( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) ) 11521 11885 ELSE … … 11525 11889 IF ( k == 1 ) THEN 11526 11890 IF ( radiation_scheme == 'clear-sky' .OR. & 11527 radiation_scheme == 'constant') THEN 11891 radiation_scheme == 'constant' .OR. & 11892 radiation_scheme == 'external' ) THEN 11528 11893 READ ( 13 ) tmp_3d2 11529 11894 rad_lw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) & … … 11571 11936 IF ( .NOT. ALLOCATED( rad_sw_in ) ) THEN 11572 11937 IF ( radiation_scheme == 'clear-sky' .OR. & 11573 radiation_scheme == 'constant') THEN 11938 radiation_scheme == 'constant' .OR. & 11939 radiation_scheme == 'external' ) THEN 11574 11940 ALLOCATE( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) ) 11575 11941 ELSE … … 11579 11945 IF ( k == 1 ) THEN 11580 11946 IF ( radiation_scheme == 'clear-sky' .OR. & 11581 radiation_scheme == 'constant') THEN 11947 radiation_scheme == 'constant' .OR. & 11948 radiation_scheme == 'external' ) THEN 11582 11949 READ ( 13 ) tmp_3d2 11583 11950 rad_sw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & … … 11593 11960 IF ( .NOT. ALLOCATED( rad_sw_in_av ) ) THEN 11594 11961 IF ( radiation_scheme == 'clear-sky' .OR. & 11595 radiation_scheme == 'constant') THEN 11962 radiation_scheme == 'constant' .OR. & 11963 radiation_scheme == 'external' ) THEN 11596 11964 ALLOCATE( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) ) 11597 11965 ELSE … … 11601 11969 IF ( k == 1 ) THEN 11602 11970 IF ( radiation_scheme == 'clear-sky' .OR. & 11603 radiation_scheme == 'constant') THEN 11971 radiation_scheme == 'constant' .OR. & 11972 radiation_scheme == 'external' ) THEN 11604 11973 READ ( 13 ) tmp_3d2 11605 11974 rad_sw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =& … … 11615 11984 IF ( .NOT. ALLOCATED( rad_sw_out ) ) THEN 11616 11985 IF ( radiation_scheme == 'clear-sky' .OR. & 11617 radiation_scheme == 'constant') THEN 11986 radiation_scheme == 'constant' .OR. & 11987 radiation_scheme == 'external' ) THEN 11618 11988 ALLOCATE( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) ) 11619 11989 ELSE … … 11623 11993 IF ( k == 1 ) THEN 11624 11994 IF ( radiation_scheme == 'clear-sky' .OR. & 11625 radiation_scheme == 'constant') THEN 11995 radiation_scheme == 'constant' .OR. & 11996 radiation_scheme == 'external' ) THEN 11626 11997 READ ( 13 ) tmp_3d2 11627 11998 rad_sw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & … … 11637 12008 IF ( .NOT. ALLOCATED( rad_sw_out_av ) ) THEN 11638 12009 IF ( radiation_scheme == 'clear-sky' .OR. & 11639 radiation_scheme == 'constant') THEN 12010 radiation_scheme == 'constant' .OR. & 12011 radiation_scheme == 'external' ) THEN 11640 12012 ALLOCATE( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) ) 11641 12013 ELSE … … 11645 12017 IF ( k == 1 ) THEN 11646 12018 IF ( radiation_scheme == 'clear-sky' .OR. & 11647 radiation_scheme == 'constant') THEN 12019 radiation_scheme == 'constant' .OR. & 12020 radiation_scheme == 'external' ) THEN 11648 12021 READ ( 13 ) tmp_3d2 11649 12022 rad_sw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) &
Note: See TracChangeset
for help on using the changeset viewer.
