Changeset 1551 for palm

Timestamp:

Mar 3, 2015 2:18:16 PM (9 years ago)

Author:

maronga

Message:

land surface model released

Location:

palm/trunk

Files:

• SOURCE/Makefile (modified) (10 diffs)
• SOURCE/average_3d_data.f90 (modified) (9 diffs)
• SOURCE/check_open.f90 (modified) (3 diffs)
• SOURCE/check_parameters.f90 (modified) (14 diffs)
• SOURCE/data_output_2d.f90 (modified) (46 diffs)
• SOURCE/data_output_3d.f90 (modified) (20 diffs)
• SOURCE/flow_statistics.f90 (modified) (12 diffs)
• SOURCE/header.f90 (modified) (20 diffs)
• SOURCE/init_3d_model.f90 (modified) (3 diffs)
• SOURCE/land_surface_model.f90 (modified) (43 diffs)
• SOURCE/modules.f90 (modified) (6 diffs)
• SOURCE/netcdf.f90 (modified) (29 diffs)
• SOURCE/package_parin.f90 (modified) (5 diffs)
• SOURCE/prandtl_fluxes.f90 (modified) (4 diffs)
• SOURCE/radiation_model.f90 (modified) (12 diffs)
• SOURCE/read_3d_binary.f90 (modified) (13 diffs)
• SOURCE/read_var_list.f90 (modified) (1 diff)
• SOURCE/sum_up_3d_data.f90 (modified) (14 diffs)
• SOURCE/swap_timelevel.f90 (modified) (5 diffs)
• SOURCE/time_integration.f90 (modified) (3 diffs)
• SOURCE/user_data_output_2d.f90 (modified) (5 diffs)
• SOURCE/user_data_output_3d.f90 (modified) (5 diffs)
• SOURCE/user_header.f90 (modified) (4 diffs)
• SOURCE/write_3d_binary.f90 (modified) (10 diffs)
• SOURCE/write_var_list.f90 (modified) (2 diffs)
• UTIL/combine_plot_fields.f90 (modified) (12 diffs)

palm/trunk/SOURCE/Makefile

@@ -20 +20 @@
 # Current revisions:
 # ------------------
-#
+# Bugfix: further adjustments for the land surface model and radiation model
 # 
 # Former revisions:
@@ -293 +293 @@
 advec_w_pw.o: modules.o mod_kinds.o
 advec_w_up.o: modules.o mod_kinds.o
-average_3d_data.o: modules.o cpulog.o mod_kinds.o
+average_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
+                   radiation_model.o
 boundary_conds.o: modules.o mod_kinds.o
 buoyancy.o: modules.o mod_kinds.o
@@ -318 +319 @@
 data_output_spectra.o: modules.o cpulog.o mod_kinds.o
 data_output_tseries.o: modules.o cpulog.o mod_kinds.o
-data_output_2d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
-data_output_3d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o
+data_output_2d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o\
+                  land_surface_model.o radiation_model.o
+data_output_3d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o\
+                  land_surface_model.o
 diffusion_e.o: modules.o mod_kinds.o
 diffusion_s.o: modules.o mod_kinds.o
@@ -326 +329 @@
 diffusion_w.o: modules.o mod_kinds.o wall_fluxes.o
 diffusivities.o: modules.o mod_kinds.o
-disturb_field.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o
+disturb_field.o: modules.o cpulog.o mod_kinds.o random_function.o\
+                 random_generator_parallel.o
 disturb_heatflux.o: modules.o cpulog.o mod_kinds.o
 eqn_state_seawater.o: modules.o mod_kinds.o
@@ -332 +336 @@
 exchange_horiz_2d.o: modules.o cpulog.o mod_kinds.o
 fft_xy.o: cuda_fft_interfaces.o modules.o mod_kinds.o singleton.o temperton_fft.o
-flow_statistics.o: modules.o cpulog.o mod_kinds.o
+flow_statistics.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
+                   radiation_model.o
 global_min_max.o: modules.o mod_kinds.o
-header.o: modules.o cpulog.o mod_kinds.o plant_canopy_model.o subsidence.o
+header.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
+          plant_canopy_model.o radiation_model.o subsidence.o 
 impact_of_latent_heat.o: modules.o mod_kinds.o
 inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
 init_1d_model.o: modules.o mod_kinds.o
-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 \
+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
 init_advec.o: modules.o mod_kinds.o
@@ -395 +401 @@
 mod_kinds.o: mod_kinds.f90
 mod_particle_attributes.o: mod_particle_attributes.f90 mod_kinds.o
-netcdf.o: modules.o mod_kinds.o
+netcdf.o: modules.o mod_kinds.o land_surface_model.o
 nudging.o: modules.o cpulog.o mod_kinds.o
-package_parin.o: modules.o mod_kinds.o land_surface_model.o plant_canopy_model.o \
-                 radiation_model.o
+package_parin.o: modules.o mod_kinds.o land_surface_model.o\
+                 plant_canopy_model.o radiation_model.o
 palm.o: modules.o cpulog.o ls_forcing.o mod_kinds.o nudging.o
 parin.o: modules.o cpulog.o mod_kinds.o progress_bar.o
@@ -404 +410 @@
 poisfft.o: modules.o cpulog.o fft_xy.o mod_kinds.o tridia_solver.o
 poismg.o: modules.o cpulog.o mod_kinds.o
-prandtl_fluxes.o: modules.o mod_kinds.o land_surface_model.o
+prandtl_fluxes.o: modules.o mod_kinds.o
 pres.o: modules.o cpulog.o mod_kinds.o poisfft.o
 print_1d.o: modules.o cpulog.o mod_kinds.o
@@ -419 +425 @@
 random_gauss.o: mod_kinds.o random_function.o random_generator_parallel.o
 random_generator_parallel.o: mod_kinds.o
-read_3d_binary.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o
+read_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
+                  radiation_model.o random_function.o\
+                  random_generator_parallel.o
 read_var_list.o: modules.o mod_kinds.o plant_canopy_model.o
 run_control.o: modules.o cpulog.o mod_kinds.o
@@ -426 +434 @@
 sor.o: modules.o mod_kinds.o
 subsidence.o: modules.o mod_kinds.o
-sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o
+sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
+                  radiation_model.o
 surface_coupler.o: modules.o cpulog.o mod_kinds.o
 swap_timelevel.o: modules.o cpulog.o mod_kinds.o land_surface_model.o
@@ -468 +477 @@
 user_statistics.o: modules.o mod_kinds.o user_module.o
 wall_fluxes.o: modules.o mod_kinds.o
-write_3d_binary.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o
+write_3d_binary.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
+                   radiation_model.o random_function.o\
+                   random_generator_parallel.o
 write_var_list.o: modules.o mod_kinds.o plant_canopy_model.o

palm/trunk/SOURCE/average_3d_data.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! Added support for land surface and radiation model parameters.
 !
 ! Former revisions:
@@ -75 +75 @@
     USE kinds
 
+    USE land_surface_model_mod,                                                &
+        ONLY:  c_liq_av, c_soil_av, c_veg_av, ghf_eb_av, lai_av, m_liq_eb_av,  &
+               m_soil_av, nzb_soil, nzt_soil, qsws_eb_av, qsws_liq_eb_av,      &
+               qsws_soil_eb_av, qsws_veg_eb_av, shf_eb_av, t_soil_av
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
 
     IMPLICIT NONE
@@ -98 +105 @@
        SELECT CASE ( TRIM( doav(ii) ) )
 
+         CASE ( 'c_liq*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_liq_av(j,i) = c_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+         CASE ( 'c_soil*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_soil_av(j,i) = c_soil_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+         CASE ( 'c_veg*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_veg_av(j,i) = c_veg_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
           CASE ( 'e' )
              DO  i = nxlg, nxrg
@@ -107 +135 @@
              ENDDO
 
+         CASE ( 'ghf_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   ghf_eb_av(j,i) = ghf_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
           CASE ( 'qsws*' )
              DO  i = nxlg, nxrg
@@ -114 +149 @@
              ENDDO
 
+         CASE ( 'lai*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   lai_av(j,i) = lai_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
           CASE ( 'lpt' )
              DO  i = nxlg, nxrg
@@ -127 +169 @@
                 DO  j = nysg, nyng
                    lwp_av(j,i) = lwp_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+         CASE ( 'm_liq_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   m_liq_eb_av(j,i) = m_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 'm_soil' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb_soil, nzt_soil
+                      m_soil_av(k,j,i) = m_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
                 ENDDO
              ENDDO
@@ -247 +305 @@
              ENDDO
 
+         CASE ( 'qsws_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_eb_av(j,i) = qsws_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+         CASE ( 'qsws_liq_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+         CASE ( 'qsws_soil_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+         CASE ( 'qsws_veg_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
           CASE ( 'qv' )
              DO  i = nxlg, nxrg
@@ -256 +342 @@
              ENDDO
 
+         CASE ( 'rad_sw_in*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   rad_sw_in_av(j,i) = rad_sw_in_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+         CASE ( 'rad_net*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
           CASE ( 'rho' )
              DO  i = nxlg, nxrg
@@ -294 +394 @@
                 DO  j = nysg, nyng
                    ts_av(j,i) = ts_av(j,i) / REAL( average_count_3d, KIND=wp )
+                ENDDO
+             ENDDO
+
+          CASE ( 't_soil' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb_soil, nzt_soil
+                      t_soil_av(k,j,i) = t_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
                 ENDDO
              ENDDO

palm/trunk/SOURCE/check_open.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Removed redundant output for combine_plot_fields
 ! 
 ! Former revisions:
@@ -99 +99 @@
 
     USE indices,                                                               &
-        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz,   &
-               nzb, nzt 
+        ONLY:  nbgp, nx, nxl, nxr, ny, nyn, nyng, nys, nysg, nz, nzb, nzt 
 
     USE kinds
@@ -462 +461 @@
 !
 !--          Specifications for combine_plot_fields
-             WRITE ( 30 )  -nbgp,nx+nbgp,-nbgp,ny+nbgp, 0 ,nz_do3d
+             WRITE ( 30 )  -nbgp,nx+nbgp,-nbgp,ny+nbgp
              WRITE ( 30 )  0,nx+1,0,ny+1,0,nz_do3d
 #endif

palm/trunk/SOURCE/check_parameters.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added various checks for land surface and radiation model. In the course of this
+! action, the length of the variable var has to be increased
 ! 
 ! Former revisions:
@@ -278 +279 @@
 
     CHARACTER (LEN=1)   ::  sq                       !: 
-    CHARACTER (LEN=6)   ::  var                      !: 
+    CHARACTER (LEN=15)  ::  var                      !: 
     CHARACTER (LEN=7)   ::  unit                     !: 
     CHARACTER (LEN=8)   ::  date                     !: 
@@ -970 +971 @@
        IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' )  THEN
           message_string = 'lsm requires setting of'//                         &
-                           'bc_pt_b = "dirichlet" and '//                        &
+                           'bc_pt_b = "dirichlet" and '//                      &
                            'bc_q_b  = "dirichlet"'
           CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
@@ -988 +989 @@
              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
           ENDIF
+ 
+          IF ( min_canopy_resistance == 9999999.9_wp)  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of min_canopy_resistance'//    &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( leaf_area_index == 9999999.9_wp)  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of leaf_area_index'//          &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( vegetation_coverage == 9999999.9_wp)  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of vegetation_coverage'//      &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of'//                          &
+                              'canopy_resistance_coefficient /= 9999999.9'
+             CALL message( 'check_parameters', 'PA0418', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( lambda_surface_stable == 9999999.9_wp)  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of lambda_surface_stable'//    &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0419', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( lambda_surface_unstable == 9999999.9_wp)  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of lambda_surface_unstable'//  &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0420', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( f_shortwave_incoming == 9999999.9_wp)  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of f_shortwave_incoming'//     &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0421', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( z0_eb == 9999999.9_wp)  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of z0_eb'//                   &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( z0h_eb == 9999999.9_wp)  THEN
+             message_string = 'veg_type = 0 (user_defined)'//                  &
+                              'requires setting of z0h_eb'//                  &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
+          ENDIF
+
+
+       ENDIF
+
+       IF ( soil_type == 0 )  THEN
+
+          IF ( alpha_vangenuchten == 9999999.9_wp)  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of alpha_vangenuchten'//       &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( l_vangenuchten == 9999999.9_wp)  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of l_vangenuchten'//           &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( n_vangenuchten == 9999999.9_wp)  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of n_vangenuchten'//           &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0424', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( hydraulic_conductivity == 9999999.9_wp)  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of hydraulic_conductivity'//   &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0425', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( saturation_moisture == 9999999.9_wp)  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of saturation_moisture'//      &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0426', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( field_capacity == 9999999.9_wp)  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of field_capacity'//           &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0427', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( wilting_point == 9999999.9_wp)  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of wilting_point'//            &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0428', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( residual_moisture == 9999999.9_wp)  THEN
+             message_string = 'soil_type = 0 (user_defined)'//                 &
+                              'requires setting of residual_moisture'//        &
+                              '/= 9999999.9'
+             CALL message( 'check_parameters', 'PA0429', 1, 2, 0, 6, 0 )
+          ENDIF
+
        ENDIF
 
@@ -996 +1122 @@
        ENDIF
 
-
     END IF 
+
+    IF ( radiation )  THEN
+       IF ( radiation_scheme == 'constant' )  THEN
+          irad_scheme = 0
+       ELSEIF ( radiation_scheme == 'clear-sky' )  THEN
+          irad_scheme = 1
+       ELSEIF ( radiation_scheme == 'rrtm' )  THEN
+          irad_scheme = 2
+       ELSE
+          message_string = 'unknown radiation_scheme = '//                     &
+                           TRIM( radiation_scheme )
+          CALL message( 'check_parameters', 'PA0430', 1, 2, 0, 6, 0 )
+       ENDIF
+    ENDIF
 
 
@@ -2862 +3001 @@
              ENDIF
 
+          CASE ( 't_soil', '#t_soil' )
+             IF ( .NOT. land_surface )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for land_surface = .FALSE.'
+                CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 89
+                dopr_unit(i)  = 'K'
+                hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
+                IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                   dopr_initial_index(i) = 90
+                   hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
+                   data_output_pr(i)     = data_output_pr(i)(2:)
+                ENDIF
+             ENDIF
+
+          CASE ( 'm_soil', '#m_soil' )
+             IF ( .NOT. land_surface )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for land_surface = .FALSE.'
+                CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 91
+                dopr_unit(i)  = 'm3/m3'
+                hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
+                IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                   dopr_initial_index(i) = 92
+                   hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
+                   data_output_pr(i)     = data_output_pr(i)(2:)
+                ENDIF
+             ENDIF
+
 
           CASE DEFAULT
@@ -2932 +3105 @@
           ENDIF
        ENDIF
+
 !
 !--    Check for allowed value and set units
@@ -2951 +3125 @@
              ENDIF
              unit = 'K'
+
+          CASE ( 'm_soil' )
+             IF ( .NOT. land_surface )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                         'land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'm3/m3'
 
           CASE ( 'nr' )
@@ -3076 +3258 @@
              unit = 'psu'
 
-          CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'qsws*', 'shf*', 'z0*', 'z0h*' )
+          CASE ( 't_soil' )
+             IF ( .NOT. land_surface )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                         'land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'K'
+
+
+          CASE ( 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'lai*', 'lwp*',     &
+                 'm_liq_eb*', 'pra*', 'prr*', 'qsws*', 'qsws_eb*',             &
+                 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*',              &
+                 'rad_net*', 'rad_sw_in*', 'shf*', 'shf_eb*', 't*', 'u*',      &
+                 'z0*', 'z0h*' )
              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
                 message_string = 'illegal value for data_output: "' //         &
@@ -3083 +3278 @@
                 CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
              ENDIF
+             IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT. land_surface )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT. land_surface )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT. land_surface )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( TRIM( var ) == 'lai*'  .AND.  .NOT. land_surface )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
+             ENDIF
              IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -3088 +3308 @@
                 CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
              ENDIF
+             IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT. land_surface )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
+             ENDIF
              IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
@@ -3108 +3333 @@
                 CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
              ENDIF
-
+             IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT. land_surface )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )  &
+             THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT. land_surface ) &
+             THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )  &
+             THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                                 'res land_surface = .TRUE.'
+                CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
+             ENDIF
+
+             IF ( TRIM( var ) == 'c_liq*' )  unit = 'none'
+             IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
+             IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
+             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 ) == 'pra*'   )  unit = 'mm'
              IF ( TRIM( var ) == 'prr*'   )  unit = 'mm/s'
              IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
+             IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
+             IF ( TRIM( var ) == 'rad_net*')       unit = 'W/m2'     
+             IF ( TRIM( var ) == 'rad_sw_in*')     unit = 'W/m2'   
              IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
+             IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
              IF ( TRIM( var ) == 't*'     )  unit = 'K'
              IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
@@ -3129 +3389 @@
 
           CASE DEFAULT
+
              CALL user_check_data_output( var, unit )
 
@@ -3137 +3398 @@
                    CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
                 ELSE
-                   message_string = 'illegal value for data_output =' //       &
+                   message_string = 'illegal value for data_output = "' //     &
                                     TRIM( data_output(i) ) // '"'
                    CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )

palm/trunk/SOURCE/data_output_2d.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added suppport for land surface model and radiation model output. In the course
+! of this action, the limits for vertical loops have been changed (from nzb and 
+! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
+! Moreover, a new vertical grid zs was introduced.
 ! 
 ! Former revisions:
@@ -127 +130 @@
                
     USE kinds
-        
+    
+    USE land_surface_model_mod,                                                &
+        ONLY:  c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb,            &
+               ghf_eb_av, lai, lai_av, m_liq_eb, m_liq_eb_av, m_soil,          &
+               m_soil_av, nzb_soil, nzt_soil, qsws_eb, qsws_eb_av,             &
+               qsws_liq_eb, qsws_liq_eb_av, qsws_soil_eb, qsws_soil_eb_av,     &
+               qsws_veg_eb, qsws_veg_eb_av, shf_eb, shf_eb_av, t_soil,         &
+               t_soil_av, zs
+    
     USE netcdf_control
 
@@ -135 +146 @@
     
     USE pegrid
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
 
     IMPLICIT NONE
@@ -157 +171 @@
     INTEGER(iwp) ::  n         !: 
     INTEGER(iwp) ::  ns        !: 
+    INTEGER(iwp) ::  nzb_do    !: lower limit of the data field (usually nzb)
+    INTEGER(iwp) ::  nzt_do    !: upper limit of the data field (usually nzt+1)
     INTEGER(iwp) ::  psi       !: 
     INTEGER(iwp) ::  s         !: 
@@ -343 +359 @@
 
        IF ( do2d_mode == mode )  THEN
+
+          nzb_do = nzb
+          nzt_do = nzt+1
 !
 !--       Store the array chosen on the temporary array.
@@ -356 +375 @@
                 IF ( mode == 'xy' )  level_z = zu
 
+             CASE ( 'c_liq*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = c_liq(j,i) * c_veg(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = c_liq_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'c_soil*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = 1.0_wp - c_veg(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = c_soil_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'c_veg*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = c_veg(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = c_veg_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'ghf_eb*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = ghf_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = ghf_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'lai*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = lai(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = lai_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
              CASE ( 'lpt_xy', 'lpt_xz', 'lpt_yz' )
                 IF ( av == 0 )  THEN
@@ -382 +491 @@
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'm_liq_eb*_xy' )        ! 2d-array
+                IF ( av == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = m_liq_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) = m_liq_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'm_soil_xy', 'm_soil_xz', 'm_soil_yz' )
+                nzb_do = nzb_soil
+                nzt_do = nzt_soil
+                IF ( av == 0 )  THEN
+                   to_be_resorted => m_soil
+                ELSE
+                   to_be_resorted => m_soil_av
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zs
 
              CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
@@ -665 +802 @@
                 level_z(nzb+1) = zu(nzb+1)
 
+             CASE ( 'qsws_eb*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  qsws_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  qsws_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'qsws_liq_eb*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  qsws_liq_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  qsws_liq_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'qsws_soil_eb*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  qsws_soil_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  qsws_soil_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'qsws_veg_eb*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  qsws_veg_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  qsws_veg_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
              CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
                 IF ( av == 0 )  THEN
@@ -680 +889 @@
                 IF ( mode == 'xy' )  level_z = zu
 
+             CASE ( 'rad_net*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  rad_net(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  rad_net_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'rad_sw_in*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  rad_sw_in(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
+                         local_pf(i,j,nzb+1) =  rad_sw_in_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
              CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )
                 IF ( av == 0 )  THEN
@@ -719 +964 @@
                 level_z(nzb+1) = zu(nzb+1)
 
+             CASE ( 'shf_eb*_xy' )        ! 2d-array
+                IF ( av == 0 ) THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  shf_eb(j,i)
+                      ENDDO
+                   ENDDO
+                ELSE
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         local_pf(i,j,nzb+1) =  shf_eb_av(j,i)
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                two_d = .TRUE.
+                level_z(nzb+1) = zu(nzb+1)
+
              CASE ( 't*_xy' )        ! 2d-array
                 IF ( av == 0 )  THEN
@@ -736 +999 @@
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 't_soil_xy', 't_soil_xz', 't_soil_yz' )
+                nzb_do = nzb_soil
+                nzt_do = nzt_soil
+                IF ( av == 0 )  THEN
+                   to_be_resorted => t_soil
+                ELSE
+                   to_be_resorted => t_soil_av
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zs
 
              CASE ( 'u_xy', 'u_xz', 'u_yz' )
@@ -839 +1112 @@
 !--             User defined quantity
                 CALL user_data_output_2d( av, do2d(av,if), found, grid,        &
-                                          local_pf, two_d )
+                                          local_pf, two_d, nzb_do, nzt_do )
                 resorted = .TRUE.
 
@@ -848 +1121 @@
                 ELSEIF ( grid == 'zu1' ) THEN
                    IF ( mode == 'xy' )  level_z(nzb+1) = zu(nzb+1)
+                ELSEIF ( grid == 'zs' ) THEN
+                   IF ( mode == 'xy' )  level_z = zs
                 ENDIF
 
@@ -863 +1138 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = to_be_resorted(k,j,i)
                    ENDDO
@@ -874 +1149 @@
 !--       section mode chosen.
           is = 1
-   loop1: DO  WHILE ( section(is,s) /= -9999  .OR.  two_d )
+   loop1: DO WHILE ( section(is,s) /= -9999  .OR.  two_d )
 
              SELECT CASE ( mode )
@@ -885 +1160 @@
                    ELSE
                       layer_xy = section(is,s)
+                   ENDIF
+
+!
+!--                Exit the loop for layers beyond the data output domain
+!--                (used for soil model)
+                   IF ( layer_xy .GT. nzt_do )  THEN 
+                      EXIT loop1
                    ENDIF
 
@@ -916 +1198 @@
 !
 !--                   Carry out the averaging (all data are on the PE)
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb_do, nzt_do
                          DO  j = nysg, nyng
                             DO  i = nxlg, nxrg
@@ -924 +1206 @@
                       ENDDO
 
-                      local_2d = local_2d / ( nzt -nzb + 2.0_wp)
+                      local_2d = local_2d / ( nzt_do - nzb_do + 1.0_wp)
 
                    ELSE
@@ -967 +1249 @@
                          DO  i = 0, io_blocks-1
                             IF ( i == io_group )  THEN
-                               WRITE ( 21 )  nxlg, nxrg, nysg, nyng
+                               WRITE ( 21 )  nxlg, nxrg, nysg, nyng, nysg, nyng
                                WRITE ( 21 )  local_2d
                             ENDIF
@@ -1103 +1385 @@
                    IF ( section(is,s) == -1 )  THEN
 
-                      ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) )
+                      ALLOCATE( local_2d_l(nxlg:nxrg,nzb_do:nzt_do) )
                       local_2d_l = 0.0_wp
-                      ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
+                      ngp = ( nxrg-nxlg + 1 ) * ( nzt_do-nzb_do + 1 )
 !
 !--                   First local averaging on the PE
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb_do, nzt_do
                          DO  j = nys, nyn
                             DO  i = nxlg, nxrg
@@ -1120 +1402 @@
 !--                   Now do the averaging over all PEs along y
                       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-                      CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb),                &
-                                          local_2d(nxlg,nzb), ngp, MPI_REAL,   &
+                      CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb_do),                &
+                                          local_2d(nxlg,nzb_do), ngp, MPI_REAL,   &
                                           MPI_SUM, comm1dy, ierr )
 #else
@@ -1136 +1418 @@
                       IF ( section(is,s) >= nys  .AND.  section(is,s) <= nyn ) &
                       THEN
-                         local_2d = local_pf(:,section(is,s),nzb:nzt+1)
+                         local_2d = local_pf(:,section(is,s),nzb_do:nzt_do)
                       ENDIF
 
@@ -1157 +1439 @@
 !--                      output file afterwards to increase the performance.
                          DO  i = nxlg, nxrg
-                            DO  k = nzb, nzt+1
+                            DO  k = nzb_do, nzt_do
                                local_2d_sections_l(i,is,k) = local_2d(i,k)
                             ENDDO
@@ -1184 +1466 @@
                                       nys-1 == -1 ) )                          &
                                THEN
-                                  WRITE (22)  nxlg, nxrg, nzb, nzt+1
+                                  WRITE (22)  nxlg, nxrg, nzb_do, nzt_do, nzb, nzt+1
                                   WRITE (22)  local_2d
                                ELSE
-                                  WRITE (22)  -1, -1, -1, -1
+                                  WRITE (22)  -1, -1, -1, -1, -1, -1
                                ENDIF
                             ENDIF
@@ -1203 +1485 @@
                          CALL MPI_BARRIER( comm2d, ierr )
 
-                         ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
+                         ngp = ( nxrg-nxlg + 1 ) * ( nzt_do-nzb_do + 1 )
                          IF ( myid == 0 )  THEN
 !
@@ -1211 +1493 @@
                                  ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
                             THEN
-                               total_2d(nxlg:nxrg,nzb:nzt+1) = local_2d
+                               total_2d(nxlg:nxrg,nzb_do:nzt_do) = local_2d
                             ENDIF
 !
@@ -1240 +1522 @@
 !--                         Relocate the local array for the next loop increment
                             DEALLOCATE( local_2d )
-                            ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
+                            ALLOCATE( local_2d(nxlg:nxrg,nzb_do:nzt_do) )
 
 #if defined( __netcdf )
                             nc_stat = NF90_PUT_VAR( id_set_xz(av),          &
                                                  id_var_do2d(av,if),        &
-                                                 total_2d(0:nx+1,nzb:nzt+1),&
+                                                 total_2d(0:nx+1,nzb_do:nzt_do),&
                             start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
-                                             count = (/ nx+2, 1, nz+2, 1 /) )
+                                             count = (/ nx+2, 1, nzt_do-nzb_do+1, 1 /) )
                             CALL handle_netcdf_error( 'data_output_2d', 58 )
 #endif
@@ -1260 +1542 @@
                             THEN
                                ind(1) = nxlg; ind(2) = nxrg
-                               ind(3) = nzb;   ind(4) = nzt+1
+                               ind(3) = nzb_do;   ind(4) = nzt_do
                             ELSE
                                ind(1) = -9999; ind(2) = -9999
@@ -1270 +1552 @@
 !--                         If applicable, send data to PE0.
                             IF ( ind(1) /= -9999 )  THEN
-                               CALL MPI_SEND( local_2d(nxlg,nzb), ngp,         &
+                               CALL MPI_SEND( local_2d(nxlg,nzb_do), ngp,         &
                                               MPI_REAL, 0, 1, comm2d, ierr )
                             ENDIF
@@ -1286 +1568 @@
                    nc_stat = NF90_PUT_VAR( id_set_xz(av),                   &
                                            id_var_do2d(av,if),              &
-                                           local_2d(nxl:nxr+1,nzb:nzt+1),   &
+                                           local_2d(nxl:nxr+1,nzb_do:nzt_do),   &
                             start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
-                                           count = (/ nx+2, 1, nz+2, 1 /) )
+                                           count = (/ nx+2, 1, nzt_do-nzb_do+1, 1 /) )
                    CALL handle_netcdf_error( 'data_output_2d', 451 )
 #endif
@@ -1322 +1604 @@
                    IF ( section(is,s) == -1 )  THEN
 
-                      ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) )
+                      ALLOCATE( local_2d_l(nysg:nyng,nzb_do:nzt_do) )
                       local_2d_l = 0.0_wp
-                      ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
+                      ngp = ( nyng-nysg+1 ) * ( nzt_do-nzb_do+1 )
 !
 !--                   First local averaging on the PE
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb_do, nzt_do
                          DO  j = nysg, nyng
                             DO  i = nxl, nxr
@@ -1339 +1621 @@
 !--                   Now do the averaging over all PEs along x
                       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-                      CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb),                &
-                                          local_2d(nysg,nzb), ngp, MPI_REAL,   &
+                      CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb_do),                &
+                                          local_2d(nysg,nzb_do), ngp, MPI_REAL,   &
                                           MPI_SUM, comm1dx, ierr )
 #else
@@ -1355 +1637 @@
                       IF ( section(is,s) >= nxl  .AND.  section(is,s) <= nxr ) &
                       THEN
-                         local_2d = local_pf(section(is,s),:,nzb:nzt+1)
+                         local_2d = local_pf(section(is,s),:,nzb_do:nzt_do)
                       ENDIF
 
@@ -1376 +1658 @@
 !--                      output file afterwards to increase the performance.
                          DO  j = nysg, nyng
-                            DO  k = nzb, nzt+1
+                            DO  k = nzb_do, nzt_do
                                local_2d_sections_l(is,j,k) = local_2d(j,k)
                             ENDDO
@@ -1403 +1685 @@
                                       nxl-1 == -1 ) )                          &
                                THEN
-                                  WRITE (23)  nysg, nyng, nzb, nzt+1
+                                  WRITE (23)  nysg, nyng, nzb_do, nzt_do, nzb, nzt+1
                                   WRITE (23)  local_2d
                                ELSE
-                                  WRITE (23)  -1, -1, -1, -1
+                                  WRITE (23)  -1, -1, -1, -1, -1, -1
                                ENDIF
                             ENDIF
@@ -1422 +1704 @@
                          CALL MPI_BARRIER( comm2d, ierr )
 
-                         ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
+                         ngp = ( nyng-nysg+1 ) * ( nzt_do-nzb_do+1 )
                          IF ( myid == 0 )  THEN
 !
@@ -1430 +1712 @@
                                  ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
                             THEN
-                               total_2d(nysg:nyng,nzb:nzt+1) = local_2d
+                               total_2d(nysg:nyng,nzb_do:nzt_do) = local_2d
                             ENDIF
 !
@@ -1459 +1741 @@
 !--                         Relocate the local array for the next loop increment
                             DEALLOCATE( local_2d )
-                            ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
+                            ALLOCATE( local_2d(nysg:nyng,nzb_do:nzt_do) )
 
 #if defined( __netcdf )
                             nc_stat = NF90_PUT_VAR( id_set_yz(av),          &
                                                  id_var_do2d(av,if),        &
-                                                 total_2d(0:ny+1,nzb:nzt+1),&
+                                                 total_2d(0:ny+1,nzb_do:nzt_do),&
                             start = (/ is, 1, 1, do2d_yz_time_count(av) /), &
-                                             count = (/ 1, ny+2, nz+2, 1 /) )
+                                             count = (/ 1, ny+2, nzt_do-nzb_do+1, 1 /) )
                             CALL handle_netcdf_error( 'data_output_2d', 61 )
 #endif
@@ -1479 +1761 @@
                             THEN
                                ind(1) = nysg; ind(2) = nyng
-                               ind(3) = nzb;   ind(4) = nzt+1
+                               ind(3) = nzb_do;   ind(4) = nzt_do
                             ELSE
                                ind(1) = -9999; ind(2) = -9999
@@ -1489 +1771 @@
 !--                         If applicable, send data to PE0.
                             IF ( ind(1) /= -9999 )  THEN
-                               CALL MPI_SEND( local_2d(nysg,nzb), ngp,         &
+                               CALL MPI_SEND( local_2d(nysg,nzb_do), ngp,         &
                                               MPI_REAL, 0, 1, comm2d, ierr )
                             ENDIF
@@ -1505 +1787 @@
                    nc_stat = NF90_PUT_VAR( id_set_yz(av),                   &
                                            id_var_do2d(av,if),              &
-                                           local_2d(nys:nyn+1,nzb:nzt+1),   &
+                                           local_2d(nys:nyn+1,nzb_do:nzt_do),   &
                             start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
-                                           count = (/ 1, ny+2, nz+2, 1 /) )
+                                           count = (/ 1, ny+2, nzt_do-nzb_do+1, 1 /) )
                    CALL handle_netcdf_error( 'data_output_2d', 452 )
 #endif
@@ -1595 +1877 @@
 !
 !--                      Distribute data over all PEs along y
-                         ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 ) * ns
+                         ngp = ( nxrg-nxlg+1 ) * ( nzt_do-nzb_do+1 ) * ns
                          IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
-                         CALL MPI_ALLREDUCE( local_2d_sections_l(nxlg,1,nzb),  &
-                                             local_2d_sections(nxlg,1,nzb),    &
+                         CALL MPI_ALLREDUCE( local_2d_sections_l(nxlg,1,nzb_do),  &
+                                             local_2d_sections(nxlg,1,nzb_do),    &
                                              ngp, MPI_REAL, MPI_SUM, comm1dy,  &
                                              ierr )
@@ -1612 +1894 @@
                                              id_var_do2d(av,if),               & 
                                              local_2d_sections(nxl:nxr+1,1:ns, &
-                                                nzb:nzt+1),                    &
+                                                nzb_do:nzt_do),                &
                                              start = (/ nxl+1, 1, 1,           &
                                                 do2d_xz_time_count(av) /),     &
-                                             count = (/ nxr-nxl+2, ns, nzt+2,  &
+                                             count = (/ nxr-nxl+2, ns, nzt_do-nzb_do+1,  &
                                                         1 /) )
                       ELSE
@@ -1621 +1903 @@
                                              id_var_do2d(av,if),               &
                                              local_2d_sections(nxl:nxr,1:ns,   &
-                                                nzb:nzt+1),                    &
+                                                nzb_do:nzt_do),                &
                                              start = (/ nxl+1, 1, 1,           &
                                                 do2d_xz_time_count(av) /),     &
-                                             count = (/ nxr-nxl+1, ns, nzt+2,  &
+                                             count = (/ nxr-nxl+1, ns, nzt_do-nzb_do+1,  &
                                                 1 /) )
                       ENDIF
@@ -1647 +1929 @@
                          ngp = ( nyng-nysg+1 ) * ( nzt-nzb + 2 ) * ns
                          IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
-                         CALL MPI_ALLREDUCE( local_2d_sections_l(1,nysg,nzb),  &
-                                             local_2d_sections(1,nysg,nzb),    &
+                         CALL MPI_ALLREDUCE( local_2d_sections_l(1,nysg,nzb_do),  &
+                                             local_2d_sections(1,nysg,nzb_do),    &
                                              ngp, MPI_REAL, MPI_SUM, comm1dx,  &
                                              ierr )
@@ -1662 +1944 @@
                                              id_var_do2d(av,if),               &
                                              local_2d_sections(1:ns,           &
-                                                nys:nyn+1,nzb:nzt+1),          &
+                                                nys:nyn+1,nzb_do:nzt_do),      &
                                              start = (/ 1, nys+1, 1,           &
                                                 do2d_yz_time_count(av) /),     &
                                              count = (/ ns, nyn-nys+2,         &
-                                                        nzt+2, 1 /) )
+                                                        nzt_do-nzb_do+1, 1 /) )
                       ELSE
                          nc_stat = NF90_PUT_VAR( id_set_yz(av),                &
                                              id_var_do2d(av,if),               &
                                              local_2d_sections(1:ns,nys:nyn,   &
-                                                nzb:nzt+1),                    &
+                                                nzb_do:nzt_do),                &
                                              start = (/ 1, nys+1, 1,           &
                                                 do2d_yz_time_count(av) /),     &
                                              count = (/ ns, nyn-nys+1,         &
-                                                        nzt+2, 1 /) )
+                                                        nzt_do-nzb_do+1, 1 /) )
                       ENDIF
 

palm/trunk/SOURCE/data_output_3d.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
+! Added suppport for land surface model and radiation model output. In the course
+! of this action, the limits for vertical loops have been changed (from nzb and 
+! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
+! Moreover, a new vertical grid zs was introduced.
 ! 
 ! Former revisions:
@@ -112 +115 @@
     USE kinds
     
+    USE land_surface_model_mod,                                                &
+        ONLY: m_soil, m_soil_av, nzb_soil, nzt_soil, t_soil, t_soil_av
+
     USE netcdf_control
         
@@ -130 +136 @@
     INTEGER(iwp) ::  k         !: 
     INTEGER(iwp) ::  n         !: 
+    INTEGER(iwp) ::  nzb_do    !: vertical lower limit for data output
+    INTEGER(iwp) ::  nzt_do    !: vertical upper limit for data output
     INTEGER(iwp) ::  pos       !: 
     INTEGER(iwp) ::  prec      !: 
@@ -183 +191 @@
 
 !
-!-- Allocate a temporary array with the desired output dimensions.
-    ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nz_do3d) )
-
-!
 !-- Update the netCDF time axis
 !-- In case of parallel output, this is only done by PE0 to increase the
@@ -209 +213 @@
 !--    Store the array chosen on the temporary array.
        resorted = .FALSE.
+       nzb_do = nzb
+       nzt_do = nz_do3d
+
+!
+!--    Allocate a temporary array with the desired output dimensions.
+       ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
+
        SELECT CASE ( TRIM( do3d(av,if) ) )
 
@@ -223 +234 @@
              ELSE
                 to_be_resorted => lpt_av
+             ENDIF
+
+          CASE ( 'm_soil' )
+             nzb_do = nzb_soil
+             nzt_do = nzt_soil
+!
+!--          For soil model quantities, it is required to re-allocate local_pf
+             DEALLOCATE ( local_pf )
+             ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
+
+             IF ( av == 0 )  THEN
+                to_be_resorted => m_soil
+             ELSE
+                to_be_resorted => m_soil_av
              ENDIF
 
@@ -251 +276 @@
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
-                      DO  k = nzb, nz_do3d
+                      DO  k = nzb_do, nzt_do
                          local_pf(i,j,k) = tend(k,j,i)
                       ENDDO
@@ -267 +292 @@
                    DO  i = nxl, nxr
                       DO  j = nys, nyn
-                         DO  k = nzb, nz_do3d
+                         DO  k = nzb_do, nzt_do
                             number_of_particles = prt_count(k,j,i)
                             IF (number_of_particles <= 0)  CYCLE
@@ -296 +321 @@
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
-                      DO  k = nzb, nz_do3d
+                      DO  k = nzb_do, nzt_do
                          local_pf(i,j,k) = tend(k,j,i)
                       ENDDO
@@ -336 +361 @@
                    DO  i = nxlg, nxrg
                       DO  j = nysg, nyng
-                         DO  k = nzb, nz_do3d
+                         DO  k = nzb_do, nzt_do
                             local_pf(i,j,k) = pt(k,j,i) + l_d_cp *             &
                                                           pt_d_t(k) *          &
@@ -389 +414 @@
                    DO  i = nxl, nxr
                       DO  j = nys, nyn
-                         DO  k = nzb, nz_do3d
+                         DO  k = nzb_do, nzt_do
                             number_of_particles = prt_count(k,j,i)
                             IF (number_of_particles <= 0)  CYCLE
@@ -409 +434 @@
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
-                      DO  k = nzb, nz_do3d
+                      DO  k = nzb_do, nzt_do
                          local_pf(i,j,k) = tend(k,j,i)
                       ENDDO
@@ -431 +456 @@
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
-                      DO  k = nzb, nz_do3d
+                      DO  k = nzb_do, nzt_do
                          local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
                       ENDDO
@@ -462 +487 @@
              ENDIF
 
+          CASE ( 't_soil' )
+             nzb_do = nzb_soil
+             nzt_do = nzt_soil
+!
+!--          For soil model quantities, it is required to re-allocate local_pf
+             DEALLOCATE ( local_pf )
+             ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
+
+             IF ( av == 0 )  THEN
+                to_be_resorted => t_soil
+             ELSE
+                to_be_resorted => t_soil_av
+             ENDIF
+
           CASE ( 'u' )
              IF ( av == 0 )  THEN
@@ -494 +533 @@
 !--          User defined quantity
              CALL user_data_output_3d( av, do3d(av,if), found, local_pf,       &
-                                       nz_do3d )
+                                       nzb_do, nzt_do )
              resorted = .TRUE.
 
@@ -510 +549 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                DO  k = nzb, nz_do3d
+                DO  k = nzb_do, nzt_do
                    local_pf(i,j,k) = to_be_resorted(k,j,i)
                 ENDDO
@@ -531 +570 @@
           DO  i = 0, io_blocks-1
              IF ( i == io_group )  THEN
-                WRITE ( 30 )  nxlg, nxrg, nysg, nyng, nzb, nz_do3d
-                WRITE ( 30 )  local_pf
+                WRITE ( 30 )  nxlg, nxrg, nysg, nyng, nzb_do, nzt_do
+                WRITE ( 30 )  local_pf(:,:,nzb_do:nzt_do)
              ENDIF
 #if defined( __parallel )
@@ -547 +586 @@
           IF ( nxr == nx  .AND.  nyn /= ny )  THEN
              nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
-                               local_pf(nxl:nxr+1,nys:nyn,nzb:nz_do3d),  &
-                start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
-                count = (/ nxr-nxl+2, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
+                               local_pf(nxl:nxr+1,nys:nyn,nzb_do:nzt_do),    &
+                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
+                count = (/ nxr-nxl+2, nyn-nys+1, nzt_do-nzb_do+1, 1 /) )
           ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
              nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
-                               local_pf(nxl:nxr,nys:nyn+1,nzb:nz_do3d),  &
-                start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
-                count = (/ nxr-nxl+1, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
+                               local_pf(nxl:nxr,nys:nyn+1,nzb_do:nzt_do),    &
+                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
+                count = (/ nxr-nxl+1, nyn-nys+2, nzt_do-nzb_do+1, 1 /) )
           ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
              nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
-                             local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d),  &
-                start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
-                count = (/ nxr-nxl+2, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
+                             local_pf(nxl:nxr+1,nys:nyn+1,nzb_do:nzt_do  ),  &
+                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
+                count = (/ nxr-nxl+2, nyn-nys+2, nzt_do-nzb_do+1, 1 /) )
           ELSE
              nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
-                                 local_pf(nxl:nxr,nys:nyn,nzb:nz_do3d),  &
-                start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
-                count = (/ nxr-nxl+1, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
+                                 local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do),    &
+                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
+                count = (/ nxr-nxl+1, nyn-nys+1, nzt_do-nzb_do+1, 1 /) )
           ENDIF
           CALL handle_netcdf_error( 'data_output_3d', 386 )
@@ -572 +611 @@
 #if defined( __netcdf )
        nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),        &
-                         local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d),      &
+                         local_pf(nxl:nxr+1,nys:nyn+1,nzb_do:nzt_do),        &
                          start = (/ 1, 1, 1, do3d_time_count(av) /),     &
-                         count = (/ nx+2, ny+2, nz_do3d-nzb+1, 1 /) )
+                         count = (/ nx+2, ny+2, nzt_do-nzb_do+1, 1 /) )
        CALL handle_netcdf_error( 'data_output_3d', 446 )
 #endif
@@ -581 +620 @@
        if = if + 1
 
+!
+!--    Deallocate temporary array
+       DEALLOCATE ( local_pf )
+
     ENDDO
-
-!
-!-- Deallocate temporary array.
-    DEALLOCATE( local_pf )
-
 
     CALL cpu_log( log_point(14), 'data_output_3d', 'stop' )

palm/trunk/SOURCE/flow_statistics.f90

@@ -21 +21 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added suppport for land surface model and radiation model output.
 !
 ! Former revisions:
@@ -132 +132 @@
         
     USE cloud_parameters,                                                      &
-        ONLY :  l_d_cp, prr, pt_d_t
+        ONLY:   l_d_cp, prr, pt_d_t
         
     USE control_parameters,                                                    &
-        ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
+        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,    &
@@ -143 +143 @@
         
     USE cpulog,                                                                &
-        ONLY :  cpu_log, log_point
+        ONLY:   cpu_log, log_point
         
     USE grid_variables,                                                        &
-        ONLY :  ddx, ddy
+        ONLY:   ddx, ddy
         
     USE indices,                                                               &
-        ONLY :  ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
+        ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
                 ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,        &
                 nzb_s_inner, nzt, nzt_diff
@@ -155 +155 @@
     USE kinds
     
+    USE land_surface_model_mod,                                                &
+        ONLY:   dots_soil, ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,   &
+                qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, shf_eb,       &
+                t_soil
+
     USE pegrid
+
+    USE radiation_model_mod,                                                   &
+        ONLY :  dots_rad, rad_net, rad_sw_in, radiation
     
     USE statistics
@@ -703 +711 @@
                                        qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
                 ENDIF
+             ENDIF
+
+             IF ( land_surface )  THEN
+                sums_l(nzb,93,tn) = sums_l(nzb,93,tn) + ghf_eb(j,i)
+                sums_l(nzb,94,tn) = sums_l(nzb,94,tn) + shf_eb(j,i)
+                sums_l(nzb,95,tn) = sums_l(nzb,95,tn) + qsws_eb(j,i)
+                sums_l(nzb,96,tn) = sums_l(nzb,96,tn) + qsws_liq_eb(j,i)
+                sums_l(nzb,97,tn) = sums_l(nzb,97,tn) + qsws_soil_eb(j,i)
+                sums_l(nzb,98,tn) = sums_l(nzb,98,tn) + qsws_veg_eb(j,i)
+             ENDIF
+
+             IF ( radiation )  THEN
+                sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  + rad_net(j,i)
+                sums_l(nzb,100,tn) = sums_l(nzb,100,tn) + rad_sw_in(j,i)
              ENDIF
 
@@ -1074 +1096 @@
        ENDIF
 
+
+       IF ( land_surface )  THEN
+          !$OMP DO
+          DO  i = nxl, nxr
+             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)
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+       
+
 !
 !--    Calculate the user-defined profiles
@@ -1133 +1169 @@
           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:pr_palm-2)    = sums(k,89:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
+          sums(k,89:100)           = sums(k,89:100)     / ngp_2dh(sr)
+          sums(k,101:pr_palm-2)    = sums(k,101:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
        ENDDO
 
@@ -1248 +1285 @@
        ENDIF
 
+       IF ( land_surface )  THEN
+          hom(:,1,89,sr) = sums(:,89)              ! t_soil
+                                                   ! 90 is initial t_soil profile
+          hom(:,1,91,sr) = sums(:,91)              ! m_soil
+                                                   ! 92 is initial m_soil profile
+          hom(:,1,93,sr) = sums(:,93)              ! ghf_eb
+          hom(:,1,94,sr) = sums(:,94)              ! shf_eb
+          hom(:,1,95,sr) = sums(:,95)              ! qsws_eb
+          hom(:,1,96,sr) = sums(:,96)              ! qsws_liq_eb
+          hom(:,1,97,sr) = sums(:,97)              ! qsws_soil_eb
+          hom(:,1,98,sr) = sums(:,98)              ! qsws_veg_eb
+       ENDIF
+
+       IF ( radiation )  THEN
+          hom(:,1,99 ,sr) = sums(:,99)             ! rad_net
+          hom(:,1,100,sr) = sums(:,100)            ! rad_sw_in
+       ENDIF
+
        hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
                                        ! upstream-parts u_x, u_y, u_z, v_x,
@@ -1392 +1447 @@
 
        ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
+
+!
+!--    Collect land surface model timeseries
+       IF ( land_surface )  THEN
+          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
+          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
+          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
+          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
+          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
+          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
+       ENDIF
+!
+!--    Collect radiation model timeseries
+       IF ( radiation )  THEN
+          ts_value(dots_rad,sr)   = hom(nzb,1,99,sr)           ! rad_net
+          ts_value(dots_rad+1,sr) = hom(nzb,1,100,sr)          ! rad_sw_in
+       ENDIF
+
 !
 !--    Calculate additional statistics provided by the user interface
@@ -2842 +2915 @@
        ENDIF
 
+
+       IF ( land_surface )  THEN
+          !$OMP DO
+          DO  i = nxl, nxr
+             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)
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+
+
 !
 !--    Calculate the user-defined profiles
@@ -2903 +2990 @@
           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:pr_palm-2)    = sums(k,89:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
+          sums(k,89:100)           = sums(k,89:100)     / ngp_2dh(sr)
+          sums(k,101:pr_palm-2)    = sums(k,101:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
        ENDDO
 
@@ -3162 +3250 @@
 
        ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
+
+!
+!--    Collect land surface model timeseries
+       IF ( land_surface )  THEN
+          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
+          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
+          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
+          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
+          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
+          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
+       ENDIF
+!
+!--    Collect radiation model timeseries
+       IF ( radiation )  THEN
+          ts_value(dots_rad,sr)   = hom(nzb,1,99,sr)           ! rad_net
+          ts_value(dots_rad+1,sr) = hom(nzb,1,100,sr)          ! rad_sw_in
+       ENDIF
+
 !
 !--    Calculate additional statistics provided by the user interface

palm/trunk/SOURCE/header.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added informal output for land surface model and radiation model. Removed typo.
 ! 
 ! Former revisions:
@@ -170 +170 @@
 ! Description:
 ! ------------
-! Writing a header with all important informations about the actual run.
+! Writing a header with all important information about the actual run.
 ! This subroutine is called three times, two times at the beginning 
 ! (writing information on files RUN_CONTROL and HEADER) and one time at the
@@ -200 +200 @@
         
     USE kinds
-    
+   
+    USE land_surface_model_mod,                                                &
+        ONLY:  conserve_water_content, dewfall, land_surface, nzb_soil,        &
+               nzt_soil, root_fraction, soil_moisture, soil_temperature,       &
+               soil_type, soil_type_name, veg_type, veg_type_name, zs
+ 
     USE model_1d,                                                              &
         ONLY:  damp_level_ind_1d, dt_pr_1d, dt_run_control_1d, end_time_1d
@@ -225 +230 @@
                lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, &
                plant_canopy
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  albedo, day_init, dt_radiation, lambda, net_radiation,          &
+               radiation, radiation_scheme, time_utc_init
     
     USE spectrum,                                                              &
@@ -263 +272 @@
     CHARACTER (LEN=86) ::  gradients           !:
     CHARACTER (LEN=86) ::  leaf_area_density   !:
+    CHARACTER (LEN=86) ::  roots               !:
     CHARACTER (LEN=86) ::  slices              !:
     CHARACTER (LEN=86) ::  temperatures        !:
@@ -311 +321 @@
 !
 !-- At the end of the run, output file (HEADER) will be rewritten with
-!-- new informations
+!-- new information
     IF ( io == 19 .AND. simulated_time_at_begin /= simulated_time ) REWIND( 19 )
 
@@ -495 +505 @@
 
 !
-!-- Runtime and timestep informations
+!-- Runtime and timestep information
     WRITE ( io, 200 )
     IF ( .NOT. dt_fixed )  THEN
@@ -850 +860 @@
 
 
+    IF ( land_surface )  THEN
+
+       temperatures = ''
+       gradients    = '' ! use for humidity here
+       coordinates  = '' ! use for height
+       roots        = '' ! use for root fraction
+       slices       = '' ! use for index
+
+       i = 1
+       DO i = nzb_soil, nzt_soil
+          WRITE (coor_chr,'(F10.2,7X)') soil_temperature(i)
+          temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F10.2,7X)') soil_moisture(i)
+          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F10.2,7X)')  - zs(i)
+          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F10.2,7X)')  root_fraction(i)
+          roots = TRIM( roots ) // ' '  // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(I10,7X)')  i
+          slices = TRIM( slices ) // ' '  // TRIM( coor_chr )
+
+
+       ENDDO
+
+!
+!--    Write land surface model header
+       WRITE( io, 419 )
+       IF ( conserve_water_content )  THEN
+          WRITE( io, 440 )
+       ELSE
+          WRITE( io, 441 )
+       ENDIF
+
+       IF ( dewfall )  THEN
+          WRITE( io, 442 )
+       ELSE
+          WRITE( io, 443 )
+       ENDIF
+
+       WRITE( io, 438 ) veg_type_name(veg_type), soil_type_name(soil_type)
+       WRITE( io, 439 ) TRIM( coordinates ), TRIM( temperatures ),             &
+                        TRIM( gradients ), TRIM( roots ), TRIM( slices )
+
+
+    ENDIF
+
+    IF ( radiation )  THEN
+!
+!--    Write land surface model header
+       WRITE( io, 444 )
+
+       IF ( radiation_scheme == "constant" )  THEN
+          WRITE( io, 445 ) net_radiation
+       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
+          WRITE( io, 446 )
+       ELSE
+          WRITE( io, 447 ) radiation_scheme
+       ENDIF
+
+       WRITE( io, 448 ) albedo
+       WRITE( io, 449 ) dt_radiation
+
+    ENDIF
+
+
 !
 !-- Boundary conditions
@@ -877 +956 @@
 
     IF ( ibc_pt_b == 0 )  THEN
-       runten = TRIM( runten ) // ' pt(0)   = pt_surface'
+       IF ( land_surface )  THEN
+          runten = TRIM( runten ) // ' pt(0)     = from soil model'
+       ELSE
+          runten = TRIM( runten ) // ' pt(0)     = pt_surface'
+       ENDIF
     ELSEIF ( ibc_pt_b == 1 )  THEN
-       runten = TRIM( runten ) // ' pt(0)   = pt(1)'
+       runten = TRIM( runten ) // ' pt(0)     = pt(1)'
     ELSEIF ( ibc_pt_b == 2 )  THEN
-       runten = TRIM( runten ) // ' pt(0) = from coupled model'
+       runten = TRIM( runten ) // ' pt(0)     = from coupled model'
     ENDIF
     IF ( ibc_pt_t == 0 )  THEN
@@ -918 +1001 @@
     IF ( humidity )  THEN
        IF ( ibc_q_b == 0 )  THEN
-          runten = 'q(0)     = q_surface'
+          IF ( land_surface )  THEN
+             runten = 'q(0)     = from soil model'
+          ELSE
+             runten = 'q(0)     = q_surface'
+          ENDIF
+
        ELSE
           runten = 'q(0)     = q(1)'
@@ -1225 +1313 @@
              coordinates = '/'
 !
-!--          Building strings with index and coordinate informations of the
+!--          Building strings with index and coordinate information of the
 !--          slices
              DO  WHILE ( section(i,1) /= -9999 )
@@ -1271 +1359 @@
              coordinates = '/'
 !
-!--          Building strings with index and coordinate informations of the
+!--          Building strings with index and coordinate information of the
 !--          slices
              DO  WHILE ( section(i,2) /= -9999 )
@@ -1313 +1401 @@
              coordinates = '/'
 !
-!--          Building strings with index and coordinate informations of the
+!--          Building strings with index and coordinate information of the
 !--          slices
              DO  WHILE ( section(i,3) /= -9999 )
@@ -1571 +1659 @@
 !
 !-- Geostrophic parameters
-    WRITE ( io, 410 )  omega, phi, f, fs
+    IF ( radiation .AND. radiation_scheme /= 'constant' )  THEN
+       WRITE ( io, 417 )  lambda
+    ENDIF
+    WRITE ( io, 410 )  phi, omega, f, fs
 
 !
 !-- Other quantities
     WRITE ( io, 411 )  g
+    IF ( radiation .AND. radiation_scheme /= 'constant' )  THEN
+       WRITE ( io, 418 )  day_init, time_utc_init
+    ENDIF
+
     WRITE ( io, 412 )  TRIM( reference_state )
     IF ( use_single_reference_value )  THEN
@@ -1732 +1827 @@
 
 !
-!-- User-defined informations
+!-- User-defined information
     CALL user_header( io )
 
@@ -1867 +1962 @@
 260 FORMAT (/' The model has a slope in x-direction. Inclination angle: ',F6.2,&
              ' degrees')
-270 FORMAT (//' Topography informations:'/ &
-              ' -----------------------'// &
+270 FORMAT (//' Topography information:'/ &
+              ' ----------------------'// &
               1X,'Topography: ',A)
 271 FORMAT (  ' Building size (x/y/z) in m: ',F5.1,' / ',F5.1,' / ',F5.1/ &
@@ -1905 +2000 @@
              ' -------------------'// &
              '                     p                    uv             ', &
-             '                   pt'// &
+             '                     pt'// &
              ' B. bound.: ',A/ &
              ' T. bound.: ',A)
@@ -2047 +2142 @@
 400 FORMAT (//' Physical quantities:'/ &
               ' -------------------'/)
-410 FORMAT ('    Angular velocity    :   omega = ',E9.3,' rad/s'/  &
-            '    Geograph. latitude  :   phi   = ',F4.1,' degr'/   &
-            '    Coriolis parameter  :   f     = ',F9.6,' 1/s'/    &
-            '                            f*    = ',F9.6,' 1/s')
-411 FORMAT (/'    Gravity             :   g     = ',F4.1,' m/s**2')
+410 FORMAT ('    Geograph. latitude  :   phi    = ',F4.1,' degr'/   &
+            '    Angular velocity    :   omega  = ',E9.3,' rad/s'/  &
+            '    Coriolis parameter  :   f      = ',F9.6,' 1/s'/    &
+            '                            f*     = ',F9.6,' 1/s')
+411 FORMAT (/'    Gravity             :   g      = ',F4.1,' m/s**2')
 412 FORMAT (/'    Reference state used in buoyancy terms: ',A)
 413 FORMAT ('       Reference density in buoyancy terms: ',F8.3,' kg/m**3')
 414 FORMAT ('       Reference temperature in buoyancy terms: ',F8.4,' K')
-415 FORMAT (/'    Cloud physics parameters:'/ &
-             '    ------------------------'/)
-416 FORMAT ('        Surface pressure   :   p_0   = ',F7.2,' hPa'/      &
-            '        Gas constant       :   R     = ',F5.1,' J/(kg K)'/ &
-            '        Density of air     :   rho_0 = ',F5.3,' kg/m**3'/  &
-            '        Specific heat cap. :   c_p   = ',F6.1,' J/(kg K)'/ &
-            '        Vapourization heat :   L_v   = ',E8.2,' J/kg')
+415 FORMAT (/' Cloud physics parameters:'/ &
+             ' ------------------------'/)
+416 FORMAT ('    Surface pressure   :   p_0   = ',F7.2,' hPa'/      &
+            '    Gas constant       :   R     = ',F5.1,' J/(kg K)'/ &
+            '    Density of air     :   rho_0 = ',F5.3,' kg/m**3'/  &
+            '    Specific heat cap. :   c_p   = ',F6.1,' J/(kg K)'/ &
+            '    Vapourization heat :   L_v   = ',E8.2,' J/kg')
+417 FORMAT ('    Geograph. longitude :   lambda = ',F4.1,' degr')
+418 FORMAT (/'    Day of the year at model start :   day_init      =     ',I3 &
+            /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
+419 FORMAT (//' Land surface model information:'/ &
+              ' ------------------------------'/)
 420 FORMAT (/'    Characteristic levels of the initial temperature profile:'// &
             '       Height:        ',A,'  m'/ &
@@ -2120 +2220 @@
                        '[0,1000] cm**2/s**3')
 437 FORMAT ('    Droplet collision is switched off')
+438 FORMAT (' --> Land surface type  : ',A,/ &
+            ' --> Soil porosity type : ',A)
+439 FORMAT (/'    Initial soil temperature and moisture profile:'// &
+            '       Height:        ',A,'  m'/ &
+            '       Temperature:   ',A,'  K'/ &
+            '       Moisture:      ',A,'  m**3/m**3'/ &
+            '       Root fraction: ',A,'  '/ &
+            '       Gridpoint:     ',A)
+440 FORMAT (/' --> Dewfall is allowed (default)')
+441 FORMAT (' --> Dewfall is inhibited')
+442 FORMAT (' --> Soil bottom is closed (water content is conserved, default)')
+443 FORMAT (' --> Soil bottom is open (water content is not conserved)')
+444 FORMAT (//' Radiation model information:'/                                 &
+              ' ----------------------------'/)
+445 FORMAT (' --> Using constant net radiation: net_radiation = ', F6.2, '  W/m**2')
+446 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,',  &
+                   ' default)')
+447 FORMAT (' --> Radiation scheme:', A)
+448 FORMAT (/'    Surface albedo: albedo = ', F5.3)
+449 FORMAT  ('    Timestep: dt_radiation = ', F5.2, '  s')
+
 450 FORMAT (//' LES / Turbulence quantities:'/ &
               ' ---------------------------'/)
@@ -2200 +2321 @@
 508 FORMAT ('    Ventilation effects on evaporation of rain drops')
 509 FORMAT ('    Slope limiter used for sedimentation process')
-510 FORMAT ('        Droplet density    :   N_c   = ',F6.1,' 1/cm**3')
-511 FORMAT ('        Sedimentation Courant number:                  '/&
+510 FORMAT ('    Droplet density    :   N_c   = ',F6.1,' 1/cm**3')
+511 FORMAT ('    Sedimentation Courant number:                  '/&
             '                               C_s   = ',F3.1,'        ')
 512 FORMAT (/' Date:                 ',A8,6X,'Run:       ',A20/      &

palm/trunk/SOURCE/init_3d_model.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
+! Allocation of land surface arrays is now done in the subroutine init_lsm_arrays,
+! which is part of land_surface_model.
 ! 
 ! Former revisions:
@@ -218 +219 @@
 
     USE land_surface_model_mod,                                                &
-        ONLY:  init_lsm, land_surface
+        ONLY:  init_lsm, init_lsm_arrays, land_surface
   
     USE ls_forcing_mod
@@ -621 +622 @@
 
 !
+!-- Allocate land surface model arrays
+    IF ( land_surface )  THEN
+       CALL init_lsm_arrays
+    ENDIF
+
+!
 !-- Allocate arrays containing the RK coefficient for calculation of 
 !-- perturbation pressure and turbulent fluxes. At this point values are

palm/trunk/SOURCE/land_surface_model.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Flux calculation is now done in prandtl_fluxes. Added support for data output.
+! Vertical indices have been replaced. Restart runs are now possible. Some
+! variables have beem renamed. Bugfix in the prognostic equation for the surface
+! temperature. Introduced z0_eb and z0h_eb, which overwrite the setting of
+! roughness_length and z0_factor. Added Clapp & Hornberger parametrization for
+! the hydraulic conductivity. Bugfix for root fraction and extraction
+! calculation
 ! 
 ! Former revisions:
@@ -43 +49 @@
 ! scheme implemented in the ECMWF IFS model, with modifications according to
 ! H-TESSEL. The implementation is based on the formulation implemented in the 
-! DALES model.
+! DALES and UCLA-LES models.
 !
 ! To do list:
 ! -----------
-! - Add support for binary I/O support
-! - Add support for lsm data output
-! - Check for time step criterion
-! - Check use with RK-2 and Euler time-stepping
-! - Adaption for use with cloud physics (liquid water potential temperature)
-! - Check reaction of plants at wilting point and at atmospheric saturation
-! - Consider partial absorption of the net shortwave radiation by the skin layer
+! - Check dewfall parametrization for fog simulations
+! - Consider partial absorption of the net shortwave radiation by the surface layer
 ! - Allow for water surfaces, check performance for bare soils 
+! - Invert indices (running from -3 to 0. Currently: nzb_soil=0, nzt_soil=3))
+! - Implement surface runoff model (required when performing long-term LES with
+!   considerable precipitation
+! Notes:
+! ------
+! - No time step criterion is required as long as the soil layers do not become
+!   too thin
 !------------------------------------------------------------------------------!
      USE arrays_3d,                                                            &
@@ -60 +68 @@
 
      USE cloud_parameters,                                                     &
-         ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
+         ONLY:  cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
 
      USE control_parameters,                                                   &
-         ONLY: dt_3d, humidity, intermediate_timestep_count,                   &
-               intermediate_timestep_count_max, precipitation, pt_surface,     &
-               rho_surface, surface_pressure, timestep_scheme, tsc
+         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
 
      USE indices,                                                              &
-         ONLY: nxlg, nxrg, nyng, nysg, nzb_s_inner 
+         ONLY:  nxlg, nxrg, nyng, nysg, nzb_s_inner 
 
      USE kinds
 
+    USE netcdf_control,                                                        &
+        ONLY:  dots_label, dots_num, dots_unit
+
      USE radiation_model_mod,                                                  &
-         ONLY: Rn, SW_in, sigma_SB
-
+         ONLY:  irad_scheme, rad_net, rad_sw_in, sigma_SB
+
+     USE statistics,                                                           &
+         ONLY:  hom, statistic_regions
 
     IMPLICIT NONE
@@ -80 +95 @@
 !
 !-- LSM model constants
-    INTEGER(iwp), PARAMETER :: soil_layers = 4 !: number of soil layers (fixed for now)
+    INTEGER(iwp), PARAMETER :: nzb_soil = 0, & !: bottom of the soil model (to be switched)
+                               nzt_soil = 3, & !: top of the soil model (to be switched)
+                               nzs = 4         !: number of soil layers (fixed for now)
+
+    INTEGER(iwp) :: dots_soil = 0  !: starting index for timeseries output 
+
+    INTEGER(iwp), DIMENSION(0:1) :: id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz,  &
+                                    id_dim_zs_3d, id_var_zs_xy,                & 
+                                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d
+                                    
+    INTEGER(iwp), DIMENSION(1:max_masks,0:1) :: id_dim_zs_mask, id_var_zs_mask
 
     REAL(wp), PARAMETER ::                     &
-              b_CH               = 6.04_wp,    & ! Clapp & Hornberger exponent
-              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil
+              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
+              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil   
               lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
               lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
               psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
-              rhoC_soil          = 2.19E6_wp,  & ! volumetric heat capacity of soil
-              rhoC_water         = 4.20E6_wp,  & ! volumetric heat capacity of water
+              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil
+              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water
               m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
 
@@ -99 +124 @@
 
     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
 
 !   value 9999999.9_wp -> generic available or user-defined value must be set
 !   otherwise -> no generic variable and user setting is optional
-    REAL(wp) :: alpha_VanGenuchten = 0.0_wp,            & !: NAMELIST alpha_VG
-                canopy_resistance_coefficient = 0.0_wp, & !: NAMELIST gD
-                C_skin   = 20000.0_wp,                  & !: Skin heat capacity
+    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !: NAMELIST alpha_vg
+                canopy_resistance_coefficient = 9999999.9_wp, & !: NAMELIST g_d
+                c_surface   = 20000.0_wp,               & !: Surface (skin) heat capacity
                 drho_l_lv,                              & !: (rho_l * l_v)**-1
                 exn,                                    & !: value of the Exner function
                 e_s = 0.0_wp,                           & !: saturation water vapour pressure
-                field_capacity = 0.0_wp,                & !: NAMELIST m_fc
-                f_shortwave_incoming = 9999999.9_wp,    & !: NAMELIST f_SW_in
-                hydraulic_conductivity = 0.0_wp,        & !: NAMELIST gamma_w_sat
-                Ke = 0.0_wp,                            & !: Kersten number
-                lambda_skin_stable = 9999999.9_wp,      & !: NAMELIST lambda_skin_s
-                lambda_skin_unstable = 9999999.9_wp,    & !: NAMELIST lambda_skin_u
-                leaf_area_index = 9999999.9_wp,         & !: NAMELIST LAI
-                l_VanGenuchten = 0.0_WP,                & !: NAMELIST l_VG
-                min_canopy_resistance = 110.0_wp,       & !: NAMELIST r_s_min
-                m_total = 0.0_wp,                       & !: weighed total water content of the soil (m3/m3)
-                n_VanGenuchten = 0.0_WP,                & !: NAMELIST n_VG
+                field_capacity = 9999999.9_wp,          & !: NAMELIST m_fc
+                f_shortwave_incoming = 9999999.9_wp,    & !: NAMELIST f_sw_in
+                hydraulic_conductivity = 9999999.9_wp,  & !: NAMELIST gamma_w_sat
+                ke = 0.0_wp,                            & !: Kersten number
+                lambda_surface_stable = 9999999.9_wp,   & !: NAMELIST lambda_surface_s
+                lambda_surface_unstable = 9999999.9_wp, & !: NAMELIST lambda_surface_u
+                leaf_area_index = 9999999.9_wp,         & !: NAMELIST lai
+                l_vangenuchten = 9999999.9_wp,          & !: NAMELIST l_vg
+                min_canopy_resistance = 9999999.9_wp,   & !: NAMELIST r_canopy_min
+                min_soil_resistance = 50.0_wp,          & !: NAMELIST r_soil_min
+                m_total = 0.0_wp,                       & !: weighted total water content of the soil (m3/m3)
+                n_vangenuchten = 9999999.9_wp,          & !: NAMELIST n_vg
                 q_s = 0.0_wp,                           & !: saturation specific humidity
-                residual_moisture = 0.0_wp,             & !: NAMELIST m_res
+                residual_moisture = 9999999.9_wp,       & !: NAMELIST m_res
                 rho_cp,                                 & !: rho_surface * cp
                 rho_lv,                                 & !: rho * l_v
                 rd_d_rv,                                & !: r_d / r_v
-                saturation_moisture = 0.0_wp,           & !: NAMELIST m_sat
+                saturation_moisture = 9999999.9_wp,     & !: NAMELIST m_sat
                 vegetation_coverage = 9999999.9_wp,     & !: NAMELIST c_veg
-                wilting_point = 0.0_wp                    !: NAMELIST m_wilt
-
-    REAL(wp), DIMENSION(0:soil_layers-1) :: &
+                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)
+
+    REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: &
               ddz_soil,                     & !: 1/dz_soil
               ddz_soil_stag,                & !: 1/dz_soil_stag
@@ -135 +164 @@
               dz_soil_stag,                 & !: soil grid spacing (edge-edge) 
               root_extr = 0.0_wp,           & !: root extraction
-              root_fraction = (/0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp/), & !: distribution of root surface area to the individual soil layers
-              soil_level = (/0.07_wp, 0.28_wp, 1.00_wp,  2.89_wp/),   & !: soil layer depths (m)
+              root_fraction = (/9999999.9_wp, 9999999.9_wp,    &
+                                9999999.9_wp, 9999999.9_wp /), & !: distribution of root surface area to the individual soil layers
+              zs = (/0.07_wp, 0.28_wp, 1.00_wp,  2.89_wp/),    & !: soil layer depths (m)
               soil_moisture = 0.0_wp          !: soil moisture content (m3/m3)
 
-    REAL(wp), DIMENSION(0:soil_layers) ::   &
+    REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) ::   &
               soil_temperature = 9999999.9_wp !: soil temperature (K)
 
 #if defined( __nopointer )
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0,    & !: skin temperature (K)
-                                                     T_0_p,  & !: progn. skin temperature (K)
-                                                     m_liq,  & !: liquid water reservoir (m)
-                                                     m_liq_p   !: progn. liquid water reservoir (m)
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface,   & !: surface temperature (K)
+                                                     t_surface_p, & !: progn. surface temperature (K)
+                                                     m_liq_eb,    & !: liquid water reservoir (m)
+                                                     m_liq_eb_av, & !: liquid water reservoir (m)
+                                                     m_liq_eb_p     !: progn. liquid water reservoir (m)
 #else
-    REAL(wp), DIMENSION(:,:), POINTER :: T_0,   &
-                                         T_0_p, & 
-                                         m_liq, & 
-                                         m_liq_p
-
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: T_0_1, T_0_2,    &
-                                                     m_liq_1, m_liq_2
+    REAL(wp), DIMENSION(:,:), POINTER :: t_surface,      &
+                                         t_surface_p,    & 
+                                         m_liq_eb,       & 
+                                         m_liq_eb_p
+
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, &
+                                                     m_liq_eb_av,              &
+                                                     m_liq_eb_1, m_liq_eb_2
 #endif
 
 !
 !-- Temporal tendencies for time stepping            
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tT_0_m,  & !: skin temperature tendency (K)
-                                             tm_liq_m   !: liquid water reservoir tendency (m)
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m,  & !: surface temperature tendency (K)
+                                             tm_liq_eb_m      !: liquid water reservoir tendency (m)
 
 !
 !-- Energy balance variables                
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::                                   &
-              alpha_VG,      & !: coef. of Van Genuchten
-              c_liq,         & !: liquid water coverage (of vegetated area)
-              c_veg,         & !: vegetation coverage   
-              f_SW_in,       & !: ?
-              G,             & !: surface soil heat flux
-              H,             & !: surface flux of sensible heat
-              gamma_w_sat,   & !: hydraulic conductivity at saturation
-              gD,            & !: coefficient for dependence of r_canopy on water vapour pressure deficit
-              LAI,           & !: leaf area index
-              LE,            & !: surface flux of latent heat (total)
-              LE_veg,        & !: surface flux of latent heat (vegetation portion)
-              LE_soil,       & !: surface flux of latent heat (soil portion)
-              LE_liq,        & !: surface flux of latent heat (liquid water portion)
-              lambda_h_sat,  & !: heat conductivity for dry soil
-              lambda_skin_s, & !: coupling between skin and soil (depends on vegetation type)
-              lambda_skin_u, & !: coupling between skin and soil (depends on vegetation type)
-              l_VG,          & !: coef. of Van Genuchten
-              m_fc,          & !: soil moisture at field capacity (m3/m3)
-              m_res,         & !: residual soil moisture
-              m_sat,         & !: saturation soil moisture (m3/m3)
-              m_wilt,        & !: soil moisture at permanent wilting point (m3/m3)
-              n_VG,          & !: coef. Van Genuchten  
-              r_a,           & !: aerodynamic resistance 
-              r_canopy,      & !: canopy resistance
-              r_soil,        & !: soil resitance
-              r_soil_min,    & !: minimum soil resistance
-              r_s,           & !: total surface resistance (combination of r_soil and r_canopy)          
-              r_s_min          !: minimum canopy resistance
+              alpha_vg,         & !: coef. of Van Genuchten
+              c_liq,            & !: liquid water coverage (of vegetated area)
+              c_liq_av,         & !: average of c_liq
+              c_soil_av,        & !: average of c_soil
+              c_veg,            & !: vegetation coverage
+              c_veg_av,         & !: average of c_veg
+              f_sw_in,          & !: fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
+              ghf_eb,           & !: ground heat flux
+              ghf_eb_av,        & !: average of ghf_eb
+              gamma_w_sat,      & !: hydraulic conductivity at saturation
+              g_d,              & !: coefficient for dependence of r_canopy on water vapour pressure deficit
+              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)
+              l_vg,             & !: coef. of Van Genuchten
+              m_fc,             & !: soil moisture at field capacity (m3/m3)
+              m_res,            & !: residual soil moisture
+              m_sat,            & !: saturation soil moisture (m3/m3)
+              m_wilt,           & !: soil moisture at permanent wilting point (m3/m3)
+              n_vg,             & !: coef. Van Genuchten  
+              qsws_eb,          & !: surface flux of latent heat (total)
+              qsws_eb_av,       & !: average of qsws_eb
+              qsws_liq_eb,      & !: surface flux of latent heat (liquid water portion)
+              qsws_liq_eb_av,   & !: average of qsws_liq_eb
+              qsws_soil_eb,     & !: surface flux of latent heat (soil portion)
+              qsws_soil_eb_av,  & !: average of qsws_soil_eb
+              qsws_veg_eb,      & !: surface flux of latent heat (vegetation portion)
+              qsws_veg_eb_av,   & !: average of qsws_veg_eb
+              r_a,              & !: aerodynamic resistance 
+              r_canopy,         & !: canopy resistance
+              r_soil,           & !: soil resitance
+              r_soil_min,       & !: minimum soil resistance
+              r_s,              & !: total surface resistance (combination of r_soil and r_canopy)          
+              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 ::                                 &
@@ -198 +240 @@
               lambda_w, &   !: hydraulic diffusivity of soil (?)
               gamma_w,  &   !: hydraulic conductivity of soil (?)
-              rhoC_total    !: volumetric heat capacity of the actual soil matrix (?) 
+              rho_c_total   !: volumetric heat capacity of the actual soil matrix (?) 
 
 #if defined( __nopointer )
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
-              T_soil,    & !: Soil temperature (K)
-              T_soil_p,  & !: Prog. soil temperature (K)
+              t_soil,    & !: Soil temperature (K)
+              t_soil_av, & !: Average of t_soil
+              t_soil_p,  & !: Prog. soil temperature (K)
               m_soil,    & !: Soil moisture (m3/m3)
+              m_soil_av, & !: Average of m_soil
               m_soil_p     !: Prog. soil moisture (m3/m3)
 #else
     REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
-              T_soil, T_soil_p, &
+              t_soil, t_soil_p, &
               m_soil, m_soil_p    
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
-              T_soil_1, T_soil_2, &
-              m_soil_1, m_soil_2
+              t_soil_av, t_soil_1, t_soil_2,                                   &
+              m_soil_av, m_soil_1, m_soil_2
 
 
@@ -220 +264 @@
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
-              tT_soil_m, & !: T_soil storage array 
+              tt_soil_m, & !: t_soil storage array 
               tm_soil_m, & !: m_soil storage array 
               root_fr      !: root fraction (sum=1) 
 
-!
-!--  Land surface parameters according to the following classes (veg_type)
-!--  (0 user defined)
-!--  1 crops, mixed farming
-!--  2 short grass
-!--  3 evergreen needleleaf trees
-!--  4 deciduous needleleaf trees
-!--  5 evergreen broadleaf trees
-!--  6 deciduous broadleaf trees
-!--  7 tall grass
-!--  8 desert
-!--  9 tundra
-!-- 10 irrigated crops
-!-- 11 semidesert 
-!-- 12 ice caps and glaciers 
-!-- 13 bogs and marshes 
-!-- 14 inland water 
-!-- 15 ocean 
-!-- 16 evergreen shrubs 
-!-- 17 deciduous shrubs 
-!-- 18 mixed forest/woodland 
-!-- 19 interrupted forest 
-
-!
-!-- Land surface parameters I     r_s_min,     LAI,   c_veg,      gD
+
+!
+!-- Predefined Land surface classes (veg_type)
+    CHARACTER(26), DIMENSION(0:19) :: veg_type_name = (/          &
+                                   'user defined',                & ! 0 
+                                   'crops, mixed farming',        & !  1
+                                   'short grass',                 & !  2
+                                   'evergreen needleleaf trees',  & !  3
+                                   'deciduous needleleaf trees',  & !  4
+                                   'evergreen broadleaf trees' ,  & !  5
+                                   'deciduous broadleaf trees',   & !  6
+                                   'tall grass',                  & !  7
+                                   'desert',                      & !  8
+                                   'tundra',                      & !  9
+                                   'irrigated crops',             & ! 10
+                                   'semidesert',                  & ! 11
+                                   'ice caps and glaciers' ,      & ! 12
+                                   'bogs and marshes',            & ! 13
+                                   'inland water',                & ! 14
+                                   'ocean',                       & ! 15
+                                   'evergreen shrubs',            & ! 16
+                                   'deciduous shrubs',            & ! 17
+                                   'mixed forest/woodland',       & ! 18
+                                   'interrupted forest'           & ! 19
+                                                       /)
+
+!
+!-- Soil model classes (soil_type)
+    CHARACTER(12), DIMENSION(0:7) :: soil_type_name = (/ &
+                                   'user defined',        & ! 0 
+                                   'coarse',              & ! 1
+                                   'medium',              & ! 2
+                                   'medium-fine',         & ! 3
+                                   'fine',                & ! 4
+                                   'very fine' ,          & ! 5
+                                   'organic',             & ! 6
+                                   'loamy (CH)'           & ! 7
+                                                        /)
+!
+!-- Land surface parameters according to the respective classes (veg_type)
+
+!
+!-- Land surface parameters I
+!--                          r_canopy_min,     lai,   c_veg,     g_d
     REAL(wp), DIMENSION(0:3,1:19) :: veg_pars = RESHAPE( (/           &
                                  180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & !  1
@@ -257 +320 @@
                                  240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, & !  6
                                  100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, & !  7
-                                 250.0_wp, 0.50_wp, 0.00_wp, 0.00_wp, & !  8
+                                 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp, & !  8
                                   80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, & !  9
                                  180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 10
@@ -296 +359 @@
 
 !
-!-- Land surface parameters III lambda_skin_s, lambda_skin_u, f_SW_in
-    REAL(wp), DIMENSION(0:2,1:19) :: skin_pars = RESHAPE( (/           &
+!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
+    REAL(wp), DIMENSION(0:2,1:19) :: surface_pars = RESHAPE( (/           &
                                       10.0_wp,       10.0_wp, 0.05_wp, & !  1
                                       10.0_wp,       10.0_wp, 0.05_wp, & !  2
@@ -345 +408 @@
 !
 !-- Soil parameters according to the following porosity classes (soil_type)
-!-- (0 user defined)
-!-- 1 coarse
-!-- 2 medium
-!-- 3 medium-fine
-!-- 4 fine
-!-- 5 very fine
-!-- 6 organic
-!
-!-- Soil parameters I           alpha_VG,      l_VG,    n_VG, gamma_w_sat
-    REAL(wp), DIMENSION(0:3,1:6) :: soil_pars = RESHAPE( (/                &
+
+!
+!-- Soil parameters I           alpha_vg,      l_vg,    n_vg, gamma_w_sat
+    REAL(wp), DIMENSION(0:3,1:7) :: soil_pars = RESHAPE( (/                &
                                  3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
                                  3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
@@ -360 +417 @@
                                  3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, & ! 4
                                  2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, & ! 5
-                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp  & ! 6
-                                 /), (/ 4, 6 /) )
+                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, & ! 6
+                                 0.00_wp,  0.00_wp,  0.00_wp,  0.57E-6_wp  & ! 7
+                                 /), (/ 4, 7 /) )
 
 !
 !-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res  
-    REAL(wp), DIMENSION(0:3,1:6) :: m_soil_pars = RESHAPE( (/            &
+    REAL(wp), DIMENSION(0:3,1:7) :: m_soil_pars = RESHAPE( (/            &
                                  0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
                                  0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
@@ -371 +429 @@
                                  0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
                                  0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
-                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp  & ! 6
-                                 /), (/ 4, 6 /) )
+                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6
+                                 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp  & ! 7
+                                 /), (/ 4, 7 /) )
 
 
@@ -381 +440 @@
 
 
-    PUBLIC alpha_VanGenuchten, C_skin, canopy_resistance_coefficient,          &
-           conserve_water_content,      field_capacity, f_shortwave_incoming,  &
-           hydraulic_conductivity, init_lsm, lambda_skin_stable,               &
-           lambda_skin_unstable, land_surface, leaf_area_index,                &
-           lsm_energy_balance, lsm_soil_model, l_VanGenuchten,                 &
-           min_canopy_resistance, n_VanGenuchten, residual_moisture,           &
-           root_fraction, saturation_moisture, soil_level, soil_moisture,      &
-           soil_temperature, soil_type, vegetation_coverage, veg_type,         &
-           wilting_point
-
-#if defined( __nopointer )
-    PUBLIC m_liq, m_liq_p, m_soil, m_soil_p, T_0, T_0_p, T_soil, T_soil_p
-#else
-    PUBLIC m_liq, m_liq_1, m_liq_2, m_liq_p, m_soil, m_soil_1, m_soil_2,       &
-           m_soil_p, T_0, T_0_1, T_0_2, T_0_p, T_soil, T_soil_1, T_soil_2,     &
-           T_soil_p
-#endif
+!
+!-- Public parameters, constants and initial values
+    PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient,       &
+           conserve_water_content, dewfall, field_capacity,                    &
+           f_shortwave_incoming, hydraulic_conductivity, init_lsm,             &
+           init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable,    &
+           land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model,  &
+           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
+
+!
+!-- Public grid and NetCDF variables
+    PUBLIC dots_soil, id_dim_zs_xy, id_dim_zs_xz, id_dim_zs_yz,                &
+           id_dim_zs_3d, id_dim_zs_mask, id_var_zs_xy, id_var_zs_xz,           &
+           id_var_zs_yz, id_var_zs_3d, id_var_zs_mask, nzb_soil, nzs, nzt_soil,&
+           zs
+
+
+!
+!-- Public 2D output variables
+    PUBLIC c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb, ghf_eb_av,     &
+           lai, lai_av, qsws_eb, qsws_eb_av, qsws_liq_eb, qsws_liq_eb_av,      &
+           qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av,         &
+           shf_eb, shf_eb_av
+
+
+!
+!-- Public prognostic variables
+    PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
 
 
@@ -412 +486 @@
     END INTERFACE lsm_soil_model
 
+    INTERFACE lsm_swap_timelevel
+       MODULE PROCEDURE lsm_swap_timelevel
+    END INTERFACE lsm_swap_timelevel
 
  CONTAINS
 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!-- Allocate land surface model arrays and define pointers
+!------------------------------------------------------------------------------!
+    SUBROUTINE init_lsm_arrays
+    
+
+       IMPLICIT NONE
+
+!
+!--    Allocate surface and soil temperature / humidity
+#if defined( __nopointer )
+       ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
+#else
+       ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
+#endif
+
+!
+!--    Allocate intermediate timestep arrays
+       ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+
+!
+!--    Allocate 2D vegetation model arrays
+       ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) )
+       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) )
+       ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
+
+#if ! defined( __nopointer )
+!
+!-- Initial assignment of the pointers
+       t_soil    => t_soil_1;    t_soil_p    => t_soil_2
+       t_surface => t_surface_1; t_surface_p => t_surface_2
+       m_soil    => m_soil_1;    m_soil_p    => m_soil_2
+       m_liq_eb  => m_liq_eb_1;  m_liq_eb_p  => m_liq_eb_2
+#endif
+
+
+    END SUBROUTINE init_lsm_arrays
 
 !------------------------------------------------------------------------------!
@@ -432 +591 @@
 
 !
+!--    Calculate Exner function
+       exn       = ( surface_pressure / 1000.0_wp )**0.286_wp
+
+
+!
+!--    If no cloud physics is used, rho_surface has not been calculated before
+       IF ( .NOT. cloud_physics )  THEN
+          rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
+       ENDIF
+
+!
 !--    Calculate frequently used parameters
        rho_cp    = cp * rho_surface
        rd_d_rv   = r_d / r_v
        rho_lv    = rho_surface * l_v
-       drho_l_lv = 1.0 / (rho_l * l_v)
-
-!
-!--    Allocate skin and soil temperature / humidity
-#if defined( __nopointer )
-       ALLOCATE ( T_0(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( T_0_p(nysg:nyng,nxlg:nxrg) )
-#else
-       ALLOCATE ( T_0_1(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( T_0_2(nysg:nyng,nxlg:nxrg) )
-#endif
-
-       ALLOCATE ( tT_0_m(nysg:nyng,nxlg:nxrg) )
-
-#if defined( __nopointer )
-       ALLOCATE ( T_soil(0:soil_layers,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( T_soil_p(0:soil_layers,nysg:nyng,nxlg:nxrg) )
-#else
-       ALLOCATE ( T_soil_1(0:soil_layers,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( T_soil_2(0:soil_layers,nysg:nyng,nxlg:nxrg) )
-#endif
-
-       ALLOCATE ( tT_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
-
-#if defined( __nopointer )
-       ALLOCATE ( m_liq(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_liq_p(nysg:nyng,nxlg:nxrg) )
-#else
-       ALLOCATE ( m_liq_1(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_liq_2(nysg:nyng,nxlg:nxrg) )
-#endif
-
-       ALLOCATE ( tm_liq_m(nysg:nyng,nxlg:nxrg) )
-
-#if defined( __nopointer )
-       ALLOCATE ( m_soil(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_soil_p(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
-#else
-       ALLOCATE ( m_soil_1(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_soil_2(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
-#endif
-
-       ALLOCATE ( tm_soil_m(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
-
-
-#if ! defined( __nopointer )
-!
-!--    Initial assignment of the pointers
-       T_soil => T_soil_1; T_soil_p => T_soil_2
-       T_0 => T_0_1; T_0_p => T_0_2
-       m_soil => m_soil_1; m_soil_p => m_soil_2
-       m_liq => m_liq_1; m_liq_p => m_liq_2
-#endif
-
-       T_0    = 0.0_wp
-       T_0_p  = 0.0_wp
-       tT_0_m = 0.0_wp
-
-       T_soil    = 0.0_wp
-       T_soil_p  = 0.0_wp
-       tT_soil_m = 0.0_wp
-
-       m_liq    = 0.0_wp
-       m_liq_p  = 0.0_wp
-       tm_liq_m = 0.0_wp
-
-       m_soil    = 0.0_wp
-       m_soil_p  = 0.0_wp
-       tm_soil_m = 0.0_wp
-
-!
-!--    Allocate 2D vegetation model arrays
-       ALLOCATE ( alpha_VG(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( f_SW_in(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( G(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( H(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( gD(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( LAI(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( LE(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( LE_veg(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( LE_soil(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( LE_liq(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( l_VG(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( lambda_skin_u(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( lambda_skin_s(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( n_VG(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( r_s_min(nysg:nyng,nxlg:nxrg) )
-
-!
-!--    Set initial and default values
-       c_liq   = 0.0_wp
-       c_veg   = 0.0_wp
-       f_SW_in = 0.05_wp
-       gD      = 0.0_wp
-       LAI     = 0.0_wp
-       lambda_skin_u = 10.0_wp
-       lambda_skin_s = 10.0_wp
-
-
-       G       = 0.0_wp
-       H       = rho_cp * shf
+       drho_l_lv = 1.0_wp / (rho_l * l_v)
+
+!
+!--    Set inital values for prognostic quantities
+       tt_surface_m = 0.0_wp
+       tt_soil_m    = 0.0_wp
+       tm_liq_eb_m  = 0.0_wp
+       c_liq        = 0.0_wp
+
+       ghf_eb = 0.0_wp
+       shf_eb = rho_cp * shf
 
        IF ( humidity )  THEN
-          LE = rho_l * l_v * qsws
+          qsws_eb = rho_l * l_v * qsws
        ELSE
-          LE = 0.0_wp
+          qsws_eb = 0.0_wp
        ENDIF
 
-       LE_veg  = 0.0_wp
-       LE_soil = LE
-       LE_liq  = 0.0_wp
+       qsws_liq_eb  = 0.0_wp
+       qsws_soil_eb = qsws_eb
+       qsws_veg_eb  = 0.0_wp
 
        r_a        = 50.0_wp
+       r_s        = 50.0_wp
        r_canopy   = 0.0_wp
        r_soil     = 0.0_wp
-       r_soil_min = 50.0_wp
-       r_s        = 110.0_wp
-       r_s_min    = min_canopy_resistance
 
 !
 !--    Allocate 3D soil model arrays
-       ALLOCATE ( root_fr(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( lambda_h(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( rhoC_total(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
 
        lambda_h = 0.0_wp
@@ -578 +643 @@
 !--    If required, allocate humidity-related variables for the soil model
        IF ( humidity )  THEN
-          ALLOCATE ( lambda_w(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )
-          ALLOCATE ( gamma_w(0:soil_layers-1,nysg:nyng,nxlg:nxrg) )   
+          ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+          ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
 
           lambda_w = 0.0_wp 
@@ -588 +653 @@
 !--    the center of the soil layers, whereas gradients/fluxes are defined
 !--    at the edges (_stag)
-       dz_soil_stag(0) = soil_level(0)
-
-       DO k = 1, soil_layers-1
-          dz_soil_stag(k) = soil_level(k) - soil_level(k-1)
+       dz_soil_stag(nzb_soil) = zs(nzb_soil)
+
+       DO k = nzb_soil+1, nzt_soil
+          dz_soil_stag(k) = zs(k) - zs(k-1)
        ENDDO
 
-       DO k = 0, soil_layers-2
+       DO k = nzb_soil, nzt_soil-1
           dz_soil(k) = 0.5 * (dz_soil_stag(k+1) + dz_soil_stag(k))
        ENDDO
-       dz_soil(soil_layers-1) = dz_soil_stag(soil_layers-1)
-
-       ddz_soil      = 1.0 / dz_soil
-       ddz_soil_stag = 1.0 / dz_soil_stag
-!
-!--    Initialize soil
-       IF ( soil_type .NE. 0 )  THEN   
-          alpha_VG    = soil_pars(0,soil_type)
-          l_VG        = soil_pars(1,soil_type)
-          n_VG        = soil_pars(2,soil_type)    
-          gamma_w_sat = soil_pars(3,soil_type) 
-          m_sat       = m_soil_pars(0,soil_type)
-          m_fc        = m_soil_pars(1,soil_type)   
-          m_wilt      = m_soil_pars(2,soil_type) 
-          m_res       = m_soil_pars(3,soil_type)
+       dz_soil(nzt_soil) = dz_soil_stag(nzt_soil)
+
+       ddz_soil      = 1.0_wp / dz_soil
+       ddz_soil_stag = 1.0_wp / dz_soil_stag
+
+!
+!--    Initialize standard soil types. It is possible to overwrite each
+!--    parameter by setting the respecticy NAMELIST variable to a 
+!--    value /= 9999999.9. 
+       IF ( soil_type /= 0 )  THEN  
+ 
+          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
+             alpha_vangenuchten = soil_pars(0,soil_type)
+          ENDIF
+
+          IF ( l_vangenuchten == 9999999.9_wp )  THEN
+             l_vangenuchten = soil_pars(1,soil_type)
+          ENDIF
+
+          IF ( n_vangenuchten == 9999999.9_wp )  THEN
+             n_vangenuchten = soil_pars(2,soil_type)            
+          ENDIF
+
+          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
+             hydraulic_conductivity = soil_pars(3,soil_type)            
+          ENDIF
+
+          IF ( saturation_moisture == 9999999.9_wp )  THEN
+             saturation_moisture = m_soil_pars(0,soil_type)           
+          ENDIF
+
+          IF ( field_capacity == 9999999.9_wp )  THEN
+             field_capacity = m_soil_pars(1,soil_type)           
+          ENDIF
+
+          IF ( wilting_point == 9999999.9_wp )  THEN
+             wilting_point = m_soil_pars(2,soil_type)            
+          ENDIF
+
+          IF ( residual_moisture == 9999999.9_wp )  THEN
+             residual_moisture = m_soil_pars(3,soil_type)       
+          ENDIF
+
+       ENDIF    
+
+       alpha_vg      = alpha_vangenuchten
+       l_vg          = l_vangenuchten
+       n_vg          = n_vangenuchten 
+       gamma_w_sat   = hydraulic_conductivity
+       m_sat         = saturation_moisture
+       m_fc          = field_capacity
+       m_wilt        = wilting_point
+       m_res         = residual_moisture
+       r_soil_min    = min_soil_resistance
+
+!
+!--    Initial run actions
+       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+
+          t_soil    = 0.0_wp
+          m_liq_eb  = 0.0_wp
+          m_soil    = 0.0_wp
+
+!
+!--       Map user settings of T and q for each soil layer
+!--       (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
+             t_soil(k,:,:)    = soil_temperature(k)
+             m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
+             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
+          ENDDO
+          t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
+
+!
+!--       Calculate surface temperature
+          t_surface = pt_surface * exn
+
+!
+!--       Set artifical values for ts and us so that r_a has its initial value for 
+!--       the first time step
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                k = nzb_s_inner(j,i)
+
+!
+!--             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)
+                shf(j,i) = - us(j,i) * ts(j,i)
+             ENDDO
+          ENDDO
+
+!
+!--    Actions for restart runs
        ELSE
-          alpha_VG    = alpha_VanGenuchten
-          l_VG        = l_VanGenuchten
-          n_VG        = n_VanGenuchten 
-          gamma_w_sat = hydraulic_conductivity
-          m_sat       = saturation_moisture
-          m_fc        = field_capacity
-          m_wilt      = wilting_point
-          m_res       = residual_moisture
-       ENDIF    
-
-!
-!--    Map user settings of T and q for each soil layer
-!--    (make sure that the soil moisture does not drop below the permanent
-!--    wilting point) -> problems with devision by zero)
-       DO k = 0, soil_layers-1
-          T_soil(k,:,:)  = soil_temperature(k)
-          m_soil(k,:,:)  = MAX(soil_moisture(k),m_wilt(:,:))
-       ENDDO
-       T_soil(soil_layers,:,:) = soil_temperature(soil_layers)
-
-
-       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
-       T_0  = pt_surface * exn
-
-       T_soil_p = T_soil
-       m_soil_p = m_soil
+
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                k = nzb_s_inner(j,i)                
+                t_surface(j,i) = pt(k,j,i) * exn
+             ENDDO
+          ENDDO
+
+       ENDIF
 
 !
@@ -645 +776 @@
                            lambda_h_water ** m_sat(:,:)
 
+
+
+
+       DO k = nzb_soil, nzt_soil
+          root_fr(k,:,:) = root_fraction(k)
+       ENDDO
+
+       IF ( veg_type /= 0 )  THEN
+          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
+             min_canopy_resistance = veg_pars(0,veg_type)
+          ENDIF
+          IF ( leaf_area_index == 9999999.9_wp )  THEN
+             leaf_area_index = veg_pars(1,veg_type)         
+          ENDIF
+          IF ( vegetation_coverage == 9999999.9_wp )  THEN
+             vegetation_coverage = veg_pars(2,veg_type)     
+          ENDIF
+          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
+              canopy_resistance_coefficient= veg_pars(3,veg_type)     
+          ENDIF
+          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
+             lambda_surface_stable = surface_pars(0,veg_type)         
+          ENDIF
+          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
+             lambda_surface_unstable = surface_pars(1,veg_type)        
+          ENDIF
+          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
+             f_shortwave_incoming = surface_pars(2,veg_type)        
+          ENDIF
+          IF ( z0_eb == 9999999.9_wp )  THEN
+             roughness_length = roughness_par(0,veg_type) 
+             z0_eb            = roughness_par(0,veg_type) 
+          ENDIF
+          IF ( z0h_eb == 9999999.9_wp )  THEN
+             z0h_eb = roughness_par(1,veg_type)
+          ENDIF
+          z0h_factor = z0h_eb / z0_eb
+
+          IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
+             DO k = nzb_soil, nzt_soil
+                root_fr(k,:,:) = root_distribution(k,veg_type)
+                root_fraction(k) = root_distribution(k,veg_type)
+             ENDDO
+          ENDIF
+
+       ENDIF
+
 !
 !--    Initialize vegetation
-       IF ( veg_type .NE. 0 )  THEN
-
-          r_s_min              = veg_pars(0,veg_type)
-          LAI                  = veg_pars(1,veg_type)
-          c_veg                = veg_pars(2,veg_type)
-          gD                   = veg_pars(3,veg_type)
-          lambda_skin_s        = skin_pars(0,veg_type)
-          lambda_skin_u        = skin_pars(1,veg_type)
-          f_SW_in              = skin_pars(2,veg_type)
-          z0                   = roughness_par(0,veg_type)
-          z0h                  = roughness_par(1,veg_type)
-
-
-          DO k = 0, soil_layers-1
-             root_fr(k,:,:) = root_distribution(k,veg_type)
-          ENDDO
-
-       ELSE
-
-          DO k = 0, soil_layers-1
-             root_fr(k,:,:) = root_fraction(k)
-          ENDDO
-
-       ENDIF
+       r_canopy_min         = min_canopy_resistance
+       lai                  = leaf_area_index
+       c_veg                = vegetation_coverage
+       g_d                  = canopy_resistance_coefficient
+       lambda_surface_s     = lambda_surface_stable
+       lambda_surface_u     = lambda_surface_unstable
+       f_sw_in              = f_shortwave_incoming
+       z0                   = z0_eb
+       z0h                  = z0h_eb
 
 !
@@ -676 +839 @@
        CALL user_init_land_surface
 
-!
-!--    Set artifical values for ts and us so that r_a has its initial value for 
-!--    the first time step
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             k = nzb_s_inner(j,i)
-!
-!--          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)
-             shf(j,i) = - us(j,i) * ts(j,i)
-          ENDDO
-       ENDDO
+
+       t_soil_p    = t_soil
+       m_soil_p    = m_soil
+       m_liq_eb_p  = m_liq_eb
+
+!--    Store initial profiles of t_soil and m_soil (assuming they are 
+!--    horizontally homogeneous on this PE)
+       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
+                                                nysg,nxlg), 2,                 &
+                                                statistic_regions+1 )
+       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
+                                                nysg,nxlg), 2,                 &
+                                                statistic_regions+1 )
 
 !
 !--    Calculate humidity at the surface
        IF ( humidity )  THEN
-          CALL calc_q0
+          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_soil = dots_num + 1
+       dots_num  = dots_num + 6
+
 
        RETURN
@@ -707 +884 @@
 ! Description:
 ! ------------
-!
+! 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
@@ -730 +913 @@
                    coef_1,      & !: coef. for prognostic equation
                    coef_2,      & !: coef. for prognostic equation
-                   f_LE,        & !: factor for LE
-                   f_LE_veg,    & !: factor for LE_veg
-                   f_LE_soil,   & !: factor for LE_soil
-                   f_LE_liq,    & !: factor for LE_liq
-                   f_H,         & !: factor for H
-                   lambda_skin, & !: Current value of lambda_skin
-                   m_liq_max      !: maxmimum value of the liquid water reservoir
+                   f_qsws,      & !: factor for qsws_eb
+                   f_qsws_veg,  & !: factor for qsws_veg_eb
+                   f_qsws_soil, & !: factor for qsws_soil_eb
+                   f_qsws_liq,  & !: factor for qsws_liq_eb
+                   f_shf,       & !: factor for shf_eb
+                   lambda_surface, & !: Current value of lambda_surface
+                   m_liq_eb_max   !: maxmimum value of the liq. water reservoir
+
 
 !
@@ -748 +932 @@
 
 !
-!--          Set lambda_skin according to stratification
+!--          Set lambda_surface according to stratification
              IF ( rif(j,i) >= 0.0_wp )  THEN
-                lambda_skin = lambda_skin_s(j,i)
+                lambda_surface = lambda_surface_s(j,i)
              ELSE
-                lambda_skin = lambda_skin_u(j,i)
+                lambda_surface = lambda_surface_u(j,i)
              ENDIF
 
@@ -760 +944 @@
 !--          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)
+             r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) +       &
+                         1.0E-20_wp)
 
 !
@@ -766 +951 @@
 !--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
  
-!--          f1: correction for incoming shortwave radiation
-             f1 = MIN(1.0_wp, ( 0.004_wp * SW_in(j,i) + 0.05_wp ) /     &
-                              (0.81_wp * (0.004_wp * SW_in(j,i) + 1.0_wp) ) )
-
-!
-!--          f2: correction for soil moisture f2=0 for very dry soil
+!--          f1: correction for incoming shortwave radiation (stomata close at 
+!--          night)
+             IF ( irad_scheme /= 0 )  THEN
+                f1 = MIN(1.0_wp, ( 0.004_wp * rad_sw_in(j,i) + 0.05_wp ) /     &
+                              (0.81_wp * (0.004_wp * rad_sw_in(j,i) + 1.0_wp) ))
+             ELSE
+                f1 = 1.0_wp
+             ENDIF
+
+!
+!--          f2: correction for soil moisture availability to plants (the 
+!--          integrated soil moisture must thus be considered here)
+!--          f2 = 0 for very dry soils
              m_total = 0.0_wp
-             DO ks = 0, soil_layers-1
-                 m_total = m_total + root_fr(ks,j,i) * m_soil(ks,j,i)
+             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 .LE. m_fc(j,i) )  THEN
+             IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LT. 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
+                f2 = 1.0_wp
              ELSE
                 f2 = 1.0E-20_wp
@@ -785 +980 @@
 !
 !--          Calculate water vapour pressure at saturation 
-!--          (T_0 should be replaced by liquid water temp?!)
-             e_s = 0.01 * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) - 273.16_wp )&
-                                           / ( T_0(j,i) - 35.86_wp ) )
+             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 ) )
 
 !
 !--          f3: correction for vapour pressure deficit
-             IF ( gD(j,i) .NE. 0.0_wp )  THEN
+             IF ( g_d(j,i) /= 0.0_wp )  THEN
 !
 !--             Calculate vapour pressure
-                e  = q_p(k+1,j,i) * surface_pressure / 0.622
-                f3 = EXP ( -gD(j,i) * (e_s - e) )
+                e  = q_p(k+1,j,i) * surface_pressure / 0.622_wp
+                f3 = EXP ( -g_d(j,i) * (e_s - e) )
              ELSE
                 f3 = 1.0_wp
@@ -801 +995 @@
 
 !
+!--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
+!--          this calculation is obsolete, as r_canopy is not used below.
 !--          To do: check for very dry soil -> r_canopy goes to infinity
-             r_canopy(j,i)  = r_s_min(j,i) / (LAI(j,i) * f1 * f2 * f3 + 1.0E-20)
-
-!
-!--          Third step: calculate bare soil resistance r_soil
-             m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *        &
-                     m_res(j,i)
-
-             f2 = ( m_soil(0,j,i) - m_min ) / ( m_fc(j,i) - m_min )
+             r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
+                                             + 1.0E-20_wp)
+
+!
+!--          Third step: calculate bare soil resistance r_soil. The Clapp & 
+!--          Hornberger parametrization does not consider c_veg.
+             IF ( soil_type /= 7 )  THEN
+                m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
+                        m_res(j,i)
+             ELSE
+                m_min = m_wilt(j,i)
+             ENDIF
+
+             f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
              f2 = MAX(f2,1.0E-20_wp)
+             f2 = MIN(f2,1.0_wp)
 
              r_soil(j,i) = r_soil_min(j,i) / f2
@@ -816 +1019 @@
 !
 !--          Calculate fraction of liquid water reservoir
-             m_liq_max = m_max_depth * LAI(j,i)
-             c_liq(j,i) = MIN(1.0_wp, m_liq(j,i)/m_liq_max)
-
+             m_liq_eb_max = m_max_depth * lai(j,i)
+             c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp))
+             
+
+!
+!--          Calculate saturation specific humidity
              q_s = 0.622_wp * e_s / surface_pressure
 
 !
-!--          In case of dew fall, set resistances to zero.
-!--          To do: what does that physically reasoning behind this?
+!--          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
+                r_canopy(j,i) = 0.0_wp
+                r_soil(j,i)   = 0.0_wp
+             ENDIF 
+
+!
+!--          Calculate coefficients for the total evapotranspiration 
+             f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
+                           (r_a(j,i) + r_canopy(j,i))
+
+             f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
+                                                          r_soil(j,i))
+             f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
+
+
+!
+!--          If soil moisture is below wilting point, plants do no longer
+!--          transpirate.
+!              IF ( m_soil(k,j,i) .LT. m_wilt(j,i) )  THEN
+!                 f_qsws_veg = 0.0_wp
+!              ENDIF
+
+!
+!--          If dewfall is deactivated, vegetation, soil and liquid water
+!--          reservoir are not allowed to take up water from the super-saturated
+!--          air.
              IF ( humidity )  THEN
                 IF ( q_s .LE. q_p(k+1,j,i) )  THEN
-                   r_canopy(j,i) = 0.0_wp
-                   r_soil(j,i) = 0.0_wp
+                   IF ( .NOT. dewfall )  THEN
+                      f_qsws_veg  = 0.0_wp
+                      f_qsws_soil = 0.0_wp
+                      f_qsws_liq  = 0.0_wp
+                   ENDIF
                 ENDIF
              ENDIF 
 
-
-!
-!--          Calculate coefficients for the total evapotranspiration 
-             f_LE_veg  = rho_lv * c_veg(j,i) * (1.0 - c_liq(j,i)) / (r_a(j,i)  &
-                                                + r_canopy(j,i))
-             f_LE_soil = rho_lv * (1.0 - c_veg(j,i)) / (r_a(j,i) + r_soil(j,i))
-             f_LE_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
-
-
-!
-!--          If soil moisture is below wilting point, plants do no longer
-!--          transpirate.
-             IF ( m_soil(k,j,i) .LT. m_wilt(j,i) )  THEN
-                f_LE_veg = 0.0
+             f_shf  = rho_cp / r_a(j,i)
+             f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
+
+!
+!--          Calculate derivative of q_s for Taylor series expansion
+             e_s_dT = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
+                              17.269_wp*(t_surface(j,i) - 273.16_wp)           &
+                              / (t_surface(j,i) - 35.86_wp)**2 )
+
+             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(j,i) = rad_net(j,i) + sigma_SB * t_surface(j,i) ** 4
+
+             IF ( humidity )  THEN
+
+!
+!--             Numerator of the prognostic equation
+                coef_1 = rad_net(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   &
+                         + 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
+
+             ELSE
+
+!
+!--             Numerator of the prognostic equation
+                coef_1 = rad_net(j,i) + 3.0_wp * sigma_SB * t_surface(j,i) ** 4&
+                         + f_shf / exn * T_1 + 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               &
+                         + lambda_surface + f_shf / exn
+
              ENDIF
 
-             f_H  = rho_cp / r_a(j,i)
-             f_LE = f_LE_veg + f_LE_soil + f_LE_liq
-       
-!
-!--          Calculate derivative of q_s for Taylor series expansion
-             e_s_dT = e_s * ( 17.269_wp / (T_0(j,i) - 35.86_wp) -              &
-                              17.269_wp*(T_0(j,i) - 273.16_wp) / (T_0(j,i)     &
-                              - 35.86_wp)**2 )
-
-             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
-             Rn(j,i) = Rn(j,i) + sigma_SB * T_0(j,i) ** 4
-
-             IF ( humidity )  THEN
-
-!
-!--             Numerator of the prognostic equation
-                coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H     &
-                         / exn * T_1 + f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT   &
-                         * T_0(j,i) ) + lambda_skin * T_soil(0,j,i)
-
-!
-!--             Denominator of the prognostic equation
-                coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 + f_LE * dq_s_dT    &
-                         + lambda_skin + f_H / exn
-
-             ELSE
-
-!
-!--             Numerator of the prognostic equation
-                coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H /   &
-                         exn * T_1 + lambda_skin * T_soil(0,j,i)
-
-!
-!--             Denominator of the prognostic equation
-                coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3                     &
-                         + lambda_skin + f_H / exn
-
-             ENDIF
-
              tend = 0.0_wp
 
 !
-!--          Implicit solution when the skin layer has no heat capacity,
+!--          Implicit solution when the surface layer has no heat capacity,
 !--          otherwise use RK3 scheme.
-             T_0_p(j,i) = ( coef_1 * dt_3d * tsc(2) + C_skin * T_0(j,i) ) /    &
-                          ( C_skin + coef_2 * dt_3d * tsc(2) ) 
-
+             t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface *        &
+                                t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d&
+                                * tsc(2) ) 
 !
 !--          Add RK3 term
-             T_0_p(j,i) = T_0_p(j,i) + dt_3d * tsc(3) * tT_soil_m(0,j,i)
+             t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)              &
+                                * tt_surface_m(j,i)
 
 !
 !--          Calculate true tendency
-             tend = (T_0_p(j,i) - T_0(j,i) - tsc(3) * tT_0_m(j,i)) / (dt_3d    &
-                      * tsc(2))
-
-!
-!--          Calculate T_0 tendencies for the next Runge-Kutta step
+             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_0_m(j,i) = tend
+                   tt_surface_m(j,i) = tend
                 ELSEIF ( intermediate_timestep_count <                         &
                          intermediate_timestep_count_max )  THEN
-                   tT_0_m(j,i) = -9.5625_wp * tend + 5.3125_wp * tT_0_m(j,i)
+                   tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp           &
+                                       * tt_surface_m(j,i)
                 ENDIF
              ENDIF
 
-             pt_p(k,j,i) = T_0_p(j,i) / exn
+             pt_p(k,j,i) = t_surface_p(j,i) / exn
 !
 !--          Calculate fluxes
-             Rn(j,i)        = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i)**4        &
-                              - 4.0_wp * sigma_SB * T_0(j,i)**3 * T_0_p(j,i)
-             G(j,i)         = lambda_skin * (T_0_p(j,i) - T_soil(0,j,i))
-             H(j,i)         = - f_H  * ( pt_p(k+1,j,i) - pt_p(k,j,i) )
+             rad_net(j,i)   = rad_net(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)
+             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) )
 
              IF ( humidity )  THEN
-                LE(j,i)        = - f_LE      * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
-                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
-
-                LE_veg(j,i)    = - f_LE_veg  * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
-                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
-                LE_soil(j,i)   = - f_LE_soil * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
-                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
-                LE_liq(j,i)    = - f_LE_liq  * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
-                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
+                qsws_eb(j,i)  = - f_qsws    * ( q_p(k+1,j,i) - 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       &
+                                    + 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       &
+                                    + 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       &
+                                    + dq_s_dT * t_surface(j,i) - dq_s_dT       &
+                                    * t_surface_p(j,i) )
              ENDIF
 
-!              IF ( i == 1 .AND. j == 1 )  THEN
-!                 PRINT*, "Rn", Rn(j,i)
-!                  PRINT*, "H", H(j,i)
-!                 PRINT*, "LE", LE(j,i)
-!                 PRINT*, "LE_liq", LE_liq(j,i)
-!                 PRINT*, "LE_veg", LE_veg(j,i)
-!                 PRINT*, "LE_soil", LE_soil(j,i)
-!                 PRINT*, "G", G(j,i)
-!              ENDIF
-
 !
 !--          Calculate the true surface resistance
-             IF ( LE(j,i) .EQ. 0.0 )  THEN
-                r_s(j,i) = 1.0E10
+             IF ( qsws_eb(j,i) .EQ. 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 * T_0(j,i)&
-                           - dq_s_dT * T_0_p(j,i) ) / LE(j,i) - r_a(j,i)
+                r_s(j,i) = - rho_lv * ( q_p(k+1,j,i) - 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)
              ENDIF
-
-!
-!--          Calculate fluxes in the atmosphere
-             shf(j,i) = H(j,i) / rho_cp
 
 !
@@ -967 +1183 @@
              IF ( humidity )  THEN
 !
-!--             If precipitation is activated, add rain water to LE_liq.
+!--             If precipitation is activated, add rain water to qsws_liq_eb.
 !--             precipitation_rate is given in mm.
                 IF ( precipitation )  THEN
-                   LE_liq(j,i) = LE_liq(j,i) + precipitation_rate(j,i)         &
-                                               * 0.001_wp * rho_l * l_v
+                   qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                         &
+                                       + precipitation_rate(j,i) * 0.001_wp    &
+                                       * rho_l * l_v
                 ENDIF
 !
@@ -977 +1194 @@
                 IF ( q_s .LE. q_p(k+1,j,i))  THEN
 !
-!--                Check if reservoir is full (avoid values > m_liq_max)
-!--                In that case, LE_liq goes to LE_soil. In this case
-!--                LE_veg is zero anyway (because c_liq = 1), so that tend is
-!--                zero and no further check is needed
-                   IF ( m_liq(j,i) .EQ. m_liq_max )  THEN
-                      LE_soil(j,i) = LE_soil(j,i) + LE_liq(j,i)
-                      LE_liq(j,i) = 0.0_wp
+!--                Check if reservoir is full (avoid values > m_liq_eb_max)
+!--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this 
+!--                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
+                      qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
+                                           + qsws_liq_eb(j,i)
+                      qsws_liq_eb(j,i)  = 0.0_wp
                    ENDIF
 
 !
-!--                In case LE_veg becomes negative (unphysical behavior), let
-!--                the water enter the liquid water reservoir as dew on the
+!--                In case qsws_veg_eb becomes negative (unphysical behavior), 
+!--                let the water enter the liquid water reservoir as dew on the
 !--                plant
-                   IF ( LE_veg(j,i) .LT. 0.0_wp )  THEN
-                      LE_liq(j,i) = LE_liq(j,i) + LE_veg(j,i)
-                      LE_veg(j,i) = 0.0_wp
+                   IF ( qsws_veg_eb(j,i) .LT. 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
                    ENDIF
                 ENDIF                    
   
-                tend = - LE_liq(j,i) * drho_l_lv
-
-                m_liq_p(j,i) = m_liq(j,i) + dt_3d * ( tsc(2) * tend            &
-                                                   + tsc(3) * tm_liq_m(j,i) )
+                tend = - qsws_liq_eb(j,i) * drho_l_lv
+
+                m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
+                                                   + tsc(3) * tm_liq_eb_m(j,i) )
 
 !
 !--             Check if reservoir is overfull -> reduce to maximum
 !--             (conservation of water is violated here)
-                m_liq_p(j,i) = MIN(m_liq_p(j,i),m_liq_max)
+                m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
 
 !
 !--             Check if reservoir is empty (avoid values < 0.0)
 !--             (conservation of water is violated here)
-                m_liq_p(j,i) = MAX(m_liq_p(j,i),0.0_wp)
-
-
-!
-!--             Calculate m_liq tendencies for the next Runge-Kutta step
+                m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
+
+
+!
+!--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
                    IF ( intermediate_timestep_count == 1 )  THEN
-                      tm_liq_m(j,i) = tend
+                      tm_liq_eb_m(j,i) = tend
                    ELSEIF ( intermediate_timestep_count <                      &
                             intermediate_timestep_count_max )  THEN
-                      tm_liq_m(j,i) = -9.5625_wp * tend + 5.3125_wp            &
-                                      * tm_liq_m(j,i)
+                      tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
+                                      * tm_liq_eb_m(j,i)
                    ENDIF
                 ENDIF
 
-!
-!--             Calculate moisture flux in the atmosphere
-                qsws(j,i) = LE(j,i) / rho_lv
-
              ENDIF
 
           ENDDO
        ENDDO
-
-
 
     END SUBROUTINE lsm_energy_balance
@@ -1042 +1254 @@
 ! ------------
 !
+! Soil model as part of the land surface model. The model predicts soil
+! temperature and water content.
 !------------------------------------------------------------------------------!
     SUBROUTINE lsm_soil_model
@@ -1052 +1266 @@
        INTEGER(iwp) ::  k   !: running index
 
-       REAL(wp)     :: h_VG !: Van Genuchten coef. h
-
-       REAL(wp), DIMENSION(0:soil_layers-1) :: gamma_temp,  & !: temp. gamma
+       REAL(wp)     :: h_vg !: Van Genuchten coef. h
+
+       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !: temp. gamma
                                                lambda_temp, & !: temp. lambda
                                                tend           !: tendency
@@ -1060 +1274 @@
        DO i = nxlg, nxrg    
           DO j = nysg, nyng
-             DO k = 0, soil_layers-1
+             DO k = nzb_soil, nzt_soil
 !
 !--             Calculate volumetric heat capacity of the soil, taking into 
 !--             account water content 
-                rhoC_total(k,j,i) = (rhoC_soil * (1.0 - m_sat(j,i))            &
-                                     + rhoC_water * m_soil(k,j,i))
+                rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i))       &
+                                     + rho_c_water * m_soil(k,j,i))
 
 !
 !--             Calculate soil heat conductivity at the center of the soil 
 !--             layers
-                Ke = 1.0 + 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) +     &
+                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_h_dry
 
@@ -1079 +1293 @@
 !--          Calculate soil heat conductivity (lambda_h) at the _stag level 
 !--          using linear interpolation
-             DO k = 0, soil_layers-2
+             DO k = nzb_soil, nzt_soil-1
                   
                 lambda_h(k,j,i) = lambda_temp(k) +                             &
@@ -1086 +1300 @@
 
              ENDDO
-             lambda_h(soil_layers-1,j,i) = lambda_temp(soil_layers-1)
-
-!
-!--          Prognostic equation for soil temperature T_soil
+             lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
+
+!
+!--          Prognostic equation for soil temperature t_soil
              tend(:) = 0.0_wp
-             tend(0) = (1.0/rhoC_total(0,j,i)) *                               &
-                       ( lambda_h(0,j,i) * ( T_soil(1,j,i) - T_soil(0,j,i) )   &
-                         * ddz_soil(0) + G(j,i) ) * ddz_soil_stag(0)
-
-             DO  k = 1, soil_layers-1
-                tend(k) = (1.0/rhoC_total(k,j,i))                              &
+             tend(0) = (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)         &
+                         + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
+
+             DO  k = 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) )            &
+                              * ( t_soil(k+1,j,i) - t_soil(k,j,i) )            &
                               * ddz_soil(k)                                    &
                               - lambda_h(k-1,j,i)                              &
-                              * ( T_soil(k,j,i) - T_soil(k-1,j,i) )            &
+                              * ( t_soil(k,j,i) - t_soil(k-1,j,i) )            &
                               * ddz_soil(k-1)                                  &
                             ) * ddz_soil_stag(k)
              ENDDO
 
-             T_soil_p(0:soil_layers-1,j,i) = T_soil(0:soil_layers-1,j,i)       &
+             t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)   &
                                              + dt_3d * ( tsc(2)                &
                                              * tend(:) + tsc(3)                &
-                                             * tT_soil_m(:,j,i) )   
-
-!
-!--          Calculate T_soil tendencies for the next Runge-Kutta step
+                                             * tt_soil_m(:,j,i) )   
+
+!
+!--          Calculate t_soil tendencies for the next Runge-Kutta step
              IF ( timestep_scheme(1:5) == 'runge' )  THEN
                 IF ( intermediate_timestep_count == 1 )  THEN
-                   DO  k = 0, soil_layers-1
-                      tT_soil_m(k,j,i) = tend(k)
+                   DO  k = nzb_soil, nzt_soil
+                      tt_soil_m(k,j,i) = tend(k)
                    ENDDO
                 ELSEIF ( intermediate_timestep_count <                         &
                          intermediate_timestep_count_max )  THEN
-                   DO  k = 0, soil_layers-1
-                      tT_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp      &
-                                         * tT_soil_m(k,j,i)
+                   DO  k = nzb_soil, nzt_soil
+                      tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp      &
+                                         * tt_soil_m(k,j,i)
                    ENDDO
                 ENDIF
@@ -1128 +1343 @@
 
 
-             DO k = 0, soil_layers-1
+             DO k = nzb_soil, nzt_soil
+
 !
 !--             Calculate soil diffusivity at the center of the soil layers
-                lambda_temp(k) = (- b_CH * gamma_w_sat(j,i) * psi_sat          &
+                lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat          &
                                   / m_sat(j,i) ) * ( MAX(m_soil(k,j,i),        &
-                                  m_wilt(j,i)) / m_sat(j,i) )**(b_CH + 2.0_wp)
-
-!
-!--             Calculate the hydraulic conductivity after Van Genuchten (1980)
-                h_VG = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -          &
-                           MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_VG(j,i)      &
-                           / (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 +               &
-                                (alpha_VG(j,i)*h_VG)**n_VG(j,i))**(1.0_wp      &
-                                -1.0_wp/n_VG(j,i)) - (alpha_VG(j,i)*h_VG       &
-                                )**(n_VG(j,i)-1.0_wp))**2 )                    &
-                                / ( (1.0_wp + (alpha_VG(j,i)*h_VG)**n_VG(j,i)  &
-                                )**((1.0_wp - 1.0_wp/n_VG(j,i))*(l_VG(j,i)     &
-                                + 2.0)) )
+                                  m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp)
+
+!
+!--             Parametrization of Van Genuchten
+                IF ( soil_type /= 7 )  THEN
+!
+!--                Calculate the hydraulic conductivity after Van Genuchten 
+!--                (1980)
+                   h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -       &
+                              MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i)   &
+                              / (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 +            &
+                                   (alpha_vg(j,i)*h_vg)**n_vg(j,i))**(1.0_wp   &
+                                   -1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg    &
+                                   )**(n_vg(j,i)-1.0_wp))**2 )                 &
+                                   / ( (1.0_wp + (alpha_vg(j,i)*h_vg           &
+                                   )**n_vg(j,i))**((1.0_wp - 1.0_wp/n_vg(j,i)) &
+                                   *(l_vg(j,i) + 2.0_wp)) )
+
+!
+!--             Parametrization of Clapp & Hornberger
+                ELSE
+                   gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i)           &
+                                   / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
+                ENDIF
 
              ENDDO
@@ -1156 +1383 @@
 !
 !--             Calculate soil diffusivity (lambda_w) at the _stag level 
-!--             using linear interpolation
-                DO k = 0, soil_layers-2
+!--             using linear interpolation. To do: replace this with
+!--             ECMWF-IFS Eq. 8.81
+                DO k = nzb_soil, nzt_soil-1
                      
                    lambda_w(k,j,i) = lambda_temp(k) +                          &
                                      ( lambda_temp(k+1) - lambda_temp(k) )     &
-                                     * 0.5 * dz_soil_stag(k) * ddz_soil(k+1)
+                                     * 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 * dz_soil_stag(k) * ddz_soil(k+1)                 
+                                     * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1)                 
 
                 ENDDO
@@ -1174 +1402 @@
 !--             in the bottom layer.
                 IF ( conserve_water_content )  THEN
-                   gamma_w(soil_layers-1,j,i) = 0.0_wp
+                   gamma_w(nzt_soil,j,i) = 0.0_wp
                 ELSE
-                   gamma_w(soil_layers-1,j,i) = lambda_temp(soil_layers-1)
+                   gamma_w(nzt_soil,j,i) = lambda_temp(nzt_soil)
                 ENDIF     
 
-!--             The root extraction (= root_extr * LE_veg / (rho_l * l_v))
+!--             The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v))
 !--             ensures the mass conservation for water. The transpiration of
 !--             plants equals the cumulative withdrawals by the roots in the 
 !--             soil. The scheme takes into account the availability of water
 !--             in the soil layers as well as the root fraction in the
-!--             respective layer
-
-!
-!--             Calculate the root extraction (ECMWF 7.69, with some
-!--             modifications)
+!--             respective layer. Layer with moisture below wilting point will
+!--             not contribute, which reflects the preference of plants to
+!--             take water from moister layers.
+
+!
+!--             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 = 0, soil_layers-1
-                    m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) *       &
-                              dz_soil_stag(k) 
-
+                DO k = nzb_soil, nzt_soil
+                    IF ( m_soil(k,j,i) .GT. m_wilt(j,i) )  THEN
+                       m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
+                    ENDIF
                 ENDDO  
 
-!
-!--             For conservation of mass, the sum of root_extr must be 1
-                DO k = 0, soil_layers-1 
-                   root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i)               &
-                                  * dz_soil_stag(k) / m_total
+                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
-
 
 !
 !--             Prognostic equation for soil water content m_soil
                 tend(:) = 0.0_wp
-                tend(0) = ( lambda_w(0,j,i) * ( m_soil(1,j,i) - m_soil(0,j,i) )&
-                            * ddz_soil(0) - gamma_w(0,j,i) - ( root_extr(0)    &
-                            * LE_veg(j,i) + LE_soil(j,i) ) * drho_l_lv         &
-                          ) * ddz_soil_stag(0)
-
-                DO  k = 1, soil_layers-2
+                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) - (   &
+                            root_extr(nzb_soil) * qsws_veg_eb(j,i)             &
+                            + qsws_soil_eb(j,i) ) * drho_l_lv )                &
+                            * ddz_soil_stag(nzb_soil)
+
+                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) * LE_veg(j,i)&
-                               * drho_l_lv)                                    &
+                               + gamma_w(k-1,j,i) - (root_extr(k)              &
+                               * qsws_veg_eb(j,i) * drho_l_lv)                 &
                              ) * ddz_soil_stag(k)
 
                 ENDDO
-                tend(soil_layers-1) = ( - gamma_w(soil_layers-1,j,i)           &
-                                        - lambda_w(soil_layers-2,j,i)          &
-                                        * (m_soil(soil_layers-1,j,i)           &
-                                        - m_soil(soil_layers-2,j,i))           &
-                                        * ddz_soil(soil_layers-2)              &
-                                        + gamma_w(soil_layers-2,j,i) - (       &
-                                          root_extr(soil_layers-1)             &
-                                        * LE_veg(j,i) * drho_l_lv      )       &
-                                      ) * ddz_soil_stag(soil_layers-1)             
-
-                m_soil_p(0:soil_layers-1,j,i) = m_soil(0:soil_layers-1,j,i)    &
+                tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                     &
+                                        - lambda_w(nzt_soil-1,j,i)             &
+                                        * (m_soil(nzt_soil,j,i)                &
+                                        - m_soil(nzt_soil-1,j,i))              &
+                                        * ddz_soil(nzt_soil-1)                 &
+                                        + gamma_w(nzt_soil-1,j,i) - (          &
+                                          root_extr(nzt_soil)                  &
+                                        * qsws_veg_eb(j,i) * drho_l_lv  )      &
+                                      ) * ddz_soil_stag(nzt_soil)             
+
+                m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
                                                 + dt_3d * ( tsc(2) * tend(:)   &
                                                 + tsc(3) * tm_soil_m(:,j,i) )   
@@ -1244 +1479 @@
                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
                    IF ( intermediate_timestep_count == 1 )  THEN
-                      DO  k = 0, soil_layers-1
+                      DO  k = nzb_soil, nzt_soil
                          tm_soil_m(k,j,i) = tend(k)
                       ENDDO
                    ELSEIF ( intermediate_timestep_count <                      &
                             intermediate_timestep_count_max )  THEN
-                      DO  k = 0, soil_layers-1
+                      DO  k = nzb_soil, nzt_soil
                          tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
                                      * tm_soil_m(k,j,i)
@@ -1264 +1499 @@
 !--    Calculate surface specific humidity
        IF ( humidity )  THEN
-          CALL calc_q0
+          CALL calc_q_surface
        ENDIF
 
@@ -1274 +1509 @@
 ! Description:
 ! ------------
-!
+! Calculation of specific humidity of the surface layer (surface)
 !------------------------------------------------------------------------------!
-    SUBROUTINE calc_q0
+    SUBROUTINE calc_q_surface
 
        IMPLICIT NONE
@@ -1288 +1523 @@
           DO j = nysg, nyng
              k = nzb_s_inner(j,i)
-!
-!--          Temporary solution as long as T_0 is prescribed
-
-             pt_p(k,j,i) = T_0(j,i) / exn
+
 !
 !--          Calculate water vapour pressure at saturation
-             e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( T_0(j,i) -         &
-                                              273.16_wp ) /  ( T_0(j,i) -      &
-                                              35.86_wp ) )
+             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 ) )
 
 !
@@ -1313 +1545 @@
        ENDDO
 
-    END SUBROUTINE calc_q0
+    END SUBROUTINE calc_q_surface
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+! Swapping of timelevels
+!------------------------------------------------------------------------------!
+    SUBROUTINE lsm_swap_timelevel ( mod_count )
+
+       IMPLICIT NONE
+
+       INTEGER, INTENT(IN) :: mod_count
+
+#if defined( __nopointer )
+
+       t_surface    = t_surface_p
+       t_soil       = t_soil_p
+       IF ( humidity )  THEN
+          m_soil    = m_soil_p
+          m_liq_eb  = m_liq_eb_p
+       ENDIF
+
+#else
+    
+       SELECT CASE ( mod_count )
+
+          CASE ( 0 )
+
+             t_surface = t_surface_p
+             t_soil    = t_soil_p
+             IF ( humidity )  THEN
+                m_soil    = m_soil_p
+                m_liq_eb  = m_liq_eb_p
+             ENDIF
+
+
+          CASE ( 1 )
+
+             t_surface  => t_surface_1; t_surface_p  => t_surface_2
+             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
+             IF ( humidity )  THEN
+                m_soil    => m_soil_1;   m_soil_p    => m_soil_2
+                m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
+             ENDIF
+
+       END SELECT
+#endif
+
+    END SUBROUTINE lsm_swap_timelevel
 
 

palm/trunk/SOURCE/modules.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
+! Increased pr_palm to 120. Increased length of dots_unit and dots_label to 13
+! digits. Increased length of domask, do2d, and do3d to 20 digits.
 ! 
 ! Former revisions:
@@ -568 +569 @@
     CHARACTER (LEN=20), DIMENSION(11)  ::  netcdf_precision = ' '
 
-    CHARACTER (LEN=10), DIMENSION(max_masks,0:1,100) ::  domask = ' '
-    CHARACTER (LEN=10), DIMENSION(0:1,100) ::  do2d = ' ', do3d = ' '
+    CHARACTER (LEN=20), DIMENSION(max_masks,0:1,100) ::  domask = ' '
+    CHARACTER (LEN=20), DIMENSION(0:1,100) ::  do2d = ' ', do3d = ' '
 
     INTEGER(iwp) ::  abort_mode = 1, average_count_pr = 0, average_count_sp = 0, &
@@ -1108 +1109 @@
     INTEGER(iwp) ::  dots_num = 23
 
-    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_names =   &
-         (/ 'wpt0  ', 'ws2   ', 'tsw2  ', 'ws3   ', 'ws2tsw', 'wstsw2', &
+    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_names =          &
+         (/ 'wpt0  ', 'ws2   ', 'tsw2  ', 'ws3   ', 'ws2tsw', 'wstsw2',        &
             'z_i   ' /)
 
-    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames = &
-         (/ 'wpt0  ', 'w*2   ', 't*w2  ', 'w*3   ', 'w*2t*w', 'w*t*w2',   &
+    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames =      &
+         (/ 'wpt0  ', 'w*2   ', 't*w2  ', 'w*3   ', 'w*2t*w', 'w*t*w2',        &
             'z_i   ' /)
 
-    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label = &
+    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label =                   &
           (/ 'tnpt   ', 'x_     ', 'y_     ', 'z_     ', 'z_abs  ', 'u      ', &
              'v      ', 'w      ', 'u"     ', 'v"     ', 'w"     ', 'npt_up ', &
@@ -1123 +1124 @@
              'w*2    ', 'u"2    ', 'v"2    ', 'w"2    ', 'npt*2  ' /)
 
-    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit = &
+    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit =                    &
           (/ 'number ', 'm      ', 'm      ', 'm      ', 'm      ', 'm/s    ', &
              'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'number ', &
@@ -1130 +1131 @@
              'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)
 
-    CHARACTER (LEN=7), DIMENSION(dots_max) :: dots_label = &
-          (/ 'E      ', 'E*     ', 'dt     ', 'u*     ', 'th*    ', 'umax   ', &
-             'vmax   ', 'wmax   ', 'div_new', 'div_old', 'z_i_wpt', 'z_i_pt ', &
-             'w*     ', 'w"pt"0 ', 'w"pt"  ', 'wpt    ', 'pt(0)  ', 'pt(zp) ', &
-             'w"u"0  ', 'w"v"0  ', 'w"q"0  ', 'mo_L   ', 'q*     ',            &
+    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
+          (/ 'E            ', 'E*           ', 'dt           ',                &
+             'u*           ', 'th*          ', 'umax         ',                &
+             'vmax         ', 'wmax         ', 'div_new      ',                &
+             'div_old      ', 'z_i_wpt      ', 'z_i_pt       ',                &
+             'w*           ', 'w"pt"0       ', 'w"pt"        ',                &
+             'wpt          ', 'pt(0)        ', 'pt(zp)       ',                &
+             'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
+             'mo_L         ', 'q*           ',                                 &
              ( 'unknown', i9 = 1, dots_max-23 ) /)
 
-    CHARACTER (LEN=7), DIMENSION(dots_max) :: dots_unit = &
-          (/ 'm2/s2  ', 'm2/s2  ', 's      ', 'm/s    ', 'K      ', 'm/s    ', &
-             'm/s    ', 'm/s    ', 's-1    ', 's-1    ', 'm      ', 'm      ', &
-             'm/s    ', 'K m/s  ', 'K m/s  ', 'K m/s  ', 'K      ', 'K      ', &
-             'm2/s2  ', 'm2/s2  ', 'kg m/s ', 'm      ', 'kg/kg  ',            &
+    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit =                     &
+          (/ 'm2/s2        ', 'm2/s2        ', 's            ',                &
+             'm/s          ', 'K            ', 'm/s          ',                &
+             'm/s          ', 'm/s          ', 's-1          ',                &
+             's-1          ', 'm            ', 'm            ',                &
+             'm/s          ', 'K m/s        ', 'K m/s        ',                &
+             'K m/s        ', 'K            ', 'K            ',                &
+             'm2/s2        ', 'm2/s2        ', 'kg m/s       ',                &
+             'm            ', 'kg/kg        ',                                 &
              ( 'unknown', i9 = 1, dots_max-23 ) /)
 
@@ -1422 +1431 @@
 
     CHARACTER (LEN=40) ::  region(0:9)
-    INTEGER(iwp) ::  pr_palm = 100, statistic_regions = 0
+    INTEGER(iwp) ::  pr_palm = 120, 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

@@ -23 +23 @@
 ! Current revisions:
 ! ------------------
-! 
+! Added support for land surface model and radiation model output. In the course
+! of this action a new vertical grid zs (soil) was introduced.
 ! 
 ! Former revisions:
@@ -97 +98 @@
 ! In case of extend = .TRUE.:
 ! Find out if dimensions and variables of an existing file match the values
-! of the actual run. If so, get all necessary informations (ids, etc.) from
+! of the actual run. If so, get all necessary information (ids, etc.) from
 ! this file.
 !
@@ -130 +131 @@
 
     USE kinds
+
+    USE land_surface_model_mod,                                                &
+        ONLY: land_surface, nzb_soil, nzt_soil, id_dim_zs_xy, id_dim_zs_xz,    &
+              id_dim_zs_yz, id_dim_zs_3d, id_dim_zs_mask, id_var_zs_xy,        &
+              id_var_zs_xz, id_var_zs_yz ,id_var_zs_3d, id_var_zs_mask,        &
+              nzs, zs
 
     USE pegrid
@@ -181 +188 @@
     INTEGER(iwp) ::  kk                                      !:
     INTEGER(iwp) ::  ns                                      !:
+    INTEGER(iwp) ::  ns_do                                   !: actual value of ns for soil model data
     INTEGER(iwp) ::  ns_old                                  !:
     INTEGER(iwp) ::  ntime_count                             !:
@@ -439 +447 @@
 !
 !--       In case of non-flat topography define 2d-arrays containing the height
-!--       informations
+!--       information
           IF ( TRIM( topography ) /= 'flat' )  THEN
 !
@@ -478 +486 @@
 
           ENDIF             
-
+ 
+          IF ( land_surface )  THEN
+!
+!--          Define vertical coordinate grid (zw grid)
+             nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'zs_3d', &
+                                     mask_size(mid,3), id_dim_zs_mask(mid,av) )
+             CALL handle_netcdf_error( 'netcdf', 536 )
+
+             nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zs_3d', NF90_DOUBLE, &
+                                     id_dim_zs_mask(mid,av), &
+                                     id_var_zs_mask(mid,av) )
+             CALL handle_netcdf_error( 'netcdf', 536 )
+
+             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
+                                  'units', 'meters' )
+             CALL handle_netcdf_error( 'netcdf', 537 )
+
+          ENDIF
 
 !
@@ -521 +546 @@
                    grid_y = 'y'
                    grid_z = 'zw'
+!
+!--             soil grid
+                CASE ( 't_soil', 'm_soil' )
+
+                   grid_x = 'x'
+                   grid_y = 'y'
+                   grid_z = 'zs'
 
                 CASE DEFAULT
@@ -548 +580 @@
              ELSEIF ( grid_z == 'zw' )  THEN
                 id_z = id_dim_zw_mask(mid,av)
+             ELSEIF ( grid_z == "zs" )  THEN
+                id_z = id_dim_zs_mask(mid,av)
              ENDIF
 
@@ -692 +726 @@
 
           ENDIF
+
+          IF ( land_surface )  THEN
+!
+!--          Write zs data (vertical axes for soil model), use negative values
+!--          to indicate soil depth
+             ALLOCATE( netcdf_data(mask_size(mid,3)) )
+
+             netcdf_data = zs( mask_k_global(mid,:mask_size(mid,3)) )
+
+             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
+                                     netcdf_data, start = (/ 1 /), &
+                                     count = (/ mask_size(mid,3) /) )
+             CALL handle_netcdf_error( 'netcdf', 538 )
+
+             DEALLOCATE( netcdf_data )
+
+          ENDIF
+
 !
 !--       restore original parameter file_id (=formal parameter av) into av
@@ -982 +1034 @@
 !
 !--       In case of non-flat topography define 2d-arrays containing the height
-!--       informations
+!--       information
           IF ( TRIM( topography ) /= 'flat' )  THEN
 !
@@ -1016 +1068 @@
           ENDIF             
 
+          IF ( land_surface )  THEN
+!
+!--          Define vertical coordinate grid (zs grid)
+             nc_stat = NF90_DEF_DIM( id_set_3d(av), 'zs_3d', nzt_soil-nzb_soil+1, &
+                                     id_dim_zs_3d(av) )
+             CALL handle_netcdf_error( 'netcdf', 70 )
+
+             nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zs_3d', NF90_DOUBLE, &
+                                     id_dim_zs_3d(av), id_var_zs_3d(av) )
+             CALL handle_netcdf_error( 'netcdf', 71 )
+
+             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zs_3d(av), 'units', &
+                                     'meters' )
+             CALL handle_netcdf_error( 'netcdf', 72 )
+
+          ENDIF
 
 !
@@ -1058 +1126 @@
                    grid_y = 'y'
                    grid_z = 'zw'
+!
+!--             soil grid
+                CASE ( 't_soil', 'm_soil' )
+
+                   grid_x = 'x'
+                   grid_y = 'y'
+                   grid_z = 'zs'
 
                 CASE DEFAULT
@@ -1085 +1160 @@
              ELSEIF ( grid_z == 'zw' )  THEN
                 id_z = id_dim_zw_3d(av)
+             ELSEIF ( grid_z == 'zs' )  THEN
+                id_z = id_dim_zs_3d(av)
              ENDIF
 
@@ -1239 +1316 @@
 
              ENDIF
+
+             IF ( land_surface )  THEN
+!
+!--             Write zs grid
+                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zs_3d(av),  &
+                                        - zs(nzb_soil:nzt_soil), start = (/ 1 /), &
+                                        count = (/ nzt_soil-nzb_soil+1 /) )
+                CALL handle_netcdf_error( 'netcdf', 86 )
+             ENDIF
+
           ENDIF
 
@@ -1496 +1583 @@
           CALL handle_netcdf_error( 'netcdf', 107 )
 
+
+          IF ( land_surface )  THEN
+
+             ns_do = 0
+             DO WHILE ( section(ns_do+1,1) < nzs )
+                ns_do = ns_do + 1
+             ENDDO
+!
+!--          Define vertical coordinate grid (zs grid)
+             nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zs_xy', ns_do, id_dim_zs_xy(av) )
+             CALL handle_netcdf_error( 'netcdf', 539 )
+
+             nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zs_xy', NF90_DOUBLE, &
+                                     id_dim_zs_xy(av), id_var_zs_xy(av) )
+             CALL handle_netcdf_error( 'netcdf', 540 )
+
+             nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zs_xy(av), 'units', &
+                                  'meters' )
+             CALL handle_netcdf_error( 'netcdf', 541 )
+
+          ENDIF
+
 !
 !--       Define a pseudo vertical coordinate grid for the surface variables
@@ -1577 +1686 @@
 !
 !--       In case of non-flat topography define 2d-arrays containing the height
-!--       informations
+!--       information
           IF ( TRIM( topography ) /= 'flat' )  THEN
 !
@@ -1611 +1720 @@
           ENDIF
 
-
 !
 !--       Define the variables
@@ -1671 +1779 @@
                          grid_y = 'y'
                          grid_z = 'zw'
+!
+!--                   soil grid
+                      CASE ( 't_soil_xy', 'm_soil_xy' )
+                         grid_x = 'x'
+                         grid_y = 'y'
+                         grid_z = 'zs'
 
                       CASE DEFAULT
@@ -1698 +1812 @@
                    ELSEIF ( grid_z == 'zw' )  THEN
                       id_z = id_dim_zw_xy(av)
+                   ELSEIF ( grid_z == 'zs' )  THEN
+                      id_z = id_dim_zs_xy(av)
                    ENDIF
 
@@ -1813 +1929 @@
 
 !
+!--             Write zs data
+             IF ( land_surface )  THEN
+                ns_do = 0
+                DO  i = 1, ns
+                   IF( section(i,1) == -1 )  THEN
+                      netcdf_data(i) = 1.0_wp  ! section averaged along z
+                      ns_do = ns_do + 1
+                   ELSEIF ( section(i,1) < nzs )  THEN
+                      netcdf_data(i) = - zs( section(i,1) )
+                      ns_do = ns_do + 1
+                   ENDIF
+                ENDDO
+
+                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zs_xy(av), &
+                                        netcdf_data(1:ns_do), start = (/ 1 /),    &
+                                        count = (/ ns_do /) )
+                CALL handle_netcdf_error( 'netcdf', 124 )
+
+             ENDIF
+
+!
 !--          Write gridpoint number data
              netcdf_data(1:ns) = section(1:ns,1)
@@ -1894 +2031 @@
 
              ENDIF
+
+
+
           ENDIF
 
@@ -2240 +2380 @@
 
 !
-!--       Define the two z-axes (zu and zw)
+!--       Define the three z-axes (zu, zw, and zs)
           nc_stat = NF90_DEF_DIM( id_set_xz(av), 'zu', nz+2, id_dim_zu_xz(av) )
           CALL handle_netcdf_error( 'netcdf', 153 )
@@ -2263 +2403 @@
           CALL handle_netcdf_error( 'netcdf', 158 )
 
+          IF ( land_surface )  THEN
+
+             nc_stat = NF90_DEF_DIM( id_set_xz(av), 'zs', nzs, id_dim_zs_xz(av) )
+             CALL handle_netcdf_error( 'netcdf', 542 )
+
+             nc_stat = NF90_DEF_VAR( id_set_xz(av), 'zs', NF90_DOUBLE, &
+                                     id_dim_zs_xz(av), id_var_zs_xz(av) )
+             CALL handle_netcdf_error( 'netcdf', 543 )
+
+             nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_zs_xz(av), 'units', &
+                                     'meters' )
+             CALL handle_netcdf_error( 'netcdf', 544 )
+
+          ENDIF
 
 !
@@ -2308 +2462 @@
                       grid_y = 'y'
                       grid_z = 'zw'
+
+!
+!--                soil grid
+                   CASE ( 't_soil_xz', 'm_soil_xz' )
+
+                      grid_x = 'x'
+                      grid_y = 'y'
+                      grid_z = 'zs'
 
                    CASE DEFAULT
@@ -2335 +2497 @@
                 ELSEIF ( grid_z == 'zw' )  THEN
                    id_z = id_dim_zw_xz(av)
+                ELSEIF ( grid_z == 'zs' )  THEN
+                   id_z = id_dim_zs_xz(av)
                 ENDIF
 
@@ -2514 +2678 @@
                                      count = (/ nz+2 /) )
              CALL handle_netcdf_error( 'netcdf', 167 )
+
+!
+!--          Write zs data
+             IF ( land_surface )  THEN
+                netcdf_data(0:nzs-1) = - zs(nzb_soil:nzt_soil)
+                nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_zs_xz(av), &
+                                        netcdf_data(0:nzs), start = (/ 1 /),    &
+                                        count = (/ nzt_soil-nzb_soil+1 /) )
+               CALL handle_netcdf_error( 'netcdf', 548 )
+             ENDIF
+
 
              DEALLOCATE( netcdf_data )
@@ -2903 +3078 @@
           CALL handle_netcdf_error( 'netcdf', 197 )
 
+          IF ( land_surface )  THEN
+
+             nc_stat = NF90_DEF_DIM( id_set_yz(av), 'zs', nzs, id_dim_zs_yz(av) )
+             CALL handle_netcdf_error( 'netcdf', 545 )
+
+             nc_stat = NF90_DEF_VAR( id_set_yz(av), 'zs', NF90_DOUBLE, &
+                                     id_dim_zs_yz(av), id_var_zs_yz(av) )
+             CALL handle_netcdf_error( 'netcdf', 546 )
+
+             nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_zs_yz(av), 'units', &
+                                     'meters' )
+             CALL handle_netcdf_error( 'netcdf', 547 )
+
+          ENDIF
+
 
 !
@@ -2948 +3138 @@
                       grid_y = 'y'
                       grid_z = 'zw'
+!
+!--                soil grid
+                   CASE ( 't_soil_yz', 'm_soil_yz' )
+
+                      grid_x = 'x'
+                      grid_y = 'y'
+                      grid_z = 'zs'
 
                    CASE DEFAULT
@@ -2975 +3172 @@
                 ELSEIF ( grid_z == 'zw' )  THEN
                    id_z = id_dim_zw_yz(av)
+                ELSEIF ( grid_z == 'zs' )  THEN
+                   id_z = id_dim_zs_yz(av)
                 ENDIF
 

palm/trunk/SOURCE/package_parin.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Several changes in the liste for land surface model and radiation model
 ! 
 ! Former revisions:
@@ -100 +100 @@
 
     USE land_surface_model_mod,                                                &
-        ONLY: alpha_VanGenuchten, C_skin, canopy_resistance_coefficient,       &
-              conserve_water_content, f_shortwave_incoming,                    &
-              hydraulic_conductivity, lambda_skin_stable, lambda_skin_unstable,&
-              leaf_area_index, land_surface, l_VanGenuchten,                   &
-              min_canopy_resistance, field_capacity, residual_moisture,        &
-              saturation_moisture, wilting_point, n_VanGenuchten,              &
-              root_fraction, soil_level, soil_moisture, soil_temperature,      &
-              soil_type, vegetation_coverage, veg_type 
+        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
 
     USE kinds
@@ -134 +135 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY: albedo, day_init, dt_radiation, lambda, radiation,               &
-              time_radiation, time_utc_init 
+        ONLY: albedo, day_init, dt_radiation, lambda, net_radiation, radiation,&
+              radiation_scheme, time_radiation, time_utc_init 
                 
 
@@ -171 +172 @@
                                   vc_size_y, vc_size_z
 
-    NAMELIST /lsm_par/            alpha_VanGenuchten, C_skin,                  &
+    NAMELIST /lsm_par/            alpha_vangenuchten, c_surface,               &
                                   canopy_resistance_coefficient,               &
-                                  conserve_water_content, f_shortwave_incoming,& 
-                                  hydraulic_conductivity, lambda_skin_stable,  &
-                                  lambda_skin_unstable, leaf_area_index,       &
-                                  l_VanGenuchten, min_canopy_resistance,       &
-                                  field_capacity, residual_moisture,           &
-                                  saturation_moisture, wilting_point,          &
-                                  n_VanGenuchten, root_fraction, soil_level,   &
+                                  conserve_water_content, dewfall,             &
+                                  f_shortwave_incoming,                        & 
+                                  hydraulic_conductivity,                      &
+                                  lambda_surface_stable,                       &
+                                  lambda_surface_unstable, leaf_area_index,    &
+                                  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 
+                                  vegetation_coverage, veg_type, zs, z0_eb,    &
+                                  z0h_eb
 
 
@@ -207 +211 @@
                                   write_particle_statistics
 
-    NAMELIST /radiation_par/      lambda, albedo, day_init, dt_radiation,      &
-                                  time_utc_init 
+    NAMELIST /radiation_par/      albedo, day_init, dt_radiation, lambda,      &
+                                  net_radiation, radiation_scheme, time_utc_init 
 
     NAMELIST /spectra_par/        averaging_interval_sp, comp_spectra_level,   &

palm/trunk/SOURCE/prandtl_fluxes.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Removed land surface model part. The surface fluxes are now always calculated
+! within prandtl_fluxes, based on the given surface temperature/humidity (which 
+! is either provided by the land surface model, by large scale forcing data, or
+! directly prescribed by the user.
 ! 
 ! Former revisions:
@@ -91 +94 @@
 
     USE kinds
-
-    USE land_surface_model_mod,                                                &
-        ONLY: land_surface
 
     IMPLICIT NONE
@@ -482 +482 @@
 !
 !-- Compute the vertical kinematic heat flux
-    IF ( .NOT. constant_heatflux .AND. .NOT. land_surface )  THEN
+    IF ( .NOT. constant_heatflux  )  THEN
        !$OMP PARALLEL DO
        !$acc kernels loop independent
@@ -495 +495 @@
 !
 !-- Compute the vertical water/scalar flux
-    IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar )       &
-         .AND. .NOT. land_surface )  THEN
+    IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar ) )  THEN
        !$OMP PARALLEL DO
        !$acc kernels loop independent

palm/trunk/SOURCE/radiation_model.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added support for data output. Various variables have been renamed. Added
+! interface for different radiation schemes (currently: clear-sky, constant, and
+! RRTM (not yet implemented).
 ! 
 ! Former revisions:
@@ -46 +48 @@
     USE kinds
 
+    USE netcdf_control,                                                        &
+        ONLY:  dots_label, dots_num, dots_unit
+
 
     IMPLICIT NONE
 
+    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtm'
+
     INTEGER(iwp) :: i, j, k
 
 
-    INTEGER(iwp) :: day_init = 1 !: day of the year at model start
+    INTEGER(iwp) :: day_init     = 172,  & !: day of the year at model start (21/06)
+                    dots_rad     = 0,    & !: starting index for timeseries output 
+                    irad_scheme  = 0
 
     LOGICAL :: radiation = .FALSE.  !: flag parameter indicating wheather the radiation model is used
@@ -61 +70 @@
  
     REAL(wp) :: albedo = 0.2_wp,             & !: NAMELIST alpha
-                dt_radiation = 9999999.9_wp, &
+                dt_radiation = 0.0_wp,       & !: radiation model timestep
                 exn,                         & !: Exner function
                 lon = 0.0_wp,                & !: longitude in radians
@@ -69 +78 @@
                 decl_3,                      & !: declination coef. 3
                 time_utc,                    & !: current time in UTC
-                time_utc_init = 0.0_wp,      & !: UTC time at model start
+                time_utc_init = 43200.0_wp,  & !: UTC time at model start (noon)
                 day,                         & !: current day of the year
                 lambda = 0.0_wp,             & !: longitude in degrees
                 declination,                 & !: solar declination angle
+                net_radiation = 0.0_wp,      & !: net radiation at surface
                 hour_angle,                  & !: solar hour angle
                 time_radiation = 0.0_wp,     &
@@ -80 +90 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
                 alpha,                       & !: surface albedo
-                Rn,                          & !: net radiation at the surface
-                LW_in,                       & !: incoming longwave radiation
-                LW_out,                      & !: outgoing longwave radiation
-                SW_in,                       & !: incoming shortwave radiation
-                SW_out                         !: outgoing shortwave radiation
+                rad_net,                     & !: net radiation at the surface
+                rad_net_av,                  & !: average of rad_net
+                rad_lw_in,                   & !: incoming longwave radiation
+                rad_lw_out,                  & !: outgoing longwave radiation
+                rad_sw_in,                   & !: incoming shortwave radiation
+                rad_sw_in_av,                & !: average of rad_sw_in
+                rad_sw_out                     !: outgoing shortwave radiation
 
 
@@ -91 +103 @@
     END INTERFACE init_radiation
 
-    INTERFACE lsm_radiation
-       MODULE PROCEDURE lsm_radiation
-    END INTERFACE lsm_radiation
+    INTERFACE radiation_clearsky
+       MODULE PROCEDURE radiation_clearsky
+    END INTERFACE radiation_clearsky
+
+    INTERFACE radiation_constant
+       MODULE PROCEDURE radiation_constant
+    END INTERFACE radiation_constant
+
+    INTERFACE radiation_rrtm
+       MODULE PROCEDURE radiation_rrtm
+    END INTERFACE radiation_rrtm
+
 
     SAVE
@@ -99 +120 @@
     PRIVATE
 
-    PUBLIC albedo, day_init, dt_radiation, init_radiation, lambda,             &
-           lsm_radiation, Rn, radiation, SW_in, sigma_SB, time_radiation,      &
-           time_utc_init 
+    PUBLIC albedo, day_init, dots_rad, dt_radiation, init_radiation,           &
+           irad_scheme, lambda, net_radiation, rad_net, rad_net_av, radiation, &
+           radiation_clearsky, radiation_constant, radiation_rrtm,             &
+           radiation_scheme, rad_sw_in, rad_sw_in_av, sigma_SB,                &
+           time_radiation, time_utc_init 
 
 
@@ -118 +141 @@
 
        ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( Rn(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( LW_in(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( LW_out(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( SW_in(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE ( SW_out(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( rad_lw_in(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( rad_lw_out(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( rad_sw_in(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE ( rad_sw_out(nysg:nyng,nxlg:nxrg) )
+
+       rad_sw_in  = 0.0_wp
+       rad_sw_out = 0.0_wp
+       rad_lw_in  = 0.0_wp
+       rad_lw_out = 0.0_wp
+       rad_net    = 0.0_wp
 
        alpha = albedo
 
 !
+!--    Fix net radiation in case of radiation_scheme = 'constant'
+       IF ( irad_scheme == 0 )  THEN
+          rad_net = net_radiation
+!
 !--    Calculate radiation scheme constants
-       decl_1 = SIN(23.45_wp * pi / 180.0_wp)
-       decl_2 = 2.0 * pi / 365.0_wp
-       decl_3 = decl_2 * 81.0_wp
-
-!
-!--    Calculate latitude and longitude angles (lat and lon, respectively)
-       lat = phi * pi / 180.0_wp
-       lon = lambda * pi / 180.0_wp
+       ELSEIF ( irad_scheme == 1 .OR. irad_scheme == 2 )  THEN
+          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
+          decl_2 = 2.0_wp * pi / 365.0_wp
+          decl_3 = decl_2 * 81.0_wp
+!
+!--       Calculate latitude and longitude angles (lat and lon, respectively)
+          lat = phi * pi / 180.0_wp
+          lon = lambda * pi / 180.0_wp
+       ENDIF
+!
+!--    Add timeseries for radiation model
+       dots_label(dots_num+1) = "rad_net"
+       dots_label(dots_num+2) = "rad_sw_in"
+       dots_unit(dots_num+1:dots_num+2) = "W/m2"
+
+       dots_rad  = dots_num + 1
+       dots_num  = dots_num + 2
 
        RETURN
@@ -147 +189 @@
 !-- A simple clear sky radiation model
 !------------------------------------------------------------------------------!
-    SUBROUTINE lsm_radiation
+    SUBROUTINE radiation_clearsky
     
 
@@ -167 +209 @@
 !
 !--    Calculate zenith angle
-       zenith = SIN(lat)*SIN(declination) + COS(lat) * COS(declination)        &
+       zenith = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)       &
                                             * COS(hour_angle)
        zenith = MAX(0.0_wp,zenith)
@@ -180 +222 @@
 
 !
-!--    Calculate radiation fluxes and net radiation (Rn) for each grid point
+!--    Calculate radiation fluxes and net radiation (rad_net) for each grid 
+!--    point
        DO i = nxlg, nxrg
           DO j = nysg, nyng
 
              k = nzb_s_inner(j,i)
-             SW_in(j,i)  = SW_0 * sky_trans * zenith
-             SW_out(j,i) = - alpha(j,i) * SW_in(j,i)
-             LW_out(j,i) = - sigma_SB * (pt(k,j,i) * exn)**4
-             LW_in(j,i)  = 0.8 * sigma_SB * (pt(k+1,j,i) * exn)**4
-             Rn(j,i)     = SW_in(j,i) + SW_out(j,i) + LW_in(j,i) + LW_out(j,i)
+             rad_sw_in(j,i)  = SW_0 * sky_trans * zenith
+             rad_sw_out(j,i) = - alpha(j,i) * rad_sw_in(j,i)
+             rad_lw_out(j,i) = - sigma_SB * (pt(k,j,i) * exn)**4
+             rad_lw_in(j,i)  = 0.8_wp * sigma_SB * (pt(k+1,j,i) * exn)**4
+             rad_net(j,i)    = rad_sw_in(j,i) + rad_sw_out(j,i)                &
+                                + rad_lw_in(j,i) + rad_lw_out(j,i)
 
           ENDDO
@@ -196 +240 @@
        RETURN
 
-    END SUBROUTINE lsm_radiation
-
-
+    END SUBROUTINE radiation_clearsky
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!-- Prescribed net radiation
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_constant
+
+       rad_net = net_radiation
+
+    END SUBROUTINE radiation_constant
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!-- Implementation of the RRTM radiation_scheme
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_rrtm
+
+
+    END SUBROUTINE radiation_rrtm
 
  END MODULE radiation_model_mod

palm/trunk/SOURCE/read_3d_binary.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added support for binary input of land surface and radiation model data. In the
+! course of this work, new temporary arrays tmp_3d_soil1, tmp_3d_soil2 were
+! required.
 ! 
 ! Former revisions:
@@ -94 +96 @@
     USE kinds
 
+    USE land_surface_model_mod,                                                &
+        ONLY:  c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb_av, lai_av,        &
+               qsws_eb_av, qsws_liq_eb_av, qsws_soil_eb_av, qsws_veg_eb_av,    & 
+               land_surface, m_liq_eb, m_liq_eb_av, m_soil, m_soil_av,         &
+               nzb_soil, nzt_soil, shf_eb_av, t_soil, t_soil_av
+
     USE particle_attributes,                                                   &
         ONLY:  iran_part
 
     USE pegrid
+
+    USE radiation_model_mod,                                                   &
+        ONLY: rad_net_av, rad_sw_in_av
 
     USE random_function_mod,                                                   &
@@ -157 +168 @@
     REAL(wp) ::  rdummy
 
-    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  tmp_2d     !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3d     !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwul  !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwun  !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwur  !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwus  !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwvl  !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwvn  !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwvr  !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwvs  !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwwl  !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwwn  !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwwr  !:
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwws  !:
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  tmp_2d      !: temporary array for storing 2D data
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3d      !: temporary array for storing 3D data
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3d_soil1!: temporary array for storing 3D soil model data
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3d_soil2!: temporary array for storing 3D soil model data
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwul   !:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwun   !:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwur   !:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwus   !:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwvl   !:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwvn   !:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwvr   !:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwvs   !:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwwl   !:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwwn   !:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwwr   !:
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp_3dwws   !:
 
     REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  tmp_4d     !:
@@ -317 +330 @@
 !--    First compare the version numbers
        READ ( 13 )  version_on_file
-       binary_version = '3.4'
+       binary_version = '4.0'
        IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
           WRITE( message_string, * ) 'version mismatch concerning data ',      &
@@ -387 +400 @@
 !
 !--    Allocate temporary arrays sized as the arrays on the restart file
-       ALLOCATE( tmp_2d(nys_on_file-nbgp:nyn_on_file+nbgp,           &
-                        nxl_on_file-nbgp:nxr_on_file+nbgp),          &
-                 tmp_3d(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp, &
+       ALLOCATE( tmp_2d(nys_on_file-nbgp:nyn_on_file+nbgp,                     &
+                        nxl_on_file-nbgp:nxr_on_file+nbgp),                    &
+                 tmp_3d(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,           &
                         nxl_on_file-nbgp:nxr_on_file+nbgp) )
 
+!
+!--    If the land surface model is used, allocate array for storing soil model
+!--    data
+       IF ( land_surface )  THEN
+          ALLOCATE( tmp_3d_soil1(nzb_soil:nzt_soil+1,                          &
+                                nys_on_file-nbgp:nyn_on_file+nbgp,             &
+                                nxl_on_file-nbgp:nxr_on_file+nbgp) )
+          ALLOCATE( tmp_3d_soil2(nzb_soil:nzt_soil,                            &
+                                nys_on_file-nbgp:nyn_on_file+nbgp,             &
+                                nxl_on_file-nbgp:nxr_on_file+nbgp) )
+       ENDIF
 !
 !--    Read arrays
@@ -423 +447 @@
              SELECT CASE ( TRIM( field_chr ) )
 
+                CASE ( 'c_liq_av' )
+                   IF ( .NOT. ALLOCATED( c_liq_av ) )  THEN
+                      ALLOCATE( c_liq_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   c_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'c_soil_av' )
+                   IF ( .NOT. ALLOCATED( c_soil_av ) )  THEN
+                      ALLOCATE( c_soil_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   c_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'c_veg_av' )
+                   IF ( .NOT. ALLOCATED( c_veg_av ) )  THEN
+                      ALLOCATE( c_veg_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   c_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
                 CASE ( 'e' )
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
@@ -436 +484 @@
                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+                CASE ( 'ghf_eb_av' )
+                   IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
+                      ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   ghf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
                 CASE ( 'iran' ) ! matching random numbers is still unresolved
                                 ! issue
                    IF ( k == 1 )  READ ( 13 )  iran, iran_part
 
+                CASE ( 'lai_av' )
+                   IF ( .NOT. ALLOCATED( lai_av ) )  THEN
+                      ALLOCATE( lai_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   lai_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
                 CASE ( 'kh' )
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
@@ -465 +529 @@
                    lwp_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'm_liq_eb' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   m_liq_eb(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =        &
+                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                CASE ( 'm_liq_eb_av' )
+                   IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
+                      ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   m_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
+                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'm_soil' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d_soil2(:,:,:)
+                   m_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
+                          tmp_3d_soil2(nzb_soil:nzt_soil,nysf-nbgp:nynf        &
+                          +nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'm_soil_av' )
+                   IF ( .NOT. ALLOCATED( m_soil_av ) )  THEN
+                      ALLOCATE( m_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d_soil2(:,:,:)
+                   m_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
+                                    tmp_3d_soil2(nzb_soil:nzt_soil,nysf        &
+                                    -nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'nr' )
@@ -679 +770 @@
                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+                CASE ( 'qsws_eb_av' )
+                   IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
+                      ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qsws_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qsws_liq_eb_av' )
+                   IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
+                      ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qsws_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                CASE ( 'qsws_soil_eb_av' )
+                   IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
+                      ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qsws_soil_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qsws_veg_eb_av' )
+                   IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
+                      ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qsws_veg_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
                 CASE ( 'qv_av' )
                    IF ( .NOT. ALLOCATED( qv_av ) )  THEN
@@ -686 +808 @@
                    qv_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_net_av' )
+                   IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
+                      ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'rad_sw_in_av' )
+                   IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
+                      ALLOCATE( rad_sw_in_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF  
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   rad_sw_in_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'random_iv' )  ! still unresolved issue
@@ -759 +897 @@
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
                    shf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'shf_eb_av' )
+                   IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
+                      ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   shf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -808 +954 @@
                    tswst(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 't_soil' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d_soil1
+                   t_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
+                                   tmp_3d_soil1(:,nysf-nbgp:nynf+nbgp,         &
+                                                nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 't_soil_av' )
+                   IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
+                      ALLOCATE( t_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d_soil2(:,:,:)
+                   t_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
+                                    tmp_3d_soil2(:,nysf-nbgp:nynf+nbgp,        &
+                                    nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'u' )
@@ -1085 +1246 @@
 
        DEALLOCATE( tmp_2d, tmp_3d )
+       IF ( ALLOCATED( tmp_3d_soil1 ) )  DEALLOCATE( tmp_3d_soil1 )
+       IF ( ALLOCATED( tmp_3d_soil2 ) )  DEALLOCATE( tmp_3d_soil2 )
 
     ENDDO  ! loop over restart files

palm/trunk/SOURCE/read_var_list.f90

@@ -724 +724 @@
 ! ------------
 ! Skipping the global control variables from restart-file (binary format)
-! except some informations needed when reading restart data from a previous
+! except some information needed when reading restart data from a previous
 ! run which used a smaller total domain or/and a different domain decomposition.
 !------------------------------------------------------------------------------!

palm/trunk/SOURCE/sum_up_3d_data.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added support for land surface model and radiation model data.
 ! 
 ! Former revisions:
@@ -94 +94 @@
     USE kinds
 
+    USE land_surface_model_mod,                                                &
+        ONLY:  c_liq, c_liq_av, c_soil_av, c_veg, c_veg_av, ghf_eb,            &
+               ghf_eb_av, lai, lai_av, m_liq_eb, m_liq_eb_av, m_soil,          &
+               m_soil_av, nzb_soil, nzt_soil, qsws_eb, qsws_eb_av,             &
+               qsws_liq_eb, qsws_liq_eb_av, qsws_soil_eb, qsws_soil_eb_av,     &
+               qsws_veg_eb, qsws_veg_eb_av, shf_eb, shf_eb_av, t_soil,         &
+               t_soil_av
+
     USE particle_attributes,                                                   &
         ONLY:  grid_particles, number_of_particles, particles, prt_count
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
 
     IMPLICIT NONE
@@ -123 +134 @@
           SELECT CASE ( TRIM( doav(ii) ) )
 
+             CASE ( 'c_liq*' )
+                IF ( .NOT. ALLOCATED( c_liq_av ) )  THEN
+                   ALLOCATE( c_liq_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                c_liq_av = 0.0_wp
+
+             CASE ( 'c_soil*' )
+                IF ( .NOT. ALLOCATED( c_soil_av ) )  THEN
+                   ALLOCATE( c_soil_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                c_soil_av = 0.0_wp
+
+             CASE ( 'c_veg*' )
+                IF ( .NOT. ALLOCATED( c_veg_av ) )  THEN
+                   ALLOCATE( c_veg_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                c_veg_av = 0.0_wp
+
              CASE ( 'e' )
                 IF ( .NOT. ALLOCATED( e_av ) )  THEN
@@ -129 +158 @@
                 e_av = 0.0_wp
 
+             CASE ( 'ghf_eb*' )
+                IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
+                   ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                ghf_eb_av = 0.0_wp
+
+             CASE ( 'lai*' )
+                IF ( .NOT. ALLOCATED( lai_av ) )  THEN
+                   ALLOCATE( lai_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                lai_av = 0.0_wp
+
              CASE ( 'lpt' )
                 IF ( .NOT. ALLOCATED( lpt_av ) )  THEN
@@ -141 +182 @@
                 lwp_av = 0.0_wp
 
+             CASE ( 'm_liq_eb*' )
+                IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
+                   ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                m_liq_eb_av = 0.0_wp
+
+             CASE ( 'm_soil' )
+                IF ( .NOT. ALLOCATED( m_soil_av ) )  THEN
+                   ALLOCATE( m_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                m_soil_av = 0.0_wp
+
              CASE ( 'nr' )
                 IF ( .NOT. ALLOCATED( nr_av ) )  THEN
@@ -231 +284 @@
                 qsws_av = 0.0_wp
 
+             CASE ( 'qsws_eb*' )
+                IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
+                   ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                qsws_eb_av = 0.0_wp
+
+             CASE ( 'qsws_liq_eb*' )
+                IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
+                   ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                qsws_liq_eb_av = 0.0_wp
+
+             CASE ( 'qsws_soil_eb*' )
+                IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
+                   ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                qsws_soil_eb_av = 0.0_wp
+
+             CASE ( 'qsws_veg_eb*' )
+                IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
+                   ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                qsws_veg_eb_av = 0.0_wp
+
              CASE ( 'qv' )
                 IF ( .NOT. ALLOCATED( qv_av ) )  THEN
@@ -237 +314 @@
                 qv_av = 0.0_wp
 
+             CASE ( 'rad_net*' )
+                IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
+                   ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_net_av = 0.0_wp
+
+             CASE ( 'rad_sw_in*' )
+                IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
+                   ALLOCATE( rad_sw_in_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                rad_sw_in_av = 0.0_wp
+
              CASE ( 'rho' )
                 IF ( .NOT. ALLOCATED( rho_av ) )  THEN
@@ -261 +350 @@
                 shf_av = 0.0_wp
 
+             CASE ( 'shf_eb*' )
+                IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
+                   ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                shf_eb_av = 0.0_wp
+
+             CASE ( 't_soil' )
+                IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
+                   ALLOCATE( t_soil_av(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                t_soil_av = 0.0_wp
+
              CASE ( 't*' )
                 IF ( .NOT. ALLOCATED( ts_av ) )  THEN
@@ -328 +429 @@
        SELECT CASE ( TRIM( doav(ii) ) )
 
+          CASE ( 'c_liq*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_liq_av(j,i) = c_liq_av(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'c_soil*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_soil_av(j,i) = c_soil_av(j,i) + (1.0_wp - c_veg(j,i))
+                ENDDO
+             ENDDO
+
+          CASE ( 'c_veg*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   c_veg_av(j,i) = c_veg_av(j,i)
+                ENDDO
+             ENDDO
+
           CASE ( 'e' )
              DO  i = nxlg, nxrg
@@ -334 +456 @@
                       e_av(k,j,i) = e_av(k,j,i) + e(k,j,i)
                    ENDDO
+                ENDDO
+             ENDDO
+
+          CASE ( 'ghf_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   ghf_eb_av(j,i) = ghf_eb_av(j,i) + ghf_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'lai*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   lai_av(j,i) = lai_av(j,i)
                 ENDDO
              ENDDO
@@ -351 +487 @@
                    lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i) * &
                                                     dzw(1:nzt+1) )
+                ENDDO
+             ENDDO
+
+          CASE ( 'm_liq_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   m_liq_eb_av(j,i) = m_liq_eb_av(j,i) + m_liq_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'm_soil' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb_soil, nzt_soil
+                      m_soil_av(k,j,i) = m_soil_av(k,j,i) + m_soil(k,j,i)
+                   ENDDO
                 ENDDO
              ENDDO
@@ -518 +670 @@
              ENDDO
 
+          CASE ( 'qsws_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_eb_av(j,i) = qsws_eb_av(j,i) + qsws_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'qsws_liq_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) + qsws_liq_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'qsws_soil_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) + qsws_soil_eb(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'qsws_veg_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) + qsws_veg_eb(j,i)
+                ENDDO
+             ENDDO
+
           CASE ( 'qv' )
              DO  i = nxlg, nxrg
@@ -527 +707 @@
              ENDDO
 
+          CASE ( 'rad_net*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   rad_net_av(j,i) = rad_net_av(j,i) + rad_net(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 'rad_sw_in*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   rad_sw_in_av(j,i) = rad_sw_in_av(j,i) + rad_sw_in(j,i)
+                ENDDO
+             ENDDO
+
           CASE ( 'rho' )
              DO  i = nxlg, nxrg
@@ -561 +755 @@
              ENDDO
 
+          CASE ( 'shf_eb*' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   shf_eb_av(j,i) = shf_eb_av(j,i) + shf_eb(j,i)
+                ENDDO
+             ENDDO
+
           CASE ( 't*' )
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
                    ts_av(j,i) = ts_av(j,i) + ts(j,i)
+                ENDDO
+             ENDDO
+
+          CASE ( 't_soil' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb_soil, nzt_soil
+                      t_soil_av(k,j,i) = t_soil_av(k,j,i) + t_soil(k,j,i)
+                   ENDDO
                 ENDDO
              ENDDO

palm/trunk/SOURCE/swap_timelevel.f90

@@ -72 +72 @@
         ONLY:  e, e_p, nr, nr_p, pt, pt_p, q, q_p, qr, qr_p, sa, sa_p, u, u_p, &
                v, v_p, w, w_p
-
-    USE land_surface_model_mod,                                                &
-        ONLY: land_surface, m_liq, m_liq_p, m_soil, m_soil_p, T_0, T_0_p,      &
-              T_soil, T_soil_p
 #else
     USE arrays_3d,                                                             &
@@ -82 +78 @@
                u_1, u_2, u_p, v, v_1, v_2, v_p, w, w_1, w_2, w_p
 
+#endif
+
     USE land_surface_model_mod,                                                &
-        ONLY: land_surface, m_liq, m_liq_1, m_liq_2, m_liq_p, m_soil,          &
-              m_soil_1, m_soil_2, m_soil_p, T_0, T_0_1, T_0_2, T_0_p, T_soil,  &
-              T_soil_1, T_soil_2, T_soil_p
-#endif
+        ONLY: land_surface, lsm_swap_timelevel
 
     USE cpulog,                                                                &
@@ -129 +124 @@
 
     IF ( land_surface )  THEN
-       T_0    = T_0_p
-       T_soil = T_soil_p
-       IF ( humidity )  THEN
-          m_soil = m_soil_p
-          m_liq  = m_liq_p
-       ENDIF
+       CALL lsm_swap_timelevel ( 0 )
     ENDIF
 
@@ -168 +158 @@
 
           IF ( land_surface )  THEN
-             T_0    => T_0_1;    T_0_p    => T_0_2
-             T_soil => T_soil_1; T_soil_p => T_soil_2
-             IF ( humidity )  THEN
-                m_soil => m_soil_1; m_soil_p  => m_soil_2
-                m_liq  => m_liq_1;  m_liq_p   => m_liq_2
-             ENDIF
+             CALL lsm_swap_timelevel ( MOD( timestep_count, 2) )
           ENDIF
 
@@ -201 +186 @@
 
           IF ( land_surface )  THEN
-             T_0    => T_0_2;    T_0_p    => T_0_1
-             T_soil => T_soil_2; T_soil_p => T_soil_1
-             IF ( humidity )  THEN
-                m_soil => m_soil_2; m_soil_p  => m_soil_1
-                m_liq  => m_liq_2;  m_liq_p   => m_liq_1
-             ENDIF
+             CALL lsm_swap_timelevel ( MOD( timestep_count, 2) )
           ENDIF
 

palm/trunk/SOURCE/time_integration.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
+! Added interface for different radiation schemes.
 ! 
 ! Former revisions:
@@ -223 +223 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY: dt_radiation, lsm_radiation, radiation, time_radiation
+        ONLY: dt_radiation, irad_scheme, radiation, radiation_clearsky,        &
+              radiation_constant, radiation_rrtm, time_radiation
 
     USE statistics,                                                            &
@@ -601 +602 @@
 
                 time_radiation = time_radiation - dt_radiation
-                CALL lsm_radiation
+                IF ( irad_scheme == 0 )  THEN
+                   CALL radiation_constant              
+                ELSEIF ( irad_scheme == 1 )  THEN
+                   CALL radiation_clearsky
+                ELSEIF ( irad_scheme == 2 )  THEN
+                   CALL radiation_rrtm
+                ENDIF
 
                 CALL cpu_log( log_point(50), 'radiation', 'stop' )

palm/trunk/SOURCE/user_data_output_2d.f90

@@ -1 @@
- SUBROUTINE user_data_output_2d( av, variable, found, grid, local_pf, two_d )
+ SUBROUTINE user_data_output_2d( av, variable, found, grid, local_pf, two_d, nzb_do, nzt_do )
 
 !--------------------------------------------------------------------------------!
@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Replaced nzb and nzt+1 with the new array bounds nzb_do and nzt_do.
 ! 
 ! Former revisions:
@@ -61 +61 @@
     INTEGER(iwp) ::  j  !: 
     INTEGER(iwp) ::  k  !: 
+    INTEGER(iwp) ::  nzb_do !: lower limit of the domain (usually nzb)
+    INTEGER(iwp) ::  nzt_do !: upper limit of the domain (usually nzt+1)
 
     LOGICAL      ::  found !: 
-    LOGICAL      ::  two_d !: 
+    LOGICAL      ::  two_d !: flag parameter that indicates 2D variables (horizontal cross sections)
 
     REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !: 
@@ -81 +83 @@
 !             DO  i = nxlg, nxrg
 !                DO  j = nysg, nyng
-!                   DO  k = nzb, nzt+1
+!                   DO  k = nzb_do, nzt_do
 !                      local_pf(i,j,k) = u2(k,j,i)
 !                   ENDDO
@@ -89 +91 @@
 !             DO  i = nxlg, nxrg
 !                DO  j = nysg, nyng
-!                   DO  k = nzb, nzt+1
+!                   DO  k = nzb_do, nzt_do
 !                      local_pf(i,j,k) = u2_av(k,j,i)
 !                   ENDDO

palm/trunk/SOURCE/user_data_output_3d.f90

@@ -1 @@
- SUBROUTINE user_data_output_3d( av, variable, found, local_pf, nz_do )
+ SUBROUTINE user_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
 
 !--------------------------------------------------------------------------------!
@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
+! Replaced nzb and nzt+1 with the new array bounds nzb_do and nzt_do.
 ! 
 ! Former revisions:
@@ -62 +62 @@
     INTEGER(iwp) ::  j     !: 
     INTEGER(iwp) ::  k     !: 
-    INTEGER(iwp) ::  nz_do !: 
+    INTEGER(iwp) ::  nzb_do !: lower limit of the data output (usually 0)
+    INTEGER(iwp) ::  nzt_do !: vertical upper limit of the data output (usually nz_do3d)
 
     LOGICAL      ::  found !: 
 
-   REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nz_do) ::  local_pf !: 
+   REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) ::  local_pf !: 
 
 
@@ -82 +83 @@
 !             DO  i = nxlg, nxrg
 !                DO  j = nysg, nyng
-!                   DO  k = nzb, nz_do
+!                   DO  k = nzb_do, nzt_do
 !                      local_pf(i,j,k) = u2(k,j,i)
 !                   ENDDO
@@ -90 +91 @@
 !             DO  i = nxlg, nxrg
 !                DO  j = nysg, nyng
-!                   DO  k = nzb, nz_do
+!                   DO  k = nzb_do, nzt_do
 !                      local_pf(i,j,k) = u2_av(k,j,i)
 !                   ENDDO

palm/trunk/SOURCE/user_header.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Typo removed
 !
 ! Former revisions:
@@ -44 +44 @@
 ! Description:
 ! ------------
-! Print a header with user-defined informations.
+! Print a header with user-defined information.
 !------------------------------------------------------------------------------!
 
@@ -61 +61 @@
 !
 !-- If no user-defined variables are read from the namelist-file, no
-!-- informations will be printed.
+!-- information will be printed.
     IF ( .NOT. user_defined_namelist_found )  THEN
        WRITE ( io, 100 )
@@ -68 +68 @@
 
 !
-!-- Printing the informations.
+!-- Printing the information.
     WRITE ( io, 110 )
 

palm/trunk/SOURCE/write_3d_binary.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added support for binary ouput of land surface and radiation model data. 
 ! 
 ! Former revisions:
@@ -86 +86 @@
     USE indices,                                                               &
         ONLY:  nxl, nxr, nys, nyn, nzb, nzt
+
+    USE land_surface_model_mod,                                                &
+        ONLY:  c_liq_av, c_soil_av, c_veg_av, ghf_eb_av, lai_av,               &
+               qsws_eb_av, qsws_liq_eb_av, qsws_soil_eb_av, qsws_veg_eb_av,    & 
+               land_surface, m_liq_eb, m_liq_eb_av, m_soil, m_soil_av,         &
+               nzb_soil, nzt_soil, shf_eb_av, t_soil, t_soil_av
         
     USE particle_attributes,                                                   &
@@ -92 +98 @@
     USE pegrid
     
+    USE radiation_model_mod,                                                   &
+        ONLY: radiation, rad_net_av, rad_sw_in_av
+
     USE random_function_mod,                                                   &
         ONLY:  random_iv, random_iy
@@ -112 +121 @@
 !
 !-- Write arrays.
-    binary_version = '3.4'
+    binary_version = '4.0'
 
     WRITE ( 14 )  binary_version
@@ -124 +133 @@
 !--            adjusted accordingly.
 
+
+    IF ( land_surface )  THEN
+       IF ( ALLOCATED( c_liq_av ) )  THEN
+          WRITE ( 14 )  'c_liq_av                 ';  WRITE ( 14 ) c_liq_av
+       ENDIF
+       IF ( ALLOCATED( c_soil_av ) )  THEN
+          WRITE ( 14 )  'c_soil_av                ';  WRITE ( 14 ) c_soil_av
+       ENDIF
+       IF ( ALLOCATED( c_veg_av ) )  THEN
+          WRITE ( 14 )  'c_veg_av                 ';  WRITE ( 14 ) c_veg_av
+       ENDIF
+    ENDIF
     WRITE ( 14 )  'e                   ';  WRITE ( 14 )  e
     IF ( ALLOCATED( e_av ) )  THEN
@@ -129 +150 @@
     ENDIF
     WRITE ( 14 )  'iran                ';  WRITE ( 14 )  iran, iran_part
+    IF ( land_surface )  THEN
+       IF ( ALLOCATED( ghf_eb_av ) )  THEN
+          WRITE ( 14 )  'ghf_eb_av                ';  WRITE ( 14 )  ghf_eb_av
+       ENDIF
+    ENDIF
     WRITE ( 14 )  'kh                  ';  WRITE ( 14 )  kh
     WRITE ( 14 )  'km                  ';  WRITE ( 14 )  km
+    IF ( land_surface )  THEN
+       IF ( ALLOCATED( lai_av ) )  THEN
+          WRITE ( 14 )  'lai_av                   ';  WRITE ( 14 )  lai_av
+       ENDIF
+    ENDIF
     IF ( ALLOCATED( lpt_av ) )  THEN
        WRITE ( 14 )  'lpt_av              ';  WRITE ( 14 )  lpt_av
@@ -136 +167 @@
     IF ( ALLOCATED( lwp_av ) )  THEN
        WRITE ( 14 )  'lwp_av              ';  WRITE ( 14 )  lwp_av
+    ENDIF
+    IF ( land_surface )  THEN
+       WRITE ( 14 )  'm_liq_eb            ';  WRITE ( 14 )  m_liq_eb
+       IF ( ALLOCATED( m_liq_eb_av ) )  THEN
+          WRITE ( 14 )  'm_liq_eb_av              ';  WRITE ( 14 )  m_liq_eb_av
+       ENDIF
+       WRITE ( 14 )  'm_soil              ';  WRITE ( 14 )  m_soil
+       IF ( ALLOCATED( m_soil_av ) )  THEN
+          WRITE ( 14 )  'm_soil_av                ';  WRITE ( 14 )  m_soil_av
+       ENDIF
     ENDIF
     WRITE ( 14 )  'p                   ';  WRITE ( 14 )  p
@@ -204 +245 @@
        WRITE ( 14 )  'qswst               ';  WRITE ( 14 ) qswst
     ENDIF
+    IF ( land_surface )  THEN
+       IF ( ALLOCATED( qsws_eb_av ) )  THEN
+          WRITE ( 14 )  'qsws_eb_av          ';  WRITE ( 14 )  qsws_eb_av
+       ENDIF   
+       IF ( ALLOCATED( qsws_liq_eb_av ) )  THEN
+          WRITE ( 14 )  'qsws_liq_eb_av      ';  WRITE ( 14 )  qsws_liq_eb_av
+       ENDIF  
+       IF ( ALLOCATED( qsws_soil_eb_av ) )  THEN
+          WRITE ( 14 )  'qsws_soil_eb_av     ';  WRITE ( 14 )  qsws_soil_eb_av
+       ENDIF 
+       IF ( ALLOCATED( qsws_veg_eb_av ) )  THEN
+          WRITE ( 14 )  'qsws_veg_eb_av      ';  WRITE ( 14 )  qsws_veg_eb_av
+       ENDIF 
+    ENDIF
+    IF ( radiation )  THEN
+       IF ( ALLOCATED( rad_net_av ) )  THEN
+          WRITE ( 14 )  'rad_net_av          ';  WRITE ( 14 )  rad_net_av  
+       ENDIF  
+       IF ( ALLOCATED( rad_sw_in_av ) )  THEN
+          WRITE ( 14 )  'rad_sw_in_av          ';  WRITE ( 14 )  rad_sw_in_av  
+       ENDIF 
+    ENDIF
     IF ( ocean )  THEN
        IF ( ALLOCATED( rho_av ) )  THEN
@@ -214 +277 @@
        WRITE ( 14 )  'saswsb              ';  WRITE ( 14 )  saswsb
        WRITE ( 14 )  'saswst              ';  WRITE ( 14 )  saswst
+    ENDIF
+    IF ( land_surface )  THEN
+       WRITE ( 14 )  't_soil              ';  WRITE ( 14 )  t_soil
+       IF ( ALLOCATED( t_soil_av ) )  THEN
+          WRITE ( 14 )  't_soil_av                ';  WRITE ( 14 )  t_soil_av
+       ENDIF
     ENDIF
     IF ( ALLOCATED( ql_c_av ) )  THEN
@@ -244 +313 @@
        WRITE ( 14 )  'shf_av              ';  WRITE ( 14 )  shf_av
     ENDIF
+    IF ( land_surface )  THEN
+       IF ( ALLOCATED( shf_eb_av ) )  THEN
+          WRITE ( 14 )  'shf_eb_av           ';  WRITE ( 14 )  shf_eb_av
+       ENDIF
+    ENDIF
     IF ( ALLOCATED( spectrum_x ) )  THEN
        WRITE ( 14 )  'spectrum_x          ';  WRITE ( 14 )  spectrum_x

palm/trunk/SOURCE/write_var_list.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Typo removed
 ! 
 ! Former revisions:
@@ -116 +116 @@
 ! ------------
 ! Writing values of control variables to restart-file (binary format).
-! These informations are only written to the file opened by PE0.
+! This information are only written to the file opened by PE0.
 !------------------------------------------------------------------------------!
 

palm/trunk/UTIL/combine_plot_fields.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adjustments for data output of soil model quantities
 ! 
 ! Former revisions:
@@ -91 +91 @@
     CHARACTER (LEN=2)    ::  modus, model_string
     CHARACTER (LEN=6)    ::  id_string
-    CHARACTER (LEN=10)   ::  dimname, var_name
+    CHARACTER (LEN=20)   ::  dimname, var_name
     CHARACTER (LEN=40)   ::  filename
 
@@ -100 +100 @@
     INTEGER(iwp) ::  av, danz, i, id, j, k, model, models, nc_stat,            &
                      nxa, nxag, nxe, nxeg, nya, nyag, nye, nyeg,               &
-                     nza, nzag, nze, nzeg, pos, time_step,                     &
-                     xa, xe, xxa, xxe, ya, ye, yya, yye, za, ze, zza, zze
+                     nza, nze, pos, time_step,                                 &
+                     xa, xe, xxa, xxe, ya, ya_do, ya_tot, ye, ye_do, ye_tot, yya, yye, za,   &
+                     ze, zza, zze
 
 #if defined( __lc ) || defined( __decalpha )
@@ -372 +373 @@
 !--                concerning the lower and upper indices of the total-domain
 !--                used in PALM (nxag, nxeg, nyag, nyeg) and the lower and
-!--                upper indices of the array to be writte by this routine
+!--                upper indices of the array to be written by this routine
 !--                (nxa, nxe, nya, nye). Usually in the horizontal directions
 !--                nxag=-1 and nxa=0 while all other variables have the same
@@ -385 +386 @@
                       READ ( id+110 )  dx, eta, hu, ho
 
+
+!
+!--                   Set actual domain bounds to total domain
+                      ya_do = nya
+                      ye_do = nye
+
                       IF ( iso2d_output )  THEN
                          OPEN ( 2, FILE='PLOT2D_'//modus//TRIM( model_string ),&
@@ -405 +412 @@
 !
 !--                Read subdomain indices
-                   READ ( id+110, END=998 )  xa, xe, ya, ye
+                   READ ( id+110, END=998 )  xa, xe, ya, ye, ya_tot, ye_tot
+
 !
 !--                IF the PE made no output (in case that no part of the
@@ -412 +420 @@
                    IF ( .NOT. ( xa == -1  .AND.  xe == -1  .AND. &
                                 ya == -1  .AND.  ye == -1 ) )  THEN
+
+
 !
 !--                   Read the subdomain grid-point values
                       ALLOCATE( pf_tmp(xa:xe,ya:ye) )
                       READ ( id+110 )  pf_tmp
-                      pf(xa:xe,ya:ye) = pf_tmp
+!
+!--                   Calculate indices on atmospheric grid (e.g. for soil model
+!--                   quantities)
+                      IF ( ya /= ya_tot .OR. ye /= ye_tot )  THEN
+                         ye_do = ye - ya
+                         ya_do = ya
+                         pf(xa:xe,ya_do:ye_do) = pf_tmp
+                      ELSE
+                         ye_do = nye
+                         ya_do = nya
+                         pf(xa:xe,ya:ye) = pf_tmp
+                      ENDIF
+
                       DEALLOCATE( pf_tmp )
                    ENDIF
@@ -459 +481 @@
                          nc_stat = NF90_PUT_VAR( id_set(av),                   &
                                            id_var(av,current_var(av)),         &
-                                           pf(nxa:nxe,nya:nye),                &
+                                           pf(nxa:nxe,ya_do:ye_do),            &
                              start = (/ 1, current_level(av), 1, time_step /), &
-                                      count = (/ nxe-nxa+1, 1, nye-nya+1, 1 /) )
+                                      count = (/ nxe-nxa+1, 1, ye_do-ya_do+1, 1 /) )
                          IF ( nc_stat /= NF90_NOERR )  THEN
                             CALL handle_netcdf_error( 9 )
@@ -469 +491 @@
                          nc_stat = NF90_PUT_VAR( id_set(av),                   &
                                            id_var(av,current_var(av)),         &
-                                           pf(nxa:nxe,nya:nye),                &
+                                           pf(nxa:nxe,ya_do:ye_do),            &
                              start = (/ current_level(av), 1, 1, time_step /), &
-                                      count = (/ 1, nxe-nxa+1, nye-nya+1, 1 /) )
+                                      count = (/ 1, nxe-nxa+1, ye_do-ya_do+1, 1 /) )
                          IF ( nc_stat /= NF90_NOERR )  THEN
                             CALL handle_netcdf_error( 10 )
@@ -743 +765 @@
 !--             File from PE0 contains special information at the beginning,
 !--             concerning the lower and upper indices of the total-domain used
-!--             in PALM (nxag, nxeg, nyag, nyeg, nzag, nzeg) and the lower and
+!--             in PALM (nxag, nxeg, nyag, nyeg) and the lower and
 !--             upper indices of the array to be written by this routine (nxa,
 !--             nxe, nya, nye, nza, nze). Usually nxag=-1 and nxa=0, nyag=-1
-!--             and nya=0, nzeg=nz and nze=nz_plot3d.
+!--             and nya=0, nzeg=nz and nze=nz_do3d.
 !--             Allocate necessary array and open the output file.
                 IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 )  THEN
-                   READ ( id+110 )  nxag, nxeg, nyag, nyeg, nzag, nzeg
+                   READ ( id+110 )  nxag, nxeg, nyag, nyeg
                    READ ( id+110 )  nxa, nxe, nya, nye, nza, nze
                    ALLOCATE ( pf3d(nxa:nxe,nya:nye,nza:nze) )
@@ -781 +803 @@
                 yya = MAX( nya, ya )
                 yye = MIN( nye, ye )
-                zza = MAX( nza, za )
-                zze = MIN( nze, ze )
-                DO  k = zza, zze
+                DO  k = za, ze
                    DO  j = yya, yye
                       DO  i = xxa, xxe
@@ -818 +838 @@
 !--             Now write the data
                 nc_stat = NF90_PUT_VAR( id_set(av), id_var(av,current_var(av)),&
-                                        pf3d, start = (/ 1, 1, 1, time_step /),&
-                              count = (/ nxe-nxa+1, nye-nya+1, nze-nza+1, 1 /) )
+                                        pf3d(nxa:nxe,nya:nye,za:ze), start = (/ 1, 1, 1, time_step /),&
+                              count = (/ nxe-nxa+1, nye-nya+1, ze-za+1, 1 /) )
                 IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 18 )