Changeset 4178 for palm/trunk/SOURCE


Ignore:
Timestamp:
Aug 21, 2019 11:13:06 AM (5 years ago)
Author:
suehring
Message:

Enable external time-dependent radiative forcing with downwelling short- and longwave radiation. Optionally, also downwelling diffuse radiation can be provided. Radiation data will be provided via dynamic input file.

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/netcdf_data_input_mod.f90

    r4150 r4178  
    2525! -----------------
    2626! $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
    2731! Some variables are given the public attribute, in order to call netcdf input
    2832! from single routines
     
    608612       LOGICAL ::  from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used
    609613    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   
    611623!
    612624!-- Define data type to read 2D real variables
     
    834846    CHARACTER(LEN=100) ::  input_file_vm      = 'PIDS_VM'      !< Name of file which comprises virtual measurement data
    835847   
    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
    837850
    838851    INTEGER(iwp)                                     ::  id_emis        !< NetCDF id of input file for chemistry emissions: TBD: It has to be removed
    839852
    840853    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
    841856
    842857    LOGICAL ::  input_pids_static  = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing static information exists
     
    949964!
    950965!-- Public data structures
    951     PUBLIC real_2d
     966    PUBLIC real_1d,                                                            &
     967           real_2d
    952968!
    953969!-- Public variables
     
    956972           chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type,   &
    957973           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,                                                  &
    959977           input_pids_static,                                                  &
    960978           input_pids_dynamic, input_pids_vm, input_file_vm,                   &
    961979           leaf_area_density_f, nest_offl,                                     &
     980           num_var_pids,                                                       &
    962981           pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f,       &
     982           pids_id,                                                            &
    963983           root_area_density_lad_f, root_area_density_lsm_f, soil_pars_f,      &
    964984           soil_type_f, street_crossing_f, street_type_f, surface_fraction_f,  &
    965985           terrain_height_f, vegetation_pars_f, vegetation_type_f,             &
     986           vars_pids,                                                          &
    966987           water_pars_f, water_type_f
    967988!
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r4168 r4178  
    2828! -----------------
    2929! $Id$
     30! External radiation forcing implemented.
     31!
     32! 4168 2019-08-16 13:50:17Z suehring
    3033! Replace function get_topography_top_index by topo_top_ind
    3134!
     
    667670        ONLY:  cloud_droplets, coupling_char,                                  &
    668671               debug_output, debug_output_timestep, debug_string,              &
     672               dt_3d,                                                          &
    669673               dz, dt_spinup, end_time,                                        &
    670674               humidity,                                                       &
     
    704708
    705709    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
    708730
    709731    USE plant_canopy_model_mod,                                                &
     
    12361258    REAL(wp), DIMENSION(:), ALLOCATABLE            :: albedo_surf        !< albedo of the surface
    12371259    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
    12381266
    12391267
     
    14141442          CASE ( 'rrtmg' )
    14151443             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
    14161453
    14171454          CASE DEFAULT
     
    18041841       IF ( radiation_scheme /= 'constant'   .AND.                             &
    18051842            radiation_scheme /= 'clear-sky'  .AND.                             &
    1806             radiation_scheme /= 'rrtmg' )  THEN
     1843            radiation_scheme /= 'rrtmg'      .AND.                             &
     1844            radiation_scheme /= 'external' )  THEN
    18071845          message_string = 'unknown radiation_scheme = '//                     &
    18081846                           TRIM( radiation_scheme )
     
    18991937       INTEGER(iwp) ::  l         !< running index for orientation of vertical surfaces
    19001938       INTEGER(iwp) ::  m         !< running index for surface elements
     1939       INTEGER(iwp) ::  ntime = 0 !< number of available external radiation timesteps
    19011940#if defined( __rrtmg )
    19021941       INTEGER(iwp) ::  ind_type  !< running index for subgrid-surface tiles
     
    21132152
    21142153       IF ( radiation_scheme == 'clear-sky'  .OR.                              &
    2115             radiation_scheme == 'constant')  THEN
    2116 
    2117 
     2154            radiation_scheme == 'constant'   .OR.                              &
     2155            radiation_scheme == 'external' )  THEN
    21182156!
    21192157!--       Allocate arrays for incoming/outgoing short/longwave radiation
     
    28962934#endif
    28972935       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
    28993061!
    29003062!--    Perform user actions if required
     
    29133075          CASE ( 'constant' )
    29143076             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
    29153086
    29163087          CASE DEFAULT
    29173088
    29183089       END SELECT
    2919 
    2920 ! readjust date and time to its initial value
     3090!
     3091!--    Readjust date and time to its initial value
    29213092       CALL init_date_and_time
    29223093
     
    29583129
    29593130
     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   
    29603308!------------------------------------------------------------------------------!
    29613309! Description:
     
    30613409               
    30623410                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
    30633413
    30643414                surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4   &
     
    31363486                i = surf%i(m)
    31373487                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)
    31393489                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)
    31413491                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
    31423492             ENDDO
     
    33553705          IF ( .NOT. lw_radiation )  WRITE( io, 10 )
    33563706          IF ( .NOT. sw_radiation )  WRITE( io, 11 )
     3707       ELSEIF ( radiation_scheme == "external" )  THEN
     3708          WRITE( io, 14 )
    33573709       ENDIF
    33583710
     
    3390374213 FORMAT (/'    Albedo is set individually for each xy-location, according ', &
    33913743                 'to given surface type.')
     374414 FORMAT ('    --> External radiation forcing is used')
    33923745
    33933746
     
    53285681         ENDIF
    53295682     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
    53345689
    53355690!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    61186473    SUBROUTINE calc_diffusion_radiation
    61196474   
    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
    61296477       
     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
     6486print*, "in diffusion_radiation"
     6487
    61306488!--     Calculate current day and time based on the initial values and simulation time
    61316489        year_angle = ( (day_of_year_init * 86400) + time_utc_init              &
     
    61386496                          0.000719_wp * cos(2.0_wp * year_angle) +             &
    61396497                          0.000077_wp * sin(2.0_wp * year_angle))
    6140         
     6498!       
    61416499!--   
    61426500!--     Under a very low angle, we keep extraterestrial radiation at
    61436501!--     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.   
    61466503        IF ( cos_zenith <= lowest_solarUp )  THEN
    61476504            corrected_solarUp = lowest_solarUp
     
    1145111808          IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
    1145211809             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11453                   radiation_scheme == 'constant')  THEN
     11810                  radiation_scheme == 'constant'   .OR.                    &
     11811                  radiation_scheme == 'external' )  THEN
    1145411812                ALLOCATE( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
    1145511813             ELSE
     
    1145911817          IF ( k == 1 )  THEN
    1146011818             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11461                   radiation_scheme == 'constant')  THEN
     11819                  radiation_scheme == 'constant'   .OR.                    &
     11820                  radiation_scheme == 'external' )  THEN
    1146211821                READ ( 13 )  tmp_3d2
    1146311822                rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
     
    1147311832          IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
    1147411833             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11475                   radiation_scheme == 'constant')  THEN
     11834                  radiation_scheme == 'constant'   .OR.                    &
     11835                  radiation_scheme == 'external' )  THEN
    1147611836                ALLOCATE( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
    1147711837             ELSE
     
    1148111841          IF ( k == 1 )  THEN
    1148211842             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11483                   radiation_scheme == 'constant')  THEN
     11843                  radiation_scheme == 'constant'   .OR.                    &
     11844                  radiation_scheme == 'external' )  THEN
    1148411845                READ ( 13 )  tmp_3d2
    1148511846                rad_lw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
     
    1149511856          IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
    1149611857             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11497                   radiation_scheme == 'constant')  THEN
     11858                  radiation_scheme == 'constant'   .OR.                    &
     11859                  radiation_scheme == 'external' )  THEN
    1149811860                ALLOCATE( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
    1149911861             ELSE
     
    1150311865          IF ( k == 1 )  THEN
    1150411866             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11505                   radiation_scheme == 'constant')  THEN
     11867                  radiation_scheme == 'constant'   .OR.                    &
     11868                  radiation_scheme == 'external' )  THEN
    1150611869                READ ( 13 )  tmp_3d2
    1150711870                rad_lw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
     
    1151711880          IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
    1151811881             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11519                   radiation_scheme == 'constant')  THEN
     11882                  radiation_scheme == 'constant'   .OR.                    &
     11883                  radiation_scheme == 'external' )  THEN
    1152011884                ALLOCATE( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
    1152111885             ELSE
     
    1152511889          IF ( k == 1 )  THEN
    1152611890             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11527                   radiation_scheme == 'constant')  THEN
     11891                  radiation_scheme == 'constant'   .OR.                    &
     11892                  radiation_scheme == 'external' )  THEN
    1152811893                READ ( 13 )  tmp_3d2
    1152911894                rad_lw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) &
     
    1157111936          IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
    1157211937             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11573                   radiation_scheme == 'constant')  THEN
     11938                  radiation_scheme == 'constant'   .OR.                    &
     11939                  radiation_scheme == 'external' )  THEN
    1157411940                ALLOCATE( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
    1157511941             ELSE
     
    1157911945          IF ( k == 1 )  THEN
    1158011946             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11581                   radiation_scheme == 'constant')  THEN
     11947                  radiation_scheme == 'constant'   .OR.                    &
     11948                  radiation_scheme == 'external' )  THEN
    1158211949                READ ( 13 )  tmp_3d2
    1158311950                rad_sw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
     
    1159311960          IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
    1159411961             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11595                   radiation_scheme == 'constant')  THEN
     11962                  radiation_scheme == 'constant'   .OR.                    &
     11963                  radiation_scheme == 'external' )  THEN
    1159611964                ALLOCATE( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
    1159711965             ELSE
     
    1160111969          IF ( k == 1 )  THEN
    1160211970             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11603                   radiation_scheme == 'constant')  THEN
     11971                  radiation_scheme == 'constant'   .OR.                    &
     11972                  radiation_scheme == 'external' )  THEN
    1160411973                READ ( 13 )  tmp_3d2
    1160511974                rad_sw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
     
    1161511984          IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
    1161611985             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11617                   radiation_scheme == 'constant')  THEN
     11986                  radiation_scheme == 'constant'   .OR.                    &
     11987                  radiation_scheme == 'external' )  THEN
    1161811988                ALLOCATE( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
    1161911989             ELSE
     
    1162311993          IF ( k == 1 )  THEN
    1162411994             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11625                   radiation_scheme == 'constant')  THEN
     11995                  radiation_scheme == 'constant'   .OR.                    &
     11996                  radiation_scheme == 'external' )  THEN
    1162611997                READ ( 13 )  tmp_3d2
    1162711998                rad_sw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
     
    1163712008          IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
    1163812009             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11639                   radiation_scheme == 'constant')  THEN
     12010                  radiation_scheme == 'constant'   .OR.                    &
     12011                  radiation_scheme == 'external' )  THEN
    1164012012                ALLOCATE( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
    1164112013             ELSE
     
    1164512017          IF ( k == 1 )  THEN
    1164612018             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    11647                   radiation_scheme == 'constant')  THEN
     12019                  radiation_scheme == 'constant'   .OR.                    &
     12020                  radiation_scheme == 'external' )  THEN
    1164812021                READ ( 13 )  tmp_3d2
    1164912022                rad_sw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) &
Note: See TracChangeset for help on using the changeset viewer.