Changeset 4190 for palm


Ignore:
Timestamp:
Aug 27, 2019 3:42:37 PM (5 years ago)
Author:
suehring
Message:

Implement external radiation forcing also for level-of-detail = 2 (horizontally 2D radiation)

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

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

    r4186 r4190  
    2525! -----------------
    2626! $Id$
     27! type real_1d changed to real_1d_3d
     28!
     29! 4186 2019-08-23 16:06:14Z suehring
    2730! Minor formatting adjustments
    2831!
     
    396399    END TYPE int_2d_32bit
    397400!
    398 !-- Define data type to read 1D real variables
    399     TYPE real_1d
     401!-- Define data type to read 1D or 3D real variables.
     402    TYPE real_1d_3d
    400403       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
    401404
     405       INTEGER(iwp) ::  lod = -1        !< level-of-detail
     406       
    402407       REAL(wp) ::  fill = -9999.9_wp                  !< fill value
    403408       
    404        REAL(wp), DIMENSION(:), ALLOCATABLE ::  var     !< respective variable
    405     END TYPE real_1d   
     409       REAL(wp), DIMENSION(:),     ALLOCATABLE ::  var1d     !< respective 1D variable
     410       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  var3d     !< respective 3D variable
     411    END TYPE real_1d_3d   
    406412!
    407413!-- Define data type to read 2D real variables
     
    409415       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
    410416
     417       INTEGER(iwp) ::  lod             !< level-of-detail
     418       
    411419       REAL(wp) ::  fill = -9999.9_wp                !< fill value
    412420       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var !< respective variable
     
    747755!
    748756!-- Public data structures
    749     PUBLIC real_1d,                                                            &
     757    PUBLIC real_1d_3d,                                                         &
    750758           real_2d
    751759!
    752760!-- Public variables
    753761    PUBLIC albedo_pars_f, albedo_type_f, basal_area_density_f, buildings_f,    &
    754            building_id_f, building_pars_f, building_type_f, char_fill,         &
     762           building_id_f, building_pars_f, building_type_f,                    &
     763           char_fill,                                                          &
     764           char_lod,                                                           &
    755765           chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type,   &
    756766           coord_ref_sys,                                                      &
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r4188 r4190  
    2828! -----------------
    2929! $Id$
     30! Implement external radiation forcing also for level-of-detail = 2
     31! (horizontally 2D radiation)
     32!
     33! 4188 2019-08-26 14:15:47Z suehring
    3034! Minor adjustment in error message
    3135!
     
    235239               water_type_f,                                                   &
    236240               char_fill,                                                      &
     241               char_lod,                                                       &
    237242               check_existence,                                                &
    238243               close_input_file,                                               &
     
    247252               pids_id,                                                        &
    248253               open_read_file,                                                 &
    249                real_1d,                                                        &
     254               real_1d_3d,                                                     &
    250255               vars_pids
    251256
     
    780785    REAL(wp), DIMENSION(:), ALLOCATABLE            :: emiss_surf         !< emissivity of the wall surface
    781786!
    782 !-- External radiation
    783     TYPE( real_1d ) ::  rad_lw_in_f     !< external incoming longwave radiation, from observation or model
    784     TYPE( real_1d ) ::  rad_sw_in_f     !< external incoming shortwave radiation, from observation or model
    785     TYPE( real_1d ) ::  rad_sw_in_dif_f !< external incoming shortwave radiation, diffuse part, from observation or model
    786     TYPE( real_1d ) ::  time_rad_f      !< time dimension for external radiation, from observation or model
    787 
     787!-- External radiation. Depending on the given level of detail either a 1D or
     788!-- a 3D array will be allocated.
     789    TYPE( real_1d_3d ) ::  rad_lw_in_f     !< external incoming longwave radiation, from observation or model
     790    TYPE( real_1d_3d ) ::  rad_sw_in_f     !< external incoming shortwave radiation, from observation or model
     791    TYPE( real_1d_3d ) ::  rad_sw_in_dif_f !< external incoming shortwave radiation, diffuse part, from observation or model
     792    TYPE( real_1d_3d ) ::  time_rad_f      !< time dimension for external radiation, from observation or model
    788793
    789794    INTERFACE radiation_check_data_output
     
    14581463       INTEGER(iwp) ::  l         !< running index for orientation of vertical surfaces
    14591464       INTEGER(iwp) ::  m         !< running index for surface elements
    1460        INTEGER(iwp) ::  ntime = 0 !< number of available external radiation timesteps 
     1465       INTEGER(iwp) ::  ntime = 0 !< number of available external radiation timesteps
    14611466#if defined( __rrtmg )
    14621467       INTEGER(iwp) ::  ind_type  !< running index for subgrid-surface tiles
     
    25012506          IF ( check_existence( vars_pids, 'time_rad' ) )  THEN
    25022507             CALL netcdf_data_input_get_dimension_length( pids_id,             &
    2503                                                           ntime, 'time_rad' )
     2508                                                          ntime,               &
     2509                                                          'time_rad' )
    25042510         
    2505              ALLOCATE( time_rad_f%var(0:ntime-1) )
     2511             ALLOCATE( time_rad_f%var1d(0:ntime-1) )
    25062512!                                                                                 
    25072513!--          Read variable                   
    2508              CALL get_variable( pids_id, 'time_rad', time_rad_f%var )
     2514             CALL get_variable( pids_id, 'time_rad', time_rad_f%var1d )
    25092515                               
    25102516             time_rad_f%from_file = .TRUE.
     
    25132519!--       Input shortwave downwelling.
    25142520          IF ( check_existence( vars_pids, 'rad_sw_in' ) )  THEN
    2515              ALLOCATE( rad_sw_in_f%var(0:ntime-1) )
    25162521!         
    2517 !--          Read _FillValue attribute
     2522!--          Get _FillValue attribute
    25182523             CALL get_attribute( pids_id, char_fill, rad_sw_in_f%fill,         &
    25192524                                 .FALSE., 'rad_sw_in' )
    2520 !                                                                                 
    2521 !--          Read variable                   
    2522              CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var )
    2523                                
    2524              rad_sw_in_f%from_file = .TRUE.
     2525!         
     2526!--          Get level-of-detail
     2527             CALL get_attribute( pids_id, char_lod, rad_sw_in_f%lod,           &
     2528                                 .FALSE., 'rad_sw_in' )
     2529!
     2530!--          Level-of-detail 1 - radiation depends only on time_rad
     2531             IF ( rad_sw_in_f%lod == 1 )  THEN
     2532                ALLOCATE( rad_sw_in_f%var1d(0:ntime-1) )
     2533                CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var1d )
     2534                rad_sw_in_f%from_file = .TRUE.
     2535!
     2536!--          Level-of-detail 2 - radiation depends on time_rad, y, x
     2537             ELSEIF ( rad_sw_in_f%lod == 2 )  THEN 
     2538                ALLOCATE( rad_sw_in_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
     2539               
     2540                CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var3d,    &
     2541                                   nxl, nxr, nys, nyn, 0, ntime-1 )
     2542                                   
     2543                rad_sw_in_f%from_file = .TRUE.
     2544             ELSE
     2545                message_string = '"rad_sw_in" has no valid lod attribute'
     2546                CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
     2547             ENDIF
    25252548          ENDIF
    25262549!         
    25272550!--       Input longwave downwelling.
    25282551          IF ( check_existence( vars_pids, 'rad_lw_in' ) )  THEN
    2529              ALLOCATE( rad_lw_in_f%var(0:ntime-1) )
    25302552!         
    2531 !--          Read _FillValue attribute
     2553!--          Get _FillValue attribute
    25322554             CALL get_attribute( pids_id, char_fill, rad_lw_in_f%fill,         &
    25332555                                 .FALSE., 'rad_lw_in' )
    2534 !                                                                                 
    2535 !--          Read variable                   
    2536              CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var )
    2537                                
    2538              rad_lw_in_f%from_file = .TRUE.
     2556!         
     2557!--          Get level-of-detail
     2558             CALL get_attribute( pids_id, char_lod, rad_lw_in_f%lod,           &
     2559                                 .FALSE., 'rad_lw_in' )
     2560!
     2561!--          Level-of-detail 1 - radiation depends only on time_rad
     2562             IF ( rad_lw_in_f%lod == 1 )  THEN
     2563                ALLOCATE( rad_lw_in_f%var1d(0:ntime-1) )
     2564                CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var1d )
     2565                rad_lw_in_f%from_file = .TRUE.
     2566!
     2567!--          Level-of-detail 2 - radiation depends on time_rad, y, x
     2568             ELSEIF ( rad_lw_in_f%lod == 2 )  THEN 
     2569                ALLOCATE( rad_lw_in_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
     2570               
     2571                CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var3d,    &
     2572                                   nxl, nxr, nys, nyn, 0, ntime-1 )
     2573                                   
     2574                rad_lw_in_f%from_file = .TRUE.
     2575             ELSE
     2576                message_string = '"rad_lw_in" has no valid lod attribute'
     2577                CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
     2578             ENDIF
    25392579          ENDIF
    25402580!         
    25412581!--       Input shortwave downwelling, diffuse part.
    25422582          IF ( check_existence( vars_pids, 'rad_sw_in_dif' ) )  THEN
    2543              ALLOCATE( rad_sw_in_dif_f%var(0:ntime-1) )
    25442583!         
    25452584!--          Read _FillValue attribute
    25462585             CALL get_attribute( pids_id, char_fill, rad_sw_in_dif_f%fill,     &
    25472586                                 .FALSE., 'rad_sw_in_dif' )
    2548 !                                                                                 
    2549 !--          Read variable                   
    2550              CALL get_variable( pids_id, 'rad_sw_in_dif', rad_sw_in_dif_f%var )
    2551                                
    2552              rad_sw_in_dif_f%from_file = .TRUE.
     2587!         
     2588!--          Get level-of-detail
     2589             CALL get_attribute( pids_id, char_lod, rad_sw_in_dif_f%lod,       &
     2590                                 .FALSE., 'rad_sw_in_dif' )
     2591!
     2592!--          Level-of-detail 1 - radiation depends only on time_rad
     2593             IF ( rad_sw_in_dif_f%lod == 1 )  THEN
     2594                ALLOCATE( rad_sw_in_dif_f%var1d(0:ntime-1) )
     2595                CALL get_variable( pids_id, 'rad_sw_in_dif',                   &
     2596                                   rad_sw_in_dif_f%var1d )
     2597                rad_sw_in_dif_f%from_file = .TRUE.
     2598!
     2599!--          Level-of-detail 2 - radiation depends on time_rad, y, x
     2600             ELSEIF ( rad_sw_in_dif_f%lod == 2 )  THEN 
     2601                ALLOCATE( rad_sw_in_dif_f%var3d(0:ntime-1,nys:nyn,nxl:nxr) )
     2602               
     2603                CALL get_variable( pids_id, 'rad_sw_in_dif',                   &
     2604                                   rad_sw_in_dif_f%var3d,                      &
     2605                                   nxl, nxr, nys, nyn, 0, ntime-1 )
     2606                                   
     2607                rad_sw_in_dif_f%from_file = .TRUE.
     2608             ELSE
     2609                message_string = '"rad_sw_in_dif" has no valid lod attribute'
     2610                CALL message( 'radiation_init', 'PA0646', 1, 2, 0, 6, 0 )
     2611             ENDIF
    25532612          ENDIF
    25542613!         
     
    25732632          ENDIF
    25742633         
    2575           IF ( ANY( rad_sw_in_f%var ==  rad_sw_in_f%fill ) )  THEN
    2576              message_string = 'External radiation forcing: rad_sw_in must ' // &
    2577                               'not contain any fill values'
    2578              CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 )
    2579           ENDIF
    2580          
    2581           IF ( ANY( rad_lw_in_f%var ==  rad_lw_in_f%fill ) )  THEN
    2582              message_string = 'External radiation forcing: rad_lw_in must ' // &
    2583                               'not contain any fill values'
    2584              CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 )
    2585           ENDIF
    2586          
    2587           IF ( rad_sw_in_dif_f%from_file )  THEN
    2588              IF ( ANY( rad_sw_in_dif_f%var ==  rad_sw_in_dif_f%fill ) )  THEN
    2589                 message_string = 'External radiation forcing: ' //             &
    2590                                  'rad_lw_in_dif must not contain any fill ' // &
    2591                                  'values'
    2592                 CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 )
    2593              ENDIF
    2594           ENDIF
    2595          
    2596           IF ( time_rad_f%var(0) /= 0.0_wp )  THEN
     2634          IF ( time_rad_f%var1d(0) /= 0.0_wp )  THEN
    25972635             message_string = 'External radiation forcing: first point in ' // &
    25982636                              'time is /= 0.0.'
     
    26002638          ENDIF
    26012639         
    2602           IF ( end_time - spinup_time > time_rad_f%var(ntime-1) )  THEN
     2640          IF ( end_time - spinup_time > time_rad_f%var1d(ntime-1) )  THEN
    26032641             message_string = 'External radiation forcing does not cover ' //  &
    26042642                              'the entire simulation time.'
    26052643             CALL message( 'radiation_init', 'PA0314', 1, 2, 0, 6, 0 )
    26062644          ENDIF
     2645!
     2646!--       Check for fill values in radiation
     2647          IF ( ALLOCATED( rad_sw_in_f%var1d ) )  THEN
     2648             IF ( ANY( rad_sw_in_f%var1d == rad_sw_in_f%fill ) )  THEN
     2649                message_string = 'External radiation array "rad_sw_in" ' //    &
     2650                                 'must not contain any fill values.'
     2651                CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 )
     2652             ENDIF
     2653          ENDIF
    26072654         
     2655          IF ( ALLOCATED( rad_lw_in_f%var1d ) )  THEN
     2656             IF ( ANY( rad_lw_in_f%var1d == rad_lw_in_f%fill ) )  THEN
     2657                message_string = 'External radiation array "rad_lw_in" ' //    &
     2658                                 'must not contain any fill values.'
     2659                CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 )
     2660             ENDIF
     2661          ENDIF
     2662         
     2663          IF ( ALLOCATED( rad_sw_in_dif_f%var1d ) )  THEN
     2664             IF ( ANY( rad_sw_in_dif_f%var1d == rad_sw_in_dif_f%fill ) )  THEN
     2665                message_string = 'External radiation array "rad_sw_in_dif" ' //&
     2666                                 'must not contain any fill values.'
     2667                CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 )
     2668             ENDIF
     2669          ENDIF
     2670         
     2671          IF ( ALLOCATED( rad_sw_in_f%var3d ) )  THEN
     2672             IF ( ANY( rad_sw_in_f%var3d == rad_sw_in_f%fill ) )  THEN
     2673                message_string = 'External radiation array "rad_sw_in" ' //    &
     2674                                 'must not contain any fill values.'
     2675                CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 )
     2676             ENDIF
     2677          ENDIF
     2678         
     2679          IF ( ALLOCATED( rad_lw_in_f%var3d ) )  THEN
     2680             IF ( ANY( rad_lw_in_f%var3d == rad_lw_in_f%fill ) )  THEN
     2681                message_string = 'External radiation array "rad_lw_in" ' //    &
     2682                                 'must not contain any fill values.'
     2683                CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 )
     2684             ENDIF
     2685          ENDIF
     2686         
     2687          IF ( ALLOCATED( rad_sw_in_dif_f%var3d ) )  THEN
     2688             IF ( ANY( rad_sw_in_dif_f%var3d == rad_sw_in_dif_f%fill ) )  THEN
     2689                message_string = 'External radiation array "rad_sw_in_dif" ' //&
     2690                                 'must not contain any fill values.'
     2691                CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 )
     2692             ENDIF
     2693          ENDIF
     2694!
     2695!--       Currently, 2D external radiation input is not possible in
     2696!--       combination with topography where average radiation is used.
     2697          IF ( ( rad_sw_in_f%lod == 2  .OR.  rad_sw_in_f%lod == 2  .OR.      &
     2698                 rad_sw_in_dif_f%lod == 2  )  .AND. average_radiation )  THEN
     2699             message_string = 'External radiation with lod = 2 is currently '//&
     2700                              'not possible with average_radiation = .T..'
     2701                CALL message( 'radiation_init', 'PA0670', 1, 2, 0, 6, 0 )
     2702          ENDIF
     2703!
     2704!--       All radiation input should have the same level of detail. The sum
     2705!--       of lods divided by the number of available radiation arrays must be
     2706!--       1 (if all are lod = 1) or 2 (if all are lod = 2).
     2707          IF ( REAL( MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) +       &
     2708                     MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) +       &
     2709                     MERGE( rad_sw_in_dif_f%lod, 0, rad_sw_in_dif_f%from_file ),&
     2710                     KIND = wp ) /                                              &
     2711                   ( MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) +           &
     2712                     MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) +           &
     2713                     MERGE( 1.0_wp, 0.0_wp, rad_sw_in_dif_f%from_file ) )       &
     2714                     /= 1.0_wp  .AND.                                           &
     2715               REAL( MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) +       &
     2716                     MERGE( rad_sw_in_f%lod, 0, rad_sw_in_f%from_file ) +       &
     2717                     MERGE( rad_sw_in_dif_f%lod, 0, rad_sw_in_dif_f%from_file ),&
     2718                     KIND = wp ) /                                              &
     2719                   ( MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) +           &
     2720                     MERGE( 1.0_wp, 0.0_wp, rad_sw_in_f%from_file ) +           &
     2721                     MERGE( 1.0_wp, 0.0_wp, rad_sw_in_dif_f%from_file ) )       &
     2722                     /= 2.0_wp )  THEN
     2723             message_string = 'External radiation input should have the same '//&
     2724                              'lod.'
     2725             CALL message( 'radiation_init', 'PA0673', 1, 2, 0, 6, 0 )
     2726          ENDIF
     2727
    26082728       ENDIF
    26092729!
     
    26902810       INTEGER(iwp) ::  tm  !< index of previous timestep
    26912811       
    2692        REAL(wp) ::  fac_dt     !< interpolation factor
    2693        REAL(wp) ::  lw_in      !< downwelling longwave radiation, interpolated value     
    2694        REAL(wp) ::  sw_in      !< downwelling shortwave radiation, interpolated value       
    2695        REAL(wp) ::  sw_in_dif  !< downwelling diffuse shortwave radiation, interpolated value   
     2812       REAL(wp) ::  fac_dt     !< interpolation factor 
    26962813
    26972814       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine   
     
    27082825       ELSE
    27092826          t = 0
    2710           DO WHILE ( time_rad_f%var(t) <= time_since_reference_point )
     2827          DO WHILE ( time_rad_f%var1d(t) <= time_since_reference_point )
    27112828             t = t + 1
    27122829          ENDDO
     
    27142831          tm = MAX( t-1, 0 )
    27152832         
    2716           fac_dt = ( time_since_reference_point - time_rad_f%var(tm) + dt_3d ) &
    2717                  / ( time_rad_f%var(t)  - time_rad_f%var(tm) )
     2833          fac_dt = ( time_since_reference_point - time_rad_f%var1d(tm) + dt_3d ) &
     2834                 / ( time_rad_f%var1d(t)  - time_rad_f%var1d(tm) )
    27182835          fac_dt = MIN( 1.0_wp, fac_dt )
    27192836       ENDIF
    2720              
    2721        sw_in = ( 1.0_wp - fac_dt ) * rad_sw_in_f%var(tm)                       &
    2722                           + fac_dt * rad_sw_in_f%var(t)
    2723                                          
    2724        lw_in = ( 1.0_wp - fac_dt ) * rad_lw_in_f%var(tm)                       &
    2725                           + fac_dt * rad_lw_in_f%var(t)
    2726 !
    2727 !--    Limit shortwave incoming radiation to positive values, in order
    2728 !--    to overcome possible observation errors.
    2729        sw_in = MAX( 0.0_wp, sw_in )
    2730        sw_in = MERGE( sw_in, 0.0_wp, sun_up )
    27312837!
    27322838!--    Call clear-sky calculation for each surface orientation.
     
    27442850          CALL radiation_external_surf
    27452851       ENDDO
    2746 !
    2747 !--    In case of average radiation check if external diffuse shortwave
    2748 !--    radiation is available and, if required, interpolate it onto the
    2749 !--    respective arrays.
    2750        IF ( average_radiation )  THEN
    2751           IF ( rad_sw_in_dif_f%from_file )  THEN
    2752              sw_in_dif= ( 1.0_wp - fac_dt ) * rad_sw_in_dif_f%var(tm)          &
    2753                                    + fac_dt * rad_sw_in_dif_f%var(t)
    2754                                          
    2755              IF ( ALLOCATED( rad_sw_in_diff ) )  rad_sw_in_diff = sw_in_dif
    2756              IF ( ALLOCATED( rad_sw_in_dir  ) )  rad_sw_in_dir  = sw_in - sw_in_dif
    2757 !
    2758 !--          Diffuse longwave radiation equals the total downwelling longwaver
    2759 !--          radiation
    2760              IF ( ALLOCATED( rad_lw_in_diff ) )  rad_lw_in_diff = lw_in
    2761           ENDIF
    2762        ENDIF
    27632852       
    27642853       CONTAINS
     
    27702859             IMPLICIT NONE
    27712860
    2772              INTEGER(iwp) ::  i         !< grid index along x-dimension   
    2773              INTEGER(iwp) ::  j         !< grid index along y-dimension 
    2774              INTEGER(iwp) ::  k         !< grid index along z-dimension   
    2775              INTEGER(iwp) ::  m         !< running index for surface elements       
     2861             INTEGER(iwp) ::  i    !< grid index along x-dimension   
     2862             INTEGER(iwp) ::  j    !< grid index along y-dimension 
     2863             INTEGER(iwp) ::  k    !< grid index along z-dimension   
     2864             INTEGER(iwp) ::  m    !< running index for surface elements
     2865             
     2866             REAL(wp) ::  lw_in     !< downwelling longwave radiation, interpolated value     
     2867             REAL(wp) ::  sw_in     !< downwelling shortwave radiation, interpolated value
     2868             REAL(wp) ::  sw_in_dif !< downwelling diffuse shortwave radiation, interpolated value   
    27762869
    27772870             IF ( surf%ns < 1 )  RETURN
     2871!
     2872!--          level-of-detail = 1. Note, here it must be distinguished between
     2873!--          averaged radiation and non-averaged radiation for the upwelling
     2874!--          fluxes.
     2875             IF ( rad_sw_in_f%lod == 1 )  THEN
     2876             
     2877                sw_in = ( 1.0_wp - fac_dt ) * rad_sw_in_f%var1d(tm)            &
     2878                                   + fac_dt * rad_sw_in_f%var1d(t)
     2879                                         
     2880                lw_in = ( 1.0_wp - fac_dt ) * rad_lw_in_f%var1d(tm)            &
     2881                                   + fac_dt * rad_lw_in_f%var1d(t)
     2882!
     2883!--             Limit shortwave incoming radiation to positive values, in order
     2884!--             to overcome possible observation errors.
     2885                sw_in = MAX( 0.0_wp, sw_in )
     2886                sw_in = MERGE( sw_in, 0.0_wp, sun_up )
    27782887                         
    2779              surf%rad_sw_in = sw_in                                         
    2780              surf%rad_lw_in = lw_in
    2781 
    2782              IF ( average_radiation )  THEN
     2888                surf%rad_sw_in = sw_in                                         
     2889                surf%rad_lw_in = lw_in
    27832890             
    2784                 surf%rad_sw_out = albedo_urb * surf%rad_sw_in
     2891                IF ( average_radiation )  THEN
     2892                   surf%rad_sw_out = albedo_urb * surf%rad_sw_in
     2893                   
     2894                   surf%rad_lw_out = emissivity_urb * sigma_sb * t_rad_urb**4  &
     2895                                  + ( 1.0_wp - emissivity_urb ) * surf%rad_lw_in
     2896                   
     2897                   surf%rad_net = surf%rad_sw_in - surf%rad_sw_out             &
     2898                                + surf%rad_lw_in - surf%rad_lw_out
     2899                   
     2900                   surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb          &
     2901                                                     * sigma_sb                &
     2902                                                     * t_rad_urb**3
     2903                ELSE
     2904                   DO  m = 1, surf%ns
     2905                      k = surf%k(m)
     2906                      surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m)  *      &
     2907                                             surf%albedo(ind_veg_wall,m)       &
     2908                                           + surf%frac(ind_pav_green,m) *      &
     2909                                             surf%albedo(ind_pav_green,m)      &
     2910                                           + surf%frac(ind_wat_win,m)   *      &
     2911                                             surf%albedo(ind_wat_win,m) )      &
     2912                                           * surf%rad_sw_in(m)
     2913                   
     2914                      surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m)  *      &
     2915                                             surf%emissivity(ind_veg_wall,m)   &
     2916                                           + surf%frac(ind_pav_green,m) *      &
     2917                                             surf%emissivity(ind_pav_green,m)  &
     2918                                           + surf%frac(ind_wat_win,m)   *      &
     2919                                             surf%emissivity(ind_wat_win,m)    &
     2920                                           )                                   &
     2921                                           * sigma_sb                          &
     2922                                           * ( surf%pt_surface(m) * exner(k) )**4
     2923                   
     2924                      surf%rad_lw_out_change_0(m) =                            &
     2925                                         ( surf%frac(ind_veg_wall,m)  *        &
     2926                                           surf%emissivity(ind_veg_wall,m)     &
     2927                                         + surf%frac(ind_pav_green,m) *        &
     2928                                           surf%emissivity(ind_pav_green,m)    &
     2929                                         + surf%frac(ind_wat_win,m)   *        &
     2930                                           surf%emissivity(ind_wat_win,m)      &
     2931                                         ) * 4.0_wp * sigma_sb                 &
     2932                                         * ( surf%pt_surface(m) * exner(k) )**3
     2933                   ENDDO
    27852934               
    2786                 surf%rad_lw_out = emissivity_urb * sigma_sb * t_rad_urb**4     &
    2787                                  + ( 1.0_wp - emissivity_urb ) * surf%rad_lw_in
    2788                
    2789                 surf%rad_net = surf%rad_sw_in - surf%rad_sw_out                &
    2790                              + surf%rad_lw_in - surf%rad_lw_out
    2791                
    2792                 surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb * sigma_sb  &
    2793                                                   * t_rad_urb**3
    2794 !
    2795 !--          Calculate upwelling radiation fluxes and net radiation for each
    2796 !--          surface element depending on its properties.
     2935                ENDIF
     2936!
     2937!--             If diffuse shortwave radiation is available, store it on
     2938!--             the respective files.
     2939                IF ( rad_sw_in_dif_f%from_file )  THEN
     2940                   sw_in_dif= ( 1.0_wp - fac_dt ) * rad_sw_in_dif_f%var1d(tm)  &
     2941                                         + fac_dt * rad_sw_in_dif_f%var1d(t)
     2942                                         
     2943                   IF ( ALLOCATED( rad_sw_in_diff ) )  rad_sw_in_diff = sw_in_dif
     2944                   IF ( ALLOCATED( rad_sw_in_dir  ) )  rad_sw_in_dir  = sw_in  &
     2945                                                                    - sw_in_dif
     2946!             
     2947!--                Diffuse longwave radiation equals the total downwelling
     2948!--                longwave radiation
     2949                   IF ( ALLOCATED( rad_lw_in_diff ) )  rad_lw_in_diff = lw_in
     2950                ENDIF
     2951!
     2952!--          level-of-detail = 2
    27972953             ELSE
    27982954
     
    28012957                   j = surf%j(m)
    28022958                   k = surf%k(m)
     2959                   
     2960                   surf%rad_sw_in(m) = ( 1.0_wp - fac_dt )                     &
     2961                                            * rad_sw_in_f%var3d(tm,j,i)        &
     2962                                   + fac_dt * rad_sw_in_f%var3d(t,j,i)                                   
     2963!
     2964!--                Limit shortwave incoming radiation to positive values, in
     2965!--                order to overcome possible observation errors.
     2966                   surf%rad_sw_in(m) = MAX( 0.0_wp, surf%rad_sw_in(m) )
     2967                   surf%rad_sw_in(m) = MERGE( surf%rad_sw_in(m), 0.0_wp, sun_up )
     2968                                         
     2969                   surf%rad_lw_in(m) = ( 1.0_wp - fac_dt )                     &
     2970                                            * rad_lw_in_f%var3d(tm,j,i)        &
     2971                                   + fac_dt * rad_lw_in_f%var3d(t,j,i)
    28032972!
    28042973!--                Weighted average according to surface fraction.
     
    28333002                   surf%rad_net(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m)    &
    28343003                                   + surf%rad_lw_in(m) - surf%rad_lw_out(m)
     3004!
     3005!--                If diffuse shortwave radiation is available, store it on
     3006!--                the respective files.
     3007                   IF ( rad_sw_in_dif_f%from_file )  THEN
     3008                      IF ( ALLOCATED( rad_sw_in_diff ) )                       &
     3009                         rad_sw_in_diff(j,i) = ( 1.0_wp - fac_dt )             &
     3010                                              * rad_sw_in_dif_f%var3d(tm,j,i)  &
     3011                                     + fac_dt * rad_sw_in_dif_f%var3d(t,j,i)
     3012!
     3013!--                   dir = sw_in - sw_in_dif.
     3014                      IF ( ALLOCATED( rad_sw_in_dir  ) )                       &
     3015                         rad_sw_in_dir(j,i)  = surf%rad_sw_in(m) -             &
     3016                                               rad_sw_in_diff(j,i)
     3017!                 
     3018!--                   Diffuse longwave radiation equals the total downwelling
     3019!--                   longwave radiation
     3020                      IF ( ALLOCATED( rad_lw_in_diff ) )                       &
     3021                         rad_lw_in_diff = surf%rad_lw_in(m)
     3022                   ENDIF
    28353023
    28363024                ENDDO
    28373025
    28383026             ENDIF
    2839                
    2840 !
    2841 !--          Fill out values in radiation arrays
     3027!
     3028!--          Store radiation also on 2D arrays, which are still used for
     3029!--          direct-diffuse splitting.
    28423030             DO  m = 1, surf%ns
    28433031                i = surf%i(m)
    28443032                j = surf%j(m)
     3033               
     3034                rad_sw_in(0,:,:)  = surf%rad_sw_in(m)
     3035                rad_lw_in(0,:,:)  = surf%rad_lw_in(m)
    28453036                rad_sw_out(0,j,i) = surf%rad_sw_out(m)
    28463037                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
    28473038             ENDDO
    2848              
    2849              rad_sw_in(0,:,:) = sw_in
    2850              rad_lw_in(0,:,:) = lw_in
    28513039 
    28523040          END SUBROUTINE radiation_external_surf
Note: See TracChangeset for help on using the changeset viewer.