Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 1857)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 1976)
@@ -19,5 +19,8 @@
 ! Current revisions:
 ! -----------------
-! 
+! Output of 2D/3D/masked data is now directly done within this module. The
+! radiation schemes have been simplified for better usability so that
+! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
+! the radiation code used.
 ! 
 ! Former revisions:
@@ -130,5 +133,4 @@
 
 #if defined ( __rrtmg )
-
     USE parrrsw,                                                               &
         ONLY:  naerec, nbndsw
@@ -276,5 +278,6 @@
 !
 !-- Flag parameters for RRTMGS (should not be changed)
-    INTEGER(iwp), PARAMETER :: rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
+    INTEGER(iwp), PARAMETER :: rrtm_idrv     = 1, & !< flag for longwave upward flux calculation option (0,1)
+                               rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
                                rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
                                rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
@@ -289,6 +292,5 @@
     INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
                     rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
-                    rrtm_iaer = 0,     & !< aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
-                    rrtm_idrv = 1        !< longwave upward flux calculation option (0,1)
+                    rrtm_iaer = 0        !< aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
 
     INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
@@ -385,4 +387,28 @@
     END INTERFACE radiation_constant
  
+    INTERFACE radiation_control
+       MODULE PROCEDURE radiation_control
+    END INTERFACE radiation_control
+
+    INTERFACE radiation_3d_data_averaging
+       MODULE PROCEDURE radiation_3d_data_averaging
+    END INTERFACE radiation_3d_data_averaging
+
+    INTERFACE radiation_data_output_2d
+       MODULE PROCEDURE radiation_data_output_2d
+    END INTERFACE radiation_data_output_2d
+
+    INTERFACE radiation_data_output_3d
+       MODULE PROCEDURE radiation_data_output_3d
+    END INTERFACE radiation_data_output_3d
+
+    INTERFACE radiation_data_output_mask
+       MODULE PROCEDURE radiation_data_output_mask
+    END INTERFACE radiation_data_output_mask
+
+    INTERFACE radiation_define_netcdf_grid
+       MODULE PROCEDURE radiation_define_netcdf_grid
+    END INTERFACE radiation_define_netcdf_grid
+
     INTERFACE radiation_header
        MODULE PROCEDURE radiation_header
@@ -406,4 +432,12 @@
     END INTERFACE radiation_tendency
 
+    INTERFACE radiation_read_restart_data
+       MODULE PROCEDURE radiation_read_restart_data
+    END INTERFACE radiation_read_restart_data
+
+    INTERFACE radiation_last_actions
+       MODULE PROCEDURE radiation_last_actions
+    END INTERFACE radiation_last_actions
+
     SAVE
 
@@ -411,12 +445,15 @@
 
 !
-!-- Public functions
+!-- Public functions / NEEDS SORTING
     PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
-           radiation_check_parameters, radiation_clearsky, radiation_constant, &
-           radiation_header, radiation_init, radiation_parin, radiation_rrtmg, &
-           radiation_tendency
+           radiation_check_parameters, radiation_control,                      &
+           radiation_header, radiation_init, radiation_parin,                  &
+           radiation_3d_data_averaging, radiation_tendency,                    &
+           radiation_data_output_2d, radiation_data_output_3d,                 &
+           radiation_define_netcdf_grid, radiation_last_actions,               &
+           radiation_read_restart_data, radiation_data_output_mask
     
 !
-!-- Public variables and constants
+!-- Public variables and constants / NEEDS SORTING
     PUBLIC dots_rad, dt_radiation, force_radiation_call,                       &
            rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
@@ -424,13 +461,43 @@
            rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
            rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
-           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb,                 &
+           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av,                           &
            skip_time_do_radiation, time_radiation, unscheduled_radiation_calls
 
 
 #if defined ( __rrtmg )
-    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv
+    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
 #endif
 
  CONTAINS
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine controls the calls of the radiation schemes
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_control
+ 
+ 
+       IMPLICIT NONE
+
+
+       SELECT CASE ( TRIM( radiation_scheme ) )
+
+          CASE ( 'constant' )
+             CALL radiation_constant
+          
+          CASE ( 'clear-sky' )
+             CALL radiation_clearsky
+       
+          CASE ( 'rrtmg' )
+             CALL radiation_rrtmg
+
+          CASE DEFAULT
+
+       END SELECT
+
+
+    END SUBROUTINE radiation_control
 
 !------------------------------------------------------------------------------!
@@ -456,6 +523,5 @@
        SELECT CASE ( TRIM( var ) )
 
-         CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_cs_hr', 'rad_lw_hr',       &
-                 'rad_sw_in', 'rad_sw_out', 'rad_sw_cs_hr', 'rad_sw_hr' )
+         CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )
              IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
                 message_string = '"output of "' // TRIM( var ) // '" requi' // &
@@ -779,14 +845,6 @@
        ENDIF
 
-
-       IF ( radiation_scheme == 'constant' )  THEN
-
-          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
-             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
-          ENDIF
-
-       ENDIF
-
-       IF ( radiation_scheme == 'clear-sky' )  THEN
+       IF ( radiation_scheme == 'clear-sky'  .OR.                              &
+            radiation_scheme == 'constant')  THEN
 
           ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
@@ -1075,4 +1133,8 @@
                             + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
 
+
+             rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emissivity         &
+                                        * (pt(k,j,i) * exn) ** 3
+
           ENDDO
        ENDDO
@@ -1093,4 +1155,5 @@
        INTEGER(iwp) :: i, j, k   !< loop indices
        REAL(wp)     :: exn,   &  !< Exner functions at surface
+                       exn1,  &  !< Exner functions at first grid level
                        pt1       !< potential temperature at first grid level
 
@@ -1099,6 +1162,5 @@
        exn = (surface_pressure / 1000.0_wp )**0.286_wp
 !
-!--    Prescribe net radiation and estimate a longwave outgoing radiative
-!--    flux (needed in land surface model)
+!--    Prescribe net radiation and estimate the remaining radiative fluxes
        DO i = nxlg, nxrg
           DO j = nysg, nyng
@@ -1106,5 +1168,19 @@
 
              rad_net(j,i)      = net_radiation
+
+             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
+
+             IF ( cloud_physics )  THEN
+                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
+                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
+             ELSE
+                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
+             ENDIF
+
              rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
+
+             rad_sw_in(0,j,i) = ( rad_net(j,i) - rad_lw_in(0,j,i)              &
+                                  + rad_lw_out(0,j,i) )                        &
+                                  / ( 1.0_wp - alpha(j,i) )
 
           ENDDO
@@ -2147,24 +2223,25 @@
 !> Cache-optimized version.
 !------------------------------------------------------------------------------!
-    SUBROUTINE radiation_tendency_ij ( i, j, tend )
-
-       USE cloud_parameters,                                                   &
-           ONLY:  pt_d_t
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) :: i, j, k !< loop indices
-
-       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
-
-#if defined ( __rrtmg )
+ SUBROUTINE radiation_tendency_ij ( i, j, tend )
+
+    USE cloud_parameters,                                                      &
+        ONLY:  pt_d_t
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) :: i, j, k !< loop indices
+
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
+
+    IF ( radiation_scheme == 'rrtmg' )  THEN
+#if defined  ( __rrtmg )
 !
 !--    Calculate tendency based on heating rate
        DO k = nzb+1, nzt+1
           tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i))    &
-                                      * pt_d_t(k) * d_seconds_hour
+                                         * pt_d_t(k) * d_seconds_hour
        ENDDO
-
 #endif
+    ENDIF
 
     END SUBROUTINE radiation_tendency_ij
@@ -2177,19 +2254,20 @@
 !> Vector-optimized version
 !------------------------------------------------------------------------------!
-    SUBROUTINE radiation_tendency ( tend )
-
-       USE cloud_parameters,                                                   &
-           ONLY:  pt_d_t
-
-       USE indices,                                                            &
-           ONLY:  nxl, nxr, nyn, nys
-
-       IMPLICIT NONE
-
-       INTEGER(iwp) :: i, j, k !< loop indices
-
-       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
-
-#if defined ( __rrtmg )
+ SUBROUTINE radiation_tendency ( tend )
+
+    USE cloud_parameters,                                                      &
+        ONLY:  pt_d_t
+
+    USE indices,                                                               &
+        ONLY:  nxl, nxr, nyn, nys
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) :: i, j, k !< loop indices
+
+    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
+
+    IF ( radiation_scheme == 'rrtmg' )  THEN
+#if defined  ( __rrtmg )
 !
 !--    Calculate tendency based on heating rate
@@ -2198,12 +2276,1280 @@
              DO k = nzb+1, nzt+1
                 tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i)                 &
-                                            +  rad_sw_hr(k,j,i) ) * pt_d_t(k)  &
-                                            * d_seconds_hour
-             ENDDO
-         ENDDO
+                                          +  rad_sw_hr(k,j,i) ) * pt_d_t(k)    &
+                                          * d_seconds_hour
+             ENDDO
+          ENDDO
        ENDDO
 #endif
-
-    END SUBROUTINE radiation_tendency
+    ENDIF
+
+
+ END SUBROUTINE radiation_tendency
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine for averaging 3D data
+!------------------------------------------------------------------------------!
+SUBROUTINE radiation_3d_data_averaging( mode, variable )
+ 
+
+    USE control_parameters
+
+    USE indices
+
+    USE kinds
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  mode    !< 
+    CHARACTER (LEN=*) :: variable !< 
+
+    INTEGER(iwp) ::  i !< 
+    INTEGER(iwp) ::  j !< 
+    INTEGER(iwp) ::  k !< 
+
+    IF ( mode == 'allocate' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+             CASE ( 'rad_net*' )
+                IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
+                   ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_net_av = 0.0_wp
+
+             CASE ( 'rad_lw_in' )
+                IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
+                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_lw_in_av = 0.0_wp
+
+             CASE ( 'rad_lw_out' )
+                IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
+                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_lw_out_av = 0.0_wp
+
+             CASE ( 'rad_lw_cs_hr' )
+                IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
+                   ALLOCATE( rad_lw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_lw_cs_hr_av = 0.0_wp
+
+             CASE ( 'rad_lw_hr' )
+                IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
+                   ALLOCATE( rad_lw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_lw_hr_av = 0.0_wp
+
+             CASE ( 'rad_sw_in' )
+                IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
+                   ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_sw_in_av = 0.0_wp
+
+             CASE ( 'rad_sw_out' )
+                IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
+                   ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_sw_out_av = 0.0_wp
+
+             CASE ( 'rad_sw_cs_hr' )
+                IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
+                   ALLOCATE( rad_sw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_sw_cs_hr_av = 0.0_wp
+
+             CASE ( 'rad_sw_hr' )
+                IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
+                   ALLOCATE( rad_sw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_sw_hr_av = 0.0_wp
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'sum' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'rad_net*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   rad_net_av(j,i) = rad_net_av(j,i) + rad_net(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_in' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) + rad_lw_in(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_out' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) + rad_lw_out(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_cs_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) + rad_lw_cs_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) + rad_lw_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_in' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) + rad_sw_in(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_out' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) + rad_sw_out(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_cs_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) + rad_sw_cs_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) + rad_sw_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'average' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+         CASE ( 'rad_net*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_in' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_out' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_cs_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_lw_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_in' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_out' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_cs_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_hr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+       END SELECT
+
+    ENDIF
+
+END SUBROUTINE radiation_3d_data_averaging
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining appropriate grid for netcdf variables.
+!> It is called out from subroutine netcdf.
+!------------------------------------------------------------------------------!
+SUBROUTINE radiation_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
+    
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*), INTENT(IN)  ::  var         !< 
+    LOGICAL, INTENT(OUT)           ::  found       !< 
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !< 
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !< 
+    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !< 
+
+    found  = .TRUE.
+
+
+!
+!-- Check for the grid
+    SELECT CASE ( TRIM( var ) )
+
+       CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr',        &
+              'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy',            &
+              'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
+              'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
+              'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' )
+          grid_x = 'x'
+          grid_y = 'y'
+          grid_z = 'zu'
+
+       CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',            &
+              'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
+              'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
+              'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
+          grid_x = 'x'
+          grid_y = 'y'
+          grid_z = 'zw'
+
+
+       CASE DEFAULT
+          found  = .FALSE.
+          grid_x = 'none'
+          grid_y = 'none'
+          grid_z = 'none'
+
+        END SELECT
+
+    END SUBROUTINE radiation_define_netcdf_grid
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining 3D output variables
+!------------------------------------------------------------------------------!
+ SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode,         &
+                                      local_pf, two_d )
+ 
+    USE indices
+
+    USE kinds
+
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  grid     !< 
+    CHARACTER (LEN=*) ::  mode     !< 
+    CHARACTER (LEN=*) ::  variable !< 
+
+    INTEGER(iwp) ::  av !< 
+    INTEGER(iwp) ::  i  !< 
+    INTEGER(iwp) ::  j  !< 
+    INTEGER(iwp) ::  k  !< 
+
+    LOGICAL      ::  found !< 
+    LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
+
+    REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !< 
+
+    found = .TRUE.
+
+    SELECT CASE ( TRIM( variable ) )
+
+       CASE ( 'rad_net*_xy' )        ! 2d-array
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   local_pf(i,j,nzb+1) = rad_net(j,i)
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   local_pf(i,j,nzb+1) = rad_net_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
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_in(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_in_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zu'
+
+       CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_out(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_out_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF   
+          IF ( mode == 'xy' )  grid = 'zu'
+
+       CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zw'
+
+       CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zw'
+
+       CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_in(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_in_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zu'
+
+       CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_out(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_out_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zu'
+
+       CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zw'
+
+       CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
+          IF ( av == 0 ) THEN
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_hr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng 
+                   DO  k = nzb, nzt+1
+                      local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+          IF ( mode == 'xy' )  grid = 'zw'
+
+       CASE DEFAULT
+          found = .FALSE.
+          grid  = 'none'
+
+    END SELECT
+ 
+ END SUBROUTINE radiation_data_output_2d
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining 3D output variables
+!------------------------------------------------------------------------------!
+ SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf )
+ 
+
+    USE indices
+
+    USE kinds
+
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  variable !< 
+
+    INTEGER(iwp) ::  av    !< 
+    INTEGER(iwp) ::  i     !< 
+    INTEGER(iwp) ::  j     !< 
+    INTEGER(iwp) ::  k     !< 
+
+    LOGICAL      ::  found !< 
+
+    REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !< 
+
+
+    found = .TRUE.
+
+
+    SELECT CASE ( TRIM( variable ) )
+
+      CASE ( 'rad_sw_in' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_in(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_in_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_sw_out' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_out(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_out_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_sw_cs_hr' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_sw_hr' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_hr(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_lw_in' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_in(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_in_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_lw_out' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_out(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_out_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_lw_cs_hr' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+      CASE ( 'rad_lw_hr' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_hr(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ELSE
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb, nzt+1
+                     local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+
+
+ END SUBROUTINE radiation_data_output_3d
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining masked data output
+!------------------------------------------------------------------------------!
+ SUBROUTINE radiation_data_output_mask( av, variable, found, local_pf )
+ 
+    USE control_parameters
+        
+    USE indices
+    
+    USE kinds
+    
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  variable   !< 
+
+    INTEGER(iwp) ::  av   !< 
+    INTEGER(iwp) ::  i    !< 
+    INTEGER(iwp) ::  j    !< 
+    INTEGER(iwp) ::  k    !< 
+
+    LOGICAL ::  found     !< 
+
+    REAL(wp),                                                                  &
+       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
+          local_pf   !< 
+
+
+    found = .TRUE.
+
+    SELECT CASE ( TRIM( variable ) )
+
+
+       CASE ( 'rad_lw_in' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_in(mask_k(mid,k),              &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_in_av(mask_k(mid,k),           &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_lw_out' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_out(mask_k(mid,k),             &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_out_av(mask_k(mid,k),          &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_lw_cs_hr' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_cs_hr(mask_k(mid,k),           &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_cs_hr_av(mask_k(mid,k),        &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_lw_hr' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_hr(mask_k(mid,k),              &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_lw_hr_av(mask_k(mid,k),           &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_sw_in' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_in(mask_k(mid,k),              &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_in_av(mask_k(mid,k),           &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_sw_out' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_out(mask_k(mid,k),             &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_out_av(mask_k(mid,k),          &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_sw_cs_hr' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_cs_hr(mask_k(mid,k),           &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_cs_hr_av(mask_k(mid,k),        &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'rad_sw_hr' )
+          IF ( av == 0 )  THEN
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_hr(mask_k(mid,k),              &
+                                            mask_j(mid,j),mask_i(mid,i))
+                    ENDDO
+                 ENDDO
+              ENDDO
+          ELSE
+             DO  i = 1, mask_size_l(mid,1)
+                DO  j = 1, mask_size_l(mid,2)
+                   DO  k = 1, mask_size_l(mid,3)
+                       local_pf(i,j,k) = rad_sw_hr_av(mask_k(mid,k),           &
+                                               mask_j(mid,j),mask_i(mid,i))
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+
+
+ END SUBROUTINE radiation_data_output_mask
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defines masked output variables
+!------------------------------------------------------------------------------!
+ SUBROUTINE radiation_last_actions
+ 
+
+    USE control_parameters
+        
+    USE kinds
+
+    IMPLICIT NONE
+
+    IF ( write_binary(1:4) == 'true' )  THEN
+       IF ( ALLOCATED( rad_net ) )  THEN
+          WRITE ( 14 )  'rad_net             ';  WRITE ( 14 )  rad_net  
+       ENDIF
+       IF ( ALLOCATED( rad_net_av ) )  THEN
+          WRITE ( 14 )  'rad_net_av          ';  WRITE ( 14 )  rad_net_av  
+       ENDIF  
+       IF ( ALLOCATED( rad_lw_in ) )  THEN
+          WRITE ( 14 )  'rad_lw_in           ';  WRITE ( 14 )  rad_lw_in  
+       ENDIF
+       IF ( ALLOCATED( rad_lw_in_av ) )  THEN
+          WRITE ( 14 )  'rad_lw_in_av        ';  WRITE ( 14 )  rad_lw_in_av  
+       ENDIF 
+       IF ( ALLOCATED( rad_lw_out ) )  THEN
+          WRITE ( 14 )  'rad_lw_out          ';  WRITE ( 14 )  rad_lw_out 
+       ENDIF 
+       IF ( ALLOCATED( rad_lw_out_av ) )  THEN
+          WRITE ( 14 )  'rad_lw_out_av       ';  WRITE ( 14 )  rad_lw_out_av  
+       ENDIF 
+       IF ( ALLOCATED( rad_lw_out_change_0 ) )  THEN
+          WRITE ( 14 )  'rad_lw_out_change_0 '
+          WRITE ( 14 )  rad_lw_out_change_0
+       ENDIF
+       IF ( ALLOCATED( rad_lw_cs_hr ) )  THEN
+          WRITE ( 14 )  'rad_lw_cs_hr        ';  WRITE ( 14 )  rad_lw_cs_hr
+       ENDIF
+       IF ( ALLOCATED( rad_lw_cs_hr_av ) )  THEN
+          WRITE ( 14 )  'rad_lw_cs_hr_av     ';  WRITE ( 14 )  rad_lw_cs_hr_av
+       ENDIF
+       IF ( ALLOCATED( rad_lw_hr ) )  THEN
+          WRITE ( 14 )  'rad_lw_hr           ';  WRITE ( 14 )  rad_lw_hr
+       ENDIF
+       IF ( ALLOCATED( rad_lw_hr_av ) )  THEN
+          WRITE ( 14 )  'rad_lw_hr_av        ';  WRITE ( 14 )  rad_lw_hr_av
+       ENDIF
+       IF ( ALLOCATED( rad_sw_in ) )  THEN
+          WRITE ( 14 )  'rad_sw_in           ';  WRITE ( 14 )  rad_sw_in  
+       ENDIF 
+       IF ( ALLOCATED( rad_sw_in_av ) )  THEN
+          WRITE ( 14 )  'rad_sw_in_av        ';  WRITE ( 14 )  rad_sw_in_av  
+       ENDIF 
+       IF ( ALLOCATED( rad_sw_out ) )  THEN
+          WRITE ( 14 )  'rad_sw_out          ';  WRITE ( 14 )  rad_sw_out  
+       ENDIF 
+       IF ( ALLOCATED( rad_sw_out_av ) )  THEN
+          WRITE ( 14 )  'rad_sw_out_av       ';  WRITE ( 14 )  rad_sw_out_av  
+       ENDIF 
+       IF ( ALLOCATED( rad_sw_cs_hr ) )  THEN
+          WRITE ( 14 )  'rad_sw_cs_hr        ';  WRITE ( 14 )  rad_sw_cs_hr
+       ENDIF
+       IF ( ALLOCATED( rad_sw_cs_hr_av ) )  THEN
+          WRITE ( 14 )  'rad_sw_cs_hr_av     ';  WRITE ( 14 )  rad_sw_cs_hr_av
+       ENDIF
+       IF ( ALLOCATED( rad_sw_hr ) )  THEN
+          WRITE ( 14 )  'rad_sw_hr           ';  WRITE ( 14 )  rad_sw_hr
+       ENDIF
+       IF ( ALLOCATED( rad_sw_hr_av ) )  THEN
+          WRITE ( 14 )  'rad_sw_hr_av        ';  WRITE ( 14 )  rad_sw_hr_av
+       ENDIF
+
+       WRITE ( 14 )  '*** end rad ***     '
+
+    ENDIF
+
+ END SUBROUTINE radiation_last_actions
+
+
+SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file,   &
+                                     nynfa, nyn_on_file, nysfa, nys_on_file,   &
+                                     offset_xa, offset_ya, overlap_count,      &
+                                     tmp_2d, tmp_3d )
+ 
+
+    USE control_parameters
+        
+    USE indices
+    
+    USE kinds
+    
+    USE pegrid
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=20) :: field_char   !< 
+
+    INTEGER(iwp) ::  i               !< 
+    INTEGER(iwp) ::  k               !< 
+    INTEGER(iwp) ::  nxlc            !< 
+    INTEGER(iwp) ::  nxlf            !< 
+    INTEGER(iwp) ::  nxl_on_file     !< 
+    INTEGER(iwp) ::  nxrc            !< 
+    INTEGER(iwp) ::  nxrf            !< 
+    INTEGER(iwp) ::  nxr_on_file     !< 
+    INTEGER(iwp) ::  nync            !< 
+    INTEGER(iwp) ::  nynf            !< 
+    INTEGER(iwp) ::  nyn_on_file     !< 
+    INTEGER(iwp) ::  nysc            !< 
+    INTEGER(iwp) ::  nysf            !< 
+    INTEGER(iwp) ::  nys_on_file     !< 
+    INTEGER(iwp) ::  overlap_count   !< 
+
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxlfa       !< 
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxrfa       !< 
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nynfa       !< 
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nysfa       !< 
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_xa   !< 
+    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !< 
+
+    REAL(wp),                                                                  &
+       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
+          tmp_2d   !< 
+
+    REAL(wp),                                                                  &
+       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
+          tmp_3d   !< 
+
+
+
+   IF ( initializing_actions == 'read_restart_data' )  THEN
+      READ ( 13 )  field_char
+
+      DO  WHILE ( TRIM( field_char ) /= '*** end rad ***' )
+
+         DO  k = 1, overlap_count
+
+            nxlf = nxlfa(i,k)
+            nxlc = nxlfa(i,k) + offset_xa(i,k)
+            nxrf = nxrfa(i,k)
+            nxrc = nxrfa(i,k) + offset_xa(i,k)
+            nysf = nysfa(i,k)
+            nysc = nysfa(i,k) + offset_ya(i,k)
+            nynf = nynfa(i,k)
+            nync = nynfa(i,k) + offset_ya(i,k)
+
+
+            SELECT CASE ( TRIM( field_char ) )
+
+                CASE ( 'rad_net' )
+                   IF ( .NOT. ALLOCATED( rad_net ) )  THEN
+                      ALLOCATE( rad_net(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_net_av' )
+                   IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
+                      ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   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' )
+                   IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
+                      ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_in_av' )
+                   IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
+                      ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_out' )
+                   IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
+                      ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_out_av' )
+                   IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
+                      ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_out_change_0' )
+                   IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) )  THEN
+                      ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
+                              = tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_cs_hr' )
+                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) )  THEN
+                      ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_cs_hr_av' )
+                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
+                      ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_hr' )
+                   IF ( .NOT. ALLOCATED( rad_lw_hr ) )  THEN
+                      ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_lw_hr_av' )
+                   IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
+                      ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_in' )
+                   IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
+                      ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_in_av' )
+                   IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
+                      ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_out' )
+                   IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
+                      ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_out_av' )
+                   IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
+                      ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_cs_hr' )
+                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) )  THEN
+                      ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_cs_hr_av' )
+                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
+                      ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_hr' )
+                   IF ( .NOT. ALLOCATED( rad_sw_hr ) )  THEN
+                      ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_hr_av' )
+                   IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
+                      ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+               CASE DEFAULT
+                  WRITE( message_string, * ) 'unknown variable named "',       &
+                                        TRIM( field_char ), '" found in',      &
+                                        '&data from prior run on PE ', myid
+                  CALL message( 'radiation_read_restart_data', 'PA0441', 1, 2, 0, 6, &
+                                 0 )
+
+            END SELECT
+
+         ENDDO
+
+         READ ( 13 )  field_char
+
+      ENDDO
+   ENDIF
+
+ END SUBROUTINE radiation_read_restart_data
+
 
  END MODULE radiation_model_mod
