Changeset 3607 for palm/trunk/SOURCE


Ignore:
Timestamp:
Dec 7, 2018 11:56:58 AM (3 years ago)
Author:
suehring
Message:

Output of radiation-related quantities migrated from urban_surface_model_mod to radiation_model_mod

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

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

    r3589 r3607  
    2828! -----------------
    2929! $Id$
     30! Output of radiation-related quantities migrated to radiation_model_mod.
     31!
     32! 3589 2018-11-30 15:09:51Z suehring
    3033! Remove erroneous UTF encoding
    3134!
     
    545548               message_string, plant_canopy, pt_surface,                       &
    546549               rho_surface, simulated_time, spinup_time, surface_pressure,     &
    547                time_since_reference_point, urban_surface
     550               time_since_reference_point, urban_surface, varnamelength
    548551
    549552    USE cpulog,                                                                &
     
    899902    INTEGER(iwp)                                   ::  nwalls           !< number of wall surfaces in local processor
    900903
     904!-- indices needed for RTM netcdf output subroutines
     905    INTEGER(iwp), PARAMETER                        :: nd = 5
     906    CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
     907    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
     908    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)
     909    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirstart
     910    INTEGER(iwp), DIMENSION(0:nd-1)                :: dirend
     911
    901912!-- indices and sizes of urban and land surface models
    902913    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_l          !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
     
    10651076    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_dist
    10661077
     1078!-- arrays for time averages
     1079    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfradnet_av    !< average of net radiation to local surface including radiation from reflections
     1080    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw_av      !< average of sw radiation falling to local surface including radiation from reflections
     1081    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw_av      !< average of lw radiation falling to local surface including radiation from reflections
     1082    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir_av   !< average of direct sw radiation falling to local surface
     1083    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif_av   !< average of diffuse sw radiation from sky and model boundary falling to local surface
     1084    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif_av   !< average of diffuse lw radiation from sky and model boundary falling to local surface
     1085    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswref_av   !< average of sw radiation falling to surface from reflections
     1086    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwref_av   !< average of lw radiation falling to surface from reflections
     1087    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
     1088    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
     1089    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
     1090    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
     1091    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw_av       !< Average of pcbinlw
     1092    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw_av       !< Average of pcbinsw
     1093    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir_av    !< Average of pcbinswdir
     1094    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif_av    !< Average of pcbinswdif
     1095    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswref_av    !< Average of pcbinswref
    10671096
    10681097
     
    12081237           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
    12091238           zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
    1210            nrefsteps, nsvfl, svf,                                              &
    1211            svfsurf, surfinsw, surfinlw, surfins, surfinl, surfinswdir,         &
    1212            surfinswdif, surfoutsw, surfoutlw, surfinlwdif, rad_sw_in_dir,      &
    1213            rad_sw_in_diff, rad_lw_in_diff, surfouts, surfoutl, surfoutsl,      &
    1214            surfoutll, idir, jdir, kdir, id, iz, iy, ix,                        &
    1215            surf, surfl, nsurfl, pcbinswdir, pcbinswdif, pcbinsw, pcbinlw,      &
     1239           idir, jdir, kdir, id, iz, iy, ix,                                   &
    12161240           iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
    12171241           iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
     
    12191243           idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
    12201244           radiation_interactions, startwall, startland, endland, endwall,     &
    1221            skyvf, skyvft, radiation_interactions_on, average_radiation, npcbl, &
    1222            pcbl
     1245           skyvf, skyvft, radiation_interactions_on, average_radiation
     1246
    12231247
    12241248#if defined ( __rrtmg )
     
    12631287!> Check data output for radiation model
    12641288!------------------------------------------------------------------------------!
    1265     SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
     1289    SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k )
    12661290 
    12671291 
     
    12711295       IMPLICIT NONE
    12721296
    1273        CHARACTER (LEN=*) ::  unit     !<
    1274        CHARACTER (LEN=*) ::  var      !<
    1275 
    1276        INTEGER(iwp) :: i
     1297       CHARACTER (LEN=*) ::  unit          !<
     1298       CHARACTER (LEN=*) ::  variable      !<
     1299
     1300       INTEGER(iwp) :: i, j, k, l
    12771301       INTEGER(iwp) :: ilen
    1278        INTEGER(iwp) :: k
    1279 
    1280        SELECT CASE ( TRIM( var ) )
    1281 
    1282           CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', &
    1283                  'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
    1284              IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
    1285                 message_string = '"output of "' // TRIM( var ) // '" requi' // &
    1286                                  'res radiation = .TRUE. and ' //              &
    1287                                  'radiation_scheme = "rrtmg"'
    1288                 CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
    1289              ENDIF
    1290              unit = 'K/h'     
    1291 
    1292           CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
    1293                  'rrtm_asdir*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',     &
    1294                  'rad_sw_out*')
    1295              IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
    1296                 ! Workaround for masked output (calls with i=ilen=k=0)
    1297                 unit = 'illegal'
    1298                 RETURN
    1299              ENDIF
    1300              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    1301                 message_string = 'illegal value for data_output: "' //         &
    1302                                  TRIM( var ) // '" & only 2d-horizontal ' //   &
    1303                                  'cross sections are allowed for this value'
    1304                 CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
    1305              ENDIF
    1306              IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
    1307                 IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
    1308                      TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
    1309                      TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
    1310                      TRIM( var ) == 'rrtm_asdir*'      )                       &
    1311                 THEN
    1312                    message_string = 'output of "' // TRIM( var ) // '" require'&
    1313                                     // 's radiation = .TRUE. and radiation_sch'&
    1314                                     // 'eme = "rrtmg"'
    1315                    CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
    1316                 ENDIF
    1317              ENDIF
    1318 
    1319              IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
    1320              IF ( TRIM( var ) == 'rad_lw_in*'    ) unit = 'W/m2'
    1321              IF ( TRIM( var ) == 'rad_lw_out*'   ) unit = 'W/m2'
    1322              IF ( TRIM( var ) == 'rad_sw_in*'    ) unit = 'W/m2'
    1323              IF ( TRIM( var ) == 'rad_sw_out*'   ) unit = 'W/m2'
    1324              IF ( TRIM( var ) == 'rad_sw_in'     ) unit = 'W/m2'
    1325              IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
    1326              IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = ''
    1327              IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = ''
    1328              IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = ''
    1329 
    1330           CASE ( 'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw'  )
    1331 
    1332              IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
    1333                 ! Workaround for masked output (calls with i=ilen=k=0)
    1334                 unit = 'illegal'
    1335                 RETURN
    1336              ENDIF
    1337 
    1338              IF ( .NOT.  radiation ) THEN
     1302       CHARACTER(LEN=varnamelength) :: var          !< TRIM(variable)
     1303
     1304       var = TRIM(variable)
     1305
     1306!--    first process diractional variables
     1307       IF ( var(1:12) == 'rtm_rad_net_'  .OR.  var(1:13) == 'rtm_rad_insw_'  .OR.        &
     1308            var(1:13) == 'rtm_rad_inlw_'  .OR.  var(1:16) == 'rtm_rad_inswdir_'  .OR.    &
     1309            var(1:16) == 'rtm_rad_inswdif_'  .OR.  var(1:16) == 'rtm_rad_inswref_'  .OR. &
     1310            var(1:16) == 'rtm_rad_inlwdif_'  .OR.  var(1:16) == 'rtm_rad_inlwref_'  .OR. &
     1311            var(1:14) == 'rtm_rad_outsw_'  .OR.  var(1:14) == 'rtm_rad_outlw_'  .OR.     &
     1312            var(1:14) == 'rtm_rad_ressw_'  .OR.  var(1:14) == 'rtm_rad_reslw_'  ) THEN
     1313          IF ( .NOT.  radiation ) THEN
    13391314                message_string = 'output of "' // TRIM( var ) // '" require'&
    13401315                                 // 's radiation = .TRUE.'
    13411316                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
    1342              ENDIF
    1343              IF ( mrt_nlevels == 0 ) THEN
     1317          ENDIF
     1318          unit = 'W/m2'
     1319       ELSE IF ( var(1:7) == 'rtm_svf'  .OR.  var(1:7) == 'rtm_dif'  .OR.                &
     1320                 var(1:9) == 'rtm_skyvf' .OR. var(1:9) == 'rtm_skyvft' )  THEN
     1321          IF ( .NOT.  radiation ) THEN
    13441322                message_string = 'output of "' // TRIM( var ) // '" require'&
    1345                                  // 's mrt_nlevels > 0'
    1346                 CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
    1347              ENDIF
    1348              IF ( TRIM( var ) == 'rad_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
    1349                 message_string = 'output of "' // TRIM( var ) // '" require'&
    1350                                  // 's rad_mrt_sw = .TRUE.'
    1351                 CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
    1352              ENDIF
    1353              IF ( TRIM( var ) == 'rad_mrt' ) THEN
    1354                 unit = 'K'
    1355              ELSE
    1356                 unit = 'W m-2'
    1357              ENDIF
    1358 
    1359           CASE DEFAULT
    1360              unit = 'illegal'
    1361 
    1362        END SELECT
    1363 
     1323                                 // 's radiation = .TRUE.'
     1324                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     1325          ENDIF
     1326          unit = '1'
     1327       ELSE
     1328!--       non-directional variables
     1329          SELECT CASE ( TRIM( var ) )
     1330             CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', &
     1331                    'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out'  )
     1332                IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
     1333                   message_string = '"output of "' // TRIM( var ) // '" requi' // &
     1334                                    'res radiation = .TRUE. and ' //              &
     1335                                    'radiation_scheme = "rrtmg"'
     1336                   CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
     1337                ENDIF
     1338                unit = 'K/h'
     1339
     1340             CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
     1341                    'rrtm_asdir*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',     &
     1342                    'rad_sw_out*')
     1343                IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
     1344                   ! Workaround for masked output (calls with i=ilen=k=0)
     1345                   unit = 'illegal'
     1346                   RETURN
     1347                ENDIF
     1348                IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
     1349                   message_string = 'illegal value for data_output: "' //         &
     1350                                    TRIM( var ) // '" & only 2d-horizontal ' //   &
     1351                                    'cross sections are allowed for this value'
     1352                   CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
     1353                ENDIF
     1354                IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
     1355                   IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
     1356                        TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
     1357                        TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
     1358                        TRIM( var ) == 'rrtm_asdir*'      )                       &
     1359                   THEN
     1360                      message_string = 'output of "' // TRIM( var ) // '" require'&
     1361                                       // 's radiation = .TRUE. and radiation_sch'&
     1362                                       // 'eme = "rrtmg"'
     1363                      CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
     1364                   ENDIF
     1365                ENDIF
     1366
     1367                IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'
     1368                IF ( TRIM( var ) == 'rad_lw_in*'    ) unit = 'W/m2'
     1369                IF ( TRIM( var ) == 'rad_lw_out*'   ) unit = 'W/m2'
     1370                IF ( TRIM( var ) == 'rad_sw_in*'    ) unit = 'W/m2'
     1371                IF ( TRIM( var ) == 'rad_sw_out*'   ) unit = 'W/m2'
     1372                IF ( TRIM( var ) == 'rad_sw_in'     ) unit = 'W/m2'
     1373                IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''
     1374                IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = ''
     1375                IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = ''
     1376                IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = ''
     1377
     1378             CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', &
     1379                    'rtm_rad_pc_inswdif', 'rtm_rad_pc_inswref')
     1380                IF ( .NOT.  radiation ) THEN
     1381                   message_string = 'output of "' // TRIM( var ) // '" require'&
     1382                                    // 's radiation = .TRUE.'
     1383                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     1384                ENDIF
     1385                unit = 'W'
     1386
     1387             CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw'  )
     1388                IF ( i == 0 .AND. ilen == 0 .AND. k == 0)  THEN
     1389                   ! Workaround for masked output (calls with i=ilen=k=0)
     1390                   unit = 'illegal'
     1391                   RETURN
     1392                ENDIF
     1393
     1394                IF ( .NOT.  radiation ) THEN
     1395                   message_string = 'output of "' // TRIM( var ) // '" require'&
     1396                                    // 's radiation = .TRUE.'
     1397                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     1398                ENDIF
     1399                IF ( mrt_nlevels == 0 ) THEN
     1400                   message_string = 'output of "' // TRIM( var ) // '" require'&
     1401                                    // 's mrt_nlevels > 0'
     1402                   CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
     1403                ENDIF
     1404                IF ( TRIM( var ) == 'rtm_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
     1405                   message_string = 'output of "' // TRIM( var ) // '" require'&
     1406                                    // 's rtm_mrt_sw = .TRUE.'
     1407                   CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
     1408                ENDIF
     1409                IF ( TRIM( var ) == 'rtm_mrt' ) THEN
     1410                   unit = 'K'
     1411                ELSE
     1412                   unit = 'W m-2'
     1413                ENDIF
     1414
     1415             CASE DEFAULT
     1416                unit = 'illegal'
     1417
     1418          END SELECT
     1419       ENDIF
    13641420
    13651421    END SUBROUTINE radiation_check_data_output
     
    58265882       endwall = nsurfl
    58275883       nwalls  = endwall - startwall + 1
     5884       dirstart = (/ startland, startwall, startwall, startwall, startwall /)
     5885       dirend = (/ endland, endwall, endwall, endwall, endwall /)
    58285886
    58295887!--    fill gridpcbl and pcbl
     
    86568714    INTEGER(iwp) ::  l, m !< index of current surface element
    86578715
     8716    INTEGER(iwp)                                       :: ids, idsint_u, idsint_l, isurf
     8717    CHARACTER(LEN=varnamelength)                       :: var
     8718
     8719!-- find the real name of the variable
     8720    ids = -1
     8721    l = -1
     8722    var = TRIM(variable)
     8723    DO i = 0, nd-1
     8724        k = len(TRIM(var))
     8725        j = len(TRIM(dirname(i)))
     8726        IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
     8727            ids = i
     8728            idsint_u = dirint_u(ids)
     8729            idsint_l = dirint_l(ids)
     8730            var = var(:k-j)
     8731            EXIT
     8732        ENDIF
     8733    ENDDO
     8734    IF ( ids == -1 )  THEN
     8735        var = TRIM(variable)
     8736    ENDIF
     8737
    86588738    IF ( mode == 'allocate' )  THEN
    86598739
    8660        SELECT CASE ( TRIM( variable ) )
    8661 
     8740       SELECT CASE ( TRIM( var ) )
     8741!--          block of large scale (e.g. RRTMG) radiation output variables
    86628742             CASE ( 'rad_net*' )
    86638743                IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
     
    87388818                rad_sw_hr_av = 0.0_wp
    87398819
    8740              CASE ( 'rad_mrt_sw' )
     8820!--          block of RTM output variables
     8821             CASE ( 'rtm_rad_net' )
     8822!--              array of complete radiation balance
     8823                 IF ( .NOT.  ALLOCATED(surfradnet_av) )  THEN
     8824                     ALLOCATE( surfradnet_av(nsurfl) )
     8825                     surfradnet_av = 0.0_wp
     8826                 ENDIF
     8827
     8828             CASE ( 'rtm_rad_insw' )
     8829!--                 array of sw radiation falling to surface after i-th reflection
     8830                 IF ( .NOT.  ALLOCATED(surfinsw_av) )  THEN
     8831                     ALLOCATE( surfinsw_av(nsurfl) )
     8832                     surfinsw_av = 0.0_wp
     8833                 ENDIF
     8834
     8835             CASE ( 'rtm_rad_inlw' )
     8836!--                 array of lw radiation falling to surface after i-th reflection
     8837                 IF ( .NOT.  ALLOCATED(surfinlw_av) )  THEN
     8838                     ALLOCATE( surfinlw_av(nsurfl) )
     8839                     surfinlw_av = 0.0_wp
     8840                 ENDIF
     8841
     8842             CASE ( 'rtm_rad_inswdir' )
     8843!--                 array of direct sw radiation falling to surface from sun
     8844                 IF ( .NOT.  ALLOCATED(surfinswdir_av) )  THEN
     8845                     ALLOCATE( surfinswdir_av(nsurfl) )
     8846                     surfinswdir_av = 0.0_wp
     8847                 ENDIF
     8848
     8849             CASE ( 'rtm_rad_inswdif' )
     8850!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
     8851                 IF ( .NOT.  ALLOCATED(surfinswdif_av) )  THEN
     8852                     ALLOCATE( surfinswdif_av(nsurfl) )
     8853                     surfinswdif_av = 0.0_wp
     8854                 ENDIF
     8855
     8856             CASE ( 'rtm_rad_inswref' )
     8857!--                 array of sw radiation falling to surface from reflections
     8858                 IF ( .NOT.  ALLOCATED(surfinswref_av) )  THEN
     8859                     ALLOCATE( surfinswref_av(nsurfl) )
     8860                     surfinswref_av = 0.0_wp
     8861                 ENDIF
     8862
     8863             CASE ( 'rtm_rad_inlwdif' )
     8864!--                 array of sw radiation falling to surface after i-th reflection
     8865                IF ( .NOT.  ALLOCATED(surfinlwdif_av) )  THEN
     8866                     ALLOCATE( surfinlwdif_av(nsurfl) )
     8867                     surfinlwdif_av = 0.0_wp
     8868                 ENDIF
     8869
     8870             CASE ( 'rtm_rad_inlwref' )
     8871!--                 array of lw radiation falling to surface from reflections
     8872                 IF ( .NOT.  ALLOCATED(surfinlwref_av) )  THEN
     8873                     ALLOCATE( surfinlwref_av(nsurfl) )
     8874                     surfinlwref_av = 0.0_wp
     8875                 ENDIF
     8876
     8877             CASE ( 'rtm_rad_outsw' )
     8878!--                 array of sw radiation emitted from surface after i-th reflection
     8879                 IF ( .NOT.  ALLOCATED(surfoutsw_av) )  THEN
     8880                     ALLOCATE( surfoutsw_av(nsurfl) )
     8881                     surfoutsw_av = 0.0_wp
     8882                 ENDIF
     8883
     8884             CASE ( 'rtm_rad_outlw' )
     8885!--                 array of lw radiation emitted from surface after i-th reflection
     8886                 IF ( .NOT.  ALLOCATED(surfoutlw_av) )  THEN
     8887                     ALLOCATE( surfoutlw_av(nsurfl) )
     8888                     surfoutlw_av = 0.0_wp
     8889                 ENDIF
     8890             CASE ( 'rtm_rad_ressw' )
     8891!--                 array of residua of sw radiation absorbed in surface after last reflection
     8892                 IF ( .NOT.  ALLOCATED(surfins_av) )  THEN
     8893                     ALLOCATE( surfins_av(nsurfl) )
     8894                     surfins_av = 0.0_wp
     8895                 ENDIF
     8896
     8897             CASE ( 'rtm_rad_reslw' )
     8898!--                 array of residua of lw radiation absorbed in surface after last reflection
     8899                 IF ( .NOT.  ALLOCATED(surfinl_av) )  THEN
     8900                     ALLOCATE( surfinl_av(nsurfl) )
     8901                     surfinl_av = 0.0_wp
     8902                 ENDIF
     8903
     8904             CASE ( 'rtm_rad_pc_inlw' )
     8905!--                 array of of lw radiation absorbed in plant canopy
     8906                 IF ( .NOT.  ALLOCATED(pcbinlw_av) )  THEN
     8907                     ALLOCATE( pcbinlw_av(1:npcbl) )
     8908                     pcbinlw_av = 0.0_wp
     8909                 ENDIF
     8910
     8911             CASE ( 'rtm_rad_pc_insw' )
     8912!--                 array of of sw radiation absorbed in plant canopy
     8913                 IF ( .NOT.  ALLOCATED(pcbinsw_av) )  THEN
     8914                     ALLOCATE( pcbinsw_av(1:npcbl) )
     8915                     pcbinsw_av = 0.0_wp
     8916                 ENDIF
     8917
     8918             CASE ( 'rtm_rad_pc_inswdir' )
     8919!--                 array of of direct sw radiation absorbed in plant canopy
     8920                 IF ( .NOT.  ALLOCATED(pcbinswdir_av) )  THEN
     8921                     ALLOCATE( pcbinswdir_av(1:npcbl) )
     8922                     pcbinswdir_av = 0.0_wp
     8923                 ENDIF
     8924
     8925             CASE ( 'rtm_rad_pc_inswdif' )
     8926!--                 array of of diffuse sw radiation absorbed in plant canopy
     8927                 IF ( .NOT.  ALLOCATED(pcbinswdif_av) )  THEN
     8928                     ALLOCATE( pcbinswdif_av(1:npcbl) )
     8929                     pcbinswdif_av = 0.0_wp
     8930                 ENDIF
     8931
     8932             CASE ( 'rtm_rad_pc_inswref' )
     8933!--                 array of of reflected sw radiation absorbed in plant canopy
     8934                 IF ( .NOT.  ALLOCATED(pcbinswref_av) )  THEN
     8935                     ALLOCATE( pcbinswref_av(1:npcbl) )
     8936                     pcbinswref_av = 0.0_wp
     8937                 ENDIF
     8938
     8939             CASE ( 'rtm_mrt_sw' )
    87418940                IF ( .NOT. ALLOCATED( mrtinsw_av ) )  THEN
    87428941                   ALLOCATE( mrtinsw_av(nmrtbl) )
     
    87448943                mrtinsw_av = 0.0_wp
    87458944
    8746              CASE ( 'rad_mrt_lw' )
     8945             CASE ( 'rtm_mrt_lw' )
    87478946                IF ( .NOT. ALLOCATED( mrtinlw_av ) )  THEN
    87488947                   ALLOCATE( mrtinlw_av(nmrtbl) )
     
    87508949                mrtinlw_av = 0.0_wp
    87518950
    8752              CASE ( 'rad_mrt' )
     8951             CASE ( 'rtm_mrt' )
    87538952                IF ( .NOT. ALLOCATED( mrt_av ) )  THEN
    87548953                   ALLOCATE( mrt_av(nmrtbl) )
     
    87638962    ELSEIF ( mode == 'sum' )  THEN
    87648963
    8765        SELECT CASE ( TRIM( variable ) )
    8766 
     8964       SELECT CASE ( TRIM( var ) )
     8965!--       block of large scale (e.g. RRTMG) radiation output variables
    87678966          CASE ( 'rad_net*' )
    87688967             IF ( ALLOCATED( rad_net_av ) ) THEN
     
    89719170             ENDIF
    89729171
     9172!--       block of RTM output variables
     9173          CASE ( 'rtm_rad_net' )
     9174!--           array of complete radiation balance
     9175              DO isurf = dirstart(ids), dirend(ids)
     9176                 IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9177                    surfradnet_av(isurf) = surfinsw(isurf) - surfoutsw(isurf) +  surfinlw(isurf) - surfoutlw(isurf)
     9178                 ENDIF
     9179              ENDDO
     9180
     9181          CASE ( 'rtm_rad_insw' )
     9182!--           array of sw radiation falling to surface after i-th reflection
     9183              DO isurf = dirstart(ids), dirend(ids)
     9184                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9185                      surfinsw_av(isurf) = surfinsw_av(isurf) + surfinsw(isurf)
     9186                  ENDIF
     9187              ENDDO
     9188
     9189          CASE ( 'rtm_rad_inlw' )
     9190!--           array of lw radiation falling to surface after i-th reflection
     9191              DO isurf = dirstart(ids), dirend(ids)
     9192                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9193                      surfinlw_av(isurf) = surfinlw_av(isurf) + surfinlw(isurf)
     9194                  ENDIF
     9195              ENDDO
     9196
     9197          CASE ( 'rtm_rad_inswdir' )
     9198!--           array of direct sw radiation falling to surface from sun
     9199              DO isurf = dirstart(ids), dirend(ids)
     9200                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9201                      surfinswdir_av(isurf) = surfinswdir_av(isurf) + surfinswdir(isurf)
     9202                  ENDIF
     9203              ENDDO
     9204
     9205          CASE ( 'rtm_rad_inswdif' )
     9206!--           array of difusion sw radiation falling to surface from sky and borders of the domain
     9207              DO isurf = dirstart(ids), dirend(ids)
     9208                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9209                      surfinswdif_av(isurf) = surfinswdif_av(isurf) + surfinswdif(isurf)
     9210                  ENDIF
     9211              ENDDO
     9212
     9213          CASE ( 'rtm_rad_inswref' )
     9214!--           array of sw radiation falling to surface from reflections
     9215              DO isurf = dirstart(ids), dirend(ids)
     9216                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9217                      surfinswref_av(isurf) = surfinswref_av(isurf) + surfinsw(isurf) - &
     9218                                          surfinswdir(isurf) - surfinswdif(isurf)
     9219                  ENDIF
     9220              ENDDO
     9221
     9222
     9223          CASE ( 'rtm_rad_inlwdif' )
     9224!--           array of sw radiation falling to surface after i-th reflection
     9225              DO isurf = dirstart(ids), dirend(ids)
     9226                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9227                      surfinlwdif_av(isurf) = surfinlwdif_av(isurf) + surfinlwdif(isurf)
     9228                  ENDIF
     9229              ENDDO
     9230!
     9231          CASE ( 'rtm_rad_inlwref' )
     9232!--           array of lw radiation falling to surface from reflections
     9233              DO isurf = dirstart(ids), dirend(ids)
     9234                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9235                      surfinlwref_av(isurf) = surfinlwref_av(isurf) + &
     9236                                          surfinlw(isurf) - surfinlwdif(isurf)
     9237                  ENDIF
     9238              ENDDO
     9239
     9240          CASE ( 'rtm_rad_outsw' )
     9241!--           array of sw radiation emitted from surface after i-th reflection
     9242              DO isurf = dirstart(ids), dirend(ids)
     9243                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9244                      surfoutsw_av(isurf) = surfoutsw_av(isurf) + surfoutsw(isurf)
     9245                  ENDIF
     9246              ENDDO
     9247
     9248          CASE ( 'rtm_rad_outlw' )
     9249!--           array of lw radiation emitted from surface after i-th reflection
     9250              DO isurf = dirstart(ids), dirend(ids)
     9251                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9252                      surfoutlw_av(isurf) = surfoutlw_av(isurf) + surfoutlw(isurf)
     9253                  ENDIF
     9254              ENDDO
     9255
     9256          CASE ( 'rtm_rad_ressw' )
     9257!--           array of residua of sw radiation absorbed in surface after last reflection
     9258              DO isurf = dirstart(ids), dirend(ids)
     9259                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9260                      surfins_av(isurf) = surfins_av(isurf) + surfins(isurf)
     9261                  ENDIF
     9262              ENDDO
     9263
     9264          CASE ( 'rtm_rad_reslw' )
     9265!--           array of residua of lw radiation absorbed in surface after last reflection
     9266              DO isurf = dirstart(ids), dirend(ids)
     9267                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9268                      surfinl_av(isurf) = surfinl_av(isurf) + surfinl(isurf)
     9269                  ENDIF
     9270              ENDDO
     9271
     9272          CASE ( 'rtm_rad_pc_inlw' )
     9273              DO l = 1, npcbl
     9274                 pcbinlw_av(l) = pcbinlw_av(l) + pcbinlw(l)
     9275              ENDDO
     9276
     9277          CASE ( 'rtm_rad_pc_insw' )
     9278              DO l = 1, npcbl
     9279                 pcbinsw_av(l) = pcbinsw_av(l) + pcbinsw(l)
     9280              ENDDO
     9281
     9282          CASE ( 'rtm_rad_pc_inswdir' )
     9283              DO l = 1, npcbl
     9284                 pcbinswdir_av(l) = pcbinswdir_av(l) + pcbinswdir(l)
     9285              ENDDO
     9286
     9287          CASE ( 'rtm_rad_pc_inswdif' )
     9288              DO l = 1, npcbl
     9289                 pcbinswdif_av(l) = pcbinswdif_av(l) + pcbinswdif(l)
     9290              ENDDO
     9291
     9292          CASE ( 'rtm_rad_pc_inswref' )
     9293              DO l = 1, npcbl
     9294                 pcbinswref_av(l) = pcbinswref_av(l) + pcbinsw(l) - pcbinswdir(l) - pcbinswdif(l)
     9295              ENDDO
     9296
    89739297          CASE ( 'rad_mrt_sw' )
    89749298             IF ( ALLOCATED( mrtinsw_av ) )  THEN
     
    89939317    ELSEIF ( mode == 'average' )  THEN
    89949318
    8995        SELECT CASE ( TRIM( variable ) )
    8996 
     9319       SELECT CASE ( TRIM( var ) )
     9320!--       block of large scale (e.g. RRTMG) radiation output variables
    89979321          CASE ( 'rad_net*' )
    89989322             IF ( ALLOCATED( rad_net_av ) ) THEN
     
    91419465             ENDIF
    91429466
    9143           CASE ( 'rad_mrt_sw' )
    9144              IF ( ALLOCATED( mrtinsw_av ) )  THEN
    9145                 mrtinsw_av(:) = mrtinsw_av(:)  / REAL( average_count_3d, KIND=wp )
    9146              ENDIF
     9467!--       block of RTM output variables
     9468          CASE ( 'rtm_rad_net' )
     9469!--           array of complete radiation balance
     9470              DO isurf = dirstart(ids), dirend(ids)
     9471                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9472                      surfradnet_av(isurf) = surfinsw_av(isurf) / REAL( average_count_3d, kind=wp )
     9473                  ENDIF
     9474              ENDDO
     9475
     9476          CASE ( 'rtm_rad_insw' )
     9477!--           array of sw radiation falling to surface after i-th reflection
     9478              DO isurf = dirstart(ids), dirend(ids)
     9479                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9480                      surfinsw_av(isurf) = surfinsw_av(isurf) / REAL( average_count_3d, kind=wp )
     9481                  ENDIF
     9482              ENDDO
     9483
     9484          CASE ( 'rtm_rad_inlw' )
     9485!--           array of lw radiation falling to surface after i-th reflection
     9486              DO isurf = dirstart(ids), dirend(ids)
     9487                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9488                      surfinlw_av(isurf) = surfinlw_av(isurf) / REAL( average_count_3d, kind=wp )
     9489                  ENDIF
     9490              ENDDO
     9491
     9492          CASE ( 'rtm_rad_inswdir' )
     9493!--           array of direct sw radiation falling to surface from sun
     9494              DO isurf = dirstart(ids), dirend(ids)
     9495                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9496                      surfinswdir_av(isurf) = surfinswdir_av(isurf) / REAL( average_count_3d, kind=wp )
     9497                  ENDIF
     9498              ENDDO
     9499
     9500          CASE ( 'rtm_rad_inswdif' )
     9501!--           array of difusion sw radiation falling to surface from sky and borders of the domain
     9502              DO isurf = dirstart(ids), dirend(ids)
     9503                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9504                      surfinswdif_av(isurf) = surfinswdif_av(isurf) / REAL( average_count_3d, kind=wp )
     9505                  ENDIF
     9506              ENDDO
     9507
     9508          CASE ( 'rtm_rad_inswref' )
     9509!--           array of sw radiation falling to surface from reflections
     9510              DO isurf = dirstart(ids), dirend(ids)
     9511                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9512                      surfinswref_av(isurf) = surfinswref_av(isurf) / REAL( average_count_3d, kind=wp )
     9513                  ENDIF
     9514              ENDDO
     9515
     9516          CASE ( 'rtm_rad_inlwdif' )
     9517!--           array of sw radiation falling to surface after i-th reflection
     9518              DO isurf = dirstart(ids), dirend(ids)
     9519                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9520                      surfinlwdif_av(isurf) = surfinlwdif_av(isurf) / REAL( average_count_3d, kind=wp )
     9521                  ENDIF
     9522              ENDDO
     9523
     9524          CASE ( 'rtm_rad_inlwref' )
     9525!--           array of lw radiation falling to surface from reflections
     9526              DO isurf = dirstart(ids), dirend(ids)
     9527                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9528                      surfinlwref_av(isurf) = surfinlwref_av(isurf) / REAL( average_count_3d, kind=wp )
     9529                  ENDIF
     9530              ENDDO
     9531
     9532          CASE ( 'rtm_rad_outsw' )
     9533!--           array of sw radiation emitted from surface after i-th reflection
     9534              DO isurf = dirstart(ids), dirend(ids)
     9535                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9536                      surfoutsw_av(isurf) = surfoutsw_av(isurf) / REAL( average_count_3d, kind=wp )
     9537                  ENDIF
     9538              ENDDO
     9539
     9540          CASE ( 'rtm_rad_outlw' )
     9541!--           array of lw radiation emitted from surface after i-th reflection
     9542              DO isurf = dirstart(ids), dirend(ids)
     9543                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9544                      surfoutlw_av(isurf) = surfoutlw_av(isurf) / REAL( average_count_3d, kind=wp )
     9545                  ENDIF
     9546              ENDDO
     9547
     9548          CASE ( 'rtm_rad_ressw' )
     9549!--           array of residua of sw radiation absorbed in surface after last reflection
     9550              DO isurf = dirstart(ids), dirend(ids)
     9551                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9552                      surfins_av(isurf) = surfins_av(isurf) / REAL( average_count_3d, kind=wp )
     9553                  ENDIF
     9554              ENDDO
     9555
     9556          CASE ( 'rtm_rad_reslw' )
     9557!--           array of residua of lw radiation absorbed in surface after last reflection
     9558              DO isurf = dirstart(ids), dirend(ids)
     9559                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     9560                      surfinl_av(isurf) = surfinl_av(isurf) / REAL( average_count_3d, kind=wp )
     9561                  ENDIF
     9562              ENDDO
     9563
     9564          CASE ( 'rtm_rad_pc_inlw' )
     9565              DO l = 1, npcbl
     9566                 pcbinlw_av(:) = pcbinlw_av(:) / REAL( average_count_3d, kind=wp )
     9567              ENDDO
     9568
     9569          CASE ( 'rtm_rad_pc_insw' )
     9570              DO l = 1, npcbl
     9571                 pcbinsw_av(:) = pcbinsw_av(:) / REAL( average_count_3d, kind=wp )
     9572              ENDDO
     9573
     9574          CASE ( 'rtm_rad_pc_inswdir' )
     9575              DO l = 1, npcbl
     9576                 pcbinswdir_av(:) = pcbinswdir_av(:) / REAL( average_count_3d, kind=wp )
     9577              ENDDO
     9578
     9579          CASE ( 'rtm_rad_pc_inswdif' )
     9580              DO l = 1, npcbl
     9581                 pcbinswdif_av(:) = pcbinswdif_av(:) / REAL( average_count_3d, kind=wp )
     9582              ENDDO
     9583
     9584          CASE ( 'rtm_rad_pc_inswref' )
     9585              DO l = 1, npcbl
     9586                 pcbinswref_av(:) = pcbinswref_av(:) / REAL( average_count_3d, kind=wp )
     9587              ENDDO
    91479588
    91489589          CASE ( 'rad_mrt_lw' )
     
    91709611!> It is called out from subroutine netcdf.
    91719612!------------------------------------------------------------------------------!
    9172 SUBROUTINE radiation_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
     9613SUBROUTINE radiation_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
    91739614   
    91749615    IMPLICIT NONE
    91759616
    9176     CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
     9617    CHARACTER (LEN=*), INTENT(IN)  ::  variable    !<
    91779618    LOGICAL, INTENT(OUT)           ::  found       !<
    91789619    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
     
    91809621    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
    91819622
     9623    CHARACTER (len=varnamelength)  :: var
     9624
    91829625    found  = .TRUE.
    91839626
    91849627!
    91859628!-- Check for the grid
    9186     SELECT CASE ( TRIM( var ) )
    9187 
    9188        CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr',        &
    9189               'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy',            &
    9190               'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
    9191               'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
    9192               'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz',               &
    9193               'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw' )
    9194           grid_x = 'x'
    9195           grid_y = 'y'
    9196           grid_z = 'zu'
    9197 
    9198        CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',            &
    9199               'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
    9200               'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
    9201               'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
    9202           grid_x = 'x'
    9203           grid_y = 'y'
    9204           grid_z = 'zw'
    9205 
    9206 
    9207        CASE DEFAULT
    9208           found  = .FALSE.
    9209           grid_x = 'none'
    9210           grid_y = 'none'
    9211           grid_z = 'none'
    9212 
    9213         END SELECT
     9629    var = TRIM(variable)
     9630!-- RTM directional variables
     9631    IF ( var(1:12) == 'rtm_rad_net_'  .OR.  var(1:13) == 'rtm_rad_insw_'  .OR.          &
     9632         var(1:13) == 'rtm_rad_inlw_'  .OR.  var(1:16) == 'rtm_rad_inswdir_'  .OR.      &
     9633         var(1:16) == 'rtm_rad_inswdif_'  .OR.  var(1:16) == 'rtm_rad_inswref_'  .OR.   &
     9634         var(1:16) == 'rtm_rad_inlwdif_'  .OR.  var(1:16) == 'rtm_rad_inlwref_'  .OR.   &
     9635         var(1:14) == 'rtm_rad_outsw_'  .OR.  var(1:14) == 'rtm_rad_outlw_'  .OR.       &
     9636         var(1:14) == 'rtm_rad_ressw_'  .OR.  var(1:14) == 'rtm_rad_reslw_'  .OR.       &
     9637         var == 'rtm_rad_pc_inlw'  .OR.                                                 &
     9638         var == 'rtm_rad_pc_insw'  .OR.  var == 'rtm_rad_pc_inswdir'  .OR.              &
     9639         var == 'rtm_rad_pc_inswdif'  .OR.  var == 'rtm_rad_pc_inswref'  .OR.           &
     9640         var(1:7) == 'rtm_svf'  .OR.  var(1:7) == 'rtm_dif'  .OR.                       &
     9641         var(1:9) == 'rtm_skyvf' .OR. var(1:10) == 'rtm_skyvft'  .OR.                   &
     9642         var == 'rtm_mrt'  .OR.  var ==  'rtm_mrt_sw'  .OR.  var == 'rtm_mrt_lw' )  THEN
     9643
     9644         found = .TRUE.
     9645         grid_x = 'x'
     9646         grid_y = 'y'
     9647         grid_z = 'zu'
     9648    ELSE
     9649
     9650       SELECT CASE ( TRIM( var ) )
     9651
     9652          CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr',        &
     9653                 'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy',            &
     9654                 'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
     9655                 'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
     9656                 'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz',               &
     9657                 'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw' )
     9658             grid_x = 'x'
     9659             grid_y = 'y'
     9660             grid_z = 'zu'
     9661
     9662          CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',            &
     9663                 'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
     9664                 'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
     9665                 'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
     9666             grid_x = 'x'
     9667             grid_y = 'y'
     9668             grid_z = 'zw'
     9669
     9670
     9671          CASE DEFAULT
     9672             found  = .FALSE.
     9673             grid_x = 'none'
     9674             grid_y = 'none'
     9675             grid_z = 'none'
     9676
     9677           END SELECT
     9678       ENDIF
    92149679
    92159680    END SUBROUTINE radiation_define_netcdf_grid
     
    964910114    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
    965010115
     10116    CHARACTER (len=varnamelength)                   :: var, surfid
     10117    INTEGER(iwp)                                    :: ids,idsint_u,idsint_l,isurf,isvf,isurfs,isurflt,ipcgb
     10118    INTEGER(iwp)                                    :: is, js, ks, istat
     10119
    965110120    found = .TRUE.
    965210121
    9653 
    9654     SELECT CASE ( TRIM( variable ) )
    9655 
     10122    ids = -1
     10123    var = TRIM(variable)
     10124    DO i = 0, nd-1
     10125        k = len(TRIM(var))
     10126        j = len(TRIM(dirname(i)))
     10127        IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
     10128            ids = i
     10129            idsint_u = dirint_u(ids)
     10130            idsint_l = dirint_l(ids)
     10131            var = var(:k-j)
     10132            EXIT
     10133        ENDIF
     10134    ENDDO
     10135    IF ( ids == -1 )  THEN
     10136        var = TRIM(variable)
     10137    ENDIF
     10138
     10139    IF ( (var(1:8) == 'rtm_svf_'  .OR.  var(1:8) == 'rtm_dif_')  .AND.  len(TRIM(var)) >= 13 )  THEN
     10140!--     svf values to particular surface
     10141        surfid = var(9:)
     10142        i = index(surfid,'_')
     10143        j = index(surfid(i+1:),'_')
     10144        READ(surfid(1:i-1),*, iostat=istat ) is
     10145        IF ( istat == 0 )  THEN
     10146            READ(surfid(i+1:i+j-1),*, iostat=istat ) js
     10147        ENDIF
     10148        IF ( istat == 0 )  THEN
     10149            READ(surfid(i+j+1:),*, iostat=istat ) ks
     10150        ENDIF
     10151        IF ( istat == 0 )  THEN
     10152            var = var(1:7)
     10153        ENDIF
     10154    ENDIF
     10155
     10156    local_pf = fill_value
     10157
     10158    SELECT CASE ( TRIM( var ) )
     10159!--   block of large scale radiation model (e.g. RRTMG) output variables
    965610160      CASE ( 'rad_sw_in' )
    965710161         IF ( av == 0 )  THEN
     
    983810342         ENDIF
    983910343
    9840       CASE ( 'rad_mrt_sw' )
     10344!--   block of RTM output variables
     10345!--   variables are intended mainly for debugging and detailed analyse purposes
     10346      CASE ( 'rtm_skyvf' )
     10347!--        sky view factor
     10348         DO isurf = dirstart(ids), dirend(ids)
     10349            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10350               local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = skyvf(isurf)
     10351            ENDIF
     10352         ENDDO
     10353
     10354      CASE ( 'rtm_skyvft' )
     10355!--      sky view factor
     10356         DO isurf = dirstart(ids), dirend(ids)
     10357            IF ( surfl(id,isurf) == ids )  THEN
     10358               local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = skyvft(isurf)
     10359            ENDIF
     10360         ENDDO
     10361
     10362      CASE ( 'rtm_svf', 'rtm_dif' )
     10363!--      shape view factors or iradiance factors to selected surface
     10364         IF ( TRIM(var)=='rtm_svf' )  THEN
     10365             k = 1
     10366         ELSE
     10367             k = 2
     10368         ENDIF
     10369         DO isvf = 1, nsvfl
     10370            isurflt = svfsurf(1, isvf)
     10371            isurfs = svfsurf(2, isvf)
     10372
     10373            IF ( surf(ix,isurfs) == is  .AND.  surf(iy,isurfs) == js  .AND. surf(iz,isurfs) == ks  .AND. &
     10374                 (surf(id,isurfs) == idsint_u .OR. surfl(id,isurf) == idsint_l ) ) THEN
     10375!--            correct source surface
     10376               local_pf(surfl(ix,isurflt),surfl(iy,isurflt),surfl(iz,isurflt)) = svf(k,isvf)
     10377            ENDIF
     10378         ENDDO
     10379
     10380      CASE ( 'rtm_rad_net' )
     10381!--     array of complete radiation balance
     10382         DO isurf = dirstart(ids), dirend(ids)
     10383            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10384               IF ( av == 0 )  THEN
     10385                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = &
     10386                         surfinsw(isurf) - surfoutsw(isurf) +  surfinlw(isurf) - surfoutlw(isurf)
     10387               ELSE
     10388                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfradnet_av(isurf)
     10389               ENDIF
     10390            ENDIF
     10391         ENDDO
     10392
     10393      CASE ( 'rtm_rad_insw' )
     10394!--      array of sw radiation falling to surface after i-th reflection
     10395         DO isurf = dirstart(ids), dirend(ids)
     10396            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10397               IF ( av == 0 )  THEN
     10398                 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinsw(isurf)
     10399               ELSE
     10400                 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinsw_av(isurf)
     10401               ENDIF
     10402            ENDIF
     10403         ENDDO
     10404
     10405      CASE ( 'rtm_rad_inlw' )
     10406!--      array of lw radiation falling to surface after i-th reflection
     10407         DO isurf = dirstart(ids), dirend(ids)
     10408            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10409               IF ( av == 0 )  THEN
     10410                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw(isurf)
     10411               ELSE
     10412                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw_av(isurf)
     10413               ENDIF
     10414             ENDIF
     10415         ENDDO
     10416
     10417      CASE ( 'rtm_rad_inswdir' )
     10418!--      array of direct sw radiation falling to surface from sun
     10419         DO isurf = dirstart(ids), dirend(ids)
     10420            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10421               IF ( av == 0 )  THEN
     10422                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdir(isurf)
     10423               ELSE
     10424                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdir_av(isurf)
     10425               ENDIF
     10426            ENDIF
     10427         ENDDO
     10428
     10429      CASE ( 'rtm_rad_inswdif' )
     10430!--      array of difusion sw radiation falling to surface from sky and borders of the domain
     10431         DO isurf = dirstart(ids), dirend(ids)
     10432            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10433               IF ( av == 0 )  THEN
     10434                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdif(isurf)
     10435               ELSE
     10436                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdif_av(isurf)
     10437               ENDIF
     10438            ENDIF
     10439         ENDDO
     10440
     10441      CASE ( 'rtm_rad_inswref' )
     10442!--      array of sw radiation falling to surface from reflections
     10443         DO isurf = dirstart(ids), dirend(ids)
     10444            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10445               IF ( av == 0 )  THEN
     10446                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = &
     10447                    surfinsw(isurf) - surfinswdir(isurf) - surfinswdif(isurf)
     10448               ELSE
     10449                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswref_av(isurf)
     10450               ENDIF
     10451            ENDIF
     10452         ENDDO
     10453
     10454      CASE ( 'rtm_rad_inlwdif' )
     10455!--      array of difusion lw radiation falling to surface from sky and borders of the domain
     10456         DO isurf = dirstart(ids), dirend(ids)
     10457            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10458               IF ( av == 0 )  THEN
     10459                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlwdif(isurf)
     10460               ELSE
     10461                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlwdif_av(isurf)
     10462               ENDIF
     10463            ENDIF
     10464         ENDDO
     10465
     10466      CASE ( 'rtm_rad_inlwref' )
     10467!--      array of lw radiation falling to surface from reflections
     10468         DO isurf = dirstart(ids), dirend(ids)
     10469            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10470               IF ( av == 0 )  THEN
     10471                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
     10472               ELSE
     10473                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlwref_av(isurf)
     10474               ENDIF
     10475            ENDIF
     10476         ENDDO
     10477
     10478      CASE ( 'rtm_rad_outsw' )
     10479!--      array of sw radiation emitted from surface after i-th reflection
     10480         DO isurf = dirstart(ids), dirend(ids)
     10481            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10482               IF ( av == 0 )  THEN
     10483                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutsw(isurf)
     10484               ELSE
     10485                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutsw_av(isurf)
     10486               ENDIF
     10487            ENDIF
     10488         ENDDO
     10489
     10490      CASE ( 'rtm_rad_outlw' )
     10491!--      array of lw radiation emitted from surface after i-th reflection
     10492         DO isurf = dirstart(ids), dirend(ids)
     10493            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10494               IF ( av == 0 )  THEN
     10495                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutlw(isurf)
     10496               ELSE
     10497                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutlw_av(isurf)
     10498               ENDIF
     10499            ENDIF
     10500         ENDDO
     10501
     10502      CASE ( 'rtm_rad_ressw' )
     10503!--      average of array of residua of sw radiation absorbed in surface after last reflection
     10504         DO isurf = dirstart(ids), dirend(ids)
     10505            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10506               IF ( av == 0 )  THEN
     10507                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfins(isurf)
     10508               ELSE
     10509                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfins_av(isurf)
     10510               ENDIF
     10511            ENDIF
     10512         ENDDO
     10513
     10514      CASE ( 'rtm_rad_reslw' )
     10515!--      average of array of residua of lw radiation absorbed in surface after last reflection
     10516         DO isurf = dirstart(ids), dirend(ids)
     10517            IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10518               IF ( av == 0 )  THEN
     10519                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinl(isurf)
     10520               ELSE
     10521                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinl_av(isurf)
     10522               ENDIF
     10523            ENDIF
     10524         ENDDO
     10525
     10526      CASE ( 'rtm_rad_pc_inlw' )
     10527!--      array of lw radiation absorbed by plant canopy
     10528         DO ipcgb = 1, npcbl
     10529            IF ( av == 0 )  THEN
     10530               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinlw(ipcgb)
     10531            ELSE
     10532               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinlw_av(ipcgb)
     10533            ENDIF
     10534         ENDDO
     10535
     10536      CASE ( 'rtm_rad_pc_insw' )
     10537!--      array of sw radiation absorbed by plant canopy
     10538         DO ipcgb = 1, npcbl
     10539            IF ( av == 0 )  THEN
     10540              local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinsw(ipcgb)
     10541            ELSE
     10542              local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinsw_av(ipcgb)
     10543            ENDIF
     10544         ENDDO
     10545
     10546      CASE ( 'rtm_rad_pc_inswdir' )
     10547!--      array of direct sw radiation absorbed by plant canopy
     10548         DO ipcgb = 1, npcbl
     10549            IF ( av == 0 )  THEN
     10550               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinswdir(ipcgb)
     10551            ELSE
     10552               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinswdir_av(ipcgb)
     10553            ENDIF
     10554         ENDDO
     10555
     10556      CASE ( 'rtm_rad_pc_inswdif' )
     10557!--      array of diffuse sw radiation absorbed by plant canopy
     10558         DO ipcgb = 1, npcbl
     10559            IF ( av == 0 )  THEN
     10560               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinswdif(ipcgb)
     10561            ELSE
     10562               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinswdif_av(ipcgb)
     10563            ENDIF
     10564         ENDDO
     10565
     10566      CASE ( 'rtm_rad_pc_inswref' )
     10567!--      array of reflected sw radiation absorbed by plant canopy
     10568         DO ipcgb = 1, npcbl
     10569            IF ( av == 0 )  THEN
     10570               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = &
     10571                                    pcbinsw(ipcgb) - pcbinswdir(ipcgb) - pcbinswdif(ipcgb)
     10572            ELSE
     10573               local_pf(pcbl(ix,ipcgb),pcbl(iy,ipcgb),pcbl(iz,ipcgb)) = pcbinswref_av(ipcgb)
     10574            ENDIF
     10575         ENDDO
     10576
     10577      CASE ( 'rtm_mrt_sw' )
    984110578         local_pf = REAL( fill_value, KIND = wp )
    984210579         IF ( av == 0 )  THEN
    984310580            DO  l = 1, nmrtbl
    9844                i = mrtbl(ix,l)
    9845                j = mrtbl(iy,l)
    9846                k = mrtbl(iz,l)
    9847                local_pf(i,j,k) = mrtinsw(l)
     10581               local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrtinsw(l)
    984810582            ENDDO
    984910583         ELSE
    985010584            IF ( ALLOCATED( mrtinsw_av ) ) THEN
    985110585               DO  l = 1, nmrtbl
    9852                   i = mrtbl(ix,l)
    9853                   j = mrtbl(iy,l)
    9854                   k = mrtbl(iz,l)
    9855                   local_pf(i,j,k) = mrtinsw_av(l)
     10586                  local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrtinsw_av(l)
    985610587               ENDDO
    985710588            ENDIF
    985810589         ENDIF
    985910590
    9860       CASE ( 'rad_mrt_lw' )
     10591      CASE ( 'rtm_mrt_lw' )
    986110592         local_pf = REAL( fill_value, KIND = wp )
    986210593         IF ( av == 0 )  THEN
    986310594            DO  l = 1, nmrtbl
    9864                i = mrtbl(ix,l)
    9865                j = mrtbl(iy,l)
    9866                k = mrtbl(iz,l)
    9867                local_pf(i,j,k) = mrtinlw(l)
     10595               local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrtinlw(l)
    986810596            ENDDO
    986910597         ELSE
    987010598            IF ( ALLOCATED( mrtinlw_av ) ) THEN
    987110599               DO  l = 1, nmrtbl
    9872                   i = mrtbl(ix,l)
    9873                   j = mrtbl(iy,l)
    9874                   k = mrtbl(iz,l)
    9875                   local_pf(i,j,k) = mrtinlw_av(l)
     10600                  local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrtinlw_av(l)
    987610601               ENDDO
    987710602            ENDIF
    987810603         ENDIF
    987910604
    9880       CASE ( 'rad_mrt' )
     10605      CASE ( 'rtm_mrt' )
    988110606         local_pf = REAL( fill_value, KIND = wp )
    988210607         IF ( av == 0 )  THEN
    988310608            DO  l = 1, nmrtbl
    9884                i = mrtbl(ix,l)
    9885                j = mrtbl(iy,l)
    9886                k = mrtbl(iz,l)
    9887                local_pf(i,j,k) = mrt(l)
     10609               local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrt(l)
    988810610            ENDDO
    988910611         ELSE
    989010612            IF ( ALLOCATED( mrt_av ) ) THEN
    989110613               DO  l = 1, nmrtbl
    9892                   i = mrtbl(ix,l)
    9893                   j = mrtbl(iy,l)
    9894                   k = mrtbl(iz,l)
    9895                   local_pf(i,j,k) = mrt_av(l)
     10614                  local_pf(mrtbl(ix,l),mrtbl(iy,l),mrtbl(iz,l)) = mrt_av(l)
    989610615               ENDDO
    989710616            ENDIF
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r3597 r3607  
    2828! -----------------
    2929! $Id$
     30! Output of radiation-related quantities migrated to radiation_model_mod.
     31!
     32! 3597 2018-12-04 08:40:18Z maronga
    3033! Fixed calculation method of near surface air potential temperature at 10 cm
    3134! and moved to surface_layer_fluxes. Removed unnecessary _eb strings.
     
    438441               radiation, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,        &
    439442               sigma_sb, sun_direction, sun_dir_lat, sun_dir_lon,              &
    440                force_radiation_call, surfinsw, surfinlw, surfinswdir,          &
    441                surfinswdif, surfoutsw, surfoutlw, surfins,nsvfl, svf, svfsurf, &
    442                surfinl, surfinlwdif, rad_sw_in_dir, rad_sw_in_diff,            &
    443                rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, &
    444                surfl, nsurfl, pcbinsw, pcbinlw, pcbinswdir,                    &
    445                pcbinswdif, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l, &
    446                inorth_l, isouth_l, ieast_l, iwest_l, id,                       &
     443               force_radiation_call, iup_u, inorth_u, isouth_u, ieast_u,       &
     444               iwest_u, iup_l, inorth_l, isouth_l, ieast_l, iwest_l, id,       &
    447445               iz, iy, ix,  nsurf, idsvf, ndsvf,                               &
    448446               idcsf, ndcsf, kdcsf, pct,                                       &
    449                startland, endland, startwall, endwall, skyvf, skyvft, nzub,    &
    450                nzut, npcbl, pcbl, unscheduled_radiation_calls
     447               nzub, nzut, unscheduled_radiation_calls
    451448
    452449    USE statistics,                                                            &
     
    461458    IMPLICIT NONE
    462459
    463     !
     460!
    464461!-- USM model constants
    465462
    466463    REAL(wp), PARAMETER ::                     &
    467464              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
    468               lambda_h_green_dry       = 0.19_wp,    & ! heat conductivity for dry soil   
    469               lambda_h_green_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
     465              lambda_h_green_dry = 0.19_wp,    & ! heat conductivity for dry soil   
     466              lambda_h_green_sm  = 3.44_wp,    & ! heat conductivity of the soil matrix
    470467              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
    471468              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
     
    965962!-- arrays for time averages
    966963!-- Attention: the variable rad_net_av is also used in the 3d field variable in radiation_model_mod.f90. It may be better to rename it
    967     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw_av      !< average of sw radiation falling to local surface including radiation from reflections
    968     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw_av      !< average of lw radiation falling to local surface including radiation from reflections
    969     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir_av   !< average of direct sw radiation falling to local surface
    970     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif_av   !< average of diffuse sw radiation from sky and model boundary falling to local surface
    971     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif_av   !< average of diffuse lw radiation from sky and model boundary falling to local surface
    972     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswref_av   !< average of sw radiation falling to surface from reflections
    973     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwref_av   !< average of lw radiation falling to surface from reflections
    974     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    975     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    976     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
    977     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
    978     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw_av       !< Average of pcbinlw
    979     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw_av       !< Average of pcbinsw
    980     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir_av    !< Average of pcbinswdir
    981     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif_av    !< Average of pcbinswdif
    982     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswref_av    !< Average of pcbinswref
    983     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface 
    984964    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wghf_eb_av       !< average of wghf_eb
    985965    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wshf_eb_av       !< average of wshf_eb
     
    18331813           
    18341814           SELECT CASE ( TRIM( var ) )
    1835                
    1836                 CASE ( 'usm_rad_net' )
    1837 !--                 array of complete radiation balance
    1838                     IF ( l == -1 ) THEN
    1839                         IF ( .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
    1840                             ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
    1841                             surf_usm_h%rad_net_av = 0.0_wp
    1842                         ENDIF
    1843                     ELSE
    1844                         IF ( .NOT.  ALLOCATED(surf_usm_v(l)%rad_net_av) )  THEN
    1845                             ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
    1846                             surf_usm_v(l)%rad_net_av = 0.0_wp
    1847                         ENDIF
    1848                     ENDIF
    1849                    
    1850                 CASE ( 'usm_rad_insw' )
    1851 !--                 array of sw radiation falling to surface after i-th reflection
    1852                     IF ( .NOT.  ALLOCATED(surfinsw_av) )  THEN
    1853                         ALLOCATE( surfinsw_av(nsurfl) )
    1854                         surfinsw_av = 0.0_wp
    1855                     ENDIF
    1856 
    1857                 CASE ( 'usm_rad_inlw' )
    1858 !--                 array of lw radiation falling to surface after i-th reflection
    1859                     IF ( .NOT.  ALLOCATED(surfinlw_av) )  THEN
    1860                         ALLOCATE( surfinlw_av(nsurfl) )
    1861                         surfinlw_av = 0.0_wp
    1862                     ENDIF
    1863 
    1864                 CASE ( 'usm_rad_inswdir' )
    1865 !--                 array of direct sw radiation falling to surface from sun
    1866                     IF ( .NOT.  ALLOCATED(surfinswdir_av) )  THEN
    1867                         ALLOCATE( surfinswdir_av(nsurfl) )
    1868                         surfinswdir_av = 0.0_wp
    1869                     ENDIF
    1870 
    1871                 CASE ( 'usm_rad_inswdif' )
    1872 !--                 array of difusion sw radiation falling to surface from sky and borders of the domain
    1873                     IF ( .NOT.  ALLOCATED(surfinswdif_av) )  THEN
    1874                         ALLOCATE( surfinswdif_av(nsurfl) )
    1875                         surfinswdif_av = 0.0_wp
    1876                     ENDIF
    1877 
    1878                 CASE ( 'usm_rad_inswref' )
    1879 !--                 array of sw radiation falling to surface from reflections
    1880                     IF ( .NOT.  ALLOCATED(surfinswref_av) )  THEN
    1881                         ALLOCATE( surfinswref_av(nsurfl) )
    1882                         surfinswref_av = 0.0_wp
    1883                     ENDIF
    1884 
    1885                 CASE ( 'usm_rad_inlwdif' )
    1886 !--                 array of sw radiation falling to surface after i-th reflection
    1887                    IF ( .NOT.  ALLOCATED(surfinlwdif_av) )  THEN
    1888                         ALLOCATE( surfinlwdif_av(nsurfl) )
    1889                         surfinlwdif_av = 0.0_wp
    1890                     ENDIF
    1891 
    1892                 CASE ( 'usm_rad_inlwref' )
    1893 !--                 array of lw radiation falling to surface from reflections
    1894                     IF ( .NOT.  ALLOCATED(surfinlwref_av) )  THEN
    1895                         ALLOCATE( surfinlwref_av(nsurfl) )
    1896                         surfinlwref_av = 0.0_wp
    1897                     ENDIF
    1898 
    1899                 CASE ( 'usm_rad_outsw' )
    1900 !--                 array of sw radiation emitted from surface after i-th reflection
    1901                     IF ( .NOT.  ALLOCATED(surfoutsw_av) )  THEN
    1902                         ALLOCATE( surfoutsw_av(nsurfl) )
    1903                         surfoutsw_av = 0.0_wp
    1904                     ENDIF
    1905 
    1906                 CASE ( 'usm_rad_outlw' )
    1907 !--                 array of lw radiation emitted from surface after i-th reflection
    1908                     IF ( .NOT.  ALLOCATED(surfoutlw_av) )  THEN
    1909                         ALLOCATE( surfoutlw_av(nsurfl) )
    1910                         surfoutlw_av = 0.0_wp
    1911                     ENDIF
    1912                 CASE ( 'usm_rad_ressw' )
    1913 !--                 array of residua of sw radiation absorbed in surface after last reflection
    1914                     IF ( .NOT.  ALLOCATED(surfins_av) )  THEN
    1915                         ALLOCATE( surfins_av(nsurfl) )
    1916                         surfins_av = 0.0_wp
    1917                     ENDIF
    1918                                    
    1919                 CASE ( 'usm_rad_reslw' )
    1920 !--                 array of residua of lw radiation absorbed in surface after last reflection
    1921                     IF ( .NOT.  ALLOCATED(surfinl_av) )  THEN
    1922                         ALLOCATE( surfinl_av(nsurfl) )
    1923                         surfinl_av = 0.0_wp
    1924                     ENDIF
    1925                                    
    1926                 CASE ( 'usm_rad_pc_inlw' )
    1927 !--                 array of of lw radiation absorbed in plant canopy
    1928                     IF ( .NOT.  ALLOCATED(pcbinlw_av) )  THEN
    1929                         ALLOCATE( pcbinlw_av(1:npcbl) )
    1930                         pcbinlw_av = 0.0_wp
    1931                     ENDIF
    1932                                    
    1933                 CASE ( 'usm_rad_pc_insw' )
    1934 !--                 array of of sw radiation absorbed in plant canopy
    1935                     IF ( .NOT.  ALLOCATED(pcbinsw_av) )  THEN
    1936                         ALLOCATE( pcbinsw_av(1:npcbl) )
    1937                         pcbinsw_av = 0.0_wp
    1938                     ENDIF
    1939                                    
    1940                 CASE ( 'usm_rad_pc_inswdir' )
    1941 !--                 array of of direct sw radiation absorbed in plant canopy
    1942                     IF ( .NOT.  ALLOCATED(pcbinswdir_av) )  THEN
    1943                         ALLOCATE( pcbinswdir_av(1:npcbl) )
    1944                         pcbinswdir_av = 0.0_wp
    1945                     ENDIF
    1946                                    
    1947                 CASE ( 'usm_rad_pc_inswdif' )
    1948 !--                 array of of diffuse sw radiation absorbed in plant canopy
    1949                     IF ( .NOT.  ALLOCATED(pcbinswdif_av) )  THEN
    1950                         ALLOCATE( pcbinswdif_av(1:npcbl) )
    1951                         pcbinswdif_av = 0.0_wp
    1952                     ENDIF
    1953                                    
    1954                 CASE ( 'usm_rad_pc_inswref' )
    1955 !--                 array of of reflected sw radiation absorbed in plant canopy
    1956                     IF ( .NOT.  ALLOCATED(pcbinswref_av) )  THEN
    1957                         ALLOCATE( pcbinswref_av(1:npcbl) )
    1958                         pcbinswref_av = 0.0_wp
    1959                     ENDIF
    1960                                    
    1961                 CASE ( 'usm_rad_hf' )
    1962 !--                 array of heat flux from radiation for surfaces after i-th reflection
    1963                     IF ( l == -1 ) THEN
    1964                         IF ( .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
    1965                            ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
    1966                            surf_usm_h%surfhf_av = 0.0_wp
    1967                         ENDIF
    1968                     ELSE
    1969                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfhf_av) )  THEN
    1970                            ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) )
    1971                            surf_usm_v(l)%surfhf_av = 0.0_wp
    1972                        ENDIF
    1973                     ENDIF
    19741815
    19751816                CASE ( 'usm_wshf' )
     
    22182059           
    22192060           SELECT CASE ( TRIM( var ) )
    2220                
    2221                 CASE ( 'usm_rad_net' )
    2222 !--                 array of complete radiation balance
    2223                     IF ( l == -1 ) THEN
    2224                        DO  m = 1, surf_usm_h%ns
    2225                           surf_usm_h%rad_net_av(m) =                              &
    2226                                              surf_usm_h%rad_net_av(m) +           &
    2227                                              surf_usm_h%rad_net_l(m)
    2228                        ENDDO
    2229                     ELSE
    2230                        DO  m = 1, surf_usm_v(l)%ns
    2231                           surf_usm_v(l)%rad_net_av(m) =                        &
    2232                                           surf_usm_v(l)%rad_net_av(m) +        &
    2233                                           surf_usm_v(l)%rad_net_l(m)
    2234                        ENDDO
    2235                     ENDIF
    2236 
    2237                 CASE ( 'usm_rad_insw' )
    2238 !--                 array of sw radiation falling to surface after i-th reflection
    2239                     DO l = 1, nsurfl
    2240                         IF ( surfl(id,l) == idsint )  THEN
    2241                             surfinsw_av(l) = surfinsw_av(l) + surfinsw(l)
    2242                         ENDIF
    2243                     ENDDO
    2244                              
    2245                 CASE ( 'usm_rad_inlw' )
    2246 !--                 array of lw radiation falling to surface after i-th reflection
    2247                     DO l = 1, nsurfl
    2248                         IF ( surfl(id,l) == idsint )  THEN
    2249                             surfinlw_av(l) = surfinlw_av(l) + surfinlw(l)
    2250                         ENDIF
    2251                     ENDDO
    2252                    
    2253                 CASE ( 'usm_rad_inswdir' )
    2254 !--                 array of direct sw radiation falling to surface from sun
    2255                     DO l = 1, nsurfl
    2256                         IF ( surfl(id,l) == idsint )  THEN
    2257                             surfinswdir_av(l) = surfinswdir_av(l) + surfinswdir(l)
    2258                         ENDIF
    2259                     ENDDO
    2260                    
    2261                 CASE ( 'usm_rad_inswdif' )
    2262 !--                 array of difusion sw radiation falling to surface from sky and borders of the domain
    2263                     DO l = 1, nsurfl
    2264                         IF ( surfl(id,l) == idsint )  THEN
    2265                             surfinswdif_av(l) = surfinswdif_av(l) + surfinswdif(l)
    2266                         ENDIF
    2267                     ENDDO
    2268                    
    2269                 CASE ( 'usm_rad_inswref' )
    2270 !--                 array of sw radiation falling to surface from reflections
    2271                     DO l = 1, nsurfl
    2272                         IF ( surfl(id,l) == idsint )  THEN
    2273                             surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
    2274                                                 surfinswdir(l) - surfinswdif(l)
    2275                         ENDIF
    2276                     ENDDO
    2277 
    2278                    
    2279                 CASE ( 'usm_rad_inlwdif' )
    2280 !--                 array of sw radiation falling to surface after i-th reflection
    2281                     DO l = 1, nsurfl
    2282                         IF ( surfl(id,l) == idsint )  THEN
    2283                             surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l)
    2284                         ENDIF
    2285                     ENDDO
    2286 !                     
    2287                 CASE ( 'usm_rad_inlwref' )
    2288 !--                 array of lw radiation falling to surface from reflections
    2289                     DO l = 1, nsurfl
    2290                         IF ( surfl(id,l) == idsint )  THEN
    2291                             surfinlwref_av(l) = surfinlwref_av(l) + &
    2292                                                 surfinlw(l) - surfinlwdif(l)
    2293                         ENDIF
    2294                     ENDDO
    2295                    
    2296                 CASE ( 'usm_rad_outsw' )
    2297 !--                 array of sw radiation emitted from surface after i-th reflection
    2298                     DO l = 1, nsurfl
    2299                         IF ( surfl(id,l) == idsint )  THEN
    2300                             surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
    2301                         ENDIF
    2302                     ENDDO
    2303                    
    2304                 CASE ( 'usm_rad_outlw' )
    2305 !--                 array of lw radiation emitted from surface after i-th reflection
    2306                     DO l = 1, nsurfl
    2307                         IF ( surfl(id,l) == idsint )  THEN
    2308                             surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
    2309                         ENDIF
    2310                     ENDDO
    2311                    
    2312                 CASE ( 'usm_rad_ressw' )
    2313 !--                 array of residua of sw radiation absorbed in surface after last reflection
    2314                     DO l = 1, nsurfl
    2315                         IF ( surfl(id,l) == idsint )  THEN
    2316                             surfins_av(l) = surfins_av(l) + surfins(l)
    2317                         ENDIF
    2318                     ENDDO
    2319                                    
    2320                 CASE ( 'usm_rad_reslw' )
    2321 !--                 array of residua of lw radiation absorbed in surface after last reflection
    2322                     DO l = 1, nsurfl
    2323                         IF ( surfl(id,l) == idsint )  THEN
    2324                             surfinl_av(l) = surfinl_av(l) + surfinl(l)
    2325                         ENDIF
    2326                     ENDDO
    2327                    
    2328                 CASE ( 'usm_rad_pc_inlw' )
    2329                     pcbinlw_av(:) = pcbinlw_av(:) + pcbinlw(:)
    2330                    
    2331                 CASE ( 'usm_rad_pc_insw' )
    2332                     pcbinsw_av(:) = pcbinsw_av(:) + pcbinsw(:)
    2333                    
    2334                 CASE ( 'usm_rad_pc_inswdir' )
    2335                     pcbinswdir_av(:) = pcbinswdir_av(:) + pcbinswdir(:)
    2336                    
    2337                 CASE ( 'usm_rad_pc_inswdif' )
    2338                     pcbinswdif_av(:) = pcbinswdif_av(:) + pcbinswdif(:)
    2339                    
    2340                 CASE ( 'usm_rad_pc_inswref' )
    2341                     pcbinswref_av(:) = pcbinswref_av(:) + pcbinsw(:)     &
    2342                                                         - pcbinswdir(:)  &
    2343                                                         - pcbinswdif(:)
    2344                    
    2345                 CASE ( 'usm_rad_hf' )
    2346 !--                 array of heat flux from radiation for surfaces after i-th reflection
    2347                     IF ( l == -1 ) THEN
    2348                        DO  m = 1, surf_usm_h%ns
    2349                           surf_usm_h%surfhf_av(m) =                               &
    2350                                              surf_usm_h%surfhf_av(m) +            &
    2351                                              surf_usm_h%surfhf(m)
    2352                        ENDDO
    2353                     ELSE
    2354                        DO  m = 1, surf_usm_v(l)%ns
    2355                           surf_usm_v(l)%surfhf_av(m) =                         &
    2356                                           surf_usm_v(l)%surfhf_av(m) +         &
    2357                                           surf_usm_v(l)%surfhf(m)
    2358                        ENDDO
    2359                     ENDIF
    2360                    
     2061
    23612062                CASE ( 'usm_wshf' )
    23622063!--                 array of sensible heat flux from surfaces (land, roof, wall)
     
    26392340           
    26402341           SELECT CASE ( TRIM( var ) )
    2641                
    2642                 CASE ( 'usm_rad_net' )
    2643 !--                 array of complete radiation balance
    2644                     IF ( l == -1 ) THEN
    2645                        DO  m = 1, surf_usm_h%ns
    2646                           surf_usm_h%rad_net_av(m) =                              &
    2647                                              surf_usm_h%rad_net_av(m) /           &
    2648                                              REAL( average_count_3d, kind=wp )
    2649                        ENDDO
    2650                     ELSE
    2651                        DO  m = 1, surf_usm_v(l)%ns
    2652                           surf_usm_v(l)%rad_net_av(m) =                        &
    2653                                           surf_usm_v(l)%rad_net_av(m) /        &
    2654                                           REAL( average_count_3d, kind=wp )
    2655                        ENDDO
    2656                     ENDIF
    2657                    
    2658                 CASE ( 'usm_rad_insw' )
    2659 !--                 array of sw radiation falling to surface after i-th reflection
    2660                     DO l = 1, nsurfl
    2661                         IF ( surfl(id,l) == idsint )  THEN
    2662                             surfinsw_av(l) = surfinsw_av(l) / REAL( average_count_3d, kind=wp )
    2663                         ENDIF
    2664                     ENDDO
    2665                              
    2666                 CASE ( 'usm_rad_inlw' )
    2667 !--                 array of lw radiation falling to surface after i-th reflection
    2668                     DO l = 1, nsurfl
    2669                         IF ( surfl(id,l) == idsint )  THEN
    2670                             surfinlw_av(l) = surfinlw_av(l) / REAL( average_count_3d, kind=wp )
    2671                         ENDIF
    2672                     ENDDO
    2673                    
    2674                 CASE ( 'usm_rad_inswdir' )
    2675 !--                 array of direct sw radiation falling to surface from sun
    2676                     DO l = 1, nsurfl
    2677                         IF ( surfl(id,l) == idsint )  THEN
    2678                             surfinswdir_av(l) = surfinswdir_av(l) / REAL( average_count_3d, kind=wp )
    2679                         ENDIF
    2680                     ENDDO
    2681                    
    2682                 CASE ( 'usm_rad_inswdif' )
    2683 !--                 array of difusion sw radiation falling to surface from sky and borders of the domain
    2684                     DO l = 1, nsurfl
    2685                         IF ( surfl(id,l) == idsint )  THEN
    2686                             surfinswdif_av(l) = surfinswdif_av(l) / REAL( average_count_3d, kind=wp )
    2687                         ENDIF
    2688                     ENDDO
    2689                    
    2690                 CASE ( 'usm_rad_inswref' )
    2691 !--                 array of sw radiation falling to surface from reflections
    2692                     DO l = 1, nsurfl
    2693                         IF ( surfl(id,l) == idsint )  THEN
    2694                             surfinswref_av(l) = surfinswref_av(l) / REAL( average_count_3d, kind=wp )
    2695                         ENDIF
    2696                     ENDDO
    2697                    
    2698                 CASE ( 'usm_rad_inlwdif' )
    2699 !--                 array of sw radiation falling to surface after i-th reflection
    2700                     DO l = 1, nsurfl
    2701                         IF ( surfl(id,l) == idsint )  THEN
    2702                             surfinlwdif_av(l) = surfinlwdif_av(l) / REAL( average_count_3d, kind=wp )
    2703                         ENDIF
    2704                     ENDDO
    2705                    
    2706                 CASE ( 'usm_rad_inlwref' )
    2707 !--                 array of lw radiation falling to surface from reflections
    2708                     DO l = 1, nsurfl
    2709                         IF ( surfl(id,l) == idsint )  THEN
    2710                             surfinlwref_av(l) = surfinlwref_av(l) / REAL( average_count_3d, kind=wp )
    2711                         ENDIF
    2712                     ENDDO
    2713                    
    2714                 CASE ( 'usm_rad_outsw' )
    2715 !--                 array of sw radiation emitted from surface after i-th reflection
    2716                     DO l = 1, nsurfl
    2717                         IF ( surfl(id,l) == idsint )  THEN
    2718                             surfoutsw_av(l) = surfoutsw_av(l) / REAL( average_count_3d, kind=wp )
    2719                         ENDIF
    2720                     ENDDO
    2721                    
    2722                 CASE ( 'usm_rad_outlw' )
    2723 !--                 array of lw radiation emitted from surface after i-th reflection
    2724                     DO l = 1, nsurfl
    2725                         IF ( surfl(id,l) == idsint )  THEN
    2726                             surfoutlw_av(l) = surfoutlw_av(l) / REAL( average_count_3d, kind=wp )
    2727                         ENDIF
    2728                     ENDDO
    2729                    
    2730                 CASE ( 'usm_rad_ressw' )
    2731 !--                 array of residua of sw radiation absorbed in surface after last reflection
    2732                     DO l = 1, nsurfl
    2733                         IF ( surfl(id,l) == idsint )  THEN
    2734                             surfins_av(l) = surfins_av(l) / REAL( average_count_3d, kind=wp )
    2735                         ENDIF
    2736                     ENDDO
    2737                                    
    2738                 CASE ( 'usm_rad_reslw' )
    2739 !--                 array of residua of lw radiation absorbed in surface after last reflection
    2740                     DO l = 1, nsurfl
    2741                         IF ( surfl(id,l) == idsint )  THEN
    2742                             surfinl_av(l) = surfinl_av(l) / REAL( average_count_3d, kind=wp )
    2743                         ENDIF
    2744                     ENDDO
    2745                    
    2746                 CASE ( 'usm_rad_pc_inlw' )
    2747                     pcbinlw_av(:) = pcbinlw_av(:) / REAL( average_count_3d, kind=wp )
    2748                    
    2749                 CASE ( 'usm_rad_pc_insw' )
    2750                     pcbinsw_av(:) = pcbinsw_av(:) / REAL( average_count_3d, kind=wp )
    2751                    
    2752                 CASE ( 'usm_rad_pc_inswdir' )
    2753                     pcbinswdir_av(:) = pcbinswdir_av(:) / REAL( average_count_3d, kind=wp )
    2754                    
    2755                 CASE ( 'usm_rad_pc_inswdif' )
    2756                     pcbinswdif_av(:) = pcbinswdif_av(:) / REAL( average_count_3d, kind=wp )
    2757                    
    2758                 CASE ( 'usm_rad_pc_inswref' )
    2759                     pcbinswref_av(:) = pcbinswref_av(:) / REAL( average_count_3d, kind=wp )
    2760                    
    2761                 CASE ( 'usm_rad_hf' )
    2762 !--                 array of heat flux from radiation for surfaces after i-th reflection
    2763                     IF ( l == -1 ) THEN
    2764                        DO  m = 1, surf_usm_h%ns
    2765                           surf_usm_h%surfhf_av(m) =                               &
    2766                                              surf_usm_h%surfhf_av(m) /            &
    2767                                              REAL( average_count_3d, kind=wp )
    2768                        ENDDO
    2769                     ELSE
    2770                        DO  m = 1, surf_usm_v(l)%ns
    2771                           surf_usm_v(l)%surfhf_av(m) =                         &
    2772                                           surf_usm_v(l)%surfhf_av(m) /         &
    2773                                           REAL( average_count_3d, kind=wp )
    2774                        ENDDO
    2775                     ENDIF
    2776                    
     2342
    27772343                CASE ( 'usm_wshf' )
    27782344!--                 array of sensible heat flux from surfaces (land, roof, wall)
     
    31152681        CHARACTER(LEN=2)                              :: ls
    31162682        CHARACTER(LEN=varnamelength)                  :: var          !< TRIM(variable)
    3117         INTEGER(iwp), PARAMETER                       :: nl1 = 31     !< number of directional usm variables
     2683        INTEGER(iwp), PARAMETER                       :: nl1 = 16     !< number of directional usm variables
    31182684        CHARACTER(LEN=varnamelength), DIMENSION(nl1)  :: varlist1 = & !< list of directional usm variables
    3119                   (/'usm_rad_net                   ', &
    3120                     'usm_rad_insw                  ', &
    3121                     'usm_rad_inlw                  ', &
    3122                     'usm_rad_inswdir               ', &
    3123                     'usm_rad_inswdif               ', &
    3124                     'usm_rad_inswref               ', &
    3125                     'usm_rad_inlwdif               ', &
    3126                     'usm_wshf                      ', &
    3127                     'usm_rad_inlwref               ', &
    3128                     'usm_rad_outsw                 ', &
    3129                     'usm_rad_outlw                 ', &
    3130                     'usm_rad_hf                    ', &
    3131                     'usm_rad_ressw                 ', &
    3132                     'usm_rad_reslw                 ', &
     2685                  (/'usm_wshf                      ', &
    31332686                    'usm_wghf                      ', &
    31342687                    'usm_wghf_window               ', &
     
    31452698                    'usm_t_surf_green              ', &
    31462699                    'usm_t_green                   ', &
    3147                     'usm_theta_10cm                ', &
    3148                     'usm_skyvf                     ', &
    3149                     'usm_skyvft                    '/)
    3150 
    3151         INTEGER(iwp), PARAMETER                       :: nl2 = 7      !< number of other variables
    3152         CHARACTER(LEN=varnamelength), DIMENSION(nl2)  :: varlist2 = & !< list of other usm variables
    3153                   (/'usm_svf                       ', &
    3154                     'usm_dif                       ', &
    3155                     'usm_rad_pc_inlw               ', &
    3156                     'usm_rad_pc_insw               ', &
    3157                     'usm_rad_pc_inswdir            ', &
    3158                     'usm_rad_pc_inswdif            ', &
    3159                     'usm_rad_pc_inswref            '/)
    3160 
    3161         INTEGER(iwp), PARAMETER                       :: nl3 = 3      !< number of directional layer usm variables
    3162         CHARACTER(LEN=varnamelength), DIMENSION(nl3)  :: varlist3 = & !< list of directional layer usm variables
     2700                    'usm_theta_10cm                '/)
     2701
     2702        INTEGER(iwp), PARAMETER                       :: nl2 = 3      !< number of directional layer usm variables
     2703        CHARACTER(LEN=varnamelength), DIMENSION(nl2)  :: varlist2 = & !< list of directional layer usm variables
    31632704                  (/'usm_t_wall                    ', &
    31642705                    'usm_t_window                  ', &
     
    31882729        IF ( lfound ) GOTO 10
    31892730!       directional layer variables
    3190         DO i = 1, nl3
     2731        DO i = 1, nl2
    31912732           DO j = 1, nd
    31922733              DO l = nzb_wall, nzt_wall
    31932734                 WRITE(ls,'(A1,I1)') '_',l
    3194                  IF ( TRIM(var) == TRIM(varlist3(i))//TRIM(ls)//TRIM(dirname(j)) ) THEN
     2735                 IF ( TRIM(var) == TRIM(varlist2(i))//TRIM(ls)//TRIM(dirname(j)) ) THEN
    31952736                    lfound = .TRUE.
    31962737                    EXIT
     
    32002741           ENDDO
    32012742        ENDDO
    3202         IF ( lfound ) GOTO 10
    3203 !       other variables
    3204         DO i = 1, nl2
    3205            IF ( TRIM(var) == TRIM(varlist2(i)) ) THEN
    3206               lfound = .TRUE.
    3207               EXIT
    3208            ENDIF
    3209         ENDDO
    32102743        IF ( .NOT.  lfound ) THEN
    32112744           unit = 'illegal'
     
    3214274710      CONTINUE
    32152748
    3216         IF ( var(1:12) == 'usm_rad_net_'  .OR.  var(1:13) == 'usm_rad_insw_'  .OR.        &
    3217              var(1:13) == 'usm_rad_inlw_'  .OR.  var(1:16) == 'usm_rad_inswdir_'  .OR.    &
    3218              var(1:16) == 'usm_rad_inswdif_'  .OR.  var(1:16) == 'usm_rad_inswref_'  .OR. &
    3219              var(1:16) == 'usm_rad_inlwdif_'  .OR.  var(1:16) == 'usm_rad_inlwref_'  .OR. &
    3220              var(1:14) == 'usm_rad_outsw_'  .OR.  var(1:14) == 'usm_rad_outlw_'  .OR.     &
    3221              var(1:14) == 'usm_rad_ressw_'  .OR.  var(1:14) == 'usm_rad_reslw_'  .OR.     &
    3222              var(1:11) == 'usm_rad_hf_'  .OR.                                             &
    3223              var(1:9)  == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_' .OR.                 &
     2749        IF ( var(1:9)  == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_' .OR.                 &
    32242750             var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR.     &
    32252751             var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_'    .OR.      &
     
    32342760                  var(1:14) == 'usm_theta_10cm' )  THEN
    32352761            unit = 'K'
    3236         ELSE IF ( var == 'usm_rad_pc_inlw'  .OR.  var == 'usm_rad_pc_insw'  .OR.          &
    3237                   var == 'usm_rad_pc_inswdir'  .OR.  var == 'usm_rad_pc_inswdif'  .OR.    &
    3238                   var == 'usm_rad_pc_inswref' )  THEN
    3239             unit = 'W'
    3240         ELSE IF ( var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.              &
    3241                   var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.           &
    3242                   var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  .OR.     &
    3243                   var(1:9) == 'usm_skyvf' .OR. var(1:9) == 'usm_skyvft' )  THEN
     2762        ELSE IF ( var(1:9) == 'usm_surfz'  .OR.  var(1:11) == 'usm_surfcat'  .OR.         &
     2763                  var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  )  THEN
    32442764            unit = '1'
    32452765        ELSE
     
    33362856        INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER             :: diridx =  (/       -1,        1,        0,        3,        2 /)
    33372857                                                                     !< index for surf_*_v: 0:3 = (North, South, East, West)
    3338         INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirstart
    3339         INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirend
    33402858        INTEGER(iwp)                                           :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt,ipcgb
    33412859        INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat, l, m
    3342 
    3343         dirstart = (/ startland, startwall, startwall, startwall, startwall /)
    3344         dirend = (/ endland, endwall, endwall, endwall, endwall /)
    33452860
    33462861        found = .TRUE.
     
    33882903            READ(var(9:9), '(I1)', iostat=istat ) iwl
    33892904            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
    3390                 var = var(1:7)
    3391             ENDIF
    3392         ENDIF
    3393         IF ( (var(1:8) == 'usm_svf_'  .OR.  var(1:8) == 'usm_dif_')  .AND.  len(TRIM(var)) >= 13 )  THEN
    3394 !--         svf values to particular surface
    3395             surfid = var(9:)
    3396             i = index(surfid,'_')
    3397             j = index(surfid(i+1:),'_')
    3398             READ(surfid(1:i-1),*, iostat=istat ) is
    3399             IF ( istat == 0 )  THEN
    3400                 READ(surfid(i+1:i+j-1),*, iostat=istat ) js
    3401             ENDIF
    3402             IF ( istat == 0 )  THEN
    3403                 READ(surfid(i+j+1:),*, iostat=istat ) ks
    3404             ENDIF
    3405             IF ( istat == 0 )  THEN
    34062905                var = var(1:7)
    34072906            ENDIF
     
    35253024              ENDIF
    35263025
    3527           CASE ( 'usm_skyvf' )
    3528 !--           sky view factor
    3529               DO isurf = dirstart(ids), dirend(ids)
    3530                  IF ( surfl(id,isurf) == idsint )  THEN
    3531                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = skyvf(isurf)
    3532                  ENDIF
    3533               ENDDO
    3534              
    3535           CASE ( 'usm_skyvft' )
    3536 !--           sky view factor
    3537               DO isurf = dirstart(ids), dirend(ids)
    3538                  IF ( surfl(id,isurf) == ids )  THEN
    3539                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = skyvft(isurf)
    3540                  ENDIF
    3541               ENDDO
    3542 
    3543 !
    3544 !-- Not adjusted so far             
    3545           CASE ( 'usm_svf', 'usm_dif' )
    3546 !--           shape view factors or iradiance factors to selected surface
    3547               IF ( TRIM(var)=='usm_svf' )  THEN
    3548                   k = 1
    3549               ELSE
    3550                   k = 2
    3551               ENDIF
    3552               DO isvf = 1, nsvfl
    3553                   isurflt = svfsurf(1, isvf)
    3554                   isurfs = svfsurf(2, isvf)
    3555                              
    3556                   IF ( surf(ix,isurfs) == is  .AND.  surf(iy,isurfs) == js  .AND.       &
    3557                        surf(iz,isurfs) == ks  .AND.  surf(id,isurfs) == idsint )  THEN
    3558   !--                 correct source surface
    3559                       temp_pf(surfl(iz,isurflt),surfl(iy,isurflt),surfl(ix,isurflt)) = svf(k,isvf)
    3560                   ENDIF
    3561               ENDDO
    3562 
    3563           CASE ( 'usm_rad_net' )
    3564 !--           array of complete radiation balance
    3565               IF ( av == 0 )  THEN
    3566                  IF ( idsint == iup_u )  THEN
    3567                     DO  m = 1, surf_usm_h%ns
    3568                        i = surf_usm_h%i(m)
    3569                        j = surf_usm_h%j(m)
    3570                        k = surf_usm_h%k(m)
    3571                        temp_pf(k,j,i) = surf_usm_h%rad_net_l(m)
    3572                     ENDDO
    3573                  ELSE
    3574                     l = idsidx
    3575                     DO  m = 1, surf_usm_v(l)%ns
    3576                        i = surf_usm_v(l)%i(m)
    3577                        j = surf_usm_v(l)%j(m)
    3578                        k = surf_usm_v(l)%k(m)
    3579                        temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m)
    3580                     ENDDO
    3581                  ENDIF
    3582               ELSE
    3583                  IF ( idsint == iup_u )  THEN
    3584                     DO  m = 1, surf_usm_h%ns
    3585                        i = surf_usm_h%i(m)
    3586                        j = surf_usm_h%j(m)
    3587                        k = surf_usm_h%k(m)
    3588                        temp_pf(k,j,i) = surf_usm_h%rad_net_av(m)
    3589                     ENDDO
    3590                  ELSE
    3591                     l = idsidx
    3592                     DO  m = 1, surf_usm_v(l)%ns
    3593                        i = surf_usm_v(l)%i(m)
    3594                        j = surf_usm_v(l)%j(m)
    3595                        k = surf_usm_v(l)%k(m)
    3596                        temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m)
    3597                     ENDDO
    3598                  ENDIF
    3599               ENDIF
    3600 
    3601           CASE ( 'usm_rad_insw' )
    3602 !--           array of sw radiation falling to surface after i-th reflection
    3603               DO isurf = dirstart(ids), dirend(ids)
    3604                  IF ( surfl(id,isurf) == idsint )  THEN
    3605                    IF ( av == 0 )  THEN
    3606                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw(isurf)
    3607                    ELSE
    3608                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw_av(isurf)
    3609                    ENDIF
    3610                  ENDIF
    3611               ENDDO
    3612 
    3613           CASE ( 'usm_rad_inlw' )
    3614 !--           array of lw radiation falling to surface after i-th reflection
    3615               DO isurf = dirstart(ids), dirend(ids)
    3616                  IF ( surfl(id,isurf) == idsint )  THEN
    3617                    IF ( av == 0 )  THEN
    3618                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf)
    3619                    ELSE
    3620                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw_av(isurf)
    3621                    ENDIF
    3622                  ENDIF
    3623               ENDDO
    3624 
    3625           CASE ( 'usm_rad_inswdir' )
    3626 !--           array of direct sw radiation falling to surface from sun
    3627               DO isurf = dirstart(ids), dirend(ids)
    3628                  IF ( surfl(id,isurf) == idsint )  THEN
    3629                    IF ( av == 0 )  THEN
    3630                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir(isurf)
    3631                    ELSE
    3632                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir_av(isurf)
    3633                    ENDIF
    3634                  ENDIF
    3635               ENDDO
    3636 
    3637           CASE ( 'usm_rad_inswdif' )
    3638 !--           array of difusion sw radiation falling to surface from sky and borders of the domain
    3639               DO isurf = dirstart(ids), dirend(ids)
    3640                  IF ( surfl(id,isurf) == idsint )  THEN
    3641                    IF ( av == 0 )  THEN
    3642                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif(isurf)
    3643                    ELSE
    3644                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif_av(isurf)
    3645                    ENDIF
    3646                  ENDIF
    3647               ENDDO
    3648 
    3649           CASE ( 'usm_rad_inswref' )
    3650 !--           array of sw radiation falling to surface from reflections
    3651               DO isurf = dirstart(ids), dirend(ids)
    3652                  IF ( surfl(id,isurf) == idsint )  THEN
    3653                    IF ( av == 0 )  THEN
    3654                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = &
    3655                        surfinsw(isurf) - surfinswdir(isurf) - surfinswdif(isurf)
    3656                    ELSE
    3657                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswref_av(isurf)
    3658                    ENDIF
    3659                  ENDIF
    3660               ENDDO
    3661 
    3662           CASE ( 'usm_rad_inlwdif' )
    3663 !--           array of difusion lw radiation falling to surface from sky and borders of the domain
    3664               DO isurf = dirstart(ids), dirend(ids)
    3665                  IF ( surfl(id,isurf) == idsint )  THEN
    3666                    IF ( av == 0 )  THEN
    3667                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif(isurf)
    3668                    ELSE
    3669                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif_av(isurf)
    3670                    ENDIF
    3671                  ENDIF
    3672               ENDDO
    3673 
    3674           CASE ( 'usm_rad_inlwref' )
    3675 !--           array of lw radiation falling to surface from reflections
    3676               DO isurf = dirstart(ids), dirend(ids)
    3677                  IF ( surfl(id,isurf) == idsint )  THEN
    3678                    IF ( av == 0 )  THEN
    3679                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
    3680                    ELSE
    3681                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwref_av(isurf)
    3682                    ENDIF
    3683                  ENDIF
    3684               ENDDO
    3685 
    3686           CASE ( 'usm_rad_outsw' )
    3687 !--           array of sw radiation emitted from surface after i-th reflection
    3688               DO isurf = dirstart(ids), dirend(ids)
    3689                  IF ( surfl(id,isurf) == idsint )  THEN
    3690                    IF ( av == 0 )  THEN
    3691                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw(isurf)
    3692                    ELSE
    3693                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw_av(isurf)
    3694                    ENDIF
    3695                  ENDIF
    3696               ENDDO
    3697 
    3698           CASE ( 'usm_rad_outlw' )
    3699 !--           array of lw radiation emitted from surface after i-th reflection
    3700               DO isurf = dirstart(ids), dirend(ids)
    3701                  IF ( surfl(id,isurf) == idsint )  THEN
    3702                    IF ( av == 0 )  THEN
    3703                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw(isurf)
    3704                    ELSE
    3705                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw_av(isurf)
    3706                    ENDIF
    3707                  ENDIF
    3708               ENDDO
    3709 
    3710           CASE ( 'usm_rad_ressw' )
    3711 !--           average of array of residua of sw radiation absorbed in surface after last reflection
    3712               DO isurf = dirstart(ids), dirend(ids)
    3713                  IF ( surfl(id,isurf) == idsint )  THEN
    3714                    IF ( av == 0 )  THEN
    3715                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins(isurf)
    3716                    ELSE
    3717                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins_av(isurf)
    3718                    ENDIF
    3719                  ENDIF
    3720               ENDDO
    3721 
    3722           CASE ( 'usm_rad_reslw' )
    3723 !--           average of array of residua of lw radiation absorbed in surface after last reflection
    3724               DO isurf = dirstart(ids), dirend(ids)
    3725                  IF ( surfl(id,isurf) == idsint )  THEN
    3726                    IF ( av == 0 )  THEN
    3727                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl(isurf)
    3728                    ELSE
    3729                      temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl_av(isurf)
    3730                    ENDIF
    3731                  ENDIF
    3732               ENDDO
    3733 
    3734           CASE ( 'usm_rad_pc_inlw' )
    3735 !--           array of lw radiation absorbed by plant canopy
    3736               DO ipcgb = 1, npcbl
    3737                   IF ( av == 0 )  THEN
    3738                       temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinlw(ipcgb)
    3739                   ELSE
    3740                       temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinlw_av(ipcgb)
    3741                   ENDIF
    3742               ENDDO
    3743 
    3744           CASE ( 'usm_rad_pc_insw' )
    3745 !--           array of sw radiation absorbed by plant canopy
    3746               DO ipcgb = 1, npcbl
    3747                   IF ( av == 0 )  THEN
    3748                       temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw(ipcgb)
    3749                   ELSE
    3750                       temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw_av(ipcgb)
    3751                   ENDIF
    3752               ENDDO
    3753 
    3754           CASE ( 'usm_rad_pc_inswdir' )
    3755 !--           array of direct sw radiation absorbed by plant canopy
    3756               DO ipcgb = 1, npcbl
    3757                   IF ( av == 0 )  THEN
    3758                       temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdir(ipcgb)
    3759                   ELSE
    3760                       temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdir_av(ipcgb)
    3761                   ENDIF
    3762               ENDDO
    3763 
    3764           CASE ( 'usm_rad_pc_inswdif' )
    3765 !--           array of diffuse sw radiation absorbed by plant canopy
    3766               DO ipcgb = 1, npcbl
    3767                   IF ( av == 0 )  THEN
    3768                       temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdif(ipcgb)
    3769                   ELSE
    3770                       temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdif_av(ipcgb)
    3771                   ENDIF
    3772               ENDDO
    3773 
    3774           CASE ( 'usm_rad_pc_inswref' )
    3775 !--           array of reflected sw radiation absorbed by plant canopy
    3776               DO ipcgb = 1, npcbl
    3777                   IF ( av == 0 )  THEN
    3778                       temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw(ipcgb)      &
    3779                                                                               - pcbinswdir(ipcgb) &
    3780                                                                               - pcbinswdif(ipcgb)
    3781                   ELSE
    3782                       temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswref_av(ipcgb)
    3783                   ENDIF
    3784               ENDDO
    3785  
    3786           CASE ( 'usm_rad_hf' )
    3787 !--           array of heat flux from radiation for surfaces after all reflections
    3788               IF ( av == 0 )  THEN
    3789                  IF ( idsint == iup_u )  THEN
    3790                     DO  m = 1, surf_usm_h%ns
    3791                        i = surf_usm_h%i(m)
    3792                        j = surf_usm_h%j(m)
    3793                        k = surf_usm_h%k(m)
    3794                        temp_pf(k,j,i) = surf_usm_h%surfhf(m)
    3795                     ENDDO
    3796                  ELSE
    3797                     l = idsidx
    3798                     DO  m = 1, surf_usm_v(l)%ns
    3799                        i = surf_usm_v(l)%i(m)
    3800                        j = surf_usm_v(l)%j(m)
    3801                        k = surf_usm_v(l)%k(m)
    3802                        temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m)
    3803                     ENDDO
    3804                  ENDIF
    3805               ELSE
    3806                  IF ( idsint == iup_u )  THEN
    3807                     DO  m = 1, surf_usm_h%ns
    3808                        i = surf_usm_h%i(m)
    3809                        j = surf_usm_h%j(m)
    3810                        k = surf_usm_h%k(m)
    3811                        temp_pf(k,j,i) = surf_usm_h%surfhf_av(m)
    3812                     ENDDO
    3813                  ELSE
    3814                     l = idsidx
    3815                     DO  m = 1, surf_usm_v(l)%ns
    3816                        i = surf_usm_v(l)%i(m)
    3817                        j = surf_usm_v(l)%j(m)
    3818                        k = surf_usm_v(l)%k(m)
    3819                        temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m)
    3820                     ENDDO
    3821                  ENDIF
    3822               ENDIF
    3823  
    38243026          CASE ( 'usm_wshf' )
    38253027!--           array of sensible heat flux from surfaces
     
    39753177              ENDIF
    39763178
    3977 
    3978 
    3979 
    39803179          CASE ( 'usm_wghf' )
    39813180!--           array of heat flux from ground (land, wall, roof)
     
    45273726
    45283727        var = TRIM(variable)
    4529         IF ( var(1:12) == 'usm_rad_net_'  .OR.  var(1:13) == 'usm_rad_insw_'  .OR.          &
    4530              var(1:13) == 'usm_rad_inlw_'  .OR.  var(1:16) == 'usm_rad_inswdir_'  .OR.      &
    4531              var(1:16) == 'usm_rad_inswdif_'  .OR.  var(1:16) == 'usm_rad_inswref_'  .OR.   &
    4532              var(1:16) == 'usm_rad_inlwdif_'  .OR.  var(1:16) == 'usm_rad_inlwref_'  .OR.   &
    4533              var(1:14) == 'usm_rad_outsw_'  .OR.  var(1:14) == 'usm_rad_outlw_'  .OR.       &
    4534              var(1:14) == 'usm_rad_ressw_'  .OR.  var(1:14) == 'usm_rad_reslw_'  .OR.       &
    4535              var(1:11) == 'usm_rad_hf_'  .OR.  var == 'usm_rad_pc_inlw'  .OR.               &
    4536              var == 'usm_rad_pc_insw'  .OR.  var == 'usm_rad_pc_inswdir'  .OR.              &
    4537              var == 'usm_rad_pc_inswdif'  .OR.  var == 'usm_rad_pc_inswref'  .OR.           &
    4538              var(1:9) == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_'  .OR.                   &
     3728        IF ( var(1:9) == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_'  .OR.                   &
    45393729             var(1:16) == 'usm_wghf_window_'  .OR. var(1:15) == 'usm_wghf_green_' .OR.      &
    45403730             var(1:10) == 'usm_iwghf_'  .OR. var(1:17) == 'usm_iwghf_window_' .OR.          &
     
    45453735             var(1:16) == 'usm_t_surf_green'  .OR. var(1:11) == 'usm_t_green' .OR.          &
    45463736             var(1:15) == 'usm_theta_10cm' .OR.                                             &
    4547              var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.                     &
    4548              var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.                  &
     3737             var(1:9) == 'usm_surfz'  .OR.  var(1:11) == 'usm_surfcat'  .OR.                  &
    45493738             var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  .OR.            &
    4550              var(1:16) == 'usm_surfwintrans'  .OR. var(1:7) == 'usm_swc' .OR.               &
    4551              var(1:9) == 'usm_skyvf' .OR. var(1:9) == 'usm_skyvft' ) THEN
     3739             var(1:16) == 'usm_surfwintrans'  .OR. var(1:7) == 'usm_swc' ) THEN
    45523740
    45533741            found = .TRUE.
     
    101769364          CALL wrd_write_string( 't_green_v(' // dum // ')' )
    101779365          WRITE ( 14 )  t_green_v(l)%t
    10178 !        
     9366       
    101799367       ENDDO
    101809368
Note: See TracChangeset for help on using the changeset viewer.