Changeset 1691 for palm/trunk/SOURCE

Timestamp:

Oct 26, 2015 4:17:44 PM (9 years ago)

Author:

maronga

Message:

various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings

Location:

palm/trunk/SOURCE

Files:

• Makefile (modified) (10 diffs)
• advec_s_bc.f90 (modified) (4 diffs)
• average_3d_data.f90 (modified) (5 diffs)
• check_parameters.f90 (modified) (20 diffs)
• data_output_2d.f90 (modified) (7 diffs)
• data_output_3d.f90 (modified) (4 diffs)
• data_output_mask.f90 (modified) (4 diffs)
• diffusion_s.f90 (modified) (3 diffs)
• diffusion_u.f90 (modified) (2 diffs)
• flow_statistics.f90 (modified) (22 diffs)
• header.f90 (modified) (3 diffs)
• init_1d_model.f90 (modified) (13 diffs)
• init_3d_model.f90 (modified) (10 diffs)
• init_cloud_physics.f90 (modified) (2 diffs)
• init_grid.f90 (modified) (4 diffs)
• land_surface_model.f90 (modified) (54 diffs)
• lpm_advec.f90 (modified) (4 diffs)
• lpm_exchange_horiz.f90 (modified) (2 diffs)
• lpm_init.f90 (modified) (3 diffs)
• microphysics.f90 (modified) (17 diffs)
• modules.f90 (modified) (12 diffs)
• netcdf.f90 (modified) (18 diffs)
• package_parin.f90 (modified) (6 diffs)
• parin.f90 (modified) (10 diffs)
• prandtl_fluxes.f90 (deleted)
• production_e.f90 (modified) (10 diffs)
• prognostic_equations.f90 (modified) (13 diffs)
• radiation_model.f90 (modified) (44 diffs)
• read_3d_binary.f90 (modified) (7 diffs)
• read_var_list.f90 (modified) (7 diffs)
• sum_up_3d_data.f90 (modified) (11 diffs)
• surface_layer_fluxes.f90 (added)
• time_integration.f90 (modified) (9 diffs)
• wall_fluxes.f90 (modified) (8 diffs)
• write_3d_binary.f90 (modified) (7 diffs)
• write_var_list.f90 (modified) (8 diffs)

palm/trunk/SOURCE/Makefile

@@ -13 +13 @@
 # PALM. If not, see <http://www.gnu.org/licenses/>.
 #
-# Copyright 1997-2014  Leibniz Universitaet Hannover
+# Copyright 1997-2015  Leibniz Universitaet Hannover
 #--------------------------------------------------------------------------------#
 #
@@ -20 +20 @@
 # Current revisions:
 # ------------------
-# 
+# Replaced prandtl_fluxes with surface_layer_fluxes. Added radiation_model to 
+# prognostic_equations
 # 
 # Former revisions:
@@ -239 +240 @@
         mod_particle_attributes.f90 netcdf.f90 nudging.f90 package_parin.f90 \
         palm.f90 parin.f90 plant_canopy_model.f90 poisfft.f90 poismg.f90 \
-        poismg_fast.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 production_e.f90 \
+        poismg_fast.f90 pres.f90 print_1d.f90 production_e.f90 \
         prognostic_equations.f90 progress_bar.f90 radiation_model.f90 \
         random_function.f90 random_gauss.f90 random_generator_parallel.f90 \
@@ -245 +246 @@
         set_slicer_attributes_dvrp.f90 singleton.f90 sor.f90 \
         subsidence.f90 sum_up_3d_data.f90 \
-        surface_coupler.f90 swap_timelevel.f90 temperton_fft.f90 \
+        surface_coupler.f90 surface_layer_fluxes.f90 swap_timelevel.f90 temperton_fft.f90 \
         time_integration.f90 time_to_string.f90 timestep.f90 \
         timestep_scheme_steering.f90 transpose.f90 tridia_solver.f90 \
@@ -356 +357 @@
 init_3d_model.o: modules.o cpulog.o mod_kinds.o random_function.o advec_ws.o\
         land_surface_model.o ls_forcing.o lpm_init.o plant_canopy_model.o\
-        radiation_model.o random_generator_parallel.o
+        radiation_model.o random_generator_parallel.o surface_layer_fluxes.o
 init_advec.o: modules.o mod_kinds.o
 init_cloud_physics.o: modules.o mod_kinds.o
@@ -421 +422 @@
 poismg.o: modules.o cpulog.o mod_kinds.o
 poismg_fast.o: modules.o cpulog.o mod_kinds.o
-prandtl_fluxes.o: modules.o mod_kinds.o
 pres.o: modules.o cpulog.o mod_kinds.o poisfft.o poismg_fast.o
 print_1d.o: modules.o cpulog.o mod_kinds.o
@@ -430 +430 @@
         cpulog.o diffusion_e.o diffusion_s.o diffusion_u.o diffusion_v.o diffusion_w.o \
         eqn_state_seawater.o impact_of_latent_heat.o mod_kinds.o microphysics.o \
-        nudging.o plant_canopy_model.o production_e.o subsidence.o user_actions.o
+        nudging.o plant_canopy_model.o production_e.o radiation_model.o \
+        subsidence.o user_actions.o
 progress_bar.o: modules.o mod_kinds.o
 radiation_model.o : modules.o mod_particle_attributes.o
@@ -446 +447 @@
 subsidence.o: modules.o mod_kinds.o
 sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
-                  radiation_model.o
+                  radiation_model.o 
 surface_coupler.o: modules.o cpulog.o mod_kinds.o
+surface_layer_fluxes.o: modules.o mod_kinds.o land_surface_model.o
 swap_timelevel.o: modules.o cpulog.o mod_kinds.o land_surface_model.o
 temperton_fft.o: modules.o mod_kinds.o
@@ -454 +456 @@
         ls_forcing.o mod_kinds.o nudging.o production_e.o \
         prognostic_equations.o progress_bar.o radiation_model.o \
-        user_actions.o
+         user_actions.o surface_layer_fluxes.o
 time_to_string.o: mod_kinds.o
 timestep.o: modules.o cpulog.o mod_kinds.o
@@ -492 +494 @@
                    radiation_model.o random_function.o\
                    random_generator_parallel.o
-write_var_list.o: modules.o mod_kinds.o plant_canopy_model.o
+write_var_list.o: modules.o mod_kinds.o plant_canopy_model.o 

palm/trunk/SOURCE/advec_s_bc.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Formatting corrections
 ! 
 ! Former revisions:
@@ -464 +464 @@
                 IF ( sw(k,i+1) == 1.0_wp )  THEN
                    snenn = sk_p(k,j,i) - sk_p(k,j,i+2)
-                   IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
+                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                    sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn
                    sterm = MIN( sterm, 0.9999_wp )
@@ -760 +760 @@
                 IF ( sw(k,j+1) == 1.0_wp )  THEN
                    snenn = sk_p(k,j,i) - sk_p(k,j+2,i)
-                   IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
+                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                    sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn
                    sterm = MIN( sterm, 0.9999_wp )
@@ -1209 +1209 @@
                 IF ( sw(k+1,j) == 1.0_wp )  THEN
                    snenn = sk_p(k,j,i) - sk_p(k+2,j,i)
-                   IF ( ABS( snenn ) .LT. 1E-9_wp )  snenn = 1E-9_wp
+                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
                    sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn
                    sterm = MIN( sterm, 0.9999_wp )

palm/trunk/SOURCE/average_3d_data.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added output of Obukhov length and radiative heating rates for RRTMG. 
 ! 
 ! Former revisions:
@@ -96 +96 @@
     USE radiation_model_mod,                                                   &
         ONLY:  rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, rad_lw_out,       &
-               rad_lw_out_av, rad_sw_in, rad_sw_in_av, rad_sw_out,             &
-               rad_sw_out_av
+               rad_lw_out_av, 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
+
 
     IMPLICIT NONE
@@ -213 +216 @@
              ENDDO
 
+         CASE ( 'ol*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   ol_av(j,i) = ol_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
           CASE ( 'p' )
              DO  i = nxlg, nxrg
@@ -383 +393 @@
              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
@@ -397 +425 @@
                    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

palm/trunk/SOURCE/check_parameters.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added output of Obukhov length (ol) and radiative heating rates for RRTMG. 
+! Added checks for use of radiation / lsm with topography.
 ! 
 ! Former revisions:
@@ -307 +308 @@
     USE transpose_indices
 
+
     IMPLICIT NONE
 
@@ -650 +652 @@
           WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
        ENDIF
-       IF ( .NOT. prandtl_layer )  THEN
-          WRITE( action, '(A)' )  'prandtl_layer = .FALSE.'
+       IF ( .NOT. constant_flux_layer )  THEN
+          WRITE( action, '(A)' )  'constant_flux_layer = .FALSE.'
        ENDIF
        IF ( action /= ' ' )  THEN
@@ -1007 +1009 @@
        ENDIF
 
-       IF ( .NOT. prandtl_layer )  THEN
+       IF ( .NOT. constant_flux_layer )  THEN
           message_string = 'lsm requires '//                                   &
-                           'prandtl_layer = .T.'
+                           'constant_flux_layer = .T.'
           CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( topography /= 'flat' )  THEN
+          message_string = 'lsm cannot be used ' //  & 
+                           'in combination with  topography /= "flat"'
+          CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
        ENDIF
 
@@ -1196 +1204 @@
           CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
        ENDIF
-
-    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
+    ENDIF
 
     IF ( .NOT. ( loop_optimization == 'cache'  .OR.                            &
@@ -1729 +1740 @@
 !-- In case of using a prandtl-layer, calculated (or prescribed) surface
 !-- fluxes have to be used in the diffusion-terms
-    IF ( prandtl_layer )  use_surface_fluxes = .TRUE.
+    IF ( constant_flux_layer )  use_surface_fluxes = .TRUE.
 
 !
@@ -1796 +1807 @@
     ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
        ibc_e_b = 2
-       IF ( .NOT. prandtl_layer )  THEN
+       IF ( .NOT. constant_flux_layer )  THEN
           bc_e_b = 'neumann'
           ibc_e_b = 1
@@ -2003 +2014 @@
        IF ( surface_waterflux == 9999999.9_wp  )  THEN
           constant_waterflux = .FALSE.
-          IF ( large_scale_forcing )  THEN
+          IF ( large_scale_forcing .OR. land_surface )  THEN
              IF ( ibc_q_b == 0 )  THEN
                 constant_waterflux = .FALSE.
@@ -2045 +2056 @@
     ELSEIF ( bc_uv_b == 'neumann' )  THEN
        ibc_uv_b = 1
-       IF ( prandtl_layer )  THEN
+       IF ( constant_flux_layer )  THEN
           message_string = 'boundary condition: bc_uv_b = "' // &
-               TRIM( bc_uv_b ) // '" is not allowed with prandtl_layer = .TRUE.'
+               TRIM( bc_uv_b ) // '" is not allowed with constant_flux_layer'  &
+               // ' = .TRUE.'
           CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
        ENDIF
@@ -2356 +2368 @@
              dopr_unit(i)  = 'm2/s2'
              hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
-             IF ( prandtl_layer )  hom(nzb,2,12,:) = zu(1)
+             IF ( constant_flux_layer )  hom(nzb,2,12,:) = zu(1)
 
           CASE ( 'w*u*' )
@@ -2367 +2379 @@
              dopr_unit(i)  = 'm2/s2'
              hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
-             IF ( prandtl_layer )  hom(nzb,2,14,:) = zu(1)
+             IF ( constant_flux_layer )  hom(nzb,2,14,:) = zu(1)
 
           CASE ( 'w*v*' )
@@ -2393 +2405 @@
              dopr_unit(i)  = 'm2/s2'
              hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
-             IF ( prandtl_layer )  hom(nzb,2,19,:) = zu(1)
+             IF ( constant_flux_layer )  hom(nzb,2,19,:) = zu(1)
 
           CASE ( 'wv' )
@@ -2399 +2411 @@
              dopr_unit(i)  = 'm2/s2'
              hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
-             IF ( prandtl_layer )  hom(nzb,2,20,:) = zu(1)
+             IF ( constant_flux_layer )  hom(nzb,2,20,:) = zu(1)
 
           CASE ( 'w*pt*BC' )
@@ -3153 +3165 @@
              ENDIF
 
+          CASE ( 'rad_lw_cs_hr' )
+             IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
+                                 'lable for radiation = .FALSE. or ' //        &
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 106
+                dopr_unit(i)  = 'K/h'
+                hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'rad_lw_hr' )
+             IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
+                                 'lable for radiation = .FALSE. or ' //        &
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 107
+                dopr_unit(i)  = 'K/h'
+                hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'rad_sw_cs_hr' )
+             IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
+                                 'lable for radiation = .FALSE. or ' //        &
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 108
+                dopr_unit(i)  = 'K/h'
+                hom(:,2,108,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'rad_sw_hr' )
+             IF ( (.NOT. radiation) .OR. radiation_scheme /= 'rrtmg' )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
+                                 'lable for radiation = .FALSE. or ' //        &
+                                 'radiation_scheme /= "rrtmg"'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 109
+                dopr_unit(i)  = 'K/h'
+                hom(:,2,109,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
           CASE DEFAULT
 
@@ -3351 +3415 @@
              unit = 'kg/kg'
 
-
-          CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out' )
+          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' //  &
+                message_string = '"output of "' // TRIM( var ) // '" requi' // &
                                  'res radiation = .TRUE. and ' //              &
                                  'radiation_scheme = "rrtmg"'
@@ -3395 +3459 @@
 
           CASE ( 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'lai*', 'lwp*',     &
-                 'm_liq_eb*', 'pra*', 'prr*', 'qsws*', 'qsws_eb*',             &
+                 'm_liq_eb*', 'ol*', 'pra*', 'prr*', 'qsws*', 'qsws_eb*',      &
                  'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', 'rad_net*',  &
                  'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*',   &
@@ -3514 +3578 @@
              IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
              IF ( TRIM( var ) == 'lai*'   )  unit = 'none'
-             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/kg*m'
+             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
+             IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
+             IF ( TRIM( var ) == 'ol*'   )   unit = 'm'
              IF ( TRIM( var ) == 'pra*'   )  unit = 'mm'
              IF ( TRIM( var ) == 'prr*'   )  unit = 'mm/s'
@@ -3788 +3854 @@
           constant_diffusion = .TRUE.
 
-          IF ( prandtl_layer )  THEN
-             message_string = 'prandtl_layer is not allowed with fixed ' //    &
-                              'value of km'
+          IF ( constant_flux_layer )  THEN
+             message_string = 'constant_flux_layer is not allowed with fixed ' &
+                              // 'value of km'
              CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
           ENDIF
@@ -3816 +3882 @@
 
 !
-!-- Check value range for rif
-    IF ( rif_min >= rif_max )  THEN
-       WRITE( message_string, * )  'rif_min = ', rif_min, ' must be less ',    &
-                                   'than rif_max = ', rif_max
+!-- Check value range for zeta = z/L
+    IF ( zeta_min >= zeta_max )  THEN
+       WRITE( message_string, * )  'zeta_min = ', zeta_min, ' must be less ',  &
+                                   'than zeta_max = ', zeta_max
        CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
     ENDIF
@@ -4235 +4301 @@
 
 !
+!-- Check for valid setting of most_method
+    IF ( TRIM( most_method ) /= 'circular'  .AND.                              &
+         TRIM( most_method ) /= 'newton'  .AND.                                &
+         TRIM( most_method ) /= 'lookup' )  THEN
+       message_string = 'most_method = "' // TRIM( most_method ) //      &
+                        '" is unknown'
+       CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
+    ENDIF
+
+!
 !-- Check &userpar parameters
     CALL user_check_parameters

palm/trunk/SOURCE/data_output_2d.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added output of Obukhov length (ol) and radiative heating rates  for RRTMG. 
+! Formatting corrections.
 ! 
 ! Former revisions:
@@ -112 +113 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  dzw, e, nr, p, pt, q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws,      &
+        ONLY:  dzw, e, nr, ol, p, pt, q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws,  &
                rho, sa, shf, tend, ts, u, us, v, vpt, w, z0, z0h, zu, zw
         
@@ -162 +163 @@
     USE radiation_model_mod,                                                   &
         ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
-               rad_sw_out_av, rad_lw_in, rad_lw_in_av, rad_lw_out,             &
-               rad_lw_out_av 
+               rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,        &
+               rad_sw_hr_av, rad_lw_in, rad_lw_in_av, rad_lw_out,              &
+               rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr,        &
+               rad_lw_hr_av
 
     IMPLICIT NONE
@@ -543 +546 @@
                 IF ( mode == 'xy' )  level_z = zu
 
+             CASE ( 'ol*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = ol(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = ol_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
              CASE ( 'p_xy', 'p_xz', 'p_yz' )
                 IF ( av == 0 )  THEN
@@ -937 +958 @@
                 ENDIF
 
+             CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => rad_lw_cs_hr
+                ELSE
+                   to_be_resorted => rad_lw_cs_hr_av
+                ENDIF
+
+             CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => rad_lw_hr
+                ELSE
+                   to_be_resorted => rad_lw_hr_av
+                ENDIF
+
              CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
                 IF ( av == 0 )  THEN
@@ -949 +984 @@
                 ELSE
                    to_be_resorted => rad_sw_out_av
+                ENDIF
+
+             CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => rad_sw_cs_hr
+                ELSE
+                   to_be_resorted => rad_sw_cs_hr_av
+                ENDIF
+
+             CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => rad_sw_hr
+                ELSE
+                   to_be_resorted => rad_sw_hr_av
                 ENDIF
 
@@ -1227 +1276 @@
 !--                Exit the loop for layers beyond the data output domain
 !--                (used for soil model)
-                   IF ( layer_xy .GT. nzt_do )  THEN 
+                   IF ( layer_xy > nzt_do )  THEN 
                       EXIT loop1
                    ENDIF

palm/trunk/SOURCE/data_output_3d.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-! 
+! Added output of radiative heating rates for RRTMG
 ! 
 ! Former revisions:
@@ -138 +138 @@
     USE radiation_model_mod,                                                   &
         ONLY:  rad_lw_in, rad_lw_in_av, rad_lw_out, rad_lw_out_av,             &
-               rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av
+               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
 
 
@@ -495 +497 @@
              ENDIF
 
+          CASE ( 'rad_sw_cs_hr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_sw_cs_hr
+             ELSE
+                to_be_resorted => rad_sw_cs_hr_av
+             ENDIF
+
+          CASE ( 'rad_sw_hr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_sw_hr
+             ELSE
+                to_be_resorted => rad_sw_hr_av
+             ENDIF
+
           CASE ( 'rad_lw_in' )
              IF ( av == 0 )  THEN
@@ -507 +523 @@
              ELSE
                 to_be_resorted => rad_lw_out_av
+             ENDIF
+
+          CASE ( 'rad_lw_cs_hr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_lw_cs_hr
+             ELSE
+                to_be_resorted => rad_lw_cs_hr_av
+             ENDIF
+
+          CASE ( 'rad_lw_hr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_lw_hr
+             ELSE
+                to_be_resorted => rad_lw_hr_av
              ENDIF
 

palm/trunk/SOURCE/data_output_mask.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added output of radiative heating rates for RRTMG
 ! 
 ! Former revisions:
@@ -119 +119 @@
     USE radiation_model_mod,                                                   &
         ONLY:  rad_lw_in, rad_lw_in_av, rad_lw_out, rad_lw_out_av,             &
-               rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av
+               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
+
 
     IMPLICIT NONE
@@ -405 +408 @@
              ENDIF
 
+          CASE ( 'rad_lw_in' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_lw_in
+             ELSE
+                to_be_resorted => rad_lw_in_av
+             ENDIF
+
+          CASE ( 'rad_lw_out' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_lw_out
+             ELSE
+                to_be_resorted => rad_lw_out_av
+             ENDIF
+
+          CASE ( 'rad_lw_cs_hr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_lw_cs_hr
+             ELSE
+                to_be_resorted => rad_lw_cs_hr_av
+             ENDIF
+
+          CASE ( 'rad_lw_hr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_lw_hr
+             ELSE
+                to_be_resorted => rad_lw_hr_av
+             ENDIF
+
           CASE ( 'rad_sw_in' )
              IF ( av == 0 )  THEN
@@ -419 +450 @@
              ENDIF
 
-          CASE ( 'rad_lw_in' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rad_lw_in
-             ELSE
-                to_be_resorted => rad_lw_in_av
-             ENDIF
-
-          CASE ( 'rad_lw_out' )
-             IF ( av == 0 )  THEN
-                to_be_resorted => rad_lw_out
-             ELSE
-                to_be_resorted => rad_lw_out_av
+          CASE ( 'rad_sw_cs_hr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_sw_cs_hr
+             ELSE
+                to_be_resorted => rad_sw_cs_hr_av
+             ENDIF
+
+          CASE ( 'rad_sw_hr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => rad_sw_hr
+             ELSE
+                to_be_resorted => rad_sw_hr_av
              ENDIF
 

palm/trunk/SOURCE/diffusion_s.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-! 
+! Formatting corrections.
 ! 
 ! Former revisions:
@@ -144 +144 @@
 !
 !--          Apply prescribed horizontal wall heatflux where necessary
-             IF ( ( wall_w_x(j,i) .NE. 0.0_wp ) .OR. ( wall_w_y(j,i) .NE. 0.0_wp ) ) &
-             THEN
+             IF ( ( wall_w_x(j,i) /= 0.0_wp ) .OR. ( wall_w_y(j,i) /= 0.0_wp ) &
+                )  THEN
                 DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
 
@@ -400 +400 @@
 !
 !--    Apply prescribed horizontal wall heatflux where necessary
-       IF ( ( wall_w_x(j,i) .NE. 0.0_wp ) .OR. ( wall_w_y(j,i) .NE. 0.0_wp ) )       &
+       IF ( ( wall_w_x(j,i) /= 0.0_wp ) .OR. ( wall_w_y(j,i) /= 0.0_wp ) )     &
        THEN
           DO  k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)

palm/trunk/SOURCE/diffusion_u.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Formatting corrections.
 ! 
 ! Former revisions:
@@ -508 +508 @@
 !
 !--    Wall functions at the north and south walls, respectively
-       IF ( wall_u(j,i) .NE. 0.0_wp )  THEN
+       IF ( wall_u(j,i) /= 0.0_wp )  THEN
 
 !

palm/trunk/SOURCE/flow_statistics.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Revised calculation of Obukhov length. Added output of radiative heating >
+! rates for RRTMG.
 ! 
 ! Former revisions:
@@ -161 +162 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
-               qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
-               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
-               uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
+        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, pt, q, qc, ql, qr, &
+               qs, qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt,      &
+               td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us,& 
+               usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
         
     USE cloud_parameters,                                                      &
@@ -172 +173 @@
         ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
                 dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
-                large_scale_subsidence, max_pr_user, message_string, ocean,    &
-                passive_scalar, precipitation, simulated_time,                 &
+                large_scale_subsidence, max_pr_user, message_string, neutral,  &
+                ocean, passive_scalar, precipitation, simulated_time,          &
                 use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
                 ws_scheme_mom, ws_scheme_sca
@@ -199 +200 @@
     USE radiation_model_mod,                                                   &
         ONLY:  dots_rad, radiation, radiation_scheme, rad_net,                 &
-               rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
+               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 defined ( __rrtmg )
@@ -207 +209 @@
   
     USE statistics
+
 
     IMPLICIT NONE
@@ -240 +243 @@
     CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
 
-    !$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, usws, v, vsws, w )
+    !$acc update host( km, kh, e, ol, pt, qs, qsws, shf, ts, u, usws, v, vsws, w )
 
 !
@@ -744 +747 @@
 !
 !--                   Formula does not work if ql(nzb) /= 0.0
-                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &   ! w"q" (w"qv")
-                                          qsws(j,i) * rmask(j,i,sr)
+                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  & 
+                                          qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
                    ENDIF
                 ENDIF
                 IF ( passive_scalar )  THEN
-                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
-                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
+                                       qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
                 ENDIF
              ENDIF
+
+             IF ( .NOT. neutral )  THEN
+                sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
+                                    ol(j,i)  * rmask(j,i,sr) ! L
+             ENDIF
+
 
              IF ( land_surface )  THEN
@@ -774 +783 @@
 #if defined ( __rrtmg )
                 IF ( radiation_scheme == 'rrtmg' )  THEN
-                   sums_l(nzb,106,tn)  = sums_l(nzb,106,tn)  + rrtm_aldif(0,j,i)
-                   sums_l(nzb,107,tn)  = sums_l(nzb,107,tn)  + rrtm_aldir(0,j,i)
-                   sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  + rrtm_asdif(0,j,i)
-                   sums_l(nzb,109,tn)  = sums_l(nzb,109,tn)  + rrtm_asdir(0,j,i)
+                   sums_l(nzb,110,tn)  = sums_l(nzb,110,tn)  + rrtm_aldif(0,j,i)
+                   sums_l(nzb,111,tn)  = sums_l(nzb,111,tn)  + rrtm_aldir(0,j,i)
+                   sums_l(nzb,112,tn)  = sums_l(nzb,112,tn)  + rrtm_asdif(0,j,i)
+                   sums_l(nzb,113,tn)  = sums_l(nzb,113,tn)  + rrtm_asdir(0,j,i)
                 ENDIF
 #endif
@@ -997 +1006 @@
 !--    First calculate the products, then the divergence.
 !--    Calculation is time consuming. Do it only, if profiles shall be plotted.
-       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
-
+       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )   &
+       THEN
           sums_ll = 0.0_wp  ! local array
 
@@ -1037 +1046 @@
 !
 !--    Divergence of vertical flux of SGS TKE and the flux itself (69)
-       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
-
+       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )   &
+       THEN
           !$OMP DO
           DO  i = nxl, nxr
@@ -1109 +1118 @@
 !--    Collect current large scale advection and subsidence tendencies for
 !--    data output
-       IF ( large_scale_forcing  .AND.  ( simulated_time .GT. 0.0_wp ) )  THEN
+       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
 !
 !--       Interpolation in time of LSF_DATA 
@@ -1156 +1165 @@
              DO  j =  nys, nyn
                 DO  k = nzb_soil, nzt_soil
-                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i) * rmask(j,i,sr)
-                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i)         &
+                                      * rmask(j,i,sr)
+                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i)         &
+                                      * rmask(j,i,sr)
                 ENDDO
              ENDDO
@@ -1168 +1179 @@
              DO  j =  nys, nyn
                 DO  k = nzb_s_inner(j,i)+1, nzt+1
-                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i) * rmask(j,i,sr)
-                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i) * rmask(j,i,sr)
-                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i) * rmask(j,i,sr)
-                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
+                                       * rmask(j,i,sr)
+                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
+                                       * rmask(j,i,sr)
+                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
+                                       * rmask(j,i,sr)
+                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
+                                       * rmask(j,i,sr)
+                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
+                                       * rmask(j,i,sr)
+                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
+                                       * rmask(j,i,sr)
+                   sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
+                                       * rmask(j,i,sr)
+                   sums_l(k,109,tn)  = sums_l(k,108,tn)  + rad_sw_hr(k,j,i)    &
+                                       * rmask(j,i,sr)
                 ENDDO
              ENDDO
@@ -1374 +1397 @@
           hom(:,1,105,sr) = sums(:,105)            ! rad_sw_out
 
+          IF ( radiation_scheme == 'rrtmg' )  THEN
+             hom(:,1,106,sr) = sums(:,106)            ! rad_lw_cs_hr
+             hom(:,1,107,sr) = sums(:,107)            ! rad_lw_hr
+             hom(:,1,108,sr) = sums(:,108)            ! rad_sw_cs_hr
+             hom(:,1,109,sr) = sums(:,109)            ! rad_sw_hr
+
 #if defined ( __rrtmg )
-          IF ( radiation_scheme == 'rrtmg' )  THEN
-             hom(:,1,106,sr) = sums(:,106)            ! rrtm_aldif
-             hom(:,1,107,sr) = sums(:,107)            ! rrtm_aldir
-             hom(:,1,108,sr) = sums(:,108)            ! rrtm_asdif
-             hom(:,1,109,sr) = sums(:,109)            ! rrtm_asdir
+             hom(:,1,110,sr) = sums(:,110)            ! rrtm_aldif
+             hom(:,1,111,sr) = sums(:,111)            ! rrtm_aldir
+             hom(:,1,112,sr) = sums(:,112)            ! rrtm_asdif
+             hom(:,1,113,sr) = sums(:,113)            ! rrtm_asdir
+#endif
           ENDIF
-#endif
-
-       ENDIF
+
+
+       ENDIF
+
+       hom(:,1,114,sr) = sums(:,114)            !: L
 
        hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
@@ -1520 +1551 @@
        ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
        ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
-
+ 
+!
+!--    Calculate obukhov length
        IF ( ts_value(5,sr) /= 0.0_wp )  THEN
-          ts_value(22,sr) = ts_value(4,sr)**2 / &
-                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
+!           IF ( most_method == 'circular' )  THEN
+!              ts_value(22,sr) = ts_value(4,sr)**2 /                             &
+!                           ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) 
+!           ELSE
+             ts_value(22,sr) = hom(nzb,1,114,sr)  
+!          ENDIF     
        ELSE
           ts_value(22,sr) = 10000.0_wp
@@ -1553 +1590 @@
 #if defined ( __rrtmg )
           IF ( radiation_scheme == 'rrtmg' )  THEN
-             ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
-             ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr)          ! rrtm_aldir
-             ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr)          ! rrtm_asdif
-             ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
+             ts_value(dots_rad+5,sr) = hom(nzb,1,110,sr)          ! rrtm_aldif
+             ts_value(dots_rad+6,sr) = hom(nzb,1,111,sr)          ! rrtm_aldir
+             ts_value(dots_rad+7,sr) = hom(nzb,1,112,sr)          ! rrtm_asdif
+             ts_value(dots_rad+8,sr) = hom(nzb,1,113,sr)          ! rrtm_asdir
           ENDIF
 #endif
@@ -1642 +1679 @@
 #if defined ( __rrtmg )
     USE radiation_model_mod,                                                   &
-        ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
+        ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rad_lw_cs_hr,   &
+               rad_lw_hr,  rad_sw_cs_hr, rad_sw_hr
 #endif
 
@@ -3063 +3101 @@
 !--    Collect current large scale advection and subsidence tendencies for
 !--    data output
-       IF ( large_scale_forcing  .AND.  ( simulated_time .GT. 0.0_wp ) )  THEN
+       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
 !
 !--       Interpolation in time of LSF_DATA 
@@ -3110 +3148 @@
              DO  j =  nys, nyn
                 DO  k = nzb_soil, nzt_soil
-                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i) * rmask(j,i,sr)
-                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i)         &
+                                      * rmask(j,i,sr)
+                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i)         &
+                                      * rmask(j,i,sr)
                 ENDDO
              ENDDO
@@ -3123 +3163 @@
              DO  j =  nys, nyn
                 DO  k = nzb_s_inner(j,i)+1, nzt+1
-                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i) * rmask(j,i,sr)
-                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i) * rmask(j,i,sr)
-                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i) * rmask(j,i,sr)
-                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i) * rmask(j,i,sr)
+                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
+                                       * rmask(j,i,sr)
+                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
+                                       * rmask(j,i,sr)
+                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
+                                       * rmask(j,i,sr)
+                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
+                                       * rmask(j,i,sr)
+#if defined ( __rrtmg )
+                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
+                                       * rmask(j,i,sr)
+                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
+                                       * rmask(j,i,sr)
+                   sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
+                                       * rmask(j,i,sr)
+                   sums_l(k,109,tn)  = sums_l(k,109,tn)  + rad_sw_hr(k,j,i)    &
+                                       * rmask(j,i,sr)
+#endif
                 ENDDO
              ENDDO
@@ -3193 +3247 @@
           sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
           sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
-          sums(k,89:105)           = sums(k,89:105)     / ngp_2dh(sr)
-          sums(k,106:pr_palm-2)    = sums(k,106:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
+          sums(k,89:114)           = sums(k,89:114)     / ngp_2dh(sr)
+          sums(k,115:pr_palm-2)    = sums(k,115:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
+
        ENDDO
 
@@ -3446 +3501 @@
 
        IF ( ts_value(5,sr) /= 0.0_wp )  THEN
-          ts_value(22,sr) = ts_value(4,sr)**2_wp / &
-                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
+          ts_value(22,sr) = hom(nzb,1,114,sr)
        ELSE
           ts_value(22,sr) = 10000.0_wp

palm/trunk/SOURCE/header.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Renamed prandtl_layer to constant_flux_layer, renames rif_min/rif_max to
+! zeta_min/zeta_max.
 ! 
 ! Former revisions:
@@ -1118 +1119 @@
     ENDIF
 
-    IF ( prandtl_layer )  THEN
-       WRITE ( io, 305 )  (zu(1)-zu(0)), roughness_length, &
-                          z0h_factor*roughness_length, kappa, &
-                          rif_min, rif_max
+    IF ( constant_flux_layer )  THEN
+       WRITE ( io, 305 )  (zu(1)-zu(0)), roughness_length,                     &
+                          z0h_factor*roughness_length, kappa,                  &
+                          zeta_min, zeta_max
        IF ( .NOT. constant_heatflux )  WRITE ( io, 308 )
        IF ( humidity  .AND.  .NOT. constant_waterflux )  THEN
@@ -1131 +1132 @@
     ELSE
        IF ( INDEX(initializing_actions, 'set_1d-model_profiles') /= 0 )  THEN
-          WRITE ( io, 310 )  rif_min, rif_max
+          WRITE ( io, 310 )  zeta_min, zeta_max
        ENDIF
     ENDIF

palm/trunk/SOURCE/init_1d_model.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Renamed prandtl_layer to constant_flux_layer. rif is replaced by ol and zeta.
 ! 
 ! Former revisions:
@@ -72 +72 @@
 !> of the wind profile is being computed.
 !> All subroutines required can be found within this file.
+!>
+!> @todo harmonize code with new surface_layer_fluxes module
 !------------------------------------------------------------------------------!
  SUBROUTINE init_1d_model
@@ -90 +92 @@
     
     USE control_parameters,                                                    &
-        ONLY:  constant_diffusion, f, humidity, kappa, km_constant,            &
-               mixing_length_1d, passive_scalar, prandtl_layer,                &
-               prandtl_number, roughness_length, simulated_time_chr,           &
-               z0h_factor
+        ONLY:  constant_diffusion, constant_flux_layer, f, humidity, kappa,    &
+               km_constant, mixing_length_1d, passive_scalar, prandtl_number,  &
+               roughness_length, simulated_time_chr, z0h_factor
 
     IMPLICIT NONE
@@ -166 +167 @@
 !
 !-- For u*, theta* and the momentum fluxes plausible values are set
-    IF ( prandtl_layer )  THEN
+    IF ( constant_flux_layer )  THEN
        us1d = 0.1_wp   ! without initial friction the flow would not change
     ELSE
@@ -213 +214 @@
         
     USE control_parameters,                                                    &
-        ONLY:  constant_diffusion, dissipation_1d, humidity,                   &
-               intermediate_timestep_count, intermediate_timestep_count_max,   &
-               f, g, ibc_e_b, kappa, mixing_length_1d, passive_scalar,         &
-               prandtl_layer, rif_max, rif_min, simulated_time_chr,            &
-               timestep_scheme, tsc
+        ONLY:  constant_diffusion, constant_flux_layer, dissipation_1d,        &
+               humidity, intermediate_timestep_count,                          &
+               intermediate_timestep_count_max, f, g, ibc_e_b, kappa,          &  
+               mixing_length_1d, passive_scalar,                               &
+               simulated_time_chr, timestep_scheme, tsc, zeta_max, zeta_min
                 
     USE indices,                                                               &
@@ -334 +335 @@
 !--       Finite differences of the momentum fluxes are computed using half the 
 !--       normal grid length (2.0*ddzw(k)) for the sake of enhanced accuracy 
-          IF ( prandtl_layer )  THEN
+          IF ( constant_flux_layer )  THEN
 
              k = nzb+1
@@ -465 +466 @@
 !
 !--          First compute the vertical fluxes in the Prandtl-layer
-             IF ( prandtl_layer )  THEN
+             IF ( constant_flux_layer )  THEN
 !
 !--             Compute theta* using Rif numbers of the previous time step
@@ -498 +499 @@
                 ENDIF
 
-             ENDIF    ! prandtl_layer
+             ENDIF    ! constant_flux_layer
 
 !
@@ -506 +507 @@
 !--          the rif-numbers of the previous time step are used.
 
-             IF ( prandtl_layer )  THEN
+             IF ( constant_flux_layer )  THEN
                 IF ( .NOT. humidity )  THEN
                    pt_0 = pt_init(nzb+1)
@@ -547 +548 @@
 !--          range. It is exceeded excessively for very small velocities
 !--          (u,v --> 0).
-             WHERE ( rif1d < rif_min )  rif1d = rif_min
-             WHERE ( rif1d > rif_max )  rif1d = rif_max
+             WHERE ( rif1d < zeta_min )  rif1d = zeta_min
+             WHERE ( rif1d > zeta_max )  rif1d = zeta_max
 
 !
 !--          Compute u* from the absolute velocity value
-             IF ( prandtl_layer )  THEN
+             IF ( constant_flux_layer )  THEN
                 uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 )
 
@@ -639 +640 @@
                 ENDIF             
 
-             ENDIF   !  prandtl_layer
+             ENDIF   !  constant_flux_layer
 
 !
@@ -673 +674 @@
 !--          computed via the adiabatic mixing length, for the unstability has
 !--          already been taken account of via the TKE (cf. also Diss.).
-             IF ( prandtl_layer )  THEN
+             IF ( constant_flux_layer )  THEN
                 IF ( rif1d(nzb+1) >= 0.0_wp )  THEN
                    km1d(nzb+1) = us1d * kappa * zu(nzb+1) /                    &
@@ -805 +806 @@
 
        uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 )
-       IF ( ABS( v1d(nzb+1) ) .LT. 1.0E-5_wp )  THEN
+       IF ( ABS( v1d(nzb+1) ) < 1.0E-5_wp )  THEN
           alpha = ACOS( SIGN( 1.0_wp , u1d(nzb+1) ) )
        ELSE

palm/trunk/SOURCE/init_3d_model.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-! 
+! Call to init_surface_layer added. rif is replaced by ol and zeta.
 ! 
 ! Former revisions:
@@ -266 +266 @@
                sums_l_l, sums_up_fraction_l, sums_wsts_bc_l, ts_value,         &
                var_d, weight_pres, weight_substep 
-    
+ 
+    USE surface_layer_fluxes_mod,                                              &
+        ONLY:  init_surface_layer_fluxes
+   
     USE transpose_indices 
 
@@ -313 +316 @@
     ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) )
 
-    ALLOCATE( rif(nysg:nyng,nxlg:nxrg), shf(nysg:nyng,nxlg:nxrg),     &
+    ALLOCATE( ol(nysg:nyng,nxlg:nxrg), shf(nysg:nyng,nxlg:nxrg),      &
               ts(nysg:nyng,nxlg:nxrg), tswst(nysg:nyng,nxlg:nxrg),    &
               us(nysg:nyng,nxlg:nxrg), usws(nysg:nyng,nxlg:nxrg),     &
@@ -714 +717 @@
              hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
 
-             IF ( prandtl_layer )  THEN
-                rif  = rif1d(nzb+1)
+             IF ( constant_flux_layer )  THEN
+                ol   = ( zu(nzb+1) - zw(nzb) ) / rif1d(nzb+1)
                 ts   = 0.0_wp  ! could actually be computed more accurately in the
                                ! 1D model. Update when opportunity arises.
@@ -723 +726 @@
              ELSE
                 ts   = 0.0_wp  ! must be set, because used in
-                rif  = 0.0_wp  ! flowste
+                ol   = ( zu(nzb+1) - zw(nzb) ) / zeta_min  ! flowste
                 us   = 0.0_wp
                 usws = 0.0_wp
@@ -731 +734 @@
           ELSE
              e    = 0.0_wp  ! must be set, because used in
-             rif  = 0.0_wp  ! flowste
+             ol   = ( zu(nzb+1) - zw(nzb) ) / zeta_min  ! flowste
              ts   = 0.0_wp
              us   = 0.0_wp
@@ -886 +889 @@
              e    = 0.0_wp
           ENDIF
-          rif   = 0.0_wp
+          ol    = ( zu(nzb+1) - zw(nzb) ) / zeta_min
           ts    = 0.0_wp
           us    = 0.0_wp
@@ -1095 +1098 @@
 !
 !--    Initialize Prandtl layer quantities
-       IF ( prandtl_layer )  THEN
+       IF ( constant_flux_layer )  THEN
 
           z0 = roughness_length
@@ -1103 +1106 @@
 !
 !--          Surface temperature is prescribed. Here the heat flux cannot be
-!--          simply estimated, because therefore rif, u* and theta* would have
+!--          simply estimated, because therefore ol, u* and theta* would have
 !--          to be computed by iteration. This is why the heat flux is assumed
 !--          to be zero before the first time step. It approaches its correct
@@ -1624 +1627 @@
        CALL location_message( 'initializing land surface model', .FALSE. )
        CALL init_lsm
+       CALL location_message( 'finished', .TRUE. )
+    ENDIF
+
+!
+!-- Initialize surface layer (done after LSM as roughness length are required
+!-- for initialization
+    IF ( constant_flux_layer )  THEN
+       CALL location_message( 'initializing surface layer', .FALSE. )
+       CALL init_surface_layer_fluxes
        CALL location_message( 'finished', .TRUE. )
     ENDIF

palm/trunk/SOURCE/init_cloud_physics.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Removed typo
 ! 
 ! Former revisions:
@@ -117 +117 @@
     dpirho_l = 1.0_wp / pirho_l
 !
-!-- constant fot the sedimentation of cloud water (2-moment cloud physics)
+!-- constant for the sedimentation of cloud water (2-moment cloud physics)
     sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                   &
                           )**( 2.0_wp / 3.0_wp ) *                             &

palm/trunk/SOURCE/init_grid.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Renamed prandtl_layer to constant_flux_layer.
 ! 
 ! Former revisions:
@@ -144 +144 @@
                building_length_y, building_wall_left, building_wall_south,     &
                canyon_height, canyon_wall_left, canyon_wall_south,             &
-               canyon_width_x, canyon_width_y, coupling_char, dp_level_ind_b,  &
-               dz, dz_max, dz_stretch_factor, dz_stretch_level,                &
-               dz_stretch_level_index, ibc_uv_b, io_blocks, io_group,          &
-               inflow_l, inflow_n, inflow_r, inflow_s, masking_method,         &
-               maximum_grid_level, message_string, momentum_advec, ocean,      &
-               outflow_l, outflow_n, outflow_r, outflow_s, prandtl_layer,      &
-               psolver, scalar_advec, topography, topography_grid_convention,  &
-               use_surface_fluxes, use_top_fluxes, wall_adjustment_factor 
+               canyon_width_x, canyon_width_y, constant_flux_layer,            &
+               coupling_char, dp_level_ind_b, dz, dz_max, dz_stretch_factor,   &
+               dz_stretch_level, dz_stretch_level_index, ibc_uv_b, io_blocks,  &
+               io_group, inflow_l, inflow_n, inflow_r, inflow_s,               &
+               masking_method, maximum_grid_level, message_string,             &
+               momentum_advec, ocean, outflow_l, outflow_n, outflow_r,         &
+               outflow_s, psolver, scalar_advec, topography,                   &
+               topography_grid_convention, use_surface_fluxes, use_top_fluxes, &
+               wall_adjustment_factor 
         
     USE grid_variables,                                                        &
@@ -477 +478 @@
 !-- Define vertical gridpoint from (or to) which on the usual finite difference
 !-- form (which does not use surface fluxes) is applied 
-    IF ( prandtl_layer  .OR.  use_surface_fluxes )  THEN
+    IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
        nzb_diff = nzb + 2
     ELSE
@@ -955 +956 @@
 !-- the usual finite difference form (which does not use surface fluxes) is
 !-- applied 
-    IF ( prandtl_layer  .OR.  use_surface_fluxes )  THEN
+    IF ( constant_flux_layer  .OR.  use_surface_fluxes )  THEN
        nzb_diff_u         = nzb_u_inner + 2
        nzb_diff_v         = nzb_v_inner + 2

palm/trunk/SOURCE/land_surface_model.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added skip_time_do_lsm to allow for spin-ups without LSM. Various bugfixes:
+! Soil temperatures are now defined at the edges of the layers, calculation of
+! shb_eb corrected, prognostic equation for skin temperature corrected. Surface
+! fluxes are now directly transfered to atmosphere
 ! 
 ! Former revisions:
@@ -73 +76 @@
 !> DALES and UCLA-LES models.
 !>
-!> @todo Check dewfall parametrization for fog simulations.
-!> @todo Consider partial absorption of the net shortwave radiation by the surface layer.
-!> @todo Allow for water surfaces, check performance for bare soils.
-!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, nzt_soil=3)).
-!> @todo Implement surface runoff model (required when performing long-term LES with
-!>       considerable precipitation.
+!> @todo Consider partial absorption of the net shortwave radiation by the 
+!>       surface layer.
+!> @todo Allow for water surfaces
+!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, 
+!>       nzt_soil=3)).
+!> @todo Implement surface runoff model (required when performing long-term LES 
+!>       with considerable precipitation.
 !>
-!> @note  No time step criterion is required as long as the soil layers do not become
-!>        too thin.
+!> @note  No time step criterion is required as long as the soil layers do not 
+!>        become too thin.
 !------------------------------------------------------------------------------!
  MODULE land_surface_model_mod
  
-     USE arrays_3d,                                                            &
-         ONLY:  pt, pt_p, q_p, qsws, rif, shf, ts, us, z0, z0h
-
-     USE cloud_parameters,                                                     &
-         ONLY:  cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
-
-     USE control_parameters,                                                   &
-         ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,   &
-                initializing_actions, intermediate_timestep_count_max,         &
-                max_masks, precipitation, pt_surface, rho_surface,             &
-                roughness_length, surface_pressure, timestep_scheme, tsc,      &
-                z0h_factor, time_since_reference_point
-
-     USE indices,                                                              &
-         ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner 
-
-     USE kinds
+    USE arrays_3d,                                                             &
+        ONLY:  hyp, pt, pt_p, q, q_p, ql, qsws, rif, shf, ts, us, vpt, z0, z0h
+
+    USE cloud_parameters,                                                      &
+        ONLY:  cp, hyrho, l_d_cp, l_d_r, l_v, prr, pt_d_t, rho_l, r_d, r_v
+
+    USE control_parameters,                                                    &
+        ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,    &
+               initializing_actions, intermediate_timestep_count_max,          &
+               max_masks, precipitation, pt_surface, rho_surface,              &
+               roughness_length, surface_pressure, timestep_scheme, tsc,       &
+               z0h_factor, time_since_reference_point
+
+    USE indices,                                                               &
+        ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner 
+
+    USE kinds
 
     USE netcdf_control,                                                        &
         ONLY:  dots_label, dots_num, dots_unit
 
-     USE radiation_model_mod,                                                  &
-         ONLY:  radiation_scheme, rad_net, rad_sw_in, sigma_sb
-
-     USE statistics,                                                           &
-         ONLY:  hom, statistic_regions
+    USE pegrid
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  force_radiation_call, radiation_scheme, rad_net, rad_sw_in,     &
+               rad_lw_out, sigma_sb
+        
+#if defined ( __rrtmg )
+    USE radiation_model_mod,                                                   &
+        ONLY:  rrtm_lwuflx_dt, rrtm_idrv
+#endif
+
+    USE statistics,                                                            &
+        ONLY:  hom, statistic_regions
 
     IMPLICIT NONE
@@ -144 +156 @@
                     soil_type = 3    !< soil type, 0: user-defined, 1-6: generic (see list)
 
-    LOGICAL :: conserve_water_content = .TRUE., & !< open or closed bottom surface for the soil model
-               dewfall = .TRUE.,                & !< allow/inhibit dewfall
-               land_surface = .FALSE.             !< flag parameter indicating wheather the lsm is used
+    LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
+               dewfall = .TRUE.,                 & !< allow/inhibit dewfall
+               force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls
+               land_surface = .FALSE.              !< flag parameter indicating wheather the lsm is used
 
 !   value 9999999.9_wp -> generic available or user-defined value must be set
@@ -160 +173 @@
                 hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_sat
                 ke = 0.0_wp,                            & !< Kersten number
+                lambda_h_sat = 0.0_wp,                  & !< heat conductivity for saturated soil
                 lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s
                 lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u
@@ -174 +188 @@
                 rd_d_rv,                                & !< r_d / r_v
                 saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
+                skip_time_do_lsm = 0.0_wp,              & !< LSM is not called before this time
                 vegetation_coverage = 9999999.9_wp,     & !< NAMELIST c_veg
                 wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
                 z0_eb  = 9999999.9_wp,                  & !< NAMELIST z0 (lsm_par)
-                z0h_eb = 9999999.9_wp                    !< NAMELIST z0h (lsm_par)
+                z0h_eb = 9999999.9_wp                     !< NAMELIST z0h (lsm_par)
 
     REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
-              ddz_soil,                     & !< 1/dz_soil
-              ddz_soil_stag,                & !< 1/dz_soil_stag
-              dz_soil,                      & !< soil grid spacing (center-center)
-              dz_soil_stag,                 & !< soil grid spacing (edge-edge) 
-              root_extr = 0.0_wp,           & !< root extraction
+              ddz_soil_stag,                  & !< 1/dz_soil_stag
+              dz_soil_stag,                   & !< soil grid spacing (center-center)
+              root_extr = 0.0_wp,             & !< root extraction
               root_fraction = (/9999999.9_wp, 9999999.9_wp,    &
                                 9999999.9_wp, 9999999.9_wp /), & !< distribution of root surface area to the individual soil layers
@@ -191 +204 @@
 
     REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) ::   &
-              soil_temperature = (/290.0_wp, 287.0_wp, 285.0_wp,  283.0_wp,    &
-                                   283.0_wp /) !< soil temperature (K)
+              soil_temperature = (/290.0_wp, 287.0_wp, 285.0_wp,  283.0_wp,    & !< soil temperature (K)
+                                   283.0_wp /),                                &                                   
+              ddz_soil,                                                        & !< 1/dz_soil
+              dz_soil                                                            !< soil grid spacing (edge-edge)
 
 #if defined( __nopointer )
@@ -218 +233 @@
 !
 !-- Energy balance variables                
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::                                   &
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
               alpha_vg,         & !< coef. of Van Genuchten
               c_liq,            & !< liquid water coverage (of vegetated area)
@@ -232 +247 @@
               lai,              & !< leaf area index
               lai_av,           & !< average of lai
-              lambda_h_sat,     & !< heat conductivity for dry soil
-              lambda_surface_s,    & !< coupling between surface and soil (depends on vegetation type)
-              lambda_surface_u,    & !< coupling between surface and soil (depends on vegetation type)
+              lambda_surface_s, & !< coupling between surface and soil (depends on vegetation type)
+              lambda_surface_u, & !< coupling between surface and soil (depends on vegetation type)
               l_vg,             & !< coef. of Van Genuchten
               m_fc,             & !< soil moisture at field capacity (m3/m3)
@@ -253 +267 @@
               r_a_av,           & !< avergae of r_a
               r_canopy,         & !< canopy resistance
-              r_soil,           & !< soil resitance
+              r_soil,           & !< soil resistance
               r_soil_min,       & !< minimum soil resistance
               r_s,              & !< total surface resistance (combination of r_soil and r_canopy)
-              r_s_av,           & !< avergae of r_s
+              r_s_av,           & !< average of r_s
               r_canopy_min,     & !< minimum canopy (stomatal) resistance
               shf_eb,           & !< surface flux of sensible heat
               shf_eb_av           !< average of shf_eb
+
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
@@ -472 +487 @@
            lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance,          &
            min_soil_resistance, n_vangenuchten, residual_moisture, rho_cp,     &
-           rho_lv, root_fraction, saturation_moisture, soil_moisture,          &
-           soil_temperature, soil_type, soil_type_name, vegetation_coverage,   &
-           veg_type, veg_type_name, wilting_point, z0_eb, z0h_eb
+           rho_lv, root_fraction, saturation_moisture, skip_time_do_lsm,       &
+           soil_moisture, soil_temperature, soil_type, soil_type_name,         &
+           vegetation_coverage, veg_type, veg_type_name, wilting_point, z0_eb, &
+           z0h_eb
 
 !
@@ -483 +499 @@
            zs
 
-
 !
 !-- Public 2D output variables
@@ -490 +505 @@
            qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av,         &
            r_a, r_a_av, r_s, r_s_av, shf_eb, shf_eb_av
-
 
 !
@@ -565 +579 @@
        ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) )
        ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) )
@@ -612 +625 @@
        INTEGER(iwp) ::  k !< running index
 
+       REAL(wp) :: pt_1   !< potential temperature at first grid level
+
 
 !
 !--    Calculate Exner function
-       exn       = ( surface_pressure / 1000.0_wp )**0.286_wp
+       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
 
 
@@ -674 +689 @@
 !
 !--    Calculate grid spacings. Temperature and moisture are defined at
-!--    the center of the soil layers, whereas gradients/fluxes are defined
-!--    at the edges (_stag)
-       dz_soil_stag(nzb_soil) = zs(nzb_soil)
-
-       DO k = nzb_soil+1, nzt_soil
-          dz_soil_stag(k) = zs(k) - zs(k-1)
+!--    the edges of the soil layers (_stag), whereas gradients/fluxes are defined
+!--    at the centers
+       dz_soil(nzb_soil) = zs(nzb_soil)
+
+       DO  k = nzb_soil+1, nzt_soil
+          dz_soil(k) = zs(k) - zs(k-1)
        ENDDO
-
-       DO k = nzb_soil, nzt_soil-1
-          dz_soil(k) = 0.5 * (dz_soil_stag(k+1) + dz_soil_stag(k))
+       dz_soil(nzt_soil+1) = dz_soil(nzt_soil)
+
+       DO  k = nzb_soil, nzt_soil-1
+          dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k))
        ENDDO
-       dz_soil(nzt_soil) = dz_soil_stag(nzt_soil)
+       dz_soil_stag(nzt_soil) = dz_soil(nzt_soil)
 
        ddz_soil      = 1.0_wp / dz_soil
@@ -752 +768 @@
 !--       (make sure that the soil moisture does not drop below the permanent
 !--       wilting point) -> problems with devision by zero)
-          DO k = nzb_soil, nzt_soil
+          DO  k = nzb_soil, nzt_soil
              t_soil(k,:,:)    = soil_temperature(k)
              m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
@@ -761 +777 @@
 !
 !--       Calculate surface temperature
-          t_surface = pt_surface * exn
+          t_surface   = pt_surface * exn
+          t_surface_p = t_surface
 
 !
@@ -770 +787 @@
                 k = nzb_s_inner(j,i)
 
+                IF ( cloud_physics )  THEN
+                   pt_1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
+                ELSE
+                   pt_1 = pt(k+1,j,i)
+                ENDIF
+
 !
 !--             Assure that r_a cannot be zero at model start
-                IF ( pt(k+1,j,i) == pt(k,j,i) )  pt(k+1,j,i) = pt(k+1,j,i)     &
-                                                 + 1.0E-10_wp
-
-                us(j,i) = 0.1_wp
-                ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i)
+                IF ( pt_1 == pt(k,j,i) )  pt_1 = pt_1 + 1.0E-10_wp
+
+                us(j,i)  = 0.1_wp
+                ts(j,i)  = (pt_1 - pt(k,j,i)) / r_a(j,i)
                 shf(j,i) = - us(j,i) * ts(j,i)
              ENDDO
@@ -794 +816 @@
        ENDIF
 
-!
-!--    Calculate saturation soil heat conductivity
-       lambda_h_sat(:,:) = lambda_h_sm ** (1.0_wp - m_sat(:,:)) *              &
-                           lambda_h_water ** m_sat(:,:)
-
-
-
-
-       DO k = nzb_soil, nzt_soil
+       DO  k = nzb_soil, nzt_soil
           root_fr(k,:,:) = root_fraction(k)
        ENDDO
@@ -838 +852 @@
 
           IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
-             DO k = nzb_soil, nzt_soil
+             DO  k = nzb_soil, nzt_soil
                 root_fr(k,:,:) = root_distribution(k,veg_type)
                 root_fraction(k) = root_distribution(k,veg_type)
@@ -886 +900 @@
 
 !
-!--    Calculate humidity at the surface
-       IF ( humidity )  THEN
-          CALL calc_q_surface
-       ENDIF
-
-
-
-!
 !--    Add timeseries for land surface model
-
-       dots_label(dots_num+1) = "ghf_eb"
-       dots_label(dots_num+2) = "shf_eb"
-       dots_label(dots_num+3) = "qsws_eb"
-       dots_label(dots_num+4) = "qsws_liq_eb"
-       dots_label(dots_num+5) = "qsws_soil_eb"
-       dots_label(dots_num+6) = "qsws_veg_eb"
-       dots_unit(dots_num+1:dots_num+6) = "W/m2"
-
-       dots_label(dots_num+7) = "r_a"
-       dots_label(dots_num+8) = "r_s"
-       dots_unit(dots_num+7:dots_num+8) = "s/m"
-
        dots_soil = dots_num + 1
        dots_num  = dots_num + 8
+
+       dots_label(dots_soil) = "ghf_eb"
+       dots_label(dots_soil+1) = "shf_eb"
+       dots_label(dots_soil+2) = "qsws_eb"
+       dots_label(dots_soil+3) = "qsws_liq_eb"
+       dots_label(dots_soil+4) = "qsws_soil_eb"
+       dots_label(dots_soil+5) = "qsws_veg_eb"
+       dots_label(dots_soil+6) = "r_a"
+       dots_label(dots_soil+7) = "r_s"
+
+       dots_unit(dots_soil:dots_soil+5) = "W/m2"
+       dots_unit(dots_soil+6:dots_soil+7) = "s/m"
 
        RETURN
@@ -921 +926 @@
 ! ------------
 !> Solver for the energy balance at the surface.
-!> @note Surface fluxes are calculated in the land surface model, but these are
-!>       not used in the atmospheric code. The fluxes are calculated afterwards in
-!>       prandtl_fluxes using the surface values of temperature and humidity as
-!>       provided by the land surface model. In this way, the fluxes in the land
-!>       surface model are not equal to the ones calculated in prandtl_fluxes
 !------------------------------------------------------------------------------!
     SUBROUTINE lsm_energy_balance
@@ -940 +940 @@
                    f3,          & !< resistance correction term 3
                    m_min,       & !< minimum soil moisture
-                   T_1,         & !< actual temperature at first grid point
                    e,           & !< water vapour pressure
                    e_s,         & !< water vapour saturation pressure
@@ -954 +953 @@
                    f_shf,       & !< factor for shf_eb
                    lambda_surface, & !< Current value of lambda_surface
-                   m_liq_eb_max   !< maxmimum value of the liq. water reservoir
-
+                   m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
+                   pt_1,        & !< potential temperature at first grid level
+                   qv_1           !< specific humidity at first grid level
 
 !
@@ -961 +961 @@
        exn = ( surface_pressure / 1000.0_wp )**0.286_wp
 
-       
-
-       DO i = nxlg, nxrg
-          DO j = nysg, nyng
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              k = nzb_s_inner(j,i)
 
@@ -980 +978 @@
 !--          time step is used here. Note that this formulation is the
 !--          equivalent to the ECMWF formulation using drag coefficients
-             r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) +       &
-                         1.0E-20_wp)
+             IF ( cloud_physics )  THEN
+                pt_1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
+                qv_1 = q(k+1,j,i) - ql(k+1,j,i)
+             ELSE
+                pt_1 = pt(k+1,j,i)
+                qv_1 = q(k+1,j,i)
+             ENDIF
+
+             r_a(j,i) = (pt_1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
 
 !
@@ -1002 +1007 @@
 !--          f2 = 0 for very dry soils
              m_total = 0.0_wp
-             DO ks = nzb_soil, nzt_soil
+             DO  ks = nzb_soil, nzt_soil
                  m_total = m_total + root_fr(ks,j,i)                           &
                            * MAX(m_soil(ks,j,i),m_wilt(j,i))
              ENDDO 
 
-             IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LT. m_fc(j,i) )  THEN
+             IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) )  THEN
                 f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
-             ELSEIF ( m_total .GE. m_fc(j,i) )  THEN
+             ELSEIF ( m_total >= m_fc(j,i) )  THEN
                 f2 = 1.0_wp
              ELSE
@@ -1025 +1030 @@
 !
 !--             Calculate vapour pressure
-                e  = q_p(k+1,j,i) * surface_pressure / 0.622_wp
+                e  = qv_1 * surface_pressure / 0.622_wp
                 f3 = EXP ( -g_d(j,i) * (e_s - e) )
              ELSE
@@ -1067 +1072 @@
 !--          In case of dewfall, set evapotranspiration to zero
 !--          All super-saturated water is then removed from the air
-             IF ( humidity .AND. dewfall .AND. q_s .LE. q_p(k+1,j,i) )  THEN
+             IF ( humidity .AND. q_s <= qv_1 )  THEN
                 r_canopy(j,i) = 0.0_wp
                 r_soil(j,i)   = 0.0_wp
-             ENDIF 
+             ENDIF
 
 !
@@ -1085 +1090 @@
 !--          If soil moisture is below wilting point, plants do no longer
 !--          transpirate.
-!              IF ( m_soil(k,j,i) .LT. m_wilt(j,i) )  THEN
+!              IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
 !                 f_qsws_veg = 0.0_wp
 !              ENDIF
@@ -1094 +1099 @@
 !--          air.
              IF ( humidity )  THEN
-                IF ( q_s .LE. q_p(k+1,j,i) )  THEN
+                IF ( q_s <= qv_1 )  THEN
                    IF ( .NOT. dewfall )  THEN
                       f_qsws_veg  = 0.0_wp
@@ -1114 +1119 @@
              dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
 
-             T_1 = pt_p(k+1,j,i) * exn
-
 !
 !--          Add LW up so that it can be removed in prognostic equation
-             rad_net_l(j,i) = rad_net(j,i) + sigma_sb * t_surface(j,i) ** 4
+             rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
+
 
              IF ( humidity )  THEN
-
+#if defined ( __rrtmg )
 !
 !--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
-                         + f_shf / exn * T_1 + f_qsws * ( q_p(k+1,j,i) - q_s   &
+                coef_1 = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb)                &
+                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
+                         + f_shf * pt_1 + f_qsws * ( qv_1 - q_s                &
                          + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
                          * t_soil(nzb_soil,j,i)
@@ -1131 +1136 @@
 !
 !--             Denominator of the prognostic equation
-                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
-                         * dq_s_dt + lambda_surface + f_shf / exn
-
+                coef_2 = rrtm_lwuflx_dt(0,nzb) + f_qsws * dq_s_dt              &
+                         + lambda_surface + f_shf / exn
+#else
+
+!
+!--             Numerator of the prognostic equation
+                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
+                         * t_surface(j,i) ** 4                                 &
+                         + f_shf * pt_1 + f_qsws * ( qv_1                      &
+                         - q_s + dq_s_dt * t_surface(j,i) )                    &
+                         + lambda_surface * t_soil(nzb_soil,j,i)
+
+!
+!--                Denominator of the prognostic equation
+                   coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws   &
+                            * dq_s_dt + lambda_surface + f_shf / exn
+ 
+#endif
              ELSE
 
+#if defined ( __rrtmg )
 !
 !--             Numerator of the prognostic equation
-                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
-                         + f_shf / exn * T_1 + lambda_surface                  &
+                coef_1 = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb)                &
+                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
+                         + f_shf * pt_1  + lambda_surface                      &
                          * t_soil(nzb_soil,j,i)
+
+!
+!--             Denominator of the prognostic equation
+                coef_2 = rrtm_lwuflx_dt(0,nzb) + lambda_surface + f_shf / exn
+#else
+
+!
+!--             Numerator of the prognostic equation
+                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
+                          * t_surface(j,i) ** 4 + f_shf * pt_1                 &
+                          + lambda_surface * t_soil(nzb_soil,j,i)
 
 !
@@ -1146 +1179 @@
                 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3               &
                          + lambda_surface + f_shf / exn
-
-             ENDIF
+#endif
+             ENDIF
+
 
              tend = 0.0_wp
@@ -1159 +1193 @@
 !
 !--          Add RK3 term
-             t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)              &
-                                * tt_surface_m(j,i)
-
-!
-!--          Calculate true tendency
-             tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)        &
-                    * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
-!
-!--          Calculate t_surface tendencies for the next Runge-Kutta step
-             IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                IF ( intermediate_timestep_count == 1 )  THEN
-                   tt_surface_m(j,i) = tend
-                ELSEIF ( intermediate_timestep_count <                         &
-                         intermediate_timestep_count_max )  THEN
-                   tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp           &
-                                       * tt_surface_m(j,i)
+             IF ( c_surface /= 0.0_wp )  THEN
+
+                t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)           &
+                                   * tt_surface_m(j,i)
+
+!
+!--             Calculate true tendency
+                tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)     &
+                       * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
+!
+!--             Calculate t_surface tendencies for the next Runge-Kutta step
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( intermediate_timestep_count == 1 )  THEN
+                      tt_surface_m(j,i) = tend
+                   ELSEIF ( intermediate_timestep_count <                      &
+                            intermediate_timestep_count_max )  THEN
+                      tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp        &
+                                          * tt_surface_m(j,i)
+                   ENDIF
                 ENDIF
-             ENDIF
-
-             pt_p(k,j,i) = t_surface_p(j,i) / exn
+
+             ENDIF
+
+!
+!--          In case of fast changes in the skin temperature, it is required to
+!--          update the radiative fluxes in order to keep the solution stable
+             IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp )  THEN
+                force_radiation_call_l = .TRUE.
+             ENDIF
+
+             pt(k,j,i) = t_surface_p(j,i) / exn
+
 !
 !--          Calculate fluxes
-             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb                 &
-                              * t_surface(j,i)**4 - 4.0_wp * sigma_sb          &
-                              * t_surface(j,i)**3 * t_surface_p(j,i)
+#if defined ( __rrtmg )
+             rad_net_l(j,i)   = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb)         &
+                                * t_surface(j,i) - rad_lw_out(nzb,j,i)         &
+                                - rrtm_lwuflx_dt(0,nzb) * t_surface_p(j,i)
+
+             IF ( rrtm_idrv == 1 )  THEN
+                rad_net(j,i) = rad_net_l(j,i)
+                rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i)                      &
+                                      + rrtm_lwuflx_dt(0,nzb)                  &
+                                      * ( t_surface_p(j,i) - t_surface(j,i) )
+             ENDIF
+#else
+             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb             &
+                                * t_surface(j,i)**4 - 4.0_wp * sigma_sb        &
+                                * t_surface(j,i)**3 * t_surface_p(j,i)
+             rad_net(j,i) = rad_net_l(j,i)
+#endif
+
              ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
                               - t_soil(nzb_soil,j,i))
-             shf_eb(j,i)    = - f_shf  * ( pt_p(k+1,j,i) - pt_p(k,j,i) )
+
+             shf_eb(j,i)    = - f_shf * ( pt_1 - pt(k,j,i) )
+
+             shf(j,i) = shf_eb(j,i) / rho_cp
 
              IF ( humidity )  THEN
-                qsws_eb(j,i)  = - f_qsws    * ( q_p(k+1,j,i) - q_s + dq_s_dt   &
+                qsws_eb(j,i)  = - f_qsws    * ( qv_1 - q_s + dq_s_dt           &
                                 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
 
-                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( q_p(k+1,j,i) - q_s       &
+                qsws(j,i) = qsws_eb(j,i) / rho_lv
+
+                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv_1 - q_s               &
                                     + dq_s_dt * t_surface(j,i) - dq_s_dt       &
                                     * t_surface_p(j,i) )
 
-                qsws_soil_eb(j,i) = - f_qsws_soil * ( q_p(k+1,j,i) - q_s       &
+                qsws_soil_eb(j,i) = - f_qsws_soil * ( qv_1 - q_s               &
                                     + dq_s_dt * t_surface(j,i) - dq_s_dt       &
                                     * t_surface_p(j,i) )
 
-                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( q_p(k+1,j,i) - q_s       &
+                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv_1 - q_s               &
                                     + dq_s_dt * t_surface(j,i) - dq_s_dt       &
                                     * t_surface_p(j,i) )
              ENDIF
-
 !
 !--          Calculate the true surface resistance
-             IF ( qsws_eb(j,i) .EQ. 0.0_wp )  THEN
+             IF ( qsws_eb(j,i) == 0.0_wp )  THEN
                 r_s(j,i) = 1.0E10_wp
              ELSE
-                r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - q_s + dq_s_dt           &
+                r_s(j,i) = - rho_lv * ( qv_1 - q_s + dq_s_dt                   &
                            * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
                            / qsws_eb(j,i) - r_a(j,i)
@@ -1228 +1293 @@
 !--                Add precipitation to liquid water reservoir, if possible.
 !--                Otherwise, add the water to bare soil fraction.
-                   IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max )  THEN
+                   IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
                       qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
-                                       + c_veg(j,i) * precipitation_rate(j,i)  &
+                                       + c_veg(j,i) * prr(k,j,i) * hyrho(k)    &
                                        * 0.001_wp * rho_l * l_v
                    ELSE
                       qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
-                                        + c_veg(j,i) * precipitation_rate(j,i) &
+                                        + c_veg(j,i) * prr(k,j,i) * hyrho(k)   &
                                         * 0.001_wp * rho_l * l_v
                    ENDIF
@@ -1241 +1306 @@
 !--                coverage.
                    qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp             &
-                                       - c_veg(j,i)) * precipitation_rate(j,i) &
+                                       - c_veg(j,i)) * prr(k,j,i) * hyrho(k)   &
                                        * 0.001_wp * rho_l * l_v
                 ENDIF
+
 !
 !--             If the air is saturated, check the reservoir water level
-                IF ( q_s .LE. q_p(k+1,j,i))  THEN
+                IF ( qsws_eb(j,i) < 0.0_wp )  THEN
+
 !
 !--                Check if reservoir is full (avoid values > m_liq_eb_max)
@@ -1252 +1319 @@
 !--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
 !--                so that tend is zero and no further check is needed
-                   IF ( m_liq_eb(j,i) .EQ. m_liq_eb_max )  THEN
+                   IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
                       qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
                                            + qsws_liq_eb(j,i)
@@ -1262 +1329 @@
 !--                let the water enter the liquid water reservoir as dew on the
 !--                plant
-                   IF ( qsws_veg_eb(j,i) .LT. 0.0_wp )  THEN
+                   IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
                       qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
                       qsws_veg_eb(j,i) = 0.0_wp
@@ -1301 +1368 @@
        ENDDO
 
+#if defined( __parallel )
+!
+!--    Make a logical OR for all processes. Force radiation call if at
+!--    least one processor
+       IF ( intermediate_timestep_count == intermediate_timestep_count_max-1 )&
+       THEN
+
+          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+          CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
+                              1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
+#else
+          force_radiation_call = force_radiation_call_l
+#endif
+          force_radiation_call_l = .FALSE.
+       ENDIF
+
+
     END SUBROUTINE lsm_energy_balance
 
@@ -1322 +1406 @@
 
        REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !< temp. gamma
-                                               lambda_temp, & !< temp. lambda
-                                               tend           !< tendency
-
-       DO i = nxlg, nxrg    
-          DO j = nysg, nyng
-             DO k = nzb_soil, nzt_soil
+                                                 lambda_temp, & !< temp. lambda
+                                                 tend           !< tendency
+
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             DO  k = nzb_soil, nzt_soil
 !
 !--             Calculate volumetric heat capacity of the soil, taking into 
@@ -1335 +1419 @@
 
 !
-!--             Calculate soil heat conductivity at the center of the soil 
+!--             Calculate soil heat conductivity at the center of the soil
 !--             layers
+                lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *          &
+                               lambda_h_water ** m_soil(k,j,i)
+
                 ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i)))
-                lambda_temp(k) = ke * (lambda_h_sat(j,i) + lambda_h_dry) +     &
+
+                lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +          &
                                  lambda_h_dry
 
@@ -1346 +1434 @@
 !--          Calculate soil heat conductivity (lambda_h) at the _stag level 
 !--          using linear interpolation
-             DO k = nzb_soil, nzt_soil-1
-                  
-                lambda_h(k,j,i) = lambda_temp(k) +                             &
-                                  ( lambda_temp(k+1) - lambda_temp(k) )        &
-                                  * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
-
+             DO  k = nzb_soil, nzt_soil-1
+                lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) * 0.5_wp
              ENDDO
              lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
@@ -1358 +1442 @@
 !--          Prognostic equation for soil temperature t_soil
              tend(:) = 0.0_wp
-             tend(0) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *                    &
+             tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *             &
                        ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)     &
-                         - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil)         &
+                         - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)       &
                          + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
 
-             DO  k = 1, nzt_soil
+
+
+
+             DO  k = nzb_soil+1, nzt_soil
                 tend(k) = (1.0_wp/rho_c_total(k,j,i))                          &
                           * (   lambda_h(k,j,i)                                &
                               * ( t_soil(k+1,j,i) - t_soil(k,j,i) )            &
-                              * ddz_soil(k)                                    &
+                              * ddz_soil(k+1)                                  &
                               - lambda_h(k-1,j,i)                              &
                               * ( t_soil(k,j,i) - t_soil(k-1,j,i) )            &
-                              * ddz_soil(k-1)                                  &
+                              * ddz_soil(k)                                    &
                             ) * ddz_soil_stag(k)
              ENDDO
@@ -1396 +1483 @@
 
 
-             DO k = nzb_soil, nzt_soil
+             DO  k = nzb_soil, nzt_soil
 
 !
@@ -1414 +1501 @@
                               / (n_vg(j,i)-1.0_wp)) - 1.0_wp                   &
                           )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i)
+
 
                    gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +            &
@@ -1438 +1526 @@
 !--             using linear interpolation. To do: replace this with
 !--             ECMWF-IFS Eq. 8.81
-                DO k = nzb_soil, nzt_soil-1
+                DO  k = nzb_soil, nzt_soil-1
                      
-                   lambda_w(k,j,i) = lambda_temp(k) +                          &
-                                     ( lambda_temp(k+1) - lambda_temp(k) )     &
-                                     * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)
-                   gamma_w(k,j,i)  = gamma_temp(k) +                           &
-                                     ( gamma_temp(k+1) - gamma_temp(k) )       &
-                                     * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)                 
+                   lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )     &
+                                     * 0.5_wp
+                   gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )       &
+                                     * 0.5_wp
 
                 ENDDO
@@ -1457 +1543 @@
                    gamma_w(nzt_soil,j,i) = 0.0_wp
                 ELSE
-                   gamma_w(nzt_soil,j,i) = lambda_temp(nzt_soil)
+                   gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
                 ENDIF     
 
@@ -1470 +1556 @@
 
 !
-!--             Calculate the root extraction (ECMWF 7.69, modified so that the
-!--             sum of root_extr = 1). The energy balance solver guarantees a 
-!--             positive transpiration, so that there is no need for an 
-!--             additional check.
-
-                m_total = 0.0_wp
-                DO k = nzb_soil, nzt_soil
-                    IF ( m_soil(k,j,i) .GT. m_wilt(j,i) )  THEN
+!--             Calculate the root extraction (ECMWF 7.69, the sum of root_extr
+!--             = 1). The energy balance solver guarantees a positive
+!--             transpiration, so that there is no need for an additional check.
+                DO  k = nzb_soil, nzt_soil
+                    IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
                        m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
                     ENDIF
                 ENDDO  
 
-                DO k = nzb_soil, nzt_soil 
-                   IF ( m_soil(k,j,i) .GT. m_wilt(j,i) )  THEN
-                      root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
-                   ELSE
-                      root_extr(k) = 0.0_wp
-                   ENDIF
-                ENDDO
+                IF ( m_total > 0.0_wp )  THEN
+                   DO  k = nzb_soil, nzt_soil
+                      IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
+                         root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total
+                      ELSE
+                         root_extr(k) = 0.0_wp
+                      ENDIF
+                   ENDDO
+                ENDIF
 
 !
@@ -1495 +1580 @@
                 tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (                  &
                             m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
-                            * ddz_soil(nzb_soil) - gamma_w(nzb_soil,j,i) - (   &
+                            * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
                             root_extr(nzb_soil) * qsws_veg_eb(j,i)             &
                             + qsws_soil_eb(j,i) ) * drho_l_lv )                &
@@ -1502 +1587 @@
                 DO  k = nzb_soil+1, nzt_soil-1
                    tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)             &
-                               - m_soil(k,j,i) ) * ddz_soil(k) - gamma_w(k,j,i)&
-                               - lambda_w(k-1,j,i) * (m_soil(k,j,i) -          &
-                               m_soil(k-1,j,i)) * ddz_soil(k-1)                &
-                               + gamma_w(k-1,j,i) - (root_extr(k)              &
-                               * qsws_veg_eb(j,i) * drho_l_lv)                 &
+                             - m_soil(k,j,i) ) * ddz_soil(k+1) - gamma_w(k,j,i)&
+                             - lambda_w(k-1,j,i) * (m_soil(k,j,i) -            &
+                             m_soil(k-1,j,i)) * ddz_soil(k)                    &
+                             + gamma_w(k-1,j,i) - (root_extr(k)                &
+                             * qsws_veg_eb(j,i) * drho_l_lv)                   &
                              ) * ddz_soil_stag(k)
 
@@ -1514 +1599 @@
                                         * (m_soil(nzt_soil,j,i)                &
                                         - m_soil(nzt_soil-1,j,i))              &
-                                        * ddz_soil(nzt_soil-1)                 &
+                                        * ddz_soil(nzt_soil)                   &
                                         + gamma_w(nzt_soil-1,j,i) - (          &
                                           root_extr(nzt_soil)                  &
@@ -1573 +1658 @@
        REAL(wp) :: resistance    !< aerodynamic and soil resistance term
 
-       DO i = nxlg, nxrg    
-          DO j = nysg, nyng
+       DO  i = nxlg, nxrg    
+          DO  j = nysg, nyng
              k = nzb_s_inner(j,i)
 
 !
 !--          Calculate water vapour pressure at saturation
-             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
-                                              - 273.16_wp ) / ( t_surface(j,i) &
-                                              - 35.86_wp ) )
+             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface_p(j,i)   &
+                   - 273.16_wp ) / ( t_surface_p(j,i) - 35.86_wp ) )
 
 !
@@ -1587 +1671 @@
              q_s = 0.622_wp * e_s / surface_pressure
 
-
              resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
 
 !
 !--          Calculate specific humidity at surface
-             q_p(k,j,i) = resistance * q_s + (1.0_wp - resistance)             &
-                          * q_p(k+1,j,i)
+             IF ( cloud_physics )  THEN
+                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
+                             * ( q(k+1,j,i) - ql(k+1,j,i) )
+             ELSE
+                q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
+                             * q(k+1,j,i)
+             ENDIF
+
+!
+!--          Update virtual potential temperature
+             vpt(k,j,i) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
 
           ENDDO

palm/trunk/SOURCE/lpm_advec.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-!
+! Renamed prandtl_layer to constant_flux_layer.
 !
 ! Former revisions:
@@ -83 +83 @@
 
     USE control_parameters,                                                    &
-        ONLY:  atmos_ocean_sign, cloud_droplets, dt_3d, dt_3d_reached_l, dz,   &
-               g, kappa, molecular_viscosity, prandtl_layer, topography,       &
+        ONLY:  atmos_ocean_sign, cloud_droplets, constant_flux_layer, dt_3d,   &
+               dt_3d_reached_l, dz, g, kappa, molecular_viscosity, topography, &
                u_gtrans, v_gtrans, simulated_time
 
@@ -224 +224 @@
 !--       First, check if particle is located below first vertical grid level 
 !--       (Prandtl-layer height)
-          IF ( prandtl_layer  .AND.  particles(n)%z < z_p )  THEN
+          IF ( constant_flux_layer  .AND.  particles(n)%z < z_p )  THEN
 !
 !--          Resolved-scale horizontal particle velocity is zero below z0. 
@@ -294 +294 @@
        DO  n = start_index(nb), end_index(nb)
 
-          IF ( prandtl_layer  .AND.  particles(n)%z < z_p )  THEN
+          IF ( constant_flux_layer  .AND.  particles(n)%z < z_p )  THEN
 
              IF ( particles(n)%z < z0_av_global )  THEN

palm/trunk/SOURCE/lpm_exchange_horiz.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-! 
+! Formatting corrections.
 !
 ! Former revisions:
@@ -1208 +1208 @@
              ENDDO
 
-             IF ( n_start .GE. 0 )  THEN
+             IF ( n_start >= 0 )  THEN
                 pack_done = .FALSE.
                 DO  n = n_start, np_old_cell

palm/trunk/SOURCE/lpm_init.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Renamed prandtl_layer to constant_flux_layer.
 !
 ! Former revisions:
@@ -100 +100 @@
 
     USE control_parameters,                                                    &
-        ONLY:  cloud_droplets, current_timestep_number, dz,                    &
-               initializing_actions, message_string, netcdf_data_format,       &
-               ocean, prandtl_layer, simulated_time
+        ONLY:  cloud_droplets, constant_flux_layer, current_timestep_number,   &
+               dz, initializing_actions, message_string, netcdf_data_format,   &
+               ocean, simulated_time
 
     USE dvrp_variables,                                                        &
@@ -295 +295 @@
 !-- To obtain exact height levels of particles, linear interpolation is applied
 !-- (see lpm_advec.f90).
-    IF ( prandtl_layer )  THEN
+    IF ( constant_flux_layer )  THEN
        
        ALLOCATE ( log_z_z0(0:number_of_sublayers) ) 

palm/trunk/SOURCE/microphysics.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-! 
+! Added new routine calc_precipitation_amount. The routine now allows to account
+! for precipitation due to sedimenation of cloud (fog) droplets
 ! 
 ! Former revisions:
@@ -130 +131 @@
     END INTERFACE sedimentation_rain
 
+    INTERFACE calc_precipitation_amount
+       MODULE PROCEDURE calc_precipitation_amount
+       MODULE PROCEDURE calc_precipitation_amount_ij
+    END INTERFACE calc_precipitation_amount
+
+
  CONTAINS
 
@@ -208 +215 @@
 
        IF ( drizzle )  CALL sedimentation_cloud
+
+       IF ( precipitation )  THEN
+          CALL calc_precipitation_amount
+       ENDIF
 
     END SUBROUTINE microphysics_control
@@ -726 +737 @@
 
        USE cloud_parameters,                                                   &
-           ONLY:  eps_sb, hyrho, l_d_cp, nc_const, pt_d_t, sed_qc_const
+           ONLY:  eps_sb, hyrho, l_d_cp, nc_const, prr, pt_d_t, sed_qc_const
 
        USE constants,                                                          &
@@ -732 +743 @@
 
        USE control_parameters,                                                 &
-           ONLY:  dt_do2d_xy, dt_micro, intermediate_timestep_count
+           ONLY:  call_microphysics_at_all_substeps, dt_micro,                 &
+                  intermediate_timestep_count, precipitation
 
        USE cpulog,                                                             &
@@ -741 +753 @@
 
        USE kinds
-       
+
+       USE statistics,                                                         &
+           ONLY:  weight_substep
+    
+
        IMPLICIT NONE
 
@@ -777 +793 @@
                                         pt_d_t(k) * dt_micro
 
+!
+!--             Compute the precipitation rate due to cloud (fog) droplets
+                IF ( precipitation )  THEN
+                   IF ( call_microphysics_at_all_substeps )  THEN
+                      prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)          &
+                                   * weight_substep(intermediate_timestep_count)
+                   ELSE
+                      prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
+                   ENDIF
+                ENDIF
+
              ENDDO
           ENDDO
@@ -799 +826 @@
        USE cloud_parameters,                                                   &
            ONLY:  a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho,        &
-                  limiter_sedimentation, l_d_cp, precipitation_amount, prr,    &
-                  pt_d_t, stp
+                  limiter_sedimentation, l_d_cp, prr,  pt_d_t, stp
 
        USE control_parameters,                                                 &
-           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro,     &
-                  dt_3d, intermediate_timestep_count,                          &
-                  intermediate_timestep_count_max,                             &
-                  precipitation_amount_interval, time_do2d_xy
-
+           ONLY:  call_microphysics_at_all_substeps, dt_micro,                 &
+                  intermediate_timestep_count
        USE cpulog,                                                             &
            ONLY:  cpu_log, log_point_s
@@ -1005 +1028 @@
 !--             Compute the rain rate
                 IF ( call_microphysics_at_all_substeps )  THEN
-                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *            &
-                                weight_substep(intermediate_timestep_count)
+                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
+                                * weight_substep(intermediate_timestep_count)
                 ELSE
-                   prr(k,j,i) = sed_qr(k) / hyrho(k)
+                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
                 ENDIF
 
@@ -1015 +1038 @@
        ENDDO
 
-!
-!--    Precipitation amount
-       IF ( intermediate_timestep_count == intermediate_timestep_count_max     &
-            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                             &
-                   precipitation_amount_interval )  THEN
+       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
+
+    END SUBROUTINE sedimentation_rain
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Computation of the precipitation amount due to gravitational settling of
+!> rain and cloud (fog) droplets
+!------------------------------------------------------------------------------!
+    SUBROUTINE calc_precipitation_amount
+
+       USE cloud_parameters,                                                   &
+           ONLY:  hyrho, precipitation_amount, prr
+
+       USE control_parameters,                                                 &
+           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
+                  intermediate_timestep_count, intermediate_timestep_count_max,&
+                  precipitation_amount_interval, time_do2d_xy
+
+       USE indices,                                                            &
+           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner
+
+       USE kinds
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                          !:
+       INTEGER(iwp) ::  j                          !:
+
+
+       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
+            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
+            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
+       THEN
+
           DO  i = nxl, nxr
              DO  j = nys, nyn
+
                 precipitation_amount(j,i) = precipitation_amount(j,i) +        &
                                             prr(nzb_s_inner(j,i)+1,j,i) *      &
                                             hyrho(nzb_s_inner(j,i)+1) * dt_3d
+
              ENDDO
           ENDDO
        ENDIF
 
-       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
-
-    END SUBROUTINE sedimentation_rain
+    END SUBROUTINE calc_precipitation_amount
 
 
@@ -1121 +1176 @@
 
        IF ( drizzle )  CALL sedimentation_cloud( i,j )
+
+       IF ( precipitation )  THEN
+          CALL calc_precipitation_amount( i,j )
+       ENDIF
 
 !
@@ -1591 +1650 @@
 
        USE cloud_parameters,                                                   &
-           ONLY:  eps_sb, hyrho, l_d_cp, pt_d_t, sed_qc_const
+           ONLY:  eps_sb, hyrho, l_d_cp, prr, pt_d_t, sed_qc_const
 
        USE constants,                                                          &
@@ -1597 +1656 @@
 
        USE control_parameters,                                                 &
-           ONLY:  dt_do2d_xy, dt_micro, intermediate_timestep_count
+           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro,     &
+                  intermediate_timestep_count, precipitation
 
        USE indices,                                                            &
@@ -1604 +1664 @@
        USE kinds
        
+       USE statistics,                                                         &
+           ONLY:  weight_substep
+
        IMPLICIT NONE
 
@@ -1632 +1695 @@
           pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
                                 hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
+
+!
+!--       Compute the precipitation rate of cloud (fog) droplets
+          IF ( precipitation )  THEN
+             IF ( call_microphysics_at_all_substeps )  THEN
+                prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *               &
+                             weight_substep(intermediate_timestep_count)
+             ELSE
+                prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)
+             ENDIF
+          ENDIF
 
        ENDDO
@@ -1837 +1911 @@
 !--       Compute the rain rate
           IF ( call_microphysics_at_all_substeps )  THEN
-             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *                  &
-                          weight_substep(intermediate_timestep_count)
+             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)                    &
+                          * weight_substep(intermediate_timestep_count)
           ELSE
-             prr(k,j,i) = sed_qr(k) / hyrho(k)
+             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
           ENDIF
 
        ENDDO
 
-!
-!--    Precipitation amount
-       IF ( intermediate_timestep_count == intermediate_timestep_count_max     &
-            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                             &
-                   precipitation_amount_interval )  THEN
+    END SUBROUTINE sedimentation_rain_ij
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This subroutine computes the precipitation amount due to gravitational
+!> settling of rain and cloud (fog) droplets
+!------------------------------------------------------------------------------!
+    SUBROUTINE calc_precipitation_amount_ij( i, j )
+
+       USE cloud_parameters,                                                   &
+           ONLY:  hyrho, precipitation_amount, prr
+
+       USE control_parameters,                                                 &
+           ONLY:  call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d,        &
+                  intermediate_timestep_count, intermediate_timestep_count_max,&
+                  precipitation_amount_interval,         &
+                  time_do2d_xy
+
+       USE indices,                                                            &
+           ONLY:  nzb_s_inner
+
+       USE kinds
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                          !:
+       INTEGER(iwp) ::  j                          !:
+
+
+       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
+            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
+            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
+       THEN
 
           precipitation_amount(j,i) = precipitation_amount(j,i) +              &
@@ -1856 +1960 @@
        ENDIF
 
-    END SUBROUTINE sedimentation_rain_ij
-
-    
+    END SUBROUTINE calc_precipitation_amount_ij
+
 !------------------------------------------------------------------------------!
 ! Description:

palm/trunk/SOURCE/modules.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-! 
+! Renamed Obukhov length. Added ol, removed rif. Increased number of profiles 
+! (pr_palm). Changed default values for dissipation_1d to 'detering' and 
+! (mixing_length_1d to 'blackadar'. Added most_method. rif_min and rif_max
+! renamed to zeta_min and zeta_max and new values assigned.
 ! 
 ! Former revisions:
@@ -300 +303 @@
 !
 !
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Definition of all variables
+!>
+!> @todo Add description for each variable
+!------------------------------------------------------------------------------!
+
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
@@ -338 +351 @@
           flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, flux_s_sa,      &
           flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg,                   &
-          mean_inflow_profiles, nrs, nrsws, nrswst, ptnudge, pt_slope_ref,     &
+          mean_inflow_profiles, nrs, nrsws, nrswst,                            &
+          ol,                                                                  & !< Obukhov length
+          ptnudge, pt_slope_ref,                                               &
           qnudge, qs, qsws, qswst, qswst_remote, qrs, qrsws, qrswst, rif,      &
           saswsb, saswst, shf, tnudge, td_lsa_lpt, td_lsa_q, td_sub_lpt,       &
@@ -390 +405 @@
     USE kinds
 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lwp_av, precipitation_rate_av,   &
-                                          qsws_av, shf_av,ts_av, us_av, z0_av, &
-                                          z0h_av
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lwp_av,                & !> Avg. liquid water path
+                                              precipitation_rate_av, & !> Avg. of precipitation rate
+                                              ol_av,                 & !> Avg. of Obukhov length
+                                              qsws_av,               & !> Avg. of surface moisture flux
+                                              shf_av,                & !> Avg. of surface heat flux
+                                              ts_av,                 & !> Avg. of characteristic temperature scale
+                                              us_av,                 & !> Avg. of friction velocity
+                                              z0_av,                 & !> Avg. of roughness length for momentum
+                                              z0h_av                   !> Avg. of roughness length for heat and moisture
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
@@ -541 +562 @@
     CHARACTER (LEN=2)    ::  coupling_char = ''
     CHARACTER (LEN=5)    ::  write_binary = 'false'
-    CHARACTER (LEN=8)    ::  run_date, run_time
+    CHARACTER (LEN=8)    ::  most_method = 'lookup', & !< NAMELIST parameter defining method to be used to calculate Okukhov length, 
+                             run_date,               & !< 
+                             run_time                  !< 
     CHARACTER (LEN=9)    ::  simulated_time_chr
     CHARACTER (LEN=11)   ::  topography_grid_convention = ' '
@@ -563 +586 @@
                              coupling_mode = 'uncoupled', &
                              coupling_mode_remote = 'uncoupled', &
-                             dissipation_1d = 'as_in_3d_model', &
+                             dissipation_1d = 'detering', &
                              fft_method = 'system-specific', &
-                             mixing_length_1d = 'as_in_3d_model', &
+                             mixing_length_1d = 'blackadar', &
                              random_generator = 'numerical-recipes', &
                              reference_state = 'initial_profile', &
@@ -653 +676 @@
                 cloud_top_radiation = .FALSE., &
                 conserve_volume_flow = .FALSE., constant_diffusion = .FALSE., &
+                constant_flux_layer = .TRUE., &
                 constant_heatflux = .TRUE., constant_top_heatflux = .TRUE., &
                 constant_top_momentumflux = .FALSE., &
@@ -679 +703 @@
                 outflow_l = .FALSE., outflow_n = .FALSE., outflow_r = .FALSE., &
                 outflow_s = .FALSE., passive_scalar = .FALSE., &
-                prandtl_layer = .TRUE., &
                 precipitation = .FALSE., &
                 random_heatflux = .FALSE., recycling_yshift = .FALSE.,&
@@ -745 +768 @@
                  recycling_width = 9999999.9_wp, residual_limit = 1.0E-4_wp, &
                  restart_time = 9999999.9_wp, rho_reference, rho_surface, &
-                 rif_max = 1.0_wp, rif_min = -5.0_wp, roughness_length = 0.1_wp, &
+                 roughness_length = 0.1_wp, &
                  sa_surface = 35.0_wp, &
                  simulated_time = 0.0_wp, simulated_time_at_begin, sin_alpha_surface, &
@@ -768 +791 @@
                  vg_surface = 0.0_wp, vpt_reference = 9999999.9_wp, &
                  v_bulk = 0.0_wp, v_gtrans = 0.0_wp, wall_adjustment_factor = 1.8_wp, &
+                 zeta_max = 20.0_wp,    & !< Maximum value of zeta = z/L
+                 zeta_min = -9990.0_wp, & !< Minimum value of zeta = z/L
                  z_max_do2d = -1.0_wp, z0h_factor = 1.0_wp
 
@@ -1149 +1174 @@
              'wpt          ', 'pt(0)        ', 'pt(zp)       ',                &
              'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
-             'mo_L         ', 'q*           ',                                 &
+             'ol           ', 'q*           ',                                 &
              ( 'unknown      ', i9 = 1, dots_max-23 ) /)
 
@@ -1429 +1454 @@
 
     CHARACTER (LEN=40) ::  region(0:9)
-    INTEGER(iwp) ::  pr_palm = 120, statistic_regions = 0
+    INTEGER(iwp) ::  pr_palm = 130, statistic_regions = 0
     INTEGER(iwp) ::  u_max_ijk(3) = -1, v_max_ijk(3) = -1, w_max_ijk(3) = -1
     LOGICAL ::  flow_statistics_called = .FALSE.

palm/trunk/SOURCE/netcdf.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-! 
+! Added output of radiative heating rates for RRTMG. Corrected output of 
+! radiative fluxes
 ! 
 ! Former revisions:
@@ -116 +117 @@
 #if defined( __netcdf )
 
-    USE arrays_3d,                                                              &
+    USE arrays_3d,                                                             &
         ONLY:  zu, zw
 
-    USE constants,                                                              &
+    USE constants,                                                             &
         ONLY:  pi
 
-    USE control_parameters,                                                     &
-        ONLY:  averaging_interval, averaging_interval_pr, averaging_interval_sp,&
-        data_output_pr,  domask,  dopr_n,dopr_time_count, dopts_time_count,     &
-        dots_time_count, dosp_time_count, do2d, do2d_xz_time_count, do3d,       &
-        do2d_yz_time_count, mask_size, do2d_xy_time_count, do3d_time_count,     &
-        domask_time_count, mask_i_global, mask_j_global,mask_k_global,          &
-        message_string, mid, netcdf_data_format, netcdf_precision, ntdim_2d_xy, &
-        ntdim_2d_xz, ntdim_2d_yz, ntdim_3d, nz_do3d, prt_time_count,            &
-        run_description_header, section, simulated_time, topography
-
-    USE grid_variables,                                                         &
+    USE control_parameters,                                                    &
+        ONLY:  averaging_interval, averaging_interval_pr,                      &
+               averaging_interval_sp, data_output_pr,  domask,  dopr_n,        &
+               dopr_time_count, dopts_time_count, dots_time_count,             &
+               dosp_time_count, do2d, do2d_xz_time_count, do3d,                &
+               do2d_yz_time_count, mask_size, do2d_xy_time_count,              &
+               do3d_time_count, domask_time_count, mask_i_global,              &
+               mask_j_global, mask_k_global, message_string, mid,              &
+               netcdf_data_format, netcdf_precision, ntdim_2d_xy,              &
+               ntdim_2d_xz, ntdim_2d_yz, ntdim_3d, nz_do3d, prt_time_count,    &
+              run_description_header, section, simulated_time, topography
+
+    USE grid_variables,                                                        &
         ONLY:  dx, dy, zu_s_inner, zw_w_inner
 
-    USE indices,                                                                &
+    USE indices,                                                               &
         ONLY:  nx, ny, nz ,nzb, nzt
 
@@ -150 +153 @@
     USE pegrid
 
-    USE particle_attributes,                                                    &
+    USE particle_attributes,                                                   &
         ONLY:  maximum_number_of_particles, number_of_particle_groups
 
-    USE profil_parameter,                                                       &
-        ONLY:  crmax, cross_profiles, dopr_index,profile_columns, profile_rows
-
-    USE radiation_model_mod,                                                    &
-        ONLY:  rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
-
-    USE spectrum,                                                               &
+    USE profil_parameter,                                                      &
+        ONLY:  crmax, cross_profiles, dopr_index, profile_columns, profile_rows
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  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
+
+
+    USE spectrum,                                                              &
         ONLY:  comp_spectra_level, data_output_sp, spectra_direction
 
-    USE statistics,                                                             &
+    USE statistics,                                                            &
         ONLY:  hom, statistic_regions
 
@@ -531 +536 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',&
-                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',    &
-                       'rho', 's', 'sa', 'vpt', 'rad_lw_in', 'rad_lw_out', &
-                       'rad_sw_in', 'rad_sw_out' )
+                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
+                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',        &
+                       'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr',            &
+                       'rad_sw_hr', 'rho', 's', 'sa', 'vpt'  )
 
                    grid_x = 'x'
@@ -555 +560 @@
 !
 !--             w grid
-                CASE ( 'w' )
+                CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',   &
+                       'w' )
 
                    grid_x = 'x'
@@ -562 +568 @@
 !
 !--             soil grid
-                CASE ( 't_soil', 'm_soil' )
+                CASE ( 'm_soil', 't_soil' )
 
                    grid_x = 'x'
@@ -1112 +1118 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',   &
-                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho',&
-                       's', 'sa', 'vpt' , 'rad_lw_in', 'rad_lw_out',          &
-                       'rad_sw_in', 'rad_sw_out' )
+                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
+                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho', &
+                       's', 'sa', 'vpt' , 'rad_lw_cs_hr', 'rad_lw_hr',         &
+                       'rad_sw_cs_hr', 'rad_sw_hr' )
 
                    grid_x = 'x'
@@ -1136 +1142 @@
 !
 !--             w grid
-                CASE ( 'w' )
+                CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',   &
+                       'w' )
 
                    grid_x = 'x'
@@ -1143 +1150 @@
 !
 !--             soil grid
-                CASE ( 't_soil', 'm_soil' )
+                CASE ( 'm_soil', 't_soil' )
 
                    grid_x = 'x'
@@ -1765 +1772 @@
 !
 !--                   Most variables are defined on the zu grid
-                      CASE ( 'e_xy', 'lpt_xy', 'nr_xy', 'p_xy', 'pc_xy', 'pr_xy',&
-                             'prr_xy', 'pt_xy', 'q_xy', 'qc_xy', 'ql_xy',        &
-                             'ql_c_xy', 'ql_v_xy', 'ql_vp_xy', 'qr_xy', 'qv_xy', &
-                             'rho_xy', 's_xy', 'sa_xy', 'vpt_xy', 'rad_lw_in_xy',&
-                             'rad_lw_out_xy', 'rad_sw_in_xy', 'rad_sw_out_xy' )
+                      CASE ( 'e_xy', 'lpt_xy', 'nr_xy', 'p_xy', 'pc_xy',       &
+                             'pr_xy', 'prr_xy', 'pt_xy', 'q_xy', 'qc_xy',      &
+                             'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy',        &
+                             'qr_xy', 'qv_xy', 'rad_lw_cs_hr_xy',              &
+                             'rad_lw_hr_xy', 'rad_sw_cs_hr_xy', 'rad_sw_hr_xy',& 
+                             'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
 
                          grid_x = 'x'
@@ -1790 +1798 @@
 !
 !--                   w grid
-                      CASE ( 'w_xy' )
+                      CASE ( 'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy',  &
+                             'rad_sw_out_xy' , 'w_xy' )
 
                          grid_x = 'x'
@@ -1797 +1806 @@
 !
 !--                   soil grid
-                      CASE ( 't_soil_xy', 'm_soil_xy' )
+                      CASE ( 'm_soil_xy', 't_soil_xy' )
                          grid_x = 'x'
                          grid_y = 'y'
@@ -2452 +2461 @@
                           'prr_xz', 'pt_xz', 'q_xz', 'qc_xz', 'ql_xz',         &
                           'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz', 'qv_xz',  &
-                          'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' , 'rad_lw_in_xz',&
-                             'rad_lw_out_xz', 'rad_sw_in_xz', 'rad_sw_out_xz' )
+                          'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',                   &
+                          'rad_sw_cs_hr_xz', 'rad_sw_hr_xz''rho_xz', 's_xz',   &
+                          'sa_xz', 'vpt_xz' )
 
                       grid_x = 'x'
@@ -2474 +2484 @@
 !
 !--                w grid
-                   CASE ( 'w_xz' )
+                   CASE ( 'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz',     &
+                          'rad_sw_out_xz', 'w_xz' )
 
                       grid_x = 'x'
@@ -2482 +2493 @@
 !
 !--                soil grid
-                   CASE ( 't_soil_xz', 'm_soil_xz' )
+                   CASE ( 'm_soil_xz', 't_soil_xz' )
 
                       grid_x = 'x'
@@ -3126 +3137 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_yz', 'lpt_yz', 'nr_yz', 'p_yz', 'pc_yz', 'pr_yz',&
-                          'prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz',        &
-                          'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz', &
-                          'rho_yz', 's_yz', 'sa_yz', 'vpt_yz', 'rad_lw_in_yz',&
-                             'rad_lw_out_yz', 'rad_sw_in_yz', 'rad_sw_out_yz' )
+                   CASE ( 'e_yz', 'lpt_yz', 'nr_yz', 'p_yz', 'pc_yz', 'pr_yz', &
+                          'prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz',         &
+                          'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz',  &
+                          'rad_lw_cs_hr_yz', 'rad_lw_hr_yz',                   &
+                          'rad_sw_cs_hr_yz', 'rad_sw_hr_yz''rho_yz', 's_yz',   &
+                          'sa_yz', 'vpt_yz' )
 
                       grid_x = 'x'
@@ -3151 +3163 @@
 !
 !--                w grid
-                   CASE ( 'w_yz' )
+                   CASE ( 'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz',     &
+                          'rad_sw_out_yz', 'w_yz' )
 
                       grid_x = 'x'
@@ -3158 +3171 @@
 !
 !--                soil grid
-                   CASE ( 't_soil_yz', 'm_soil_yz' )
+                   CASE ( 'm_soil_yz', 't_soil_yz' )
 
                       grid_x = 'x'

palm/trunk/SOURCE/package_parin.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added skip_time_do_lsm, skip_time_do_radiation, and emissivity
 ! 
 ! Former revisions:
@@ -92 +92 @@
 !> This subroutine reads from the NAMELIST file variables controling model
 !> software packages which are used optionally in the run.
+!>
+!> @todo Perform all actions in the respective submodules and remove 
+!>       package_parin
 !------------------------------------------------------------------------------!
  SUBROUTINE package_parin
@@ -117 +120 @@
     USE land_surface_model_mod,                                                &
         ONLY: alpha_vangenuchten, c_surface, canopy_resistance_coefficient,    &
-              conserve_water_content, dewfall, f_shortwave_incoming,           &
-              hydraulic_conductivity, lambda_surface_stable,                   &
-              lambda_surface_unstable,                                         &
-              leaf_area_index, land_surface, l_vangenuchten,                   &
-              min_canopy_resistance, min_soil_resistance, field_capacity,      &
-              residual_moisture, saturation_moisture, wilting_point,           &
-              n_vangenuchten, root_fraction, soil_moisture, soil_temperature,  &
-              soil_type, vegetation_coverage, veg_type, z0_eb, z0h_eb, zs
+              conserve_water_content, dewfall, field_capacity,                 &
+              f_shortwave_incoming, hydraulic_conductivity,                    &
+              lambda_surface_stable, lambda_surface_unstable, leaf_area_index, &
+              land_surface, l_vangenuchten, min_canopy_resistance,             &
+              min_soil_resistance, n_vangenuchten, residual_moisture,          &
+              root_fraction, saturation_moisture, wilting_point,               &
+              skip_time_do_lsm, soil_moisture, soil_temperature, soil_type,    &
+              vegetation_coverage, veg_type, zs, z0_eb, z0h_eb
 
     USE kinds
@@ -152 +155 @@
     USE radiation_model_mod,                                                   &
         ONLY: albedo, albedo_type, albedo_lw_dif, albedo_lw_dir, albedo_sw_dif,&
-              albedo_sw_dir, day_init, constant_albedo, dt_radiation, lambda,  &
-              lw_radiation, net_radiation, radiation, radiation_scheme,        &
-              sw_radiation, time_utc_init
+              albedo_sw_dir, constant_albedo, day_init, dt_radiation,          &
+              emissivity, lambda, lw_radiation, net_radiation, radiation,      &
+              radiation_scheme, skip_time_do_radiation, sw_radiation,          &
+              time_utc_init
                 
 
@@ -200 +204 @@
                                   min_soil_resistance, n_vangenuchten,         &
                                   residual_moisture, root_fraction,            &
-                                  saturation_moisture, soil_moisture,          &
-                                  soil_temperature, soil_type,                 &
+                                  saturation_moisture, skip_time_do_lsm,       &
+                                  soil_moisture, soil_temperature, soil_type,  &
                                   vegetation_coverage, veg_type, wilting_point,& 
                                   zs, z0_eb, z0h_eb
@@ -234 +238 @@
                                   constant_albedo, day_init, dt_radiation,     &
                                   lambda, lw_radiation, net_radiation,         &
-                                  radiation_scheme, sw_radiation, time_utc_init 
+                                  radiation_scheme, skip_time_do_radiation,    &
+                                  sw_radiation, time_utc_init 
 
     NAMELIST /spectra_par/        averaging_interval_sp, comp_spectra_level,   &

palm/trunk/SOURCE/parin.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added parameter most_method. Renamed prandtl_layer to constant_flux_layer.
 ! 
 ! Former revisions:
@@ -185 +185 @@
                cloud_droplets, cloud_physics, cloud_scheme,                    &
                cloud_top_radiation, conserve_volume_flow,                      &
-               conserve_volume_flow_mode, coupling_start_time,                 &
+               conserve_volume_flow_mode, constant_flux_layer,                 &
+               coupling_start_time,                                            &
                create_disturbances, cycle_mg, data_output, data_output_masks,  &
                data_output_pr, data_output_2d_on_each_pe,                      &
@@ -210 +211 @@
                maximum_parallel_io_streams, max_pr_user, message_string,       &
                mg_cycles, mg_switch_to_pe0_level, mixing_length_1d,            &
-               momentum_advec, netcdf_data_format, netcdf_precision, neutral,  &
-               ngsrb, normalizing_region, nsor, nsor_ini, nudging, ocean,      &
+               momentum_advec, most_method, netcdf_data_format,                &
+               netcdf_precision, neutral, ngsrb, normalizing_region, nsor,     &
+               nsor_ini, nudging, ocean,                                       &
                omega, omega_sor, passive_scalar, phi, nz_do3d,                 &
-               prandtl_layer, prandtl_number, precipitation,                   &
+               prandtl_number, precipitation,                                  &
                precipitation_amount_interval, psolver, pt_damping_factor,      &
                pt_damping_width, pt_reference, pt_surface,                     &
@@ -224 +226 @@
                reference_state, residual_limit,                                &
                restart_time, return_addres, return_username,                   &
-               revision, rif_max, rif_min, roughness_length, runnr,            &
+               revision, roughness_length, runnr,                              &
                run_identifier, sa_surface, sa_vertical_gradient,               &
                sa_vertical_gradient_level, scalar_advec,                       &
@@ -246 +248 @@
                vg_vertical_gradient_level, v_bulk, v_profile,                  &
                wall_adjustment, wall_heatflux, wall_humidityflux,              &
-               wall_scalarflux, write_binary, z0h_factor, z_max_do2d
+               wall_scalarflux, write_binary, zeta_max, zeta_min, z0h_factor,  &
+               z_max_do2d
 
     USE cpulog,                                                                &
@@ -274 +277 @@
     USE statistics,                                                            &
         ONLY:  hom, hom_sum, pr_palm, region, statistic_regions
+
 
     IMPLICIT NONE
@@ -291 +295 @@
              cloud_physics, cloud_scheme, cloud_top_radiation,                 &
              collective_wait, conserve_volume_flow,                            &
-             conserve_volume_flow_mode, coupling_start_time,                   &
+             conserve_volume_flow_mode, constant_flux_layer,                   &
+             coupling_start_time,                                              &
              curvature_solution_effects, cycle_mg, damp_level_1d,              &
              dissipation_1d,                                                   &
@@ -306 +311 @@
              loop_optimization, masking_method, mg_cycles,                     &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
-             nc_const, netcdf_precision, neutral, ngsrb,                       &
+             most_method, nc_const, netcdf_precision, neutral, ngsrb,          &
              nsor, nsor_ini, nudging, nx, ny, nz, ocean, omega, omega_sor,     &
-             passive_scalar, phi, prandtl_layer,                               &
+             passive_scalar, phi,                                              &
              prandtl_number, precipitation, psolver, pt_damping_factor,        &
              pt_damping_width, pt_reference, pt_surface,                       &
@@ -318 +323 @@
              recycling_width, recycling_yshift,                                &
              reference_state, residual_limit,                                  &
-             rif_max, rif_min, roughness_length, sa_surface,                   &
+             roughness_length, sa_surface,                                     &
              sa_vertical_gradient, sa_vertical_gradient_level, scalar_advec,   &
              scalar_rayleigh_damping,                                          &
@@ -335 +340 @@
              vg_vertical_gradient_level, v_bulk, v_profile, ventilation_effect,&
              wall_adjustment, wall_heatflux, wall_humidityflux,                &
-             wall_scalarflux, z0h_factor
+             wall_scalarflux, zeta_max, zeta_min, z0h_factor
       
     NAMELIST /d3par/  averaging_interval, averaging_interval_pr,               &

palm/trunk/SOURCE/production_e.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Renamed prandtl_layer to constant_flux_layer.
 ! 
 ! Former revisions:
@@ -78 +78 @@
 ! ------------
 !> Production terms (shear + buoyancy) of the TKE.
-!> @warning The case with prandtl_layer = F and use_surface_fluxes = T is
+!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
 !>          not considered well!
 !------------------------------------------------------------------------------!
@@ -128 +128 @@
 
        USE control_parameters,                                                 &
-           ONLY:  cloud_droplets, cloud_physics, g, humidity, kappa, neutral,  &
-                  ocean, prandtl_layer, pt_reference, rho_reference,           &
-                  use_single_reference_value, use_surface_fluxes,              &
-                  use_top_fluxes
+           ONLY:  cloud_droplets, cloud_physics, constant_flux_layer, g,       &
+                  humidity, kappa, neutral, ocean, pt_reference,               &
+                  rho_reference, use_single_reference_value,                   &
+                  use_surface_fluxes, use_top_fluxes
 
        USE grid_variables,                                                     &
@@ -217 +217 @@
           ENDDO
 
-          IF ( prandtl_layer )  THEN
+          IF ( constant_flux_layer )  THEN
 
 !
@@ -734 +734 @@
 
        USE control_parameters,                                                 &
-           ONLY:  cloud_droplets, cloud_physics, g, humidity, kappa, neutral,  &
-                  ocean, prandtl_layer, pt_reference, rho_reference,           &
-                  topography, use_single_reference_value, use_surface_fluxes,  &
-                  use_top_fluxes
+           ONLY:  cloud_droplets, cloud_physics, constant_flux_layer, g,       &
+                  humidity, kappa, neutral, ocean, pt_reference,               &
+                  rho_reference, topography, use_single_reference_value,       &
+                  use_surface_fluxes, use_top_fluxes
 
        USE grid_variables,                                                     &
@@ -832 +832 @@
        ENDDO
 
-       IF ( prandtl_layer )  THEN
+       IF ( constant_flux_layer )  THEN
 
 !
@@ -1402 +1402 @@
 
        USE control_parameters,                                                 &
-           ONLY:  cloud_droplets, cloud_physics, g, humidity, kappa, neutral,  &
-                  ocean, prandtl_layer, pt_reference, rho_reference,           &
-                  use_single_reference_value, use_surface_fluxes,              &
-                  use_top_fluxes
+           ONLY:  cloud_droplets, cloud_physics, constant_flux_layer, g,       &
+                  humidity, kappa, neutral, ocean, pt_reference,               &
+                  rho_reference, use_single_reference_value,                   &
+                  use_surface_fluxes, use_top_fluxes
 
        USE grid_variables,                                                     &
@@ -1473 +1473 @@
        ENDDO
 
-       IF ( prandtl_layer )  THEN
+       IF ( constant_flux_layer )  THEN
 
           IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) )  THEN
@@ -1924 +1924 @@
 
        USE control_parameters,                                                 &
-           ONLY:  kappa, prandtl_layer
+           ONLY:  constant_flux_layer, kappa
 
        USE indices,                                                            &
@@ -1937 +1937 @@
        INTEGER(iwp) ::  kv  !<
 
-       IF ( prandtl_layer )  THEN
+       IF ( constant_flux_layer )  THEN
 
           IF ( first_call )  THEN

palm/trunk/SOURCE/prognostic_equations.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-! 
+! Added optional model spin-up without radiation / land surface model calls.
+! Formatting corrections.
 ! 
 ! Former revisions:
@@ -286 +287 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  radiation, radiation_scheme, radiation_tendency
+        ONLY:  radiation, radiation_scheme, radiation_tendency,                &
+               skip_time_do_radiation
 
     USE statistics,                                                            &
@@ -410 +412 @@
 !
 !--          Nudging
-             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' ) 
+             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' )
 
              CALL user_actions( i, j, 'u-tendency' )
@@ -469 +471 @@
 !
 !--          Nudging
-             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' ) 
+             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' )
 
              CALL user_actions( i, j, 'v-tendency' )
@@ -610 +612 @@
 !
 !--          If required, add tendency due to radiative heating/cooling
-             IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
+             IF ( radiation .AND. radiation_scheme == 'rrtmg' .AND.            &
+                  simulated_time > skip_time_do_radiation )  THEN
                 CALL radiation_tendency ( i, j, tend )
              ENDIF
@@ -718 +721 @@
              ENDIF
              CALL diffusion_s( i, j, q, qsws, qswst, wall_qflux )
-     
+
 !
 !--          If required compute decrease of total water content due to
@@ -724 +727 @@
              IF ( cloud_physics  .AND.  icloud_scheme == 1  .AND.              &
                   precipitation )  THEN
-                CALL calc_precipitation( i, j )             
-             ENDIF
+                CALL calc_precipitation( i, j )
+             ENDIF
+
 !
 !--          Sink or source of scalar concentration due to canopy elements
@@ -1009 +1013 @@
 !
 !-- Nudging
-    IF ( nudging )  CALL nudge( simulated_time, 'u' ) 
+    IF ( nudging )  CALL nudge( simulated_time, 'u' )
 
     CALL user_actions( 'u-tendency' )
@@ -1085 +1089 @@
 !
 !-- Nudging
-    IF ( nudging )  CALL nudge( simulated_time, 'v' ) 
+    IF ( nudging )  CALL nudge( simulated_time, 'v' )
 
     CALL user_actions( 'v-tendency' )
@@ -1274 +1278 @@
 !
 !--    If required, add tendency due to radiative heating/cooling
-       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
-          CALL radiation_tendency ( tend )
+       IF ( radiation .AND. radiation_scheme == 'rrtmg' .AND.                  &
+            simulated_time > skip_time_do_radiation )  THEN
+            CALL radiation_tendency ( tend )
        ENDIF
 
@@ -1859 +1864 @@
 !
 !-- Nudging
-    IF ( nudging )  CALL nudge( simulated_time, 'u' ) 
+    IF ( nudging )  CALL nudge( simulated_time, 'u' )
 
     CALL user_actions( 'u-tendency' )
@@ -1926 +1931 @@
 !
 !-- Nudging
-    IF ( nudging )  CALL nudge( simulated_time, 'v' ) 
+    IF ( nudging )  CALL nudge( simulated_time, 'v' )
 
     CALL user_actions( 'v-tendency' )
@@ -2096 +2101 @@
        ENDIF
 
+       IF ( radiation .AND. radiation_scheme == 'rrtmg' .AND.                  &
+            simulated_time > skip_time_do_radiation )  THEN
+            CALL radiation_tendency ( tend )
+       ENDIF
+
        CALL user_actions( 'pt-tendency' )
 

palm/trunk/SOURCE/radiation_model.f90

@@ -19 +19 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
+! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
+! Added output of radiative heating rates.
 ! 
 ! Former revisions:
@@ -61 +63 @@
 !> @todo Output of other rrtm arrays (such as volume mixing ratios)
 !> @todo Adapt for use with topography
-!> @todo Remove 3D dummy arrays (such as clear-sky output)
 !>
 !> @note Many variables have a leading dummy dimension (0:0) in order to
@@ -70 +71 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  hyp, pt, q, ql, zw
+        ONLY:  dzw, hyp, pt, q, ql, zw
 
     USE cloud_parameters,                                                      &
-        ONLY:  cp, l_v, nc_const, rho_l, sigma_gc  
+        ONLY:  cp, l_d_cp, nc_const, rho_l, sigma_gc  
 
     USE constants,                                                             &
@@ -80 +81 @@
     USE control_parameters,                                                    &
         ONLY:  cloud_droplets, cloud_physics, g, initializing_actions,         &
-               large_scale_forcing, lsf_surf, phi, pt_surface,                 &
+               large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface,    &
                surface_pressure, time_since_reference_point
 
@@ -153 +154 @@
 
 
-    LOGICAL ::  constant_albedo = .FALSE., & !< flag parameter indicating whether the albedo may change depending on zenith
-                lw_radiation = .TRUE.,     & !< flag parameter indicing whether longwave radiation shall be calculated
-                radiation = .FALSE.,       & !< flag parameter indicating whether the radiation model is used
-                sun_up    = .TRUE.,        & !< flag parameter indicating whether the sun is up or down
-                sw_radiation = .TRUE.        !< flag parameter indicing whether shortwave radiation shall be calculated
-
-    REAL(wp), PARAMETER :: sigma_sb       = 5.67E-8_wp, & !< Stefan-Boltzmann constant
-                           solar_constant = 1368.0_wp     !< solar constant at top of atmosphere
- 
-    REAL(wp) :: albedo = 9999999.9_wp,        & !< NAMELIST alpha
-                albedo_lw_dif = 9999999.9_wp, & !< NAMELIST aldif
-                albedo_lw_dir = 9999999.9_wp, & !< NAMELIST aldir
-                albedo_sw_dif = 9999999.9_wp, & !< NAMELIST asdif
-                albedo_sw_dir = 9999999.9_wp, & !< NAMELIST asdir
-                dt_radiation = 0.0_wp,        & !< radiation model timestep
-                emissivity = 0.95_wp,         & !< NAMELIST surface emissivity
-                exn,                          & !< Exner function
-                lon = 0.0_wp,                 & !< longitude in radians
-                lat = 0.0_wp,                 & !< latitude in radians
-                decl_1,                       & !< declination coef. 1
-                decl_2,                       & !< declination coef. 2
-                decl_3,                       & !< declination coef. 3
-                time_utc,                     & !< current time in UTC
-                time_utc_init = 43200.0_wp,   & !< UTC time at model start (noon)
-                lambda = 0.0_wp,              & !< longitude in degrees
-                net_radiation = 0.0_wp,       & !< net radiation at surface
-                time_radiation = 0.0_wp,      & !< time since last call of radiation code
-                sky_trans                       !< sky transmissivity
-
+    LOGICAL ::  constant_albedo = .FALSE.,       & !< flag parameter indicating whether the albedo may change depending on zenith
+                force_radiation_call = .FALSE.,  & !< flag parameter for unscheduled radiation calls
+                lw_radiation = .TRUE.,           & !< flag parameter indicing whether longwave radiation shall be calculated
+                radiation = .FALSE.,             & !< flag parameter indicating whether the radiation model is used
+                sun_up    = .TRUE.,              & !< flag parameter indicating whether the sun is up or down
+                sw_radiation = .TRUE.              !< flag parameter indicing whether shortwave radiation shall be calculated
+
+
+    REAL(wp), PARAMETER :: d_seconds_hour  = 0.000277777777778_wp,  & !< inverse of seconds per hour (1/3600)
+                           d_hours_day    = 0.0416666666667_wp,     & !< inverse of hours per day (1/24)
+                           sigma_sb       = 5.67037321E-8_wp,       & !< Stefan-Boltzmann constant
+                           solar_constant = 1368.0_wp                 !< solar constant at top of atmosphere
+
+    REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
+                albedo_lw_dif = 9999999.9_wp,    & !< NAMELIST aldif
+                albedo_lw_dir = 9999999.9_wp,    & !< NAMELIST aldir
+                albedo_sw_dif = 9999999.9_wp,    & !< NAMELIST asdif
+                albedo_sw_dir = 9999999.9_wp,    & !< NAMELIST asdir
+                decl_1,                          & !< declination coef. 1
+                decl_2,                          & !< declination coef. 2
+                decl_3,                          & !< declination coef. 3
+                dt_radiation = 0.0_wp,           & !< radiation model timestep
+                emissivity = 0.98_wp,            & !< NAMELIST surface emissivity
+                lambda = 0.0_wp,                 & !< longitude in degrees
+                lon = 0.0_wp,                    & !< longitude in radians
+                lat = 0.0_wp,                    & !< latitude in radians
+                net_radiation = 0.0_wp,          & !< net radiation at surface
+                skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
+                sky_trans,                       & !< sky transmissivity
+                time_radiation = 0.0_wp,         & !< time since last call of radiation code
+                time_utc,                        & !< current time in UTC
+                time_utc_init = 43200.0_wp         !< UTC time at model start (noon)
 
     REAL(wp), DIMENSION(0:0) ::  zenith        !< solar zenith angle
@@ -213 +217 @@
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
+                        rad_lw_cs_hr,                  & !< longwave clear sky radiation heating rate (K/s)
+                        rad_lw_cs_hr_av,               & !< average of rad_lw_cs_hr
+                        rad_lw_hr,                     & !< longwave radiation heating rate (K/s)
+                        rad_lw_hr_av,                  & !< average of rad_sw_hr
+                        rad_lw_in,                     & !< incoming longwave radiation (W/m2)
+                        rad_lw_in_av,                  & !< average of rad_lw_in
+                        rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
+                        rad_lw_out_av,                 & !< average of rad_lw_out
+                        rad_sw_cs_hr,                  & !< shortwave clear sky radiation heating rate (K/s)
+                        rad_sw_cs_hr_av,               & !< average of rad_sw_cs_hr
+                        rad_sw_hr,                     & !< shortwave radiation heating rate (K/s)
+                        rad_sw_hr_av,                  & !< average of rad_sw_hr
                         rad_sw_in,                     & !< incoming shortwave radiation (W/m2)
                         rad_sw_in_av,                  & !< average of rad_sw_in
                         rad_sw_out,                    & !< outgoing shortwave radiation (W/m2)
-                        rad_sw_out_av,                 & !< average of rad_sw_out
-                        rad_lw_in,                     & !< incoming longwave radiation (W/m2)
-                        rad_lw_in_av,                  & !< average of rad_lw_in
-                        rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
-                        rad_lw_out_av                    !< average of rad_lw_out
+                        rad_sw_out_av                    !< average of rad_sw_out
+
 
 !
@@ -244 +257 @@
                     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 = 0        !< longwave upward flux calculation option (0,1)
+                    rrtm_idrv = 1        !< longwave upward flux calculation option (0,1)
 
     LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
 
-    REAL(wp), PARAMETER :: mol_weight_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
+    REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
 
     REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
@@ -255 +268 @@
                                            t_snd          !< actual temperature from sounding data (hPa)
 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif,         & !< longwave diffuse albedo solar angle of 60°
-                                             aldir,         & !< longwave direct albedo solar angle of 60°
-                                             asdif,         & !< shortwave diffuse albedo solar angle of 60°
-                                             asdir,         & !< shortwave direct albedo solar angle of 60°
-                                             rrtm_ccl4vmr,  & !< CCL4 volume mixing ratio (g/mol)
-                                             rrtm_cfc11vmr, & !< CFC11 volume mixing ratio (g/mol)
-                                             rrtm_cfc12vmr, & !< CFC12 volume mixing ratio (g/mol)
-                                             rrtm_cfc22vmr, & !< CFC22 volume mixing ratio (g/mol)
-                                             rrtm_ch4vmr,   & !< CH4 volume mixing ratio
-                                             rrtm_cicewp,   & !< in-cloud ice water path (g/m²)
-                                             rrtm_cldfr,    & !< cloud fraction (0,1)
-                                             rrtm_cliqwp,   & !< in-cloud liquid water path (g/m²)
-                                             rrtm_co2vmr,   & !< CO2 volume mixing ratio (g/mol)
-                                             rrtm_emis,     & !< surface emissivity (0-1)    
-                                             rrtm_h2ovmr,   & !< H2O volume mixing ratio
-                                             rrtm_n2ovmr,   & !< N2O volume mixing ratio
-                                             rrtm_o2vmr,    & !< O2 volume mixing ratio
-                                             rrtm_o3vmr,    & !< O3 volume mixing ratio
-                                             rrtm_play,     & !< pressure layers (hPa, zu-grid)
-                                             rrtm_plev,     & !< pressure layers (hPa, zw-grid)
-                                             rrtm_reice,    & !< cloud ice effective radius (microns)
-                                             rrtm_reliq,    & !< cloud water drop effective radius (microns)
-                                             rrtm_tlay,     & !< actual temperature (K, zu-grid)
-                                             rrtm_tlev,     & !< actual temperature (K, zw-grid)
-                                             rrtm_lwdflx,   & !< RRTM output of incoming longwave radiation flux (W/m2)
-                                             rrtm_lwuflx,   & !< RRTM output of outgoing longwave radiation flux (W/m2)
-                                             rrtm_lwhr,     & !< RRTM output of longwave radiation heating rate (K/d)
-                                             rrtm_lwuflxc,  & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
-                                             rrtm_lwdflxc,  & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2) 
-                                             rrtm_lwhrc,    & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
-                                             rrtm_swdflx,   & !< RRTM output of incoming shortwave radiation flux (W/m2)
-                                             rrtm_swuflx,   & !< RRTM output of outgoing shortwave radiation flux (W/m2)
-                                             rrtm_swhr,     & !< RRTM output of shortwave radiation heating rate (K/d)
-                                             rrtm_swuflxc,  & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
-                                             rrtm_swdflxc,  & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2) 
-                                             rrtm_swhrc       !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif,          & !< longwave diffuse albedo solar angle of 60°
+                                             aldir,          & !< longwave direct albedo solar angle of 60°
+                                             asdif,          & !< shortwave diffuse albedo solar angle of 60°
+                                             asdir,          & !< shortwave direct albedo solar angle of 60°
+                                             rrtm_ccl4vmr,   & !< CCL4 volume mixing ratio (g/mol)
+                                             rrtm_cfc11vmr,  & !< CFC11 volume mixing ratio (g/mol)
+                                             rrtm_cfc12vmr,  & !< CFC12 volume mixing ratio (g/mol)
+                                             rrtm_cfc22vmr,  & !< CFC22 volume mixing ratio (g/mol)
+                                             rrtm_ch4vmr,    & !< CH4 volume mixing ratio
+                                             rrtm_cicewp,    & !< in-cloud ice water path (g/m²)
+                                             rrtm_cldfr,     & !< cloud fraction (0,1)
+                                             rrtm_cliqwp,    & !< in-cloud liquid water path (g/m²)
+                                             rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
+                                             rrtm_emis,      & !< surface emissivity (0-1)    
+                                             rrtm_h2ovmr,    & !< H2O volume mixing ratio
+                                             rrtm_n2ovmr,    & !< N2O volume mixing ratio
+                                             rrtm_o2vmr,     & !< O2 volume mixing ratio
+                                             rrtm_o3vmr,     & !< O3 volume mixing ratio
+                                             rrtm_play,      & !< pressure layers (hPa, zu-grid)
+                                             rrtm_plev,      & !< pressure layers (hPa, zw-grid)
+                                             rrtm_reice,     & !< cloud ice effective radius (microns)
+                                             rrtm_reliq,     & !< cloud water drop effective radius (microns)
+                                             rrtm_tlay,      & !< actual temperature (K, zu-grid)
+                                             rrtm_tlev,      & !< actual temperature (K, zw-grid)
+                                             rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
+                                             rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2) 
+                                             rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
+                                             rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
+                                             rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
+                                             rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
+                                             rrtm_lwhr,      & !< RRTM output of longwave radiation heating rate (K/d)
+                                             rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
+                                             rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
+                                             rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2) 
+                                             rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
+                                             rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
+                                             rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
+                                             rrtm_swhrc        !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 
 
 !
@@ -296 +311 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !< incoming clear sky longwave radiation (W/m2) (not used)
                                                 rad_lw_cs_out,  & !< outgoing clear sky longwave radiation (W/m2) (not used)
-                                                rad_lw_cs_hr,   & !< longwave clear sky radiation heating rate (K/d) (not used)
-                                                rad_lw_hr,      & !< longwave radiation heating rate (K/d)
                                                 rad_sw_cs_in,   & !< incoming clear sky shortwave radiation (W/m2) (not used)
                                                 rad_sw_cs_out,  & !< outgoing clear sky shortwave radiation (W/m2) (not used)
-                                                rad_sw_cs_hr,   & !< shortwave clear sky radiation heating rate (K/d) (not used)
-                                                rad_sw_hr,      & !< shortwave radiation heating rate (K/d)
                                                 rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
                                                 rrtm_aldir,     & !< surface albedo for longwave direct radiation
@@ -316 +327 @@
                                                 rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
                                                 rrtm_sw_ecaer     !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
+
 #endif
 
@@ -341 +353 @@
     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, 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_sw_in, rad_sw_in_av,   &
-           rad_sw_out, rad_sw_out_av, sigma_sb, sw_radiation, time_radiation,  &
-           time_utc_init 
+           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_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
+
 
 #if defined ( __rrtmg )
-    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
+    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv,          &
+           rrtm_lwuflx_dt
 #endif
 
@@ -483 +499 @@
 !
 !--       Allocate 3d arrays of radiative fluxes and heating rates
-          ALLOCATE ( rad_sw_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
-          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-          ALLOCATE ( rad_sw_cs_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
-
           IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
              ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
              rad_sw_in = 0.0_wp
-
           ENDIF
 
           IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
              ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-             rad_sw_out = 0.0_wp
           ENDIF
 
           IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
              ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             rad_sw_out = 0.0_wp
           ENDIF
 
@@ -507 +517 @@
           ENDIF
 
-          ALLOCATE ( rad_lw_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
-          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-          ALLOCATE ( rad_lw_cs_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
+          IF ( .NOT. ALLOCATED ( rad_sw_hr ) )  THEN
+             ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             rad_sw_hr = 0.0_wp
+          ENDIF
+
+          IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) )  THEN
+             ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             rad_sw_hr_av = 0.0_wp
+          ENDIF
+
+          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) )  THEN
+             ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             rad_sw_cs_hr = 0.0_wp
+          ENDIF
+
+          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) )  THEN
+             ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             rad_sw_cs_hr_av = 0.0_wp
+          ENDIF
 
           IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
@@ -530 +555 @@
           ENDIF
 
-          rad_sw_hr     = 0.0_wp
+          IF ( .NOT. ALLOCATED ( rad_lw_hr ) )  THEN
+             ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             rad_lw_hr = 0.0_wp
+          ENDIF
+
+          IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) )  THEN
+             ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             rad_lw_hr_av = 0.0_wp
+          ENDIF
+
+          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) )  THEN
+             ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             rad_lw_cs_hr = 0.0_wp
+          ENDIF
+
+          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) )  THEN
+             ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             rad_lw_cs_hr_av = 0.0_wp
+          ENDIF
+
+          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
           rad_sw_cs_in  = 0.0_wp
           rad_sw_cs_out = 0.0_wp
-          rad_sw_cs_hr  = 0.0_wp
-
-          rad_lw_hr     = 0.0_wp
+
+          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
           rad_lw_cs_in  = 0.0_wp
           rad_lw_cs_out = 0.0_wp
-          rad_lw_cs_hr  = 0.0_wp
 
 !
@@ -578 +623 @@
 !--    Add timeseries for radiation model
        dots_rad = dots_num + 1
-       dots_label(dots_num+1) = "rad_net"
-       dots_label(dots_num+2) = "rad_lw_in"
-       dots_label(dots_num+3) = "rad_lw_out"
-       dots_label(dots_num+4) = "rad_sw_in"
-       dots_label(dots_num+5) = "rad_sw_out"
-       dots_unit(dots_num+1:dots_num+5) = "W/m2"
-       dots_num  = dots_num + 5
+       dots_num = dots_num + 5
+
+       dots_label(dots_rad+1) = "rad_net"
+       dots_label(dots_rad+2) = "rad_lw_in"
+       dots_label(dots_rad+3) = "rad_lw_out"
+       dots_label(dots_rad+4) = "rad_sw_in"
+       dots_label(dots_rad+5) = "rad_sw_out"
+       dots_unit(dots_rad:dots_rad+4) = "W/m2"
 
 !
 !--    Output of albedos is only required for RRTMG
        IF ( radiation_scheme == 'rrtmg' )  THEN
-          dots_label(dots_num+1) = "rrtm_aldif"
-          dots_label(dots_num+2) = "rrtm_aldir"
-          dots_label(dots_num+3) = "rrtm_asdif"
-          dots_label(dots_num+4) = "rrtm_asdir"
-          dots_unit(dots_num+1:dots_num+4) = ""
           dots_num  = dots_num + 4
+          dots_label(dots_rad+5) = "rrtm_aldif"
+          dots_label(dots_rad+6) = "rrtm_aldir"
+          dots_label(dots_rad+7) = "rrtm_asdif"
+          dots_label(dots_rad+8) = "rrtm_asdir"
+          dots_unit(dots_num+5:dots_num+8) = ""
+
        ENDIF
 
@@ -624 +671 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) :: i, j, k
+       INTEGER(iwp) :: i, j, k   !< loop indices
+       REAL(wp)     :: exn,   &  !< Exner functions at surface
+                       exn_1, &  !< Exner functions at first grid level
+                       pt_1      !< potential temperature at first grid level
 
 !
@@ -643 +693 @@
           DO j = nys, nyn
              k = nzb_s_inner(j,i)
+
+             exn_1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
+
              rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
              rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
-             rad_lw_out(0,j,i) = sigma_sb * (pt(k,j,i) * exn)**4
-             rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn)**4
-             rad_net(j,i)      = rad_sw_in(0,j,i) - rad_sw_out(0,j,i)          &
-                                 + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
+             rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
+
+             IF ( cloud_physics )  THEN
+                pt_1 = pt(k+1,j,i) + l_d_cp / exn_1 * ql(k+1,j,i)
+                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt_1 * exn_1)**4
+             ELSE
+                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn_1)**4
+             ENDIF
+
+             rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i)               &
+                            + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
 
           ENDDO
@@ -682 +742 @@
 #if defined ( __rrtmg )
 
-       INTEGER(iwp) :: i, j, k, n          !< 
-
-       REAL(wp)     ::  s_r2      !< weighted sum over all droplets with r^2
-       REAL(wp)     ::  s_r3      !< weighted sum over all droplets with r^3
+       INTEGER(iwp) :: i, j, k, n !< loop indices
+
+       REAL(wp)     ::  s_r2, &   !< weighted sum over all droplets with r^2
+                        s_r3      !< weighted sum over all droplets with r^3
 
 !
@@ -703 +763 @@
 !--    profile. nzt_rad might be modified by these calls and all required arrays
 !--    will then be re-allocated
-       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
+       IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
           CALL read_sounding_data
           CALL read_trace_gas_data
@@ -714 +774 @@
 !
 !--          Prepare profiles of temperature and H2O volume mixing ratio
-             rrtm_tlev(0,nzb+1) = pt(nzb,j,i)
+             rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure             &
+                                                  / 1000.0_wp )**0.286_wp
 
              DO k = nzb+1, nzt+1
                 rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp         &
-                                 )**0.286_wp + ( l_v / cp ) * ql(k,j,i)
-                rrtm_h2ovmr(0,k) = mol_weight_air_d_wv * (q(k,j,i) - ql(k,j,i))
+                                 )**0.286_wp + l_d_cp * ql(k,j,i)
+                rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
 
              ENDDO
@@ -754 +815 @@
              rrtm_reliq  = 0.0_wp
              rrtm_cliqwp = 0.0_wp
+             rrtm_icld   = 0
 
              DO k = nzb+1, nzt+1
-                rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *                      &
-                                  (rrtm_plev(0,k) - rrtm_plev(0,k+1)) * 100.0_wp / g 
-
-                IF ( rrtm_cliqwp(0,k) .GT. 0 )  THEN
+                rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *                    &
+                                    (rrtm_plev(0,k) - rrtm_plev(0,k+1))        &
+                                    * 100.0_wp / g 
+
+                IF ( rrtm_cliqwp(0,k) > 0.0_wp )  THEN
                    rrtm_cldfr(0,k) = 1.0_wp
-                   rrtm_icld = 1
+                   IF ( rrtm_icld == 0 )  rrtm_icld = 1
 
 !
 !--                Calculate cloud droplet effective radius
                    IF ( cloud_physics )  THEN
-                      rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)          &
-                                      / ( 4.0_wp * pi * nc_const * rho_l )     &
-                                      )**0.33333333333333_wp                   &
-                                     * EXP( LOG( sigma_gc )**2 )
+                      rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)        &
+                                        * rho_surface                          &
+                                        / ( 4.0_wp * pi * nc_const * rho_l )   &
+                                        )**0.33333333333333_wp                 &
+                                        * EXP( LOG( sigma_gc )**2 )
 
                    ELSEIF ( cloud_droplets )  THEN
@@ -794 +858 @@
 !
 !--                Limit effective radius
-                   IF ( rrtm_reliq(0,k) .GT. 0.0_wp )  THEN
+                   IF ( rrtm_reliq(0,k) > 0.0_wp )  THEN
                       rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
                       rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
                   ENDIF
-                ELSE
-                   rrtm_icld = 0
                 ENDIF
              ENDDO
@@ -817 +879 @@
                rrtm_reliq      , rrtm_lw_tauaer,                               &
                rrtm_lwuflx     , rrtm_lwdflx  , rrtm_lwhr  ,                   &
-               rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc )
-
+               rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc ,                   &
+               rrtm_lwuflx_dt  ,  rrtm_lwuflxc_dt )
+
+!
+!--             Save fluxes
                 DO k = nzb, nzt+1
                    rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k)
@@ -824 +889 @@
                 ENDDO
 
+!
+!--             Save heating rates (convert from K/d to K/h)
+                DO k = nzb+1, nzt+1
+                   rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k)  * d_hours_day
+                   rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k) * d_hours_day
+                ENDDO
 
              ENDIF
@@ -841 +912 @@
                rrtm_swuflxc    , rrtm_swdflxc , rrtm_swhrc )
   
+!
+!--             Save fluxes
                 DO k = nzb, nzt+1
                    rad_sw_in(k,j,i)  = rrtm_swdflx(0,k)
                    rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
                 ENDDO
+
+!
+!--             Save heating rates (convert from K/d to K/s)
+                DO k = nzb+1, nzt+1
+                   rad_sw_hr(k,j,i)     = rrtm_swhr(0,k)  * d_hours_day
+                   rad_sw_cs_hr(k,j,i)  = rrtm_swhrc(0,k) * d_hours_day
+                ENDDO
+
              ENDIF
 
@@ -857 +938 @@
        CALL exchange_horiz( rad_lw_in,  nbgp )
        CALL exchange_horiz( rad_lw_out, nbgp )
+       CALL exchange_horiz( rad_lw_hr,    nbgp )
+       CALL exchange_horiz( rad_lw_cs_hr, nbgp )
+
        CALL exchange_horiz( rad_sw_in,  nbgp )
        CALL exchange_horiz( rad_sw_out, nbgp ) 
+       CALL exchange_horiz( rad_sw_hr,    nbgp )
+       CALL exchange_horiz( rad_sw_cs_hr, nbgp )
+
        CALL exchange_horiz_2d( rad_net, nbgp )
 #endif
@@ -897 +984 @@
 !
 !--    Check if the sun is up (otheriwse shortwave calculations can be skipped)
-       IF ( zenith(0) .GT. 0.0_wp )  THEN
+       IF ( zenith(0) > 0.0_wp )  THEN
           sun_up = .TRUE.
        ELSE
@@ -996 +1083 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) :: id, id_dim_zrad, id_var, k, nz_snd, nz_snd_start, nz_snd_end
-
-       REAL(wp) :: t_surface
-
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp_snd_tmp, t_snd_tmp
+       INTEGER(iwp) :: id,           & !< NetCDF id of input file
+                       id_dim_zrad,  & !< pressure level id in the NetCDF file
+                       id_var,       & !< NetCDF variable id
+                       k,            & !< loop index
+                       nz_snd,       & !< number of vertical levels in the sounding data
+                       nz_snd_start, & !< start vertical index for sounding data to be used
+                       nz_snd_end      !< end vertical index for souding data to be used
+
+       REAL(wp) :: t_surface           !< actual surface temperature
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
+                                               t_snd_tmp      !< temporary temperature profile (sounding)
 
 !
@@ -1014 +1108 @@
           DEALLOCATE ( rrtm_tlay )
           DEALLOCATE ( rrtm_tlev )
+
           DEALLOCATE ( rrtm_h2ovmr )
           DEALLOCATE ( rrtm_cicewp )
@@ -1022 +1117 @@
           DEALLOCATE ( rrtm_lw_taucld )
           DEALLOCATE ( rrtm_lw_tauaer )
+
           DEALLOCATE ( rrtm_lwdflx  )
+          DEALLOCATE ( rrtm_lwdflxc )
           DEALLOCATE ( rrtm_lwuflx  )
+          DEALLOCATE ( rrtm_lwuflxc )
+          DEALLOCATE ( rrtm_lwuflx_dt )
+          DEALLOCATE ( rrtm_lwuflxc_dt )
           DEALLOCATE ( rrtm_lwhr  )
-          DEALLOCATE ( rrtm_lwuflxc )
-          DEALLOCATE ( rrtm_lwdflxc )
           DEALLOCATE ( rrtm_lwhrc )
+
           DEALLOCATE ( rrtm_sw_taucld )
           DEALLOCATE ( rrtm_sw_ssacld )
@@ -1035 +1134 @@
           DEALLOCATE ( rrtm_sw_ssaaer )
           DEALLOCATE ( rrtm_sw_asmaer ) 
-          DEALLOCATE ( rrtm_sw_ecaer )    
+          DEALLOCATE ( rrtm_sw_ecaer )   
+ 
           DEALLOCATE ( rrtm_swdflx  )
+          DEALLOCATE ( rrtm_swdflxc )
           DEALLOCATE ( rrtm_swuflx  )
+          DEALLOCATE ( rrtm_swuflxc )
           DEALLOCATE ( rrtm_swhr  )
-          DEALLOCATE ( rrtm_swuflxc )
-          DEALLOCATE ( rrtm_swdflxc )
           DEALLOCATE ( rrtm_swhrc )
+
        ENDIF
 
@@ -1063 +1164 @@
 !--    Read pressure from file
        nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
-       nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/),    &
+       nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/),      &
                                count = (/nz_snd/) )
        CALL handle_netcdf_error( 'netcdf', 552 )
@@ -1075 +1176 @@
 !--    Read temperature from file
        nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
-       nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/),      &
+       nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/),        &
                                count = (/nz_snd/) )
        CALL handle_netcdf_error( 'netcdf', 553 )
@@ -1088 +1189 @@
 !--    in Pa, hyp_snd in hPa).
        DO  k = 1, nz_snd
-          IF ( hyp_snd_tmp(k) .LT. (hyp(nzt+1) - 1000.0_wp) * 0.01_wp )  THEN
+          IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp )  THEN
              nz_snd_start = k
              EXIT
@@ -1094 +1195 @@
        END DO
 
-       IF ( nz_snd_start .LE. nz_snd )  THEN
+       IF ( nz_snd_start <= nz_snd )  THEN
           nz_snd_end = nz_snd - 1
        END IF
@@ -1128 +1229 @@
        ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2)   )
 
-       t_surface = pt_surface * (surface_pressure / 1000.0_wp )**0.286_wp
+       t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
        DO k = nzb+1, nzt+1
           rrtm_play(0,k) = hyp(k) * 0.01_wp
@@ -1160 +1261 @@
           rrtm_h2ovmr(0,k) = q_snd(k)
        ENDDO
-       rrtm_tlay(0,nzt_rad+1)   = 2.0_wp * rrtm_tlay(0,nzt_rad)                &
-                                  - rrtm_tlay(0,nzt_rad-1)
+       rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                 &
+                                - rrtm_tlay(0,nzt_rad-1)
        DO k = nzt+9, nzt_rad+1
           rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k)                &
@@ -1237 +1338 @@
        rrtm_lwhrc   = 0.0_wp
 
+       ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
+       ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
+
+       rrtm_lwuflxc_dt = 0.0_wp
+       rrtm_lwuflxc_dt = 0.0_wp
 
     END SUBROUTINE read_sounding_data
@@ -1253 +1359 @@
        IMPLICIT NONE
 
-       INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases
-
-       CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              &
+       INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
+
+       CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              & !< trace gas names
            trace_names = (/'O3   ', 'CO2  ', 'CH4  ', 'N2O  ', 'O2   ',        &
                            'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
 
-       INTEGER(iwp) :: id, k, m, n, nabs, np, id_abs, id_dim, id_var
+       INTEGER(iwp) :: id,     & !< NetCDF id
+                       k,      & !< loop index
+                       m,      & !< loop index
+                       n,      & !< loop index
+                       nabs,   & !< number of absorbers
+                       np,     & !< number of pressure levels
+                       id_abs, & !< NetCDF id of the respective absorber
+                       id_dim, & !< NetCDF id of asborber's dimension
+                       id_var    !< NetCDf id ot the absorber
 
        REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
@@ -1380 +1494 @@
 !--          When the pressure level is higher than the trace gas pressure
 !--          level, assume that 
-             IF ( rrtm_plev_tmp(k-1) .GT. p_mls(1) )  THEN             
+             IF ( rrtm_plev_tmp(k-1) > p_mls(1) )  THEN             
                 
                 trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1)     &
@@ -1399 +1513 @@
                           MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
 
-                IF ( p_mls_l .GT. p_mls_u )  THEN
+                IF ( p_mls_l > p_mls_u )  THEN
 
 !
@@ -1412 +1526 @@
                                          +  ( p_wgt_u * trace_mls(m,n)         &
                                             + p_wgt_l * trace_mls(m,n-1) )     &
-                                         * (p_mls_l + p_mls_u) / g
+                                         * (p_mls_l - p_mls_u) / g
                 ENDIF
              ENDDO
 
-             IF ( rrtm_plev_tmp(k) .LT. p_mls(np) )  THEN 
+             IF ( rrtm_plev_tmp(k) < p_mls(np) )  THEN 
                 trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np)    &
                                       * ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
@@ -1503 +1617 @@
     SUBROUTINE radiation_tendency_ij ( i, j, tend )
 
-       USE arrays_3d,                                                          &
-           ONLY:  dzw
-
        USE cloud_parameters,                                                   &
-           ONLY:  pt_d_t, cp
-
-       USE control_parameters,                                                 &
-           ONLY:  rho_surface
+           ONLY:  pt_d_t
 
        IMPLICIT NONE
 
-       INTEGER(iwp) :: i, j, k
-
-       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend
+       INTEGER(iwp) :: i, j, k !< loop indices
+
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
 
 #if defined ( __rrtmg )
-
-       REAL(wp) :: rad_sw_net_l, rad_sw_net_u, rad_lw_net_l, rad_lw_net_u
-
-       rad_sw_net_l = 0.0_wp
-       rad_sw_net_u = 0.0_wp
-       rad_lw_net_l = 0.0_wp
-       rad_lw_net_u = 0.0_wp
-
-!
-!--    Calculate radiative flux divergence
+!
+!--    Calculate tendency based on heating rate
        DO k = nzb+1, nzt+1
-
-          rad_sw_net_l = rad_sw_in(k-1,j,i) - rad_sw_out(k-1,j,i)
-          rad_sw_net_u = rad_sw_in(k,j,i)   - rad_sw_out(k,j,i)
-          rad_lw_net_l = rad_lw_in(k-1,j,i) - rad_lw_out(k-1,j,i)
-          rad_lw_net_u = rad_lw_in(k,j,i)   - rad_lw_out(k,j,i)
-
-          tend(k,j,i) = tend(k,j,i) + ( rad_sw_net_u - rad_sw_net_l            &
-                                      + rad_lw_net_u - rad_lw_net_l ) /        &
-                                     ( dzw(k) * cp * rho_surface ) * pt_d_t(k) 
+          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
 
@@ -1554 +1647 @@
     SUBROUTINE radiation_tendency ( tend )
 
-       USE arrays_3d,                                                          &
-           ONLY:  dzw
-
        USE cloud_parameters,                                                   &
-           ONLY:  pt_d_t, cp
+           ONLY:  pt_d_t
 
        USE indices,                                                            &
            ONLY:  nxl, nxr, nyn, nys
 
-       USE control_parameters,                                                 &
-           ONLY:  rho_surface
-
        IMPLICIT NONE
 
-       INTEGER(iwp) :: i, j, k
-
-       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend
+       INTEGER(iwp) :: i, j, k !< loop indices
+
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
 
 #if defined ( __rrtmg )
-
-       REAL(wp) :: rad_sw_net_l, rad_sw_net_u, rad_lw_net_l, rad_lw_net_u
-
-       rad_sw_net_l = 0.0_wp
-       rad_sw_net_u = 0.0_wp
-       rad_lw_net_l = 0.0_wp
-       rad_lw_net_u = 0.0_wp
-
+!
+!--    Calculate tendency based on heating rate
        DO  i = nxl, nxr
           DO  j = nys, nyn
-
-!
-!--          Calculate radiative flux divergence
              DO k = nzb+1, nzt+1
-
-                rad_sw_net_l = rad_sw_in(k-1,j,i) - rad_sw_out(k-1,j,i)
-                rad_sw_net_u = rad_sw_in(k  ,j,i) - rad_sw_out(k  ,j,i)
-                rad_lw_net_l = rad_lw_in(k-1,j,i) - rad_lw_out(k-1,j,i)
-                rad_lw_net_u = rad_lw_in(k  ,j,i) - rad_lw_out(k  ,j,i)
-
-                tend(k,j,i) = tend(k,j,i) + ( rad_sw_net_u - rad_sw_net_l      &
-                                      + rad_lw_net_u - rad_lw_net_l ) /        &
-                                      ( dzw(k) * cp * rho_surface ) * pt_d_t(k) 
+                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

palm/trunk/SOURCE/read_3d_binary.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added output of radiative heating rates and Obukhov length. Removed output of
+! rif.
 ! 
 ! Former revisions:
@@ -83 +84 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, kh, km, p, pt, q, ql, qc, nr, nrs, nrsws, nrswst, qr, qrs,   &
-               qrsws, qrswst, qs, qsws, qswst, sa, saswsb, saswst, rif,        &
+        ONLY:  e, kh, km, ol, p, pt, q, ql, qc, nr, nrs, nrsws, nrswst, qr,    &
+               qrs, qrsws, qrswst, qs, qsws, qswst, sa, saswsb, saswst,        &
                rif_wall, shf, ts, tswst, u, u_m_l, u_m_n, u_m_r, u_m_s, us,    &
                usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, vsws, vswst,   &
@@ -119 +120 @@
     USE radiation_model_mod,                                                   &
         ONLY: rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, rad_lw_out,        &
-              rad_lw_out_av, rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av
+              rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av,          &
+              rad_lw_out_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
 
     USE random_function_mod,                                                   &
@@ -595 +599 @@
                    nrswst(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                CASE ( 'ol' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   ol(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'p' )
@@ -867 +875 @@
                              tmp_3d(:,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_out_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
@@ -898 +938 @@
                    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_out_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 ( 'random_iv' )  ! still unresolved issue
@@ -914 +986 @@
                    rho_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
                                 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
-                CASE ( 'rif' )
-                   IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   rif(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
-                                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'rif_wall' )

palm/trunk/SOURCE/read_var_list.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-! 
+! Added output of most_method, constant_flux_layer, zeta_min, zeta_max. Removed 
+! output of prandtl_layer and rif_min, rif_max.
 ! 
 ! Former revisions:
@@ -179 +180 @@
         ONLY:  statistic_regions, hom, hom_sum, pr_palm, u_max, u_max_ijk,     &
                v_max, v_max_ijk, w_max, w_max_ijk
+
 
     IMPLICIT NONE
@@ -356 +358 @@
           CASE ( 'conserve_volume_flow_mode' )
              READ ( 13 )  conserve_volume_flow_mode
+          CASE ( 'constant_flux_layer' )
+             READ ( 13 )  constant_flux_layer
           CASE ( 'coupling_start_time' )
              READ ( 13 )  coupling_start_time
@@ -460 +464 @@
           CASE ( 'momentum_advec' )
              READ ( 13 )  momentum_advec
+          CASE ( 'most_method' )
+             READ ( 13 )  most_method
           CASE ( 'nc_const' )
              READ ( 13 )  nc_const
@@ -494 +500 @@
           CASE ( 'phi' )
              READ ( 13 )  phi
-          CASE ( 'prandtl_layer' )
-             READ ( 13 )  prandtl_layer
           CASE ( 'prandtl_number' )
              READ ( 13 )  prandtl_number
@@ -552 +556 @@
           CASE ( 'residual_limit' )
              READ ( 13 )  residual_limit
-          CASE ( 'rif_max' )
-             READ ( 13 )  rif_max
-          CASE ( 'rif_min' )
-             READ ( 13 )  rif_min
           CASE ( 'roughness_length' )
              READ ( 13 )  roughness_length
@@ -714 +714 @@
           CASE ( 'w_max_ijk' )
              READ ( 13 )  w_max_ijk
+          CASE ( 'zeta_max' )
+             READ ( 13 )  zeta_max
+          CASE ( 'zeta_min' )
+             READ ( 13 )  zeta_min
           CASE ( 'z0h_factor' )
              READ ( 13 )  z0h_factor

palm/trunk/SOURCE/sum_up_3d_data.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added output of Obukhov length and radiative heating rates for RRTMG. 
+! Corrected output of LWC.
 ! 
 ! Former revisions:
@@ -84 +85 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  dzw, e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, qsws, rho, sa,    &
+        ONLY:  dzw, e, nr, ol, p, pt, q, qc, ql, ql_c, ql_v, qr, qsws, rho, sa,&
                shf, ts, u, us, v, vpt, w, z0, z0h
 
     USE averaging,                                                             &
-        ONLY:  e_av, lpt_av, lwp_av, nr_av, p_av, pc_av, pr_av, prr_av,        &
+        ONLY:  e_av, lpt_av, lwp_av, nr_av, ol_av, p_av, pc_av, pr_av, prr_av, &
                precipitation_rate_av, pt_av, q_av, qc_av, ql_av, ql_c_av,      &
                ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_av, s_av, sa_av,  &
@@ -97 +98 @@
 
     USE control_parameters,                                                    &
-        ONLY:  average_count_3d, cloud_physics, doav, doav_n
+        ONLY:  average_count_3d, cloud_physics, doav, doav_n, rho_surface
 
     USE cpulog,                                                                &
@@ -120 +121 @@
     USE radiation_model_mod,                                                   &
         ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
-               rad_sw_out_av, rad_lw_in, rad_lw_in_av, rad_lw_out,             &
-               rad_lw_out_av
+               rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,        &
+               rad_sw_hr_av, rad_lw_in, rad_lw_in_av, rad_lw_out,              &
+               rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr,        &
+               rad_lw_hr_av
+
 
     IMPLICIT NONE
@@ -215 +219 @@
                 nr_av = 0.0_wp
 
+             CASE ( 'ol*' )
+                IF ( .NOT. ALLOCATED( ol_av ) )  THEN
+                   ALLOCATE( ol_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                ol_av = 0.0_wp
+
              CASE ( 'p' )
                 IF ( .NOT. ALLOCATED( p_av ) )  THEN
@@ -347 +357 @@
                 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
@@ -359 +381 @@
                 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 ( 'rho' )
                 IF ( .NOT. ALLOCATED( rho_av ) )  THEN
@@ -530 +564 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i) * &
-                                                    dzw(1:nzt+1) )
+                   lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i)            &
+                                               * dzw(1:nzt+1) ) * rho_surface
                 ENDDO
              ENDDO
@@ -557 +591 @@
                       nr_av(k,j,i) = nr_av(k,j,i) + nr(k,j,i)
                    ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'ol*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   ol_av(j,i) = ol_av(j,i) + ol(j,i)
                 ENDDO
              ENDDO
@@ -777 +818 @@
              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' )
@@ -792 +850 @@
                    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

palm/trunk/SOURCE/time_integration.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! ------------------
-! 
+! Added option for spin-ups without land surface and radiation models. Moved calls
+! for radiation and lan surface schemes.
 ! 
 ! Former revisions:
@@ -170 +171 @@
                averaging_interval_sp, bc_lr_cyc, bc_ns_cyc, bc_pt_t_val,       &
                bc_q_t_val, call_psolver_at_all_substeps, cloud_droplets,       &
-               cloud_physics, constant_heatflux, create_disturbances, dopr_n,  &
-               constant_diffusion, coupling_mode, coupling_start_time,         &
-               current_timestep_number, disturbance_created,                   &
-               disturbance_energy_limit, dist_range, do_sum, dt_3d,            &
-               dt_averaging_input, dt_averaging_input_pr, dt_coupling,         &
-               dt_data_output_av, dt_disturb, dt_do2d_xy, dt_do2d_xz,          &
-               dt_do2d_yz, dt_do3d, dt_domask,dt_dopts, dt_dopr,               &
+               cloud_physics, constant_flux_layer, constant_heatflux,          &
+               create_disturbances, dopr_n, constant_diffusion, coupling_mode, &
+               coupling_start_time, current_timestep_number,                   &
+               disturbance_created, disturbance_energy_limit, dist_range,      &
+               do_sum, dt_3d, dt_averaging_input, dt_averaging_input_pr,       &
+               dt_coupling, dt_data_output_av, dt_disturb, dt_do2d_xy,         &
+               dt_do2d_xz, dt_do2d_yz, dt_do3d, dt_domask,dt_dopts, dt_dopr,   &
                dt_dopr_listing, dt_dosp, dt_dots, dt_dvrp, dt_run_control,     &
                end_time, first_call_lpm, galilei_transformation, humidity,     &
@@ -183 +184 @@
                loop_optimization, lsf_surf, lsf_vert, masks, mid,              &
                netcdf_data_format, neutral, nr_timesteps_this_run, nudging,    &
-               ocean, on_device, passive_scalar, prandtl_layer, precipitation, &
+               ocean, on_device, passive_scalar, precipitation,                &
                prho_reference, pt_reference, pt_slope_offset, random_heatflux, &
                run_coupled, simulated_time, simulated_time_chr,                &
@@ -213 +214 @@
 
     USE land_surface_model_mod,                                                &
-        ONLY:  land_surface, lsm_energy_balance, lsm_soil_model
+        ONLY:  land_surface, lsm_energy_balance, lsm_soil_model,               &
+               skip_time_do_lsm
 
     USE ls_forcing_mod,                                                        &
@@ -238 +240 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY: dt_radiation, radiation, radiation_clearsky,                     &
-              radiation_rrtmg, radiation_scheme, time_radiation
+        ONLY: dt_radiation, force_radiation_call, radiation,                   &
+              radiation_clearsky, radiation_rrtmg, radiation_scheme,           &
+              skip_time_do_radiation, time_radiation
 
     USE statistics,                                                            &
         ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l
+
+    USE surface_layer_fluxes_mod,                                              &
+        ONLY:  surface_layer_fluxes
 
     USE user_actions_mod,                                                      &
@@ -374 +380 @@
           ENDIF
 
-          IF ( radiation .AND. intermediate_timestep_count                     &
-               == intermediate_timestep_count_max  )  THEN
-
-               time_radiation = time_radiation + dt_3d
-
-             IF ( time_radiation >= dt_radiation )  THEN
-
-                CALL cpu_log( log_point(50), 'radiation', 'start' )
-
-                time_radiation = time_radiation - dt_radiation
-                IF ( radiation_scheme == 'clear-sky' )  THEN
-                   CALL radiation_clearsky
-                ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
-                   CALL radiation_rrtmg
-                ENDIF
-
-                CALL cpu_log( log_point(50), 'radiation', 'stop' )
-             ENDIF
-          ENDIF
 !
 !--       Solve the prognostic equations. A fast cache optimized version with
@@ -631 +618 @@
 !--       velocities at the outflow in case of a non-cyclic lateral wall)
           CALL boundary_conds
-
-!
-!--       When using the land surface model: 
-!--       1) solve energy balance equation to calculate new skin temperature
-!--       2) run soil model 
-          IF ( land_surface )  THEN
-
-             CALL cpu_log( log_point(54), 'land_surface', 'start' )
-
-             CALL lsm_energy_balance
-             CALL lsm_soil_model
-
-             CALL cpu_log( log_point(54), 'land_surface', 'stop' )
-
-          ENDIF
 
 !
@@ -731 +703 @@
 
 !
-!--          First the vertical fluxes in the Prandtl layer are being computed
-             IF ( prandtl_layer )  THEN
-                CALL cpu_log( log_point(19), 'prandtl_fluxes', 'start' )
-                CALL prandtl_fluxes
-                CALL cpu_log( log_point(19), 'prandtl_fluxes', 'stop' )
-             ENDIF
-
+!--          First the vertical fluxes in the surface (constant flux) layer are computed
+             IF ( constant_flux_layer )  THEN
+                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
+                CALL surface_layer_fluxes
+                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
+             ENDIF
+
+!
+!--          If required, solve the energy balance for the surface and run soil 
+!--          model
+             IF ( land_surface .AND. simulated_time > skip_time_do_lsm)  THEN
+
+                CALL cpu_log( log_point(54), 'land_surface', 'start' )
+                CALL lsm_energy_balance
+                CALL lsm_soil_model
+                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
+             ENDIF
 !
 !--          Compute the diffusion coefficients
@@ -752 +734 @@
              CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
 
+          ENDIF
+
+!
+!--       If required, calculate radiative fluxes and heating rates
+          IF ( radiation .AND. intermediate_timestep_count                     &
+               == intermediate_timestep_count_max .AND. simulated_time >    &
+               skip_time_do_radiation )  THEN
+
+               time_radiation = time_radiation + dt_3d
+
+             IF ( time_radiation >= dt_radiation .OR. force_radiation_call )   &
+             THEN
+
+                CALL cpu_log( log_point(50), 'radiation', 'start' )
+
+                IF ( .NOT. force_radiation_call )  THEN
+                   time_radiation = time_radiation - dt_radiation
+                ELSE
+                   WRITE(9,*) "Unscheduled radiation call at ", simulated_time
+                   CALL LOCAL_FLUSH ( 9 )
+                ENDIF
+
+                IF ( radiation_scheme == 'clear-sky' )  THEN
+                   CALL radiation_clearsky
+                ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
+                   CALL radiation_rrtmg
+                ENDIF
+
+                CALL cpu_log( log_point(50), 'radiation', 'stop' )
+             ENDIF
           ENDIF
 

palm/trunk/SOURCE/wall_fluxes.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Renamed rif_min and rif_max with zeta_min and zeta_max, respectively.
 ! 
 ! Former revisions:
@@ -69 +69 @@
 !> it gives slightly different results from the ij-version for some unknown
 !> reason.
+!>
+!> @todo Rename rif to zeta throughout the routine
 !------------------------------------------------------------------------------!
  MODULE wall_fluxes_mod
@@ -107 +109 @@
        
        USE control_parameters,                                                 &
-           ONLY:  g, kappa, rif_max, rif_min
+           ONLY:  g, kappa, zeta_max, zeta_min
        
        USE grid_variables,                                                     &
@@ -234 +236 @@
 !--                shear stresses and very small momentum fluxes (both are
 !--                generally unrealistic).
-                   IF ( rifs < rif_min )  rifs = rif_min
-                   IF ( rifs > rif_max )  rifs = rif_max
+                   IF ( rifs < zeta_min )  rifs = zeta_min
+                   IF ( rifs > zeta_max )  rifs = zeta_max
 
 !
@@ -291 +293 @@
        
        USE control_parameters,                                                 &
-           ONLY:  g, kappa, rif_max, rif_min
+           ONLY:  g, kappa, zeta_max, zeta_min
        
        USE grid_variables,                                                     &
@@ -428 +430 @@
 !--                shear stresses and very small momentum fluxes (both are
 !--                generally unrealistic).
-                   IF ( rifs < rif_min )  rifs = rif_min
-                   IF ( rifs > rif_max )  rifs = rif_max
+                   IF ( rifs < zeta_min )  rifs = zeta_min
+                   IF ( rifs > zeta_max )  rifs = zeta_max
 
 !
@@ -485 +487 @@
        
        USE control_parameters,                                                 &
-           ONLY:  g, kappa, rif_max, rif_min
+           ONLY:  g, kappa, zeta_max, zeta_min
        
        USE grid_variables,                                                     &
@@ -596 +598 @@
 !--       consequence would result in very large shear stresses and very
 !--       small momentum fluxes (both are generally unrealistic).
-          IF ( rifs < rif_min )  rifs = rif_min
-          IF ( rifs > rif_max )  rifs = rif_max
+          IF ( rifs < zeta_min )  rifs = zeta_min
+          IF ( rifs > zeta_max )  rifs = zeta_max
 
 !

palm/trunk/SOURCE/write_3d_binary.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added output of radiative heating rates for RRTMG. Added output of ol. Removed
+! output of rif.
 ! 
 ! Former revisions:
@@ -79 +80 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, kh, km, p, pt, q, ql, qc, nr, nrs, nrsws, nrswst, qr, qrs,   &
-               qrsws, qrswst, qs, qsws, qswst, sa, saswsb, saswst, rif,        &
+        ONLY:  e, kh, km, ol, p, pt, q, ql, qc, nr, nrs, nrsws, nrswst, qr,    &
+               qrs, qrsws, qrswst, qs, qsws, qswst, sa, saswsb, saswst,        &
                rif_wall, shf, ts, tswst, u, u_m_l, u_m_n, u_m_r, u_m_s, us,    &
                usws, uswst, v, v_m_l, v_m_n, v_m_r, v_m_s, vpt, vsws, vswst,   &
@@ -110 +111 @@
     USE radiation_model_mod,                                                   &
         ONLY: radiation, rad_net, rad_net_av, rad_lw_in, rad_lw_in_av,         &
-              rad_lw_out, rad_lw_out_av, rad_sw_in, rad_sw_in_av, rad_sw_out,  &
-              rad_sw_out_av
+              rad_lw_out, rad_lw_out_av, 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
 
     USE random_function_mod,                                                   &
@@ -190 +193 @@
        ENDIF
     ENDIF
+    WRITE ( 14 )  'ol                  ';  WRITE ( 14 )  ol
     WRITE ( 14 )  'p                   ';  WRITE ( 14 )  p
     IF ( ALLOCATED( p_av ) )  THEN
@@ -290 +294 @@
           WRITE ( 14 )  'rad_lw_out_av       ';  WRITE ( 14 )  rad_lw_out_av  
        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  
@@ -302 +318 @@
           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
     ENDIF
     IF ( ocean )  THEN
@@ -338 +366 @@
                                            WRITE ( 14 )  seq_random_array
     ENDIF
-    WRITE ( 14 )  'rif                 ';  WRITE ( 14 )  rif
     IF ( topography /= 'flat' )  THEN
        WRITE ( 14 )  'rif_wall            ';  WRITE ( 14 )  rif_wall

palm/trunk/SOURCE/write_var_list.f90

@@ -14 +14 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2014 Leibniz Universitaet Hannover
+! Copyright 1997-2015 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !
 ! Current revisions:
 ! -----------------
-! 
+! Added output of most_method, constant_flux_layer, zeta_min, zeta_max. Removed 
+! output of prandtl_layer and rif_min, rif_max.
 ! 
 ! Former revisions:
@@ -161 +162 @@
         ONLY:  statistic_regions, hom, hom_sum, u_max, u_max_ijk, v_max,       &
                v_max_ijk, w_max, w_max_ijk
+
     
     IMPLICIT NONE
@@ -167 +169 @@
 
 
-    binary_version = '3.9b'
+    binary_version = '4.0'
 
     WRITE ( 14 )  binary_version
@@ -273 +275 @@
     WRITE ( 14 )  conserve_volume_flow_mode
     WRITE ( 14 )  'coupling_start_time           '
+    WRITE ( 14 )  'constant_flux_layer           '
+    WRITE ( 14 )  constant_flux_layer
     WRITE ( 14 )  coupling_start_time
     WRITE ( 14 )  'current_timestep_number       '
@@ -374 +378 @@
     WRITE ( 14 )  'momentum_advec                '
     WRITE ( 14 )  momentum_advec
+    WRITE ( 14 )  'most_method                   '
+    WRITE ( 14 )  most_method
     WRITE ( 14 )  'nc_const                      '
     WRITE ( 14 )  nc_const
@@ -406 +412 @@
     WRITE ( 14 )  'phi                           '
     WRITE ( 14 )  phi
-    WRITE ( 14 )  'prandtl_layer                 '
-    WRITE ( 14 )  prandtl_layer
     WRITE ( 14 )  'prandtl_number                '
     WRITE ( 14 )  prandtl_number
@@ -462 +466 @@
     WRITE ( 14 )  'residual_limit                '
     WRITE ( 14 )  residual_limit
-    WRITE ( 14 )  'rif_max                       '
-    WRITE ( 14 )  rif_max
-    WRITE ( 14 )  'rif_min                       '
-    WRITE ( 14 )  rif_min
     WRITE ( 14 )  'roughness_length              '
     WRITE ( 14 )  roughness_length
@@ -624 +624 @@
     WRITE ( 14 )  'w_max_ijk                     '
     WRITE ( 14 )  w_max_ijk
+    WRITE ( 14 )  'zeta_max                      '
+    WRITE ( 14 )  zeta_max
+    WRITE ( 14 )  'zeta_min                      '
+    WRITE ( 14 )  zeta_min
     WRITE ( 14 )  'z0h_factor                    '
     WRITE ( 14 )  z0h_factor