Index: palm/trunk/SOURCE/model_1d_mod.f90
===================================================================
--- palm/trunk/SOURCE/init_1d_model.f90	(revision 2337)
+++ palm/trunk/SOURCE/model_1d_mod.f90	(revision 2338)
@@ -1,3 +1,3 @@
-!> @file init_1d_model.f90
+!> @file model_1d_mod.f90
 !------------------------------------------------------------------------------!
 ! This file is part of PALM.
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! renamed init_1d_model to model_1d_mod and and formatted it as a module;
+! reformatted output of profiles
+! 
+! 
 ! revised calculation of mixing length
 ! removed rounding of time step
@@ -113,25 +117,120 @@
 !> @bug 1D model crashes when using small grid spacings in the order of 1 m
 !------------------------------------------------------------------------------!
+ MODULE model_1d_mod
+
+    USE arrays_3d,                                                             &
+        ONLY:  dd2zu, ddzu, ddzw, dzu, l_grid, pt_init, q_init, ug, u_init,    &
+               vg, v_init, zu
+    
+    USE control_parameters,                                                    &
+        ONLY:  constant_diffusion, constant_flux_layer, dissipation_1d, f, g,  &
+               humidity, ibc_e_b, intermediate_timestep_count,                 &
+               intermediate_timestep_count_max, kappa, km_constant,            &
+               message_string, mixing_length_1d, prandtl_number,               &
+               roughness_length, run_description_header, simulated_time_chr, timestep_scheme, tsc, z0h_factor
+
+    USE indices,                                                               &
+        ONLY:  nzb, nzb_diff, nzt
+    
+    USE kinds
+
+    USE pegrid
+        
+
+    IMPLICIT NONE
+
+
+    INTEGER(iwp) ::  current_timestep_number_1d = 0  !< current timestep number (1d-model)
+    INTEGER(iwp) ::  damp_level_ind_1d               !< lower grid index of damping layer (1d-model) !#
+
+    LOGICAL ::  run_control_header_1d = .FALSE.  !< flag for output of run control header (1d-model)
+    LOGICAL ::  stop_dt_1d = .FALSE.             !< termination flag, used in case of too small timestep (1d-model)
+
+    REAL(wp) ::  c_m = 0.4_wp                  !< model constant, 0.4 according to Detering and Etling (1985)
+    REAL(wp) ::  damp_level_1d = -1.0_wp       !< namelist parameter    !#
+    REAL(wp) ::  dt_1d = 60.0_wp               !< dynamic timestep (1d-model)
+    REAL(wp) ::  dt_max_1d = 300.0_wp          !< timestep limit (1d-model)
+    REAL(wp) ::  dt_pr_1d = 9999999.9_wp       !< namelist parameter     !#
+    REAL(wp) ::  dt_run_control_1d = 60.0_wp   !< namelist parameter     !#
+    REAL(wp) ::  end_time_1d = 864000.0_wp     !< namelist parameter     !#
+    REAL(wp) ::  qs1d                          !< characteristic humidity scale (1d-model)
+    REAL(wp) ::  simulated_time_1d = 0.0_wp    !< updated simulated time (1d-model)
+    REAL(wp) ::  time_pr_1d = 0.0_wp           !< updated simulated time for profile output (1d-model)
+    REAL(wp) ::  time_run_control_1d = 0.0_wp  !< updated simulated time for run-control output (1d-model)
+    REAL(wp) ::  ts1d                          !< characteristic temperature scale (1d-model)
+    REAL(wp) ::  us1d                          !< friction velocity (1d-model)                 !#
+    REAL(wp) ::  usws1d                        !< u-component of the momentum flux (1d-model)  !#
+    REAL(wp) ::  vsws1d                        !< v-component of the momentum flux (1d-model)  !#
+    REAL(wp) ::  z01d                          !< roughness length for momentum (1d-model)
+    REAL(wp) ::  z0h1d                         !< roughness length for scalars (1d-model)
+
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  e1d      !< tke (1d-model)                         !#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  e1d_p    !< prognostic value of tke (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  kh1d     !< turbulent diffusion coefficient for heat (1d-model)    !#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  km1d     !< turbulent diffusion coefficient for momentum (1d-model)!#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black  !< mixing length Blackadar (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d      !< mixing length for turbulent diffusion coefficients (1d-model) !#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l1d_diss !< mixing length for dissipation (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  rif1d    !< Richardson flux number (1d-model)      !#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_e     !< tendency of e (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_em    !< weighted tendency of e for previous sub-timestep (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_u     !< tendency of u (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_um    !< weighted tendency of u for previous sub-timestep (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_v     !< tendency of v (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  te_vm    !< weighted tendency of v for previous sub-timestep (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u1d      !< u-velocity component (1d-model)        !#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u1d_p    !< prognostic value of u-velocity component (1d-model)
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v1d      !< v-velocity component (1d-model)        !#
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v1d_p    !< prognostic value of v-velocity component (1d-model)
+
+!
+!-- Initialize 1D model
+    INTERFACE init_1d_model
+       MODULE PROCEDURE init_1d_model
+    END INTERFACE init_1d_model
+
+!
+!-- Print profiles
+    INTERFACE print_1d_model
+       MODULE PROCEDURE print_1d_model
+    END INTERFACE print_1d_model
+
+!
+!-- Print run control information
+    INTERFACE run_control_1d
+       MODULE PROCEDURE run_control_1d
+    END INTERFACE run_control_1d
+
+!
+!-- Main procedure
+    INTERFACE time_integration_1d
+       MODULE PROCEDURE time_integration_1d
+    END INTERFACE time_integration_1d
+
+!
+!-- Calculate time step
+    INTERFACE timestep_1d
+       MODULE PROCEDURE timestep_1d
+    END INTERFACE timestep_1d
+
+    SAVE
+
+    PRIVATE
+!
+!-- Public interfaces
+    PUBLIC  init_1d_model
+
+!
+!-- Public variables
+    PUBLIC  damp_level_1d, damp_level_ind_1d, dt_pr_1d, dt_run_control_1d,     &
+            e1d, end_time_1d, kh1d, km1d, l1d, rif1d, u1d, us1d, usws1d, v1d,  &
+            vsws1d
+
+
+    CONTAINS
+
  SUBROUTINE init_1d_model
  
-
-    USE arrays_3d,                                                             &
-        ONLY:  l_grid, ug, u_init, vg, v_init, zu
-    
-    USE indices,                                                               &
-        ONLY:  nzb, nzt
-    
-    USE kinds
-    
-    USE model_1d,                                                              &
-        ONLY:  e1d, e1d_p, kh1d, km1d, l1d, l1d_diss, l_black, qs1d, rif1d,    &
-               simulated_time_1d, te_e, te_em, te_u, te_um, te_v, te_vm, ts1d, &
-               u1d, u1d_p, us1d, usws1d, v1d, v1d_p, vsws1d, z01d, z0h1d
-    
-    USE control_parameters,                                                    &
-        ONLY:  constant_diffusion, constant_flux_layer, f, humidity, kappa,    &
-               km_constant, mixing_length_1d, prandtl_number,                  &
-               roughness_length, simulated_time_chr, z0h_factor
-
     IMPLICIT NONE
 
@@ -144,12 +243,10 @@
 !
 !-- Allocate required 1D-arrays
-    ALLOCATE( e1d(nzb:nzt+1),    e1d_p(nzb:nzt+1),                             &
-              kh1d(nzb:nzt+1),   km1d(nzb:nzt+1),                              &
-              l_black(nzb:nzt+1), l1d(nzb:nzt+1), l1d_diss(nzb:nzt+1),         &
-              rif1d(nzb:nzt+1),   te_e(nzb:nzt+1),                             &
-              te_em(nzb:nzt+1),  te_u(nzb:nzt+1),    te_um(nzb:nzt+1),         &
-              te_v(nzb:nzt+1),   te_vm(nzb:nzt+1),    u1d(nzb:nzt+1),          &
-              u1d_p(nzb:nzt+1),  v1d(nzb:nzt+1),                               &
-              v1d_p(nzb:nzt+1) )
+    ALLOCATE( e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1),               &
+              km1d(nzb:nzt+1), l_black(nzb:nzt+1), l1d(nzb:nzt+1),             &
+              l1d_diss(nzb:nzt+1), rif1d(nzb:nzt+1), te_e(nzb:nzt+1),          &
+              te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1),             &
+              te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1),               &
+              u1d_p(nzb:nzt+1),  v1d(nzb:nzt+1), v1d_p(nzb:nzt+1) )
 
 !
@@ -228,4 +325,8 @@
 
 !
+!-- Set model constant
+    IF ( dissipation_1d == 'as_in_3d_model' )  c_m = 0.1_wp
+
+!
 !-- Set start time in hh:mm:ss - format
     simulated_time_chr = time_to_string( simulated_time_1d )
@@ -243,52 +344,25 @@
 ! Description:
 ! ------------
-!> Leap-frog time differencing scheme for the 1D-model.
+!> Runge-Kutta time differencing scheme for the 1D-model.
 !------------------------------------------------------------------------------!
  
  SUBROUTINE time_integration_1d
 
-
-    USE arrays_3d,                                                             &
-        ONLY:  dd2zu, ddzu, ddzw, l_grid, pt_init, q_init, ug, vg, zu
-        
-    USE control_parameters,                                                    &
-        ONLY:  constant_diffusion, constant_flux_layer, dissipation_1d,        &
-               humidity, intermediate_timestep_count,                          &
-               intermediate_timestep_count_max, f, g, ibc_e_b, kappa,          &  
-               mixing_length_1d,                                               &
-               simulated_time_chr, timestep_scheme, tsc
-                
-    USE indices,                                                               &
-        ONLY:  nzb, nzb_diff, nzt
-        
-    USE kinds
-    
-    USE model_1d,                                                              &
-        ONLY:  current_timestep_number_1d, damp_level_ind_1d, dt_1d,           &
-               dt_pr_1d, dt_run_control_1d, e1d, e1d_p, end_time_1d,           &
-               kh1d, km1d, l1d, l1d_diss, l_black, qs1d, rif1d, simulated_time_1d, &
-               stop_dt_1d, te_e, te_em, te_u, te_um, te_v, te_vm, time_pr_1d,  &
-               ts1d, time_run_control_1d, u1d, u1d_p, us1d, usws1d, v1d,       &
-               v1d_p, vsws1d, z01d, z0h1d
-        
-    USE pegrid
-
     IMPLICIT NONE
 
     CHARACTER (LEN=9) ::  time_to_string  !<
     
-    INTEGER(iwp) ::  k  !<
+    INTEGER(iwp) ::  k  !< loop index
     
-    REAL(wp) ::  a            !<
-    REAL(wp) ::  b            !<
-    REAL(wp) ::  c_m = 0.4_wp !< model constant, 0.4 according to Detering and Etling (1985)
-    REAL(wp) ::  dissipation  !<
-    REAL(wp) ::  dpt_dz       !<
-    REAL(wp) ::  flux         !<
-    REAL(wp) ::  kmzm         !<
-    REAL(wp) ::  kmzp         !<
-    REAL(wp) ::  l_stable     !<
-    REAL(wp) ::  pt_0         !<
-    REAL(wp) ::  uv_total     !<
+    REAL(wp) ::  a            !< auxiliary variable
+    REAL(wp) ::  b            !< auxiliary variable
+    REAL(wp) ::  dissipation  !< dissipation of TKE
+    REAL(wp) ::  dpt_dz       !< vertical temperature gradient
+    REAL(wp) ::  flux         !< vertical temperature gradient
+    REAL(wp) ::  kmzm         !< Km(z-dz/2)
+    REAL(wp) ::  kmzp         !< Km(z+dz/2)
+    REAL(wp) ::  l_stable     !< mixing length for stable case
+    REAL(wp) ::  pt_0         !< reference temperature
+    REAL(wp) ::  uv_total     !< horizontal wind speed
 
 !
@@ -759,28 +833,14 @@
     USE constants,                                                             &
         ONLY:  pi
-        
-    USE indices,                                                               &
-        ONLY:  nzb, nzt
-        
-    USE kinds
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  k     !< loop index
     
-    USE model_1d,                                                              &
-        ONLY:  current_timestep_number_1d, dt_1d, run_control_header_1d, u1d,  &
-               us1d, v1d
-    
-    USE pegrid
-    
-    USE control_parameters,                                                    &
-        ONLY:  simulated_time_chr
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  k  !<
-    
-    REAL(wp) ::  alpha 
-    REAL(wp) ::  energy 
-    REAL(wp) ::  umax
-    REAL(wp) ::  uv_total 
-    REAL(wp) ::  vmax
+    REAL(wp) ::  alpha     !< angle of wind vector at top of constant-flux layer
+    REAL(wp) ::  energy    !< kinetic energy
+    REAL(wp) ::  umax      !< maximum of u
+    REAL(wp) ::  uv_total  !< horizontal wind speed
+    REAL(wp) ::  vmax      !< maximum of v
 
 !
@@ -844,30 +904,11 @@
  SUBROUTINE timestep_1d
 
-
-    USE arrays_3d,                                                             &
-        ONLY:  dzu, zu
-        
-    USE indices,                                                               &
-        ONLY:  nzb, nzt
-    
-    USE kinds
-    
-    USE model_1d,                                                              &
-        ONLY:  dt_1d, dt_max_1d, km1d, stop_dt_1d
-    
-    USE pegrid
-    
-    USE control_parameters,                                                    &
-        ONLY:  message_string
-
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  k !<
-    
-    REAL(wp) ::  div      !<
-    REAL(wp) ::  dt_diff  !<
-    REAL(wp) ::  fac      !<
-    REAL(wp) ::  value    !<
-
+    INTEGER(iwp) ::  k    !< loop index
+
+    REAL(wp) ::  dt_diff  !< time step accorind to diffusion criterion
+    REAL(wp) ::  fac      !< factor of criterion
+    REAL(wp) ::  value    !< auxiliary variable
 
 !
@@ -906,25 +947,9 @@
  SUBROUTINE print_1d_model
 
-
-    USE arrays_3d,                                                             &
-        ONLY:  pt_init, zu
-        
-    USE indices,                                                               &
-        ONLY:  nzb, nzt
-        
-    USE kinds
-    
-    USE model_1d,                                                              &
-        ONLY:  e1d, kh1d, km1d, l1d, rif1d, u1d, v1d
-    
-    USE pegrid
-    
-    USE control_parameters,                                                    &
-        ONLY:  run_description_header, simulated_time_chr
-
     IMPLICIT NONE
 
-
-    INTEGER(iwp) ::  k  !<
+    INTEGER(iwp) ::  k  !< loop parameter
+
+    LOGICAL, SAVE :: write_first = .TRUE. !< flag for writing header
 
 
@@ -936,10 +961,13 @@
 !
 !--    Write Header
-       WRITE ( 17, 100 )  TRIM( run_description_header ), &
-                          TRIM( simulated_time_chr )
+       IF ( write_first )  THEN
+          WRITE ( 17, 100 )  TRIM( run_description_header )
+          write_first = .FALSE.
+       ENDIF
+
+!
+!--    Write the values
+       WRITE ( 17, 104 )  TRIM( simulated_time_chr )
        WRITE ( 17, 101 )
-
-!
-!--    Write the values
        WRITE ( 17, 102 )
        WRITE ( 17, 101 )
@@ -960,8 +988,8 @@
 !
 !-- Formats
-100 FORMAT (//1X,A/1X,10('-')/' 1d-model profiles'/ &
-            ' Time: ',A)
-101 FORMAT (1X,79('-'))
-102 FORMAT ('   k     zu      u      v     pt      e    rif    Km    Kh     ', &
+100 FORMAT ('# ',A/'#',10('-')/'# 1d-model profiles')
+104 FORMAT (//'# Time: ',A)
+101 FORMAT ('#',79('-'))
+102 FORMAT ('#  k     zu      u      v     pt      e    rif    Km    Kh     ', &
             'l      zu      k')
 103 FORMAT (1X,I4,1X,F7.1,1X,F6.2,1X,F6.2,1X,F6.2,1X,F6.2,1X,F5.2,1X,F5.2, &
@@ -970,2 +998,5 @@
 
  END SUBROUTINE print_1d_model
+
+
+ END MODULE
