Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3115)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3116)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Output of long/shortwave radiation at surface
+! 
+! 3107 2018-07-06 15:55:51Z suehring
 ! Bugfix, missing index for dz
 ! 
@@ -521,5 +524,9 @@
                                  sun_dir_lon       !< solar directional vector in longitudes
 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_net_av   !< average of rad_net
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_net_av       !< average of net radiation (rad_net) at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_in_xy_av  !< average of incoming longwave radiation at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_lw_out_xy_av !< average of outgoing longwave radiation at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_xy_av  !< average of incoming shortwave radiation at surface
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_out_xy_av !< average of outgoing shortwave radiation at surface
 !
 !-- Land surface albedos for solar zenith angle of 60Â° after Briegleb (1992)     
@@ -1089,5 +1096,6 @@
 
           CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
-                 'rrtm_asdir*' )
+                 'rrtm_asdir*', 'rad_lw_in*', 'rad_lw_out*', 'rad_sw_in*',     &
+                 'rad_sw_out*')
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
                 message_string = 'illegal value for data_output: "' //         &
@@ -1110,4 +1118,8 @@
 
              IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
+             IF ( TRIM( var ) == 'rad_lw_in*'    ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'rad_lw_out*'   ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'rad_sw_in*'    ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'rad_sw_out*'   ) unit = 'W/m2'
              IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
              IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = '' 
@@ -7362,4 +7374,28 @@
                 ENDIF
                 rad_net_av = 0.0_wp
+             
+             CASE ( 'rad_lw_in*' )
+                IF ( .NOT. ALLOCATED( rad_lw_in_xy_av ) )  THEN
+                   ALLOCATE( rad_lw_in_xy_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_lw_in_xy_av = 0.0_wp
+                
+             CASE ( 'rad_lw_out*' )
+                IF ( .NOT. ALLOCATED( rad_lw_out_xy_av ) )  THEN
+                   ALLOCATE( rad_lw_out_xy_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_lw_out_xy_av = 0.0_wp
+                
+             CASE ( 'rad_sw_in*' )
+                IF ( .NOT. ALLOCATED( rad_sw_in_xy_av ) )  THEN
+                   ALLOCATE( rad_sw_in_xy_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_sw_in_xy_av = 0.0_wp
+                
+             CASE ( 'rad_sw_out*' )
+                IF ( .NOT. ALLOCATED( rad_sw_out_xy_av ) )  THEN
+                   ALLOCATE( rad_sw_out_xy_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_sw_out_xy_av = 0.0_wp                
 
              CASE ( 'rad_lw_in' )
@@ -7424,9 +7460,13 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      DO m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
-                         rad_net_av(j,i) = rad_net_av(j,i) + surf_lsm_h%rad_net(m)
+                      DO m = surf_lsm_h%start_index(j,i),                      &
+                             surf_lsm_h%end_index(j,i)
+                         rad_net_av(j,i) = rad_net_av(j,i) +                   &
+                                           surf_lsm_h%rad_net(m)
                       ENDDO
-                      DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
-                         rad_net_av(j,i) = rad_net_av(j,i) + surf_usm_h%rad_net(m)
+                      DO m = surf_usm_h%start_index(j,i),                      &
+                             surf_usm_h%end_index(j,i)
+                         rad_net_av(j,i) = rad_net_av(j,i) +                   &
+                                           surf_usm_h%rad_net(m)
                       ENDDO
                    ENDDO
@@ -7434,4 +7474,76 @@
              ENDIF
 
+          CASE ( 'rad_lw_in*' )
+             IF ( ALLOCATED( rad_lw_in_xy_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO m = surf_lsm_h%start_index(j,i),                      &
+                             surf_lsm_h%end_index(j,i)
+                         rad_lw_in_xy_av(j,i) = rad_lw_in_xy_av(j,i) +         &
+                                                surf_lsm_h%rad_lw_in(m)
+                      ENDDO
+                      DO m = surf_usm_h%start_index(j,i),                      &
+                             surf_usm_h%end_index(j,i)
+                         rad_lw_in_xy_av(j,i) = rad_lw_in_xy_av(j,i) +         &
+                                                surf_usm_h%rad_lw_in(m)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'rad_lw_out*' )
+             IF ( ALLOCATED( rad_lw_out_xy_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO m = surf_lsm_h%start_index(j,i),                      &
+                             surf_lsm_h%end_index(j,i)
+                         rad_lw_out_xy_av(j,i) = rad_lw_out_xy_av(j,i) +       &
+                                                   surf_lsm_h%rad_lw_out(m)
+                      ENDDO
+                      DO m = surf_usm_h%start_index(j,i),                      &
+                             surf_usm_h%end_index(j,i)
+                         rad_lw_out_xy_av(j,i) = rad_lw_out_xy_av(j,i) +       &
+                                                  surf_usm_h%rad_lw_out(m)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'rad_sw_in*' )
+             IF ( ALLOCATED( rad_sw_in_xy_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO m = surf_lsm_h%start_index(j,i),                      &
+                             surf_lsm_h%end_index(j,i)
+                         rad_sw_in_xy_av(j,i) = rad_sw_in_xy_av(j,i) +         &
+                                                  surf_lsm_h%rad_sw_in(m)
+                      ENDDO
+                      DO m = surf_usm_h%start_index(j,i),                      &
+                             surf_usm_h%end_index(j,i)
+                         rad_sw_in_xy_av(j,i) = rad_sw_in_xy_av(j,i) +         &
+                                                surf_usm_h%rad_sw_in(m)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'rad_sw_out*' )
+             IF ( ALLOCATED( rad_sw_out_xy_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO m = surf_lsm_h%start_index(j,i),                      &
+                             surf_lsm_h%end_index(j,i)
+                         rad_sw_out_xy_av(j,i) = rad_sw_out_xy_av(j,i) +       &
+                                                   surf_lsm_h%rad_sw_out(m)
+                      ENDDO
+                      DO m = surf_usm_h%start_index(j,i),                      &
+                             surf_usm_h%end_index(j,i)
+                         rad_sw_out_xy_av(j,i) = rad_sw_out_xy_av(j,i) +       &
+                                                   surf_usm_h%rad_sw_out(m)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             
           CASE ( 'rad_lw_in' )
              IF ( ALLOCATED( rad_lw_in_av ) ) THEN
@@ -7548,4 +7660,44 @@
                 ENDDO
              ENDIF
+             
+          CASE ( 'rad_lw_in*' )
+             IF ( ALLOCATED( rad_lw_in_xy_av ) ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      rad_lw_in_xy_av(j,i) = rad_lw_in_xy_av(j,i)              &
+                                        / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'rad_lw_out*' )
+             IF ( ALLOCATED( rad_lw_out_xy_av ) ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      rad_lw_out_xy_av(j,i) = rad_lw_out_xy_av(j,i)            &
+                                        / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'rad_sw_in*' )
+             IF ( ALLOCATED( rad_sw_in_xy_av ) ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      rad_sw_in_xy_av(j,i) = rad_sw_in_xy_av(j,i)              &
+                                        / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDIF
+             
+          CASE ( 'rad_sw_out*' )
+             IF ( ALLOCATED( rad_sw_out_xy_av ) ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      rad_sw_out_xy_av(j,i) = rad_sw_out_xy_av(j,i)             &
+                                        / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDIF
 
           CASE ( 'rad_lw_in' )
@@ -7763,8 +7915,8 @@
              ENDDO
           ELSE
-            IF ( .NOT. ALLOCATED( rad_net_av ) ) THEN
-               ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
-               rad_net_av = REAL( fill_value, KIND = wp )
-            ENDIF
+             IF ( .NOT. ALLOCATED( rad_net_av ) ) THEN
+                ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
+                rad_net_av = REAL( fill_value, KIND = wp )
+             ENDIF
              DO  i = nxl, nxr
                 DO  j = nys, nyn 
@@ -7775,6 +7927,137 @@
           two_d = .TRUE.
           grid = 'zu1'
-
- 
+          
+       CASE ( 'rad_lw_in*_xy' )        ! 2d-array
+          IF ( av == 0 ) THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+!
+!--                Obtain rad_net from its respective surface type
+!--                Natural-type surfaces
+                   DO  m = surf_lsm_h%start_index(j,i),                        &
+                           surf_lsm_h%end_index(j,i)  
+                      local_pf(i,j,nzb+1) = surf_lsm_h%rad_lw_in(m)
+                   ENDDO
+!
+!--                Urban-type surfaces
+                   DO  m = surf_usm_h%start_index(j,i),                        &
+                           surf_usm_h%end_index(j,i)  
+                      local_pf(i,j,nzb+1) = surf_usm_h%rad_lw_in(m)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             IF ( .NOT. ALLOCATED( rad_lw_in_xy_av ) ) THEN
+                ALLOCATE( rad_lw_in_xy_av(nysg:nyng,nxlg:nxrg) )
+                rad_lw_in_xy_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             DO  i = nxl, nxr
+                DO  j = nys, nyn 
+                   local_pf(i,j,nzb+1) = rad_lw_in_xy_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+          two_d = .TRUE.
+          grid = 'zu1'
+          
+       CASE ( 'rad_lw_out*_xy' )        ! 2d-array
+          IF ( av == 0 ) THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+!
+!--                Obtain rad_net from its respective surface type
+!--                Natural-type surfaces
+                   DO  m = surf_lsm_h%start_index(j,i),                        &
+                           surf_lsm_h%end_index(j,i)  
+                      local_pf(i,j,nzb+1) = surf_lsm_h%rad_lw_out(m)
+                   ENDDO
+!
+!--                Urban-type surfaces
+                   DO  m = surf_usm_h%start_index(j,i),                        &
+                           surf_usm_h%end_index(j,i)  
+                      local_pf(i,j,nzb+1) = surf_usm_h%rad_lw_out(m)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             IF ( .NOT. ALLOCATED( rad_lw_out_xy_av ) ) THEN
+                ALLOCATE( rad_lw_out_xy_av(nysg:nyng,nxlg:nxrg) )
+                rad_lw_out_xy_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             DO  i = nxl, nxr
+                DO  j = nys, nyn 
+                   local_pf(i,j,nzb+1) = rad_lw_out_xy_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+          two_d = .TRUE.
+          grid = 'zu1'
+          
+       CASE ( 'rad_sw_in*_xy' )        ! 2d-array
+          IF ( av == 0 ) THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+!
+!--                Obtain rad_net from its respective surface type
+!--                Natural-type surfaces
+                   DO  m = surf_lsm_h%start_index(j,i),                        &
+                           surf_lsm_h%end_index(j,i)  
+                      local_pf(i,j,nzb+1) = surf_lsm_h%rad_sw_in(m)
+                   ENDDO
+!
+!--                Urban-type surfaces
+                   DO  m = surf_usm_h%start_index(j,i),                        &
+                           surf_usm_h%end_index(j,i)  
+                      local_pf(i,j,nzb+1) = surf_usm_h%rad_sw_in(m)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             IF ( .NOT. ALLOCATED( rad_sw_in_xy_av ) ) THEN
+                ALLOCATE( rad_sw_in_xy_av(nysg:nyng,nxlg:nxrg) )
+                rad_sw_in_xy_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             DO  i = nxl, nxr
+                DO  j = nys, nyn 
+                   local_pf(i,j,nzb+1) = rad_sw_in_xy_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+          two_d = .TRUE.
+          grid = 'zu1'
+          
+       CASE ( 'rad_sw_out*_xy' )        ! 2d-array
+          IF ( av == 0 ) THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+!
+!--                Obtain rad_net from its respective surface type
+!--                Natural-type surfaces
+                   DO  m = surf_lsm_h%start_index(j,i),                        &
+                           surf_lsm_h%end_index(j,i)  
+                      local_pf(i,j,nzb+1) = surf_lsm_h%rad_sw_out(m)
+                   ENDDO
+!
+!--                Urban-type surfaces
+                   DO  m = surf_usm_h%start_index(j,i),                        &
+                           surf_usm_h%end_index(j,i)  
+                      local_pf(i,j,nzb+1) = surf_usm_h%rad_sw_out(m)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             IF ( .NOT. ALLOCATED( rad_sw_out_xy_av ) ) THEN
+                ALLOCATE( rad_sw_out_xy_av(nysg:nyng,nxlg:nxrg) )
+                rad_sw_out_xy_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             DO  i = nxl, nxr
+                DO  j = nys, nyn 
+                   local_pf(i,j,nzb+1) = rad_sw_out_xy_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+          two_d = .TRUE.
+          grid = 'zu1'          
+          
        CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
           IF ( av == 0 ) THEN
@@ -8435,4 +8718,24 @@
        WRITE ( 14 )  rad_net_av
     ENDIF
+    
+    IF ( ALLOCATED( rad_lw_in_xy_av ) )  THEN
+       CALL wrd_write_string( 'rad_lw_in_xy_av' )
+       WRITE ( 14 )  rad_lw_in_xy_av
+    ENDIF
+    
+    IF ( ALLOCATED( rad_lw_out_xy_av ) )  THEN
+       CALL wrd_write_string( 'rad_lw_out_xy_av' )
+       WRITE ( 14 )  rad_lw_out_xy_av
+    ENDIF
+    
+    IF ( ALLOCATED( rad_sw_in_xy_av ) )  THEN
+       CALL wrd_write_string( 'rad_sw_in_xy_av' )
+       WRITE ( 14 )  rad_sw_in_xy_av
+    ENDIF
+    
+    IF ( ALLOCATED( rad_sw_out_xy_av ) )  THEN
+       CALL wrd_write_string( 'rad_sw_out_xy_av' )
+       WRITE ( 14 )  rad_sw_out_xy_av
+    ENDIF
 
     IF ( ALLOCATED( rad_lw_in ) )  THEN
@@ -8576,4 +8879,37 @@
           rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =           &
                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                        
+       CASE ( 'rad_lw_in_xy_av' )
+          IF ( .NOT. ALLOCATED( rad_lw_in_xy_av ) )  THEN
+             ALLOCATE( rad_lw_in_xy_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF  
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
+          rad_lw_in_xy_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =      &
+                        tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                        
+       CASE ( 'rad_lw_out_xy_av' )
+          IF ( .NOT. ALLOCATED( rad_lw_out_xy_av ) )  THEN
+             ALLOCATE( rad_lw_out_xy_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF  
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
+          rad_lw_out_xy_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =     &
+                        tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                        
+       CASE ( 'rad_sw_in_xy_av' )
+          IF ( .NOT. ALLOCATED( rad_sw_in_xy_av ) )  THEN
+             ALLOCATE( rad_sw_in_xy_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF  
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
+          rad_sw_in_xy_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =      &
+                        tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                        
+       CASE ( 'rad_sw_out_xy_av' )
+          IF ( .NOT. ALLOCATED( rad_sw_out_xy_av ) )  THEN
+             ALLOCATE( rad_sw_out_xy_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF  
+          IF ( k == 1 )  READ ( 13 )  tmp_2d
+          rad_sw_out_xy_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =     &
+                        tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                        
        CASE ( 'rad_lw_in' )
           IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
