Changeset 1691 for palm/trunk/SOURCE/init_1d_model.f90

Timestamp:

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

Author:

maronga

Message:

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

File:

• palm/trunk/SOURCE/init_1d_model.f90 (modified) (13 diffs)

palm/trunk/SOURCE/init_1d_model.f90

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