SUBROUTINE init_1d_model !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: init_1d_model.f90 979 2012-08-09 08:50:11Z hoffmann $ ! ! 978 2012-08-09 08:28:32Z fricke ! roughness length for scalar quantities z0h1d added ! ! 667 2010-12-23 12:06:00Z suehring/gryschka ! replaced mirror boundary conditions for u and v at the ground ! by dirichlet boundary conditions ! ! 254 2009-03-05 15:33:42Z heinze ! Output of messages replaced by message handling routine. ! ! 184 2008-08-04 15:53:39Z letzel ! provisional solution for run_control_1d output: add 'CALL check_open( 15 )' ! ! 135 2007-11-22 12:24:23Z raasch ! Bugfix: absolute value of f must be used when calculating the Blackadar ! mixing length ! ! 82 2007-04-16 15:40:52Z raasch ! Preprocessor strings for different linux clusters changed to "lc", ! routine local_flush is used for buffer flushing ! ! 75 2007-03-22 09:54:05Z raasch ! Bugfix: preset of tendencies te_em, te_um, te_vm, ! moisture renamed humidity ! ! RCS Log replace by Id keyword, revision history cleaned up ! ! Revision 1.21 2006/06/02 15:19:57 raasch ! cpp-directives extended for lctit ! ! Revision 1.1 1998/03/09 16:22:10 raasch ! Initial revision ! ! ! Description: ! ------------ ! 1D-model to initialize the 3D-arrays. ! The temperature profile is set as steady and a corresponding steady solution ! of the wind profile is being computed. ! All subroutines required can be found within this file. !------------------------------------------------------------------------------! USE arrays_3d USE indices USE model_1d USE control_parameters IMPLICIT NONE CHARACTER (LEN=9) :: time_to_string INTEGER :: k REAL :: lambda ! !-- Allocate required 1D-arrays ALLOCATE( e1d(nzb:nzt+1), e1d_m(nzb:nzt+1), e1d_p(nzb:nzt+1), & kh1d(nzb:nzt+1), kh1d_m(nzb:nzt+1), km1d(nzb:nzt+1), & km1d_m(nzb:nzt+1), l_black(nzb:nzt+1), l1d(nzb:nzt+1), & l1d_m(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_m(nzb:nzt+1), u1d_p(nzb:nzt+1), v1d(nzb:nzt+1), & v1d_m(nzb:nzt+1), v1d_p(nzb:nzt+1) ) ! !-- Initialize arrays IF ( constant_diffusion ) THEN km1d = km_constant km1d_m = km_constant kh1d = km_constant / prandtl_number kh1d_m = km_constant / prandtl_number ELSE e1d = 0.0; e1d_m = 0.0; e1d_p = 0.0 kh1d = 0.0; kh1d_m = 0.0; km1d = 0.0; km1d_m = 0.0 rif1d = 0.0 ! !-- Compute the mixing length l_black(nzb) = 0.0 IF ( TRIM( mixing_length_1d ) == 'blackadar' ) THEN ! !-- Blackadar mixing length IF ( f /= 0.0 ) THEN lambda = 2.7E-4 * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) / & ABS( f ) + 1E-10 ELSE lambda = 30.0 ENDIF DO k = nzb+1, nzt+1 l_black(k) = kappa * zu(k) / ( 1.0 + kappa * zu(k) / lambda ) ENDDO ELSEIF ( TRIM( mixing_length_1d ) == 'as_in_3d_model' ) THEN ! !-- Use the same mixing length as in 3D model l_black(1:nzt) = l_grid l_black(nzt+1) = l_black(nzt) ENDIF ! !-- Adjust mixing length to the prandtl mixing length (within the prandtl !-- layer) IF ( adjust_mixing_length .AND. prandtl_layer ) THEN k = nzb+1 l_black(k) = MIN( l_black(k), kappa * zu(k) ) ENDIF ENDIF l1d = l_black l1d_m = l_black u1d = u_init u1d_m = u_init u1d_p = u_init v1d = v_init v1d_m = v_init v1d_p = v_init ! !-- Set initial horizontal velocities at the lowest grid levels to a very small !-- value in order to avoid too small time steps caused by the diffusion limit !-- in the initial phase of a run (at k=1, dz/2 occurs in the limiting formula!) u1d(0:1) = 0.1 u1d_m(0:1) = 0.1 u1d_p(0:1) = 0.1 v1d(0:1) = 0.1 v1d_m(0:1) = 0.1 v1d_p(0:1) = 0.1 ! !-- For u*, theta* and the momentum fluxes plausible values are set IF ( prandtl_layer ) THEN us1d = 0.1 ! without initial friction the flow would not change ELSE e1d(nzb+1) = 1.0 km1d(nzb+1) = 1.0 us1d = 0.0 ENDIF ts1d = 0.0 usws1d = 0.0; usws1d_m = 0.0 vsws1d = 0.0; vsws1d_m = 0.0 z01d = roughness_length z0h1d = z0h_factor * z01d IF ( humidity .OR. passive_scalar ) qs1d = 0.0 ! !-- Tendencies must be preset in order to avoid runtime errors within the !-- first Runge-Kutta step te_em = 0.0 te_um = 0.0 te_vm = 0.0 ! !-- Set start time in hh:mm:ss - format simulated_time_chr = time_to_string( simulated_time_1d ) ! !-- Integrate the 1D-model equations using the leap-frog scheme CALL time_integration_1d END SUBROUTINE init_1d_model SUBROUTINE time_integration_1d !------------------------------------------------------------------------------! ! Description: ! ------------ ! Leap-frog time differencing scheme for the 1D-model. !------------------------------------------------------------------------------! USE arrays_3d USE control_parameters USE indices USE model_1d USE pegrid IMPLICIT NONE CHARACTER (LEN=9) :: time_to_string INTEGER :: k REAL :: a, b, dissipation, dpt_dz, flux, kmzm, kmzp, l_stable, pt_0, & uv_total ! !-- Determine the time step at the start of a 1D-simulation and !-- determine and printout quantities used for run control CALL timestep_1d CALL run_control_1d ! !-- Start of time loop DO WHILE ( simulated_time_1d < end_time_1d .AND. .NOT. stop_dt_1d ) ! !-- Depending on the timestep scheme, carry out one or more intermediate !-- timesteps intermediate_timestep_count = 0 DO WHILE ( intermediate_timestep_count < & intermediate_timestep_count_max ) intermediate_timestep_count = intermediate_timestep_count + 1 CALL timestep_scheme_steering ! !-- Compute all tendency terms. If a Prandtl-layer is simulated, k starts !-- at nzb+2. DO k = nzb_diff, nzt kmzm = 0.5 * ( km1d_m(k-1) + km1d_m(k) ) kmzp = 0.5 * ( km1d_m(k) + km1d_m(k+1) ) ! !-- u-component te_u(k) = f * ( v1d(k) - vg(k) ) + ( & kmzp * ( u1d_m(k+1) - u1d_m(k) ) * ddzu(k+1) & - kmzm * ( u1d_m(k) - u1d_m(k-1) ) * ddzu(k) & ) * ddzw(k) ! !-- v-component te_v(k) = -f * ( u1d(k) - ug(k) ) + ( & kmzp * ( v1d_m(k+1) - v1d_m(k) ) * ddzu(k+1) & - kmzm * ( v1d_m(k) - v1d_m(k-1) ) * ddzu(k) & ) * ddzw(k) ENDDO IF ( .NOT. constant_diffusion ) THEN DO k = nzb_diff, nzt ! !-- TKE kmzm = 0.5 * ( km1d_m(k-1) + km1d_m(k) ) kmzp = 0.5 * ( km1d_m(k) + km1d_m(k+1) ) IF ( .NOT. humidity ) THEN pt_0 = pt_init(k) flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k) ELSE pt_0 = pt_init(k) * ( 1.0 + 0.61 * q_init(k) ) flux = ( ( pt_init(k+1) - pt_init(k-1) ) + & 0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) ) & ) * dd2zu(k) ENDIF IF ( dissipation_1d == 'detering' ) THEN ! !-- According to Detering, c_e=0.064 dissipation = 0.064 * e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k) ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN dissipation = ( 0.19 + 0.74 * l1d_m(k) / l_grid(k) ) & * e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k) ENDIF te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2& + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2& ) & - g / pt_0 * kh1d(k) * flux & + ( & kmzp * ( e1d_m(k+1) - e1d_m(k) ) * ddzu(k+1) & - kmzm * ( e1d_m(k) - e1d_m(k-1) ) * ddzu(k) & ) * ddzw(k) & - dissipation ENDDO ENDIF ! !-- Tendency terms at the top of the Prandtl-layer. !-- 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 k = nzb+1 kmzm = 0.5 * ( km1d_m(k-1) + km1d_m(k) ) kmzp = 0.5 * ( km1d_m(k) + km1d_m(k+1) ) IF ( .NOT. humidity ) THEN pt_0 = pt_init(k) flux = ( pt_init(k+1)-pt_init(k-1) ) * dd2zu(k) ELSE pt_0 = pt_init(k) * ( 1.0 + 0.61 * q_init(k) ) flux = ( ( pt_init(k+1) - pt_init(k-1) ) + & 0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) ) & ) * dd2zu(k) ENDIF IF ( dissipation_1d == 'detering' ) THEN ! !-- According to Detering, c_e=0.064 dissipation = 0.064 * e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k) ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN dissipation = ( 0.19 + 0.74 * l1d_m(k) / l_grid(k) ) & * e1d_m(k) * SQRT( e1d_m(k) ) / l1d_m(k) ENDIF ! !-- u-component te_u(k) = f * ( v1d(k) - vg(k) ) + ( & kmzp * ( u1d_m(k+1) - u1d_m(k) ) * ddzu(k+1) + usws1d_m & ) * 2.0 * ddzw(k) ! !-- v-component te_v(k) = -f * ( u1d(k) - ug(k) ) + ( & kmzp * ( v1d_m(k+1) - v1d_m(k) ) * ddzu(k+1) + vsws1d_m & ) * 2.0 * ddzw(k) ! !-- TKE te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & ) & - g / pt_0 * kh1d(k) * flux & + ( & kmzp * ( e1d_m(k+1) - e1d_m(k) ) * ddzu(k+1) & - kmzm * ( e1d_m(k) - e1d_m(k-1) ) * ddzu(k) & ) * ddzw(k) & - dissipation ENDIF ! !-- Prognostic equations for all 1D variables DO k = nzb+1, nzt u1d_p(k) = ( 1. - tsc(1) ) * u1d_m(k) + & tsc(1) * u1d(k) + dt_1d * ( tsc(2) * te_u(k) + & tsc(3) * te_um(k) ) v1d_p(k) = ( 1. - tsc(1) ) * v1d_m(k) + & tsc(1) * v1d(k) + dt_1d * ( tsc(2) * te_v(k) + & tsc(3) * te_vm(k) ) ENDDO IF ( .NOT. constant_diffusion ) THEN DO k = nzb+1, nzt e1d_p(k) = ( 1. - tsc(1) ) * e1d_m(k) + & tsc(1) * e1d(k) + dt_1d * ( tsc(2) * te_e(k) + & tsc(3) * te_em(k) ) ENDDO ! !-- Eliminate negative TKE values, which can result from the !-- integration due to numerical inaccuracies. In such cases the TKE !-- value is reduced to 10 percent of its old value. WHERE ( e1d_p < 0.0 ) e1d_p = 0.1 * e1d ENDIF ! !-- Calculate tendencies for the next Runge-Kutta step IF ( timestep_scheme(1:5) == 'runge' ) THEN IF ( intermediate_timestep_count == 1 ) THEN DO k = nzb+1, nzt te_um(k) = te_u(k) te_vm(k) = te_v(k) ENDDO IF ( .NOT. constant_diffusion ) THEN DO k = nzb+1, nzt te_em(k) = te_e(k) ENDDO ENDIF ELSEIF ( intermediate_timestep_count < & intermediate_timestep_count_max ) THEN DO k = nzb+1, nzt te_um(k) = -9.5625 * te_u(k) + 5.3125 * te_um(k) te_vm(k) = -9.5625 * te_v(k) + 5.3125 * te_vm(k) ENDDO IF ( .NOT. constant_diffusion ) THEN DO k = nzb+1, nzt te_em(k) = -9.5625 * te_e(k) + 5.3125 * te_em(k) ENDDO ENDIF ENDIF ENDIF ! !-- Boundary conditions for the prognostic variables. !-- At the top boundary (nzt+1) u,v and e keep their initial values !-- (ug(nzt+1), vg(nzt+1), 0), at the bottom boundary the mirror !-- boundary condition applies to u and v. !-- The boundary condition for e is set further below ( (u*/cm)**2 ). ! u1d_p(nzb) = -u1d_p(nzb+1) ! v1d_p(nzb) = -v1d_p(nzb+1) u1d_p(nzb) = 0.0 v1d_p(nzb) = 0.0 ! !-- If necessary, apply the time filter IF ( asselin_filter_factor /= 0.0 .AND. & timestep_scheme(1:5) /= 'runge' ) THEN u1d = u1d + asselin_filter_factor * ( u1d_p - 2.0 * u1d + u1d_m ) v1d = v1d + asselin_filter_factor * ( v1d_p - 2.0 * v1d + v1d_m ) IF ( .NOT. constant_diffusion ) THEN e1d = e1d + asselin_filter_factor * & ( e1d_p - 2.0 * e1d + e1d_m ) ENDIF ENDIF ! !-- Swap the time levels in preparation for the next time step. IF ( timestep_scheme(1:4) == 'leap' ) THEN u1d_m = u1d v1d_m = v1d IF ( .NOT. constant_diffusion ) THEN e1d_m = e1d kh1d_m = kh1d ! The old diffusion quantities are required for km1d_m = km1d ! explicit diffusion in the leap-frog scheme. l1d_m = l1d IF ( prandtl_layer ) THEN usws1d_m = usws1d vsws1d_m = vsws1d ENDIF ENDIF ENDIF u1d = u1d_p v1d = v1d_p IF ( .NOT. constant_diffusion ) THEN e1d = e1d_p ENDIF ! !-- Compute diffusion quantities IF ( .NOT. constant_diffusion ) THEN ! !-- First compute the vertical fluxes in the Prandtl-layer IF ( prandtl_layer ) THEN ! !-- Compute theta* using Rif numbers of the previous time step IF ( rif1d(1) >= 0.0 ) THEN ! !-- Stable stratification ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * & ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & ) ELSE ! !-- Unstable stratification a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) ) b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0h1d ) ! !-- In the borderline case the formula for stable stratification !-- must be applied, because otherwise a zero division would !-- occur in the argument of the logarithm. IF ( a == 0.0 .OR. b == 0.0 ) THEN ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * & ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & ) ELSE ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & LOG( (a-1.0) / (a+1.0) * (b+1.0) / (b-1.0) ) ENDIF ENDIF ENDIF ! prandtl_layer ! !-- Compute the Richardson-flux numbers, !-- first at the top of the Prandtl-layer using u* of the previous !-- time step (+1E-30, if u* = 0), then in the remaining area. There !-- the rif-numbers of the previous time step are used. IF ( prandtl_layer ) THEN IF ( .NOT. humidity ) THEN pt_0 = pt_init(nzb+1) flux = ts1d ELSE pt_0 = pt_init(nzb+1) * ( 1.0 + 0.61 * q_init(nzb+1) ) flux = ts1d + 0.61 * pt_init(k) * qs1d ENDIF rif1d(nzb+1) = zu(nzb+1) * kappa * g * flux / & ( pt_0 * ( us1d**2 + 1E-30 ) ) ENDIF DO k = nzb_diff, nzt IF ( .NOT. humidity ) THEN pt_0 = pt_init(k) flux = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k) ELSE pt_0 = pt_init(k) * ( 1.0 + 0.61 * q_init(k) ) flux = ( ( pt_init(k+1) - pt_init(k-1) ) & + 0.61 * pt_init(k) * ( q_init(k+1) - q_init(k-1) )& ) * dd2zu(k) ENDIF IF ( rif1d(k) >= 0.0 ) THEN rif1d(k) = g / pt_0 * flux / & ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & + 1E-30 & ) ELSE rif1d(k) = g / pt_0 * flux / & ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & + 1E-30 & ) * ( 1.0 - 16.0 * rif1d(k) )**0.25 ENDIF ENDDO ! !-- Richardson-numbers must remain restricted to a realistic value !-- 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 ! !-- Compute u* from the absolute velocity value IF ( prandtl_layer ) THEN uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 ) IF ( rif1d(nzb+1) >= 0.0 ) THEN ! !-- Stable stratification us1d = kappa * uv_total / ( & LOG( zu(nzb+1) / z01d ) + 5.0 * rif1d(nzb+1) * & ( zu(nzb+1) - z01d ) / zu(nzb+1) & ) ELSE ! !-- Unstable stratification a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif1d(nzb+1) ) ) b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) & * z01d ) ) ! !-- In the borderline case the formula for stable stratification !-- must be applied, because otherwise a zero division would !-- occur in the argument of the logarithm. IF ( a == 1.0 .OR. b == 1.0 ) THEN us1d = kappa * uv_total / ( & LOG( zu(nzb+1) / z01d ) + & 5.0 * rif1d(nzb+1) * ( zu(nzb+1) - z01d ) / & zu(nzb+1) ) ELSE us1d = kappa * uv_total / ( & LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) +& 2.0 * ( ATAN( b ) - ATAN( a ) ) & ) ENDIF ENDIF ! !-- Compute the momentum fluxes for the diffusion terms usws1d = - u1d(nzb+1) / uv_total * us1d**2 vsws1d = - v1d(nzb+1) / uv_total * us1d**2 ! !-- Boundary condition for the turbulent kinetic energy at the top !-- of the Prandtl-layer. c_m = 0.4 according to Detering. !-- Additional Neumann condition de/dz = 0 at nzb is set to ensure !-- compatibility with the 3D model. IF ( ibc_e_b == 2 ) THEN e1d(nzb+1) = ( us1d / 0.1 )**2 ! e1d(nzb+1) = ( us1d / 0.4 )**2 !not used so far, see also !prandtl_fluxes ENDIF e1d(nzb) = e1d(nzb+1) IF ( humidity .OR. passive_scalar ) THEN ! !-- Compute q* IF ( rif1d(1) >= 0.0 ) THEN ! !-- Stable stratification qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * & ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & ) ELSE ! !-- Unstable stratification a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) ) b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0h1d ) ! !-- In the borderline case the formula for stable stratification !-- must be applied, because otherwise a zero division would !-- occur in the argument of the logarithm. IF ( a == 1.0 .OR. b == 1.0 ) THEN qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * & ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & ) ELSE qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & LOG( (a-1.0) / (a+1.0) * (b+1.0) / (b-1.0) ) ENDIF ENDIF ELSE qs1d = 0.0 ENDIF ENDIF ! prandtl_layer ! !-- Compute the diabatic mixing length IF ( mixing_length_1d == 'blackadar' ) THEN DO k = nzb+1, nzt IF ( rif1d(k) >= 0.0 ) THEN l1d(k) = l_black(k) / ( 1.0 + 5.0 * rif1d(k) ) ELSE l1d(k) = l_black(k) * ( 1.0 - 16.0 * rif1d(k) )**0.25 ENDIF l1d(k) = l_black(k) ENDDO ELSEIF ( mixing_length_1d == 'as_in_3d_model' ) THEN DO k = nzb+1, nzt dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k) IF ( dpt_dz > 0.0 ) THEN l_stable = 0.76 * SQRT( e1d(k) ) / & SQRT( g / pt_init(k) * dpt_dz ) + 1E-5 ELSE l_stable = l_grid(k) ENDIF l1d(k) = MIN( l_grid(k), l_stable ) ENDDO ENDIF ! !-- Adjust mixing length to the prandtl mixing length IF ( adjust_mixing_length .AND. prandtl_layer ) THEN k = nzb+1 IF ( rif1d(k) >= 0.0 ) THEN l1d(k) = MIN( l1d(k), kappa * zu(k) / ( 1.0 + 5.0 * & rif1d(k) ) ) ELSE l1d(k) = MIN( l1d(k), kappa * zu(k) * & SQRT( SQRT( 1.0 - 16.0 * rif1d(k) ) ) ) ENDIF ENDIF ! !-- Compute the diffusion coefficients for momentum via the !-- corresponding Prandtl-layer relationship and according to !-- Prandtl-Kolmogorov, respectively. The unstable stratification is !-- 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 ( rif1d(nzb+1) >= 0.0 ) THEN km1d(nzb+1) = us1d * kappa * zu(nzb+1) / & ( 1.0 + 5.0 * rif1d(nzb+1) ) ELSE km1d(nzb+1) = us1d * kappa * zu(nzb+1) * & ( 1.0 - 16.0 * rif1d(nzb+1) )**0.25 ENDIF ENDIF DO k = nzb_diff, nzt ! km1d(k) = 0.4 * SQRT( e1d(k) ) !changed: adjustment to 3D-model km1d(k) = 0.1 * SQRT( e1d(k) ) IF ( rif1d(k) >= 0.0 ) THEN km1d(k) = km1d(k) * l1d(k) ELSE km1d(k) = km1d(k) * l_black(k) ENDIF ENDDO ! !-- Add damping layer DO k = damp_level_ind_1d+1, nzt+1 km1d(k) = 1.1 * km1d(k-1) km1d(k) = MIN( km1d(k), 10.0 ) ENDDO ! !-- Compute the diffusion coefficient for heat via the relationship !-- kh = phim / phih * km DO k = nzb+1, nzt IF ( rif1d(k) >= 0.0 ) THEN kh1d(k) = km1d(k) ELSE kh1d(k) = km1d(k) * ( 1.0 - 16.0 * rif1d(k) )**0.25 ENDIF ENDDO ENDIF ! .NOT. constant_diffusion ! !-- The Runge-Kutta scheme needs the recent diffusion quantities IF ( timestep_scheme(1:5) == 'runge' ) THEN u1d_m = u1d v1d_m = v1d IF ( .NOT. constant_diffusion ) THEN e1d_m = e1d kh1d_m = kh1d km1d_m = km1d l1d_m = l1d IF ( prandtl_layer ) THEN usws1d_m = usws1d vsws1d_m = vsws1d ENDIF ENDIF ENDIF ENDDO ! intermediate step loop ! !-- Increment simulated time and output times current_timestep_number_1d = current_timestep_number_1d + 1 simulated_time_1d = simulated_time_1d + dt_1d simulated_time_chr = time_to_string( simulated_time_1d ) time_pr_1d = time_pr_1d + dt_1d time_run_control_1d = time_run_control_1d + dt_1d ! !-- Determine and print out quantities for run control IF ( time_run_control_1d >= dt_run_control_1d ) THEN CALL run_control_1d time_run_control_1d = time_run_control_1d - dt_run_control_1d ENDIF ! !-- Profile output on file IF ( time_pr_1d >= dt_pr_1d ) THEN CALL print_1d_model time_pr_1d = time_pr_1d - dt_pr_1d ENDIF ! !-- Determine size of next time step CALL timestep_1d ENDDO ! time loop END SUBROUTINE time_integration_1d SUBROUTINE run_control_1d !------------------------------------------------------------------------------! ! Description: ! ------------ ! Compute and print out quantities for run control of the 1D model. !------------------------------------------------------------------------------! USE constants USE indices USE model_1d USE pegrid USE control_parameters IMPLICIT NONE INTEGER :: k REAL :: alpha, energy, umax, uv_total, vmax ! !-- Output IF ( myid == 0 ) THEN ! !-- If necessary, write header IF ( .NOT. run_control_header_1d ) THEN CALL check_open( 15 ) WRITE ( 15, 100 ) run_control_header_1d = .TRUE. ENDIF ! !-- Compute control quantities !-- grid level nzp is excluded due to mirror boundary condition umax = 0.0; vmax = 0.0; energy = 0.0 DO k = nzb+1, nzt+1 umax = MAX( ABS( umax ), ABS( u1d(k) ) ) vmax = MAX( ABS( vmax ), ABS( v1d(k) ) ) energy = energy + 0.5 * ( u1d(k)**2 + v1d(k)**2 ) ENDDO energy = energy / REAL( nzt - nzb + 1 ) uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 ) IF ( ABS( v1d(nzb+1) ) .LT. 1.0E-5 ) THEN alpha = ACOS( SIGN( 1.0 , u1d(nzb+1) ) ) ELSE alpha = ACOS( u1d(nzb+1) / uv_total ) IF ( v1d(nzb+1) <= 0.0 ) alpha = 2.0 * pi - alpha ENDIF alpha = alpha / ( 2.0 * pi ) * 360.0 WRITE ( 15, 101 ) current_timestep_number_1d, simulated_time_chr, & dt_1d, umax, vmax, us1d, alpha, energy ! !-- Write buffer contents to disc immediately CALL local_flush( 15 ) ENDIF ! !-- formats 100 FORMAT (///'1D-Zeitschrittkontrollausgaben:'/ & &'------------------------------'// & &'ITER. HH:MM:SS DT UMAX VMAX U* ALPHA ENERG.'/ & &'-------------------------------------------------------------') 101 FORMAT (I5,2X,A9,1X,F6.2,2X,F6.2,1X,F6.2,2X,F5.3,2X,F5.1,2X,F7.2) END SUBROUTINE run_control_1d SUBROUTINE timestep_1d !------------------------------------------------------------------------------! ! Description: ! ------------ ! Compute the time step w.r.t. the diffusion criterion !------------------------------------------------------------------------------! USE arrays_3d USE indices USE model_1d USE pegrid USE control_parameters IMPLICIT NONE INTEGER :: k REAL :: div, dt_diff, fac, percent_change, value ! !-- Compute the currently feasible time step according to the diffusion !-- criterion. At nzb+1 the half grid length is used. IF ( timestep_scheme(1:4) == 'leap' ) THEN fac = 0.25 ELSE fac = 0.35 ENDIF dt_diff = dt_max_1d DO k = nzb+2, nzt value = fac * dzu(k) * dzu(k) / ( km1d(k) + 1E-20 ) dt_diff = MIN( value, dt_diff ) ENDDO value = fac * zu(nzb+1) * zu(nzb+1) / ( km1d(nzb+1) + 1E-20 ) dt_1d = MIN( value, dt_diff ) ! !-- Set flag when the time step becomes too small IF ( dt_1d < ( 0.00001 * dt_max_1d ) ) THEN stop_dt_1d = .TRUE. WRITE( message_string, * ) 'timestep has exceeded the lower limit &', & 'dt_1d = ',dt_1d,' s simulation stopped!' CALL message( 'timestep_1d', 'PA0192', 1, 2, 0, 6, 0 ) ENDIF IF ( timestep_scheme(1:4) == 'leap' ) THEN ! !-- The current time step will only be changed if the new time step exceeds !-- its previous value by 5 % or falls 2 % below. After a time step !-- reduction at least 30 iterations must be done with this value before a !-- new reduction will be allowed again. !-- The control parameters for application of Euler- or leap-frog schemes are !-- set accordingly. percent_change = dt_1d / old_dt_1d - 1.0 IF ( percent_change > 0.05 .OR. percent_change < -0.02 ) THEN ! !-- Each time step increase is by at most 2 % IF ( percent_change > 0.0 .AND. simulated_time_1d /= 0.0 ) THEN dt_1d = 1.02 * old_dt_1d ENDIF ! !-- A more or less simple new time step value is obtained taking only the !-- first two significant digits div = 1000.0 DO WHILE ( dt_1d < div ) div = div / 10.0 ENDDO dt_1d = NINT( dt_1d * 100.0 / div ) * div / 100.0 ! !-- Now the time step can be changed. IF ( percent_change < 0.0 ) THEN ! !-- Time step reduction old_dt_1d = dt_1d last_dt_change_1d = current_timestep_number_1d ELSE ! !-- Time step will only be increased if at least 30 iterations have !-- been done since the previous time step change, and of course at !-- simulation start, respectively. IF ( current_timestep_number_1d >= last_dt_change_1d + 30 .OR. & simulated_time_1d == 0.0 ) THEN old_dt_1d = dt_1d last_dt_change_1d = current_timestep_number_1d ELSE dt_1d = old_dt_1d ENDIF ENDIF ELSE ! !-- No time step change since the difference is too small dt_1d = old_dt_1d ENDIF ELSE ! Runge-Kutta !-- A more or less simple new time step value is obtained taking only the !-- first two significant digits div = 1000.0 DO WHILE ( dt_1d < div ) div = div / 10.0 ENDDO dt_1d = NINT( dt_1d * 100.0 / div ) * div / 100.0 old_dt_1d = dt_1d last_dt_change_1d = current_timestep_number_1d ENDIF END SUBROUTINE timestep_1d SUBROUTINE print_1d_model !------------------------------------------------------------------------------! ! Description: ! ------------ ! List output of profiles from the 1D-model !------------------------------------------------------------------------------! USE arrays_3d USE indices USE model_1d USE pegrid USE control_parameters IMPLICIT NONE INTEGER :: k IF ( myid == 0 ) THEN ! !-- Open list output file for profiles from the 1D-model CALL check_open( 17 ) ! !-- Write Header WRITE ( 17, 100 ) TRIM( run_description_header ), & TRIM( simulated_time_chr ) WRITE ( 17, 101 ) ! !-- Write the values WRITE ( 17, 102 ) WRITE ( 17, 101 ) DO k = nzt+1, nzb, -1 WRITE ( 17, 103) k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), & rif1d(k), km1d(k), kh1d(k), l1d(k), zu(k), k ENDDO WRITE ( 17, 101 ) WRITE ( 17, 102 ) WRITE ( 17, 101 ) ! !-- Write buffer contents to disc immediately CALL local_flush( 17 ) ENDIF ! !-- 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 ', & '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, & 1X,F5.2,1X,F6.2,1X,F7.1,2X,I4) END SUBROUTINE print_1d_model