Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model.f90	(revision 1818)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 1826)
@@ -1,3 +1,3 @@
-!> @file radiation_model.f90
+!> @file radiation_model_mod.f90
 !--------------------------------------------------------------------------------!
 ! This file is part of PALM.
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Further modularization.
 ! 
 ! Former revisions:
@@ -80,5 +80,5 @@
 ! ------------
 !> Radiation models and interfaces
-!> @todo move variable definitions used in init_radiation only to the subroutine
+!> @todo move variable definitions used in radiation_init only to the subroutine
 !>       as they are no longer required after initialization.
 !> @todo Output of full column vertical profiles used in RRTMG
@@ -115,7 +115,4 @@
 
 #if defined ( __rrtmg )
-
-!     USE netcdf_interface,                                                      &
-!         ONLY:  nc_stat, netcdf_handle_error
 
     USE parrrsw,                                                               &
@@ -353,12 +350,32 @@
 #endif
 
-    INTERFACE init_radiation
-       MODULE PROCEDURE init_radiation
-    END INTERFACE init_radiation
-
+    INTERFACE radiation_check_data_output
+       MODULE PROCEDURE radiation_check_data_output
+    END INTERFACE radiation_check_data_output
+
+    INTERFACE radiation_check_data_output_pr
+       MODULE PROCEDURE radiation_check_data_output_pr
+    END INTERFACE radiation_check_data_output_pr
+  
+    INTERFACE radiation_check_parameters
+       MODULE PROCEDURE radiation_check_parameters
+    END INTERFACE radiation_check_parameters
+  
     INTERFACE radiation_clearsky
        MODULE PROCEDURE radiation_clearsky
     END INTERFACE radiation_clearsky
-
+  
+    INTERFACE radiation_header
+       MODULE PROCEDURE radiation_header
+    END INTERFACE radiation_header 
+  
+    INTERFACE radiation_init
+       MODULE PROCEDURE radiation_init
+    END INTERFACE radiation_init
+
+    INTERFACE radiation_parin
+       MODULE PROCEDURE radiation_parin
+    END INTERFACE radiation_parin
+    
     INTERFACE radiation_rrtmg
        MODULE PROCEDURE radiation_rrtmg
@@ -374,15 +391,19 @@
     PRIVATE
 
-    PUBLIC albedo, albedo_type, albedo_type_name, albedo_lw_dif, albedo_lw_dir,&
-           albedo_sw_dif, albedo_sw_dir, constant_albedo, day_init, dots_rad,  &
-           dt_radiation, emissivity, force_radiation_call, init_radiation,     &
-           lambda, lw_radiation, net_radiation, rad_net, rad_net_av, radiation,&
-           radiation_clearsky, radiation_rrtmg, radiation_scheme,              &
-           radiation_tendency, rad_lw_in, rad_lw_in_av, rad_lw_out,            &
-           rad_lw_out_av, rad_lw_out_change_0, 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, skip_time_do_radiation, sw_radiation,       &
-           time_radiation, time_utc_init, unscheduled_radiation_calls
+!
+!-- Public functions
+    PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
+           radiation_check_parameters, radiation_clearsky, radiation_header,   &
+           radiation_init, radiation_parin, radiation_rrtmg, radiation_tendency
+    
+!
+!-- Public variables and constants
+    PUBLIC dots_rad, dt_radiation, force_radiation_call,                       &
+           rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
+           rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0,       &
+           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,                 &
+           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls
 
 
@@ -396,7 +417,314 @@
 ! Description:
 ! ------------
+!> Check data output for radiation model
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
+ 
+ 
+       USE control_parameters,                                                 &
+           ONLY: data_output, message_string
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=*) ::  unit     !< 
+       CHARACTER (LEN=*) ::  var      !<
+
+       INTEGER(iwp) :: i
+       INTEGER(iwp) :: ilen
+       INTEGER(iwp) :: k
+
+       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' )
+             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
+                message_string = '"output of "' // TRIM( var ) // '" requi' // &
+                                 'res radiation = .TRUE. and ' //              &
+                                 'radiation_scheme = "rrtmg"'
+                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'W/m2'     
+
+          CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
+                 'rrtm_asdir*' )
+             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+                message_string = 'illegal value for data_output: "' //         &
+                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
+                                 'cross sections are allowed for this value'
+                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
+                IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
+                     TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
+                     TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
+                     TRIM( var ) == 'rrtm_asdir*'      )                       &
+                THEN
+                   message_string = 'output of "' // TRIM( var ) // '" require'&
+                                    // 's radiation = .TRUE. and radiation_sch'&
+                                    // 'eme = "rrtmg"'
+                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
+                ENDIF
+             ENDIF
+
+             IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
+             IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
+             IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = '' 
+             IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = '' 
+             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = '' 
+
+          CASE DEFAULT
+             unit = 'illegal'
+
+       END SELECT
+
+
+    END SUBROUTINE radiation_check_data_output
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check data output of profiles for radiation model
+!------------------------------------------------------------------------------!  
+    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
+ 
+       USE arrays_3d,                                                          &
+           ONLY: zu
+
+       USE control_parameters,                                                 &
+           ONLY: data_output_pr, message_string
+
+       USE indices
+
+       USE profil_parameter
+
+       USE statistics
+
+       IMPLICIT NONE
+   
+       CHARACTER (LEN=*) ::  unit      !< 
+       CHARACTER (LEN=*) ::  variable  !< 
+       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
+ 
+       INTEGER(iwp) ::  user_pr_index !< 
+       INTEGER(iwp) ::  var_count     !< 
+
+       SELECT CASE ( TRIM( variable ) )
+       
+         CASE ( 'rad_net' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 101
+                dopr_unit  = 'W/m2'
+                hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+
+          CASE ( 'rad_lw_in' )
+             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 102
+                dopr_unit  = 'W/m2'
+                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit 
+             ENDIF
+
+          CASE ( 'rad_lw_out' )
+             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 103
+                dopr_unit  = 'W/m2'
+                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit   
+             ENDIF
+
+          CASE ( 'rad_sw_in' )
+             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 104
+                dopr_unit  = 'W/m2'
+                hom(:,2,104,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+
+          CASE ( 'rad_sw_out')
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme = "constant"'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 105
+                dopr_unit  = 'W/m2'
+                hom(:,2,105,:)  = SPREAD( zw, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+
+          CASE ( 'rad_lw_cs_hr' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 106
+                dopr_unit  = 'K/h'
+                hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+
+          CASE ( 'rad_lw_hr' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 107
+                dopr_unit  = 'K/h'
+                hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+
+          CASE ( 'rad_sw_cs_hr' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 108
+                dopr_unit  = 'K/h'
+                hom(:,2,108,:)  = SPREAD( zu, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+
+          CASE ( 'rad_sw_hr' )
+             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
+             THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) // ' is' // &
+                                 'not available for radiation = .FALSE. or ' //&
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(var_count) = 109
+                dopr_unit  = 'K/h'
+                hom(:,2,109,:)  = SPREAD( zu, 2, statistic_regions+1 )
+                unit = dopr_unit
+             ENDIF
+
+
+          CASE DEFAULT
+             unit = 'illegal'
+
+       END SELECT
+
+
+    END SUBROUTINE radiation_check_data_output_pr
+ 
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check parameters routine for radiation model
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_check_parameters
+
+       USE control_parameters,                                                 &
+           ONLY: message_string, topography
+                 
+    
+       IMPLICIT NONE
+
+       IF ( radiation_scheme /= 'constant'   .AND.                             &
+            radiation_scheme /= 'clear-sky'  .AND.                             &
+            radiation_scheme /= 'rrtmg' )  THEN
+          message_string = 'unknown radiation_scheme = '//                     &
+                           TRIM( radiation_scheme )
+          CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
+       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
+#if ! defined ( __rrtmg )
+          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
+                           'compilation of PALM with pre-processor ' //        &
+                           'directive -D__rrtmg'
+          CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
+#endif
+#if defined ( __rrtmg ) && ! defined( __netcdf )
+          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
+                           'the use of NetCDF (preprocessor directive ' //     &
+                           '-D__netcdf'
+          CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
+#endif
+
+       ENDIF
+
+       IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.             &
+            radiation_scheme == 'clear-sky')  THEN
+          message_string = 'radiation_scheme = "clear-sky" in combination' //  & 
+                           'with albedo_type = 0 requires setting of albedo'// &
+                           ' /= 9999999.9'
+          CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.        &
+          (    albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
+          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 
+          ) ) THEN
+          message_string = 'radiation_scheme = "rrtmg" in combination' //      & 
+                           'with albedo_type = 0 requires setting of ' //      &
+                           'albedo_lw_dif /= 9999999.9' //                     &
+                           'albedo_lw_dir /= 9999999.9' //                     &
+                           'albedo_sw_dif /= 9999999.9 and' //                 &
+                           'albedo_sw_dir /= 9999999.9'
+          CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( topography /= 'flat' )  THEN
+          message_string = 'radiation scheme cannot be used ' //               & 
+                           'in combination with  topography /= "flat"'
+          CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
+       ENDIF
+ 
+    END SUBROUTINE radiation_check_parameters 
+ 
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Initialization of the radiation model
 !------------------------------------------------------------------------------!
-    SUBROUTINE init_radiation
+    SUBROUTINE radiation_init
     
        IMPLICIT NONE
@@ -661,5 +989,5 @@
        RETURN
 
-    END SUBROUTINE init_radiation
+    END SUBROUTINE radiation_init
 
 
@@ -719,4 +1047,115 @@
 
     END SUBROUTINE radiation_clearsky
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Header output for radiation model
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_header ( io )
+
+
+       IMPLICIT NONE
+ 
+       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
+    
+
+       
+!
+!--    Write radiation model header
+       WRITE( io, 3 )
+
+       IF ( radiation_scheme == "constant" )  THEN
+          WRITE( io, 4 ) net_radiation
+       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
+          WRITE( io, 5 )
+       ELSEIF ( radiation_scheme == "rrtmg" )  THEN
+          WRITE( io, 6 )
+          IF ( .NOT. lw_radiation )  WRITE( io, 10 )
+          IF ( .NOT. sw_radiation )  WRITE( io, 11 )
+       ENDIF 
+
+       IF ( albedo_type == 0 )  THEN
+          WRITE( io, 7 ) albedo
+       ELSE
+          WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
+       ENDIF
+       IF ( constant_albedo )  THEN
+          WRITE( io, 9 )
+       ENDIF
+       
+       IF ( radiation .AND. radiation_scheme /= 'constant' )  THEN
+          WRITE ( io, 1 )  lambda
+          WRITE ( io, 2 )  day_init, time_utc_init
+       ENDIF
+
+       WRITE( io, 12 ) dt_radiation
+ 
+
+ 1 FORMAT ('    Geograph. longitude            :   lambda = ',F4.1,' degr')
+ 2 FORMAT ('    Day of the year at model start :   day_init = ',I3   &
+            /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
+ 3 FORMAT (//' Radiation model information:'/                                  &
+              ' ----------------------------'/)
+ 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,        &
+           // 'W/m**2')
+ 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',   &
+                   ' default)')
+ 6 FORMAT ('    --> RRTMG scheme is used')
+ 7 FORMAT (/'    User-specific surface albedo: albedo =', F6.3)
+ 8 FORMAT (/'    Albedo is set for land surface type: ', A)
+ 9 FORMAT (/'    --> Albedo is fixed during the run')
+10 FORMAT (/'    --> Longwave radiation is disabled')
+11 FORMAT (/'    --> Shortwave radiation is disabled.')
+12 FORMAT  ('    Timestep: dt_radiation = ', F6.2, '  s')
+
+
+    END SUBROUTINE radiation_header
+   
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Parin for &radiation_par for radiation model
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_parin
+
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
+       
+       NAMELIST /radiation_par/   albedo, albedo_type, albedo_lw_dir,          &
+                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
+                                  constant_albedo, day_init, dt_radiation,     &
+                                  lambda, lw_radiation, net_radiation,         &
+                                  radiation_scheme, skip_time_do_radiation,    &
+                                  sw_radiation, time_utc_init,                 &
+                                  unscheduled_radiation_calls
+       
+       line = ' '
+       
+!
+!--    Try to find radiation model package
+       REWIND ( 11 )
+       line = ' '
+       DO   WHILE ( INDEX( line, '&radiation_par' ) == 0 )
+          READ ( 11, '(A)', END=10 )  line
+       ENDDO
+       BACKSPACE ( 11 )
+
+!
+!--    Read user-defined namelist
+       READ ( 11, radiation_par )
+
+!
+!--    Set flag that indicates that the radiation model is switched on
+       radiation = .TRUE.
+
+ 10    CONTINUE
+       
+
+    END SUBROUTINE radiation_parin
 
 
@@ -1608,4 +2047,5 @@
     END SUBROUTINE read_trace_gas_data
 
+
     SUBROUTINE netcdf_handle_error_rad( routine_name, errno )
 
