source: palm/trunk/SOURCE/init_3d_model.f90 @ 80

Last change on this file since 80 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:


General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:


Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

  • Property svn:keywords set to Id
File size: 33.4 KB
RevLine 
[1]1#if defined( __ibmy_special )
2@PROCESS NOOPTimize
3#endif
4 SUBROUTINE init_3d_model
5
6!------------------------------------------------------------------------------!
7! Actual revisions:
8! -----------------
[77]9!
10!
11! Former revisions:
12! -----------------
13! $Id: init_3d_model.f90 77 2007-03-29 04:26:56Z raasch $
14!
15! 75 2007-03-22 09:54:05Z raasch
[73]16! Arrays for radiation boundary conditions are allocated (u_m_l, u_m_r, etc.),
17! bugfix for cases with the outflow damping layer extending over more than one
[75]18! subdomain, moisture renamed humidity,
19! new initializing action "by_user" calls user_init_3d_model,
[72]20! precipitation_amount/rate, ts_value are allocated, +module netcdf_control,
[51]21! initial velocities at nzb+1 are regarded for volume
22! flow control in case they have been set zero before (to avoid small timesteps)
[75]23! -uvmean_outflow, uxrp, vynp eliminated
[1]24!
[39]25! 19 2007-02-23 04:53:48Z raasch
26! +handling of top fluxes
27!
[3]28! RCS Log replace by Id keyword, revision history cleaned up
29!
[1]30! Revision 1.49  2006/08/22 15:59:07  raasch
31! No optimization of this file on the ibmy (Yonsei Univ.)
32!
33! Revision 1.1  1998/03/09 16:22:22  raasch
34! Initial revision
35!
36!
37! Description:
38! ------------
39! Allocation of arrays and initialization of the 3D model via
40! a) pre-run the 1D model
41! or
42! b) pre-set constant linear profiles
43! or
44! c) read values of a previous run
45!------------------------------------------------------------------------------!
46
47    USE arrays_3d
48    USE averaging
[72]49    USE cloud_parameters
[1]50    USE constants
51    USE control_parameters
52    USE cpulog
53    USE indices
54    USE interfaces
55    USE model_1d
[51]56    USE netcdf_control
[1]57    USE particle_attributes
58    USE pegrid
59    USE profil_parameter
60    USE random_function_mod
61    USE statistics
62
63    IMPLICIT NONE
64
65    INTEGER ::  i, j, k, sr
66
67    INTEGER, DIMENSION(:), ALLOCATABLE ::  ngp_2dh_l, ngp_3d_inner_l
68
69    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l
70
71    REAL, DIMENSION(1:2) ::  volume_flow_area_l, volume_flow_initial_l
72
73
74!
75!-- Allocate arrays
76    ALLOCATE( ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
77              ngp_3d(0:statistic_regions),                                  &
78              ngp_3d_inner(0:statistic_regions),                            &
79              ngp_3d_inner_l(0:statistic_regions),                          &
80              sums_divnew_l(0:statistic_regions),                           &
81              sums_divold_l(0:statistic_regions) )
[75]82    ALLOCATE( rdf(nzb+1:nzt) )
[1]83    ALLOCATE( hom_sum(nzb:nzt+1,var_hom,0:statistic_regions),               &
84              ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                 &
85              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),               &
86              rmask(nys-1:nyn+1,nxl-1:nxr+1,0:statistic_regions),           &
87              sums(nzb:nzt+1,var_sum),                                      &
88              sums_l(nzb:nzt+1,var_sum,0:threads_per_task-1),               &
89              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
90              sums_up_fraction_l(10,3,0:statistic_regions),                 &
[48]91              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                &
92              ts_value(var_ts,0:statistic_regions) )
[1]93    ALLOCATE( km_damp_x(nxl-1:nxr+1), km_damp_y(nys-1:nyn+1) )
94
[19]95    ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), &
96              ts(nys-1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1),  &
97              us(nys-1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1),   &
98              vsws_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) )
[1]99
100    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
101!
102!--    Leapfrog scheme needs two timelevels of diffusion quantities
[19]103       ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1),   &
104                 shf_2(nys-1:nyn+1,nxl-1:nxr+1),   &
105                 tswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
106                 usws_2(nys-1:nyn+1,nxl-1:nxr+1),  &
[1]107                 vsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
108    ENDIF
109
[75]110    ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra),         &
111              e_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
112              e_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
113              e_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
114              kh_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
115              km_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
116              p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),    &
117              pt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
118              pt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
119              pt_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
120              tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
121              u_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
122              u_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
123              u_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
124              v_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
125              v_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
126              v_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
127              w_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
128              w_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
[1]129              w_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
130
131    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
132       ALLOCATE( kh_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
133                 km_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
134    ENDIF
135
[75]136    IF ( humidity  .OR.  passive_scalar ) THEN
[1]137!
[75]138!--    2D-humidity/scalar arrays
[1]139       ALLOCATE ( qs(nys-1:nyn+1,nxl-1:nxr+1),     &
[19]140                  qsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
141                  qswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
[1]142
143       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
[19]144          ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1), &
145                    qswst_2(nys-1:nyn+1,nxl-1:nxr+1) )
[1]146       ENDIF
147!
[75]148!--    3D-humidity/scalar arrays
[1]149       ALLOCATE( q_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
150                 q_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
151                 q_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
152
153!
[75]154!--    3D-arrays needed for humidity only
155       IF ( humidity )  THEN
[1]156          ALLOCATE( vpt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
157
158          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
159             ALLOCATE( vpt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
160          ENDIF
161
162          IF ( cloud_physics ) THEN
163!
164!--          Liquid water content
165             ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
[72]166!
167!--          Precipitation amount and rate (only needed if output is switched)
168             ALLOCATE( precipitation_amount(nys-1:nyn+1,nxl-1:nxr+1), &
169                       precipitation_rate(nys-1:nyn+1,nxl-1:nxr+1) )
[1]170          ENDIF
171
172          IF ( cloud_droplets )  THEN
173!
174!--          Liquid water content, change in liquid water content,
175!--          real volume of particles (with weighting), volume of particles
176             ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
177                        ql_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
178                        ql_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
179                        ql_vp(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
180          ENDIF
181
182       ENDIF
183
184    ENDIF
185
186!
187!-- 3D-array for storing the dissipation, needed for calculating the sgs
188!-- particle velocities
189    IF ( use_sgs_for_particles )  THEN
190       ALLOCATE ( diss(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
191    ENDIF
192
193    IF ( dt_dosp /= 9999999.9 )  THEN
194       ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &
195                 spectrum_y( 1:ny/2, 1:10, 1:10 ) )
196    ENDIF
197
198!
[51]199!-- 4D-array for storing the Rif-values at vertical walls
200    IF ( topography /= 'flat' )  THEN
201       ALLOCATE( rif_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1,1:4) )
202       rif_wall = 0.0
203    ENDIF
204
205!
206!-- Velocities at nzb+1 needed for volume flow control
207    IF ( conserve_volume_flow )  THEN
208       ALLOCATE( u_nzb_p1_for_vfc(nys:nyn), v_nzb_p1_for_vfc(nxl:nxr) )
209       u_nzb_p1_for_vfc = 0.0
210       v_nzb_p1_for_vfc = 0.0
211    ENDIF
212
213!
[73]214!-- Arrays to store velocity data from t-dt needed for radiation boundary
215!-- conditions
216    IF ( outflow_l )  THEN
[75]217       ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1), &
218                 v_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1), &
219                 w_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1) )
[73]220    ENDIF
221    IF ( outflow_r )  THEN
[75]222       ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1), &
223                 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1), &
224                 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1) )
[73]225    ENDIF
226    IF ( outflow_s )  THEN
[75]227       ALLOCATE( u_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1), &
228                 v_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1), &
229                 w_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1) )
[73]230    ENDIF
231    IF ( outflow_n )  THEN
[75]232       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1), &
233                 v_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1), &
234                 w_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1) )
[73]235    ENDIF
236
237!
[1]238!-- Initial assignment of the pointers
239    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
240
[19]241       rif_m   => rif_1;    rif   => rif_2
242       shf_m   => shf_1;    shf   => shf_2
243       tswst_m => tswst_1;  tswst => tswst_2
244       usws_m  => usws_1;   usws  => usws_2
245       vsws_m  => vsws_1;   vsws  => vsws_2
[1]246       e_m  => e_1;   e  => e_2;   e_p  => e_3;   te_m  => e_3
247       kh_m => kh_1;  kh => kh_2
248       km_m => km_1;  km => km_2
249       pt_m => pt_1;  pt => pt_2;  pt_p => pt_3;  tpt_m => pt_3
250       u_m  => u_1;   u  => u_2;   u_p  => u_3;   tu_m  => u_3
251       v_m  => v_1;   v  => v_2;   v_p  => v_3;   tv_m  => v_3
252       w_m  => w_1;   w  => w_2;   w_p  => w_3;   tw_m  => w_3
253
[75]254       IF ( humidity  .OR.  passive_scalar )  THEN
[19]255          qsws_m  => qsws_1;   qsws  => qsws_2
256          qswst_m => qswst_1;  qswst => qswst_2
[1]257          q_m    => q_1;     q    => q_2;     q_p => q_3;     tq_m => q_3
[75]258          IF ( humidity )        vpt_m  => vpt_1;   vpt  => vpt_2
[1]259          IF ( cloud_physics )   ql   => ql_1
260          IF ( cloud_droplets )  THEN
261             ql   => ql_1
262             ql_c => ql_2
263          ENDIF
264       ENDIF
265
266    ELSE
267
[19]268       rif   => rif_1
269       shf   => shf_1
270       tswst => tswst_1
271       usws  => usws_1
272       vsws  => vsws_1
273       e     => e_1;   e_p  => e_2;   te_m  => e_3;   e_m  => e_3
274       kh    => kh_1
275       km    => km_1
276       pt    => pt_1;  pt_p => pt_2;  tpt_m => pt_3;  pt_m => pt_3
277       u     => u_1;   u_p  => u_2;   tu_m  => u_3;   u_m  => u_3
278       v     => v_1;   v_p  => v_2;   tv_m  => v_3;   v_m  => v_3
279       w     => w_1;   w_p  => w_2;   tw_m  => w_3;   w_m  => w_3
[1]280
[75]281       IF ( humidity  .OR.  passive_scalar )  THEN
[1]282          qsws   => qsws_1
[19]283          qswst  => qswst_1
[1]284          q      => q_1;     q_p  => q_2;     tq_m => q_3;    q_m => q_3
[75]285          IF ( humidity )        vpt  => vpt_1
[1]286          IF ( cloud_physics )   ql   => ql_1
287          IF ( cloud_droplets )  THEN
288             ql   => ql_1
289             ql_c => ql_2
290          ENDIF
291       ENDIF
292
293    ENDIF
294
295!
296!-- Initialize model variables
297    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
298!
299!--    First model run of a possible job queue.
300!--    Initial profiles of the variables must be computes.
301       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
302!
303!--       Use solutions of the 1D model as initial profiles,
304!--       start 1D model
305          CALL init_1d_model
306!
307!--       Transfer initial profiles to the arrays of the 3D model
308          DO  i = nxl-1, nxr+1
309             DO  j = nys-1, nyn+1
310                e(:,j,i)  = e1d
311                kh(:,j,i) = kh1d
312                km(:,j,i) = km1d
313                pt(:,j,i) = pt_init
314                u(:,j,i)  = u1d
315                v(:,j,i)  = v1d
316             ENDDO
317          ENDDO
318
[75]319          IF ( humidity  .OR.  passive_scalar )  THEN
[1]320             DO  i = nxl-1, nxr+1
321                DO  j = nys-1, nyn+1
322                   q(:,j,i) = q_init
323                ENDDO
324             ENDDO
325          ENDIF
326
327          IF ( .NOT. constant_diffusion )  THEN
328             DO  i = nxl-1, nxr+1
329                DO  j = nys-1, nyn+1
330                   e(:,j,i)  = e1d
331                ENDDO
332             ENDDO
333!
334!--          Store initial profiles for output purposes etc.
335             hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
336
337             IF ( prandtl_layer )  THEN
338                rif  = rif1d(nzb+1)
339                ts   = 0.0  ! could actually be computed more accurately in the
340                            ! 1D model. Update when opportunity arises.
341                us   = us1d
342                usws = usws1d
343                vsws = vsws1d
344             ELSE
345                ts   = 0.0  ! must be set, because used in
346                rif  = 0.0  ! flowste
347                us   = 0.0
348                usws = 0.0
349                vsws = 0.0
350             ENDIF
351
352          ELSE
353             e    = 0.0  ! must be set, because used in
354             rif  = 0.0  ! flowste
355             ts   = 0.0
356             us   = 0.0
357             usws = 0.0
358             vsws = 0.0
359          ENDIF
360
361!
362!--       In every case qs = 0.0 (see also pt)
363!--       This could actually be computed more accurately in the 1D model.
364!--       Update when opportunity arises!
[75]365          IF ( humidity  .OR.  passive_scalar )  qs = 0.0
[1]366
367!
368!--       inside buildings set velocities back to zero
369          IF ( topography /= 'flat' )  THEN
370             DO  i = nxl-1, nxr+1
371                DO  j = nys-1, nyn+1
372                   u(nzb:nzb_u_inner(j,i),j,i) = 0.0
373                   v(nzb:nzb_v_inner(j,i),j,i) = 0.0
374                ENDDO
375             ENDDO
376!
377!--          WARNING: The extra boundary conditions set after running the
378!--          -------  1D model impose an error on the divergence one layer
379!--                   below the topography; need to correct later
380!--          ATTENTION: Provisional correction for Piacsek & Williams
381!--          ---------  advection scheme: keep u and v zero one layer below
382!--                     the topography.
383             IF ( ibc_uv_b == 0 )  THEN
384!
385!--             Satisfying the Dirichlet condition with an extra layer below
386!--             the surface where the u and v component change their sign.
387                DO  i = nxl-1, nxr+1
388                   DO  j = nys-1, nyn+1
389                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)
390                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)
391                   ENDDO
392                ENDDO
393
394             ELSE
395!
396!--             Neumann condition
397                DO  i = nxl-1, nxr+1
398                   DO  j = nys-1, nyn+1
399                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = u(1,j,i)
400                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = v(1,j,i)
401                   ENDDO
402                ENDDO
403
404             ENDIF
405
406          ENDIF
407
408       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
409       THEN
410!
411!--       Use constructed initial profiles (velocity constant with height,
412!--       temperature profile with constant gradient)
413          DO  i = nxl-1, nxr+1
414             DO  j = nys-1, nyn+1
415                pt(:,j,i) = pt_init
416                u(:,j,i)  = u_init
417                v(:,j,i)  = v_init
418             ENDDO
419          ENDDO
[75]420
[1]421!
[51]422!--       Set initial horizontal velocities at the lowest computational grid levels
423!--       to zero in order to avoid too small time steps caused by the diffusion
[1]424!--       limit in the initial phase of a run (at k=1, dz/2 occurs in the
[51]425!--       limiting formula!). The original values are stored to be later used for
426!--       volume flow control.
[1]427          DO  i = nxl-1, nxr+1
428             DO  j = nys-1, nyn+1
429                u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0
430                v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0
431             ENDDO
432          ENDDO
[51]433          IF ( conserve_volume_flow )  THEN
434             IF ( nxr == nx )  THEN
435                DO  j = nys, nyn
436                   k = nzb_u_inner(j,nx) + 1
437                   u_nzb_p1_for_vfc(j) = u_init(k) * dzu(k)
438                ENDDO
439             ENDIF
440             IF ( nyn == ny )  THEN
441                DO  i = nxl, nxr
442                   k = nzb_v_inner(ny,i) + 1
443                   v_nzb_p1_for_vfc(i) = v_init(k) * dzu(k)
444                ENDDO
445             ENDIF
446          ENDIF
[1]447
[75]448          IF ( humidity  .OR.  passive_scalar )  THEN
[1]449             DO  i = nxl-1, nxr+1
450                DO  j = nys-1, nyn+1
451                   q(:,j,i) = q_init
452                ENDDO
453             ENDDO
454          ENDIF
455
456         
457          IF ( constant_diffusion )  THEN
458             km   = km_constant
459             kh   = km / prandtl_number
460          ELSE
461             kh   = 0.01   ! there must exist an initial diffusion, because
462             km   = 0.01   ! otherwise no TKE would be produced by the
463                           ! production terms, as long as not yet
464                           ! e = (u*/cm)**2 at k=nzb+1
465          ENDIF
466          e    = 0.0
467          rif  = 0.0
468          ts   = 0.0
469          us   = 0.0
470          usws = 0.0
471          vsws = 0.0
[75]472          IF ( humidity  .OR.  passive_scalar ) qs = 0.0
[1]473
474!
475!--       Compute initial temperature field and other constants used in case
476!--       of a sloping surface
477          IF ( sloping_surface )  CALL init_slope
478
[46]479       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
480       THEN
481!
482!--       Initialization will completely be done by the user
483          CALL user_init_3d_model
484
[1]485       ENDIF
486
487!
488!--    Calculate virtual potential temperature
[75]489       IF ( humidity ) vpt = pt * ( 1.0 + 0.61 * q )
[1]490
491!
492!--    Store initial profiles for output purposes etc.
493       hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
494       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
495       IF ( ibc_uv_b == 0 )  THEN
496          hom(nzb,1,5,:) = -hom(nzb+1,1,5,:)  ! due to satisfying the Dirichlet
497          hom(nzb,1,6,:) = -hom(nzb+1,1,6,:)  ! condition with an extra layer
498              ! below the surface where the u and v component change their sign
499       ENDIF
500       hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
501       hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
502       hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
503
504
[75]505       IF ( humidity )  THEN
[1]506!
507!--       Store initial profile of total water content, virtual potential
508!--       temperature
509          hom(:,1,26,:) = SPREAD(   q(:,nys,nxl), 2, statistic_regions+1 )
510          hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
511          IF ( cloud_physics  .OR.  cloud_droplets ) THEN
512!
513!--          Store initial profile of specific humidity and potential
514!--          temperature
515             hom(:,1,27,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
516             hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
517          ENDIF
518       ENDIF
519
520       IF ( passive_scalar )  THEN
521!
522!--       Store initial scalar profile
523          hom(:,1,26,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
524       ENDIF
525
526!
[19]527!--    Initialize fluxes at bottom surface
[1]528       IF ( use_surface_fluxes )  THEN
529
530          IF ( constant_heatflux )  THEN
531!
532!--          Heat flux is prescribed
533             IF ( random_heatflux )  THEN
534                CALL disturb_heatflux
535             ELSE
536                shf = surface_heatflux
537!
538!--             Over topography surface_heatflux is replaced by wall_heatflux(0)
539                IF ( TRIM( topography ) /= 'flat' )  THEN
540                   DO  i = nxl-1, nxr+1
541                      DO  j = nys-1, nyn+1
542                         IF ( nzb_s_inner(j,i) /= 0 )  THEN
543                            shf(j,i) = wall_heatflux(0)
544                         ENDIF
545                      ENDDO
546                   ENDDO
547                ENDIF
548             ENDIF
549             IF ( ASSOCIATED( shf_m ) )  shf_m = shf
550          ENDIF
551
552!
553!--       Determine the near-surface water flux
[75]554          IF ( humidity  .OR.  passive_scalar )  THEN
[1]555             IF ( constant_waterflux )  THEN
556                qsws   = surface_waterflux
557                IF ( ASSOCIATED( qsws_m ) )  qsws_m = qsws
558             ENDIF
559          ENDIF
560
561       ENDIF
562
563!
[19]564!--    Initialize fluxes at top surface
565!--    Currently, only the heatflux can be prescribed. The latent flux is
[40]566!--    zero in this case!
[19]567       IF ( use_top_fluxes )  THEN
568
569          IF ( constant_top_heatflux )  THEN
570!
571!--          Heat flux is prescribed
572             tswst = top_heatflux
573             IF ( ASSOCIATED( tswst_m ) )  tswst_m = tswst
574
[75]575             IF ( humidity  .OR.  passive_scalar )  THEN
[19]576                qswst = 0.0
577                IF ( ASSOCIATED( qswst_m ) )  qswst_m = qswst
578             ENDIF
579         ENDIF
580
581       ENDIF
582
583!
[1]584!--    Initialize Prandtl layer quantities
585       IF ( prandtl_layer )  THEN
586
587          z0 = roughness_length
588
589          IF ( .NOT. constant_heatflux )  THEN 
590!
591!--          Surface temperature is prescribed. Here the heat flux cannot be
592!--          simply estimated, because therefore rif, u* and theta* would have
593!--          to be computed by iteration. This is why the heat flux is assumed
594!--          to be zero before the first time step. It approaches its correct
595!--          value in the course of the first few time steps.
596             shf   = 0.0
597             IF ( ASSOCIATED( shf_m ) )  shf_m = 0.0
598          ENDIF
599
[75]600          IF ( humidity  .OR.  passive_scalar )  THEN
[1]601             IF ( .NOT. constant_waterflux )  THEN
602                qsws   = 0.0
603                IF ( ASSOCIATED( qsws_m ) )   qsws_m = 0.0
604             ENDIF
605          ENDIF
606
607       ENDIF
608
609!
610!--    Calculate the initial volume flow at the right and north boundary
611       IF ( conserve_volume_flow )  THEN
612
613          volume_flow_initial_l = 0.0
614          volume_flow_area_l    = 0.0
615 
616          IF ( nxr == nx )  THEN
617             DO  j = nys, nyn
618                DO  k = nzb_2d(j,nx) + 1, nzt
619                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
620                                              u(k,j,nx) * dzu(k)
621                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzu(k)
622                ENDDO
[51]623!
624!--             Correction if velocity at nzb+1 has been set zero further above
625                volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
626                                           u_nzb_p1_for_vfc(j)
[1]627             ENDDO
628          ENDIF
629
630          IF ( nyn == ny )  THEN
631             DO  i = nxl, nxr
632                DO  k = nzb_2d(ny,i) + 1, nzt 
633                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
634                                              v(k,ny,i) * dzu(k)
635                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzu(k)
636                ENDDO
[51]637!
638!--             Correction if velocity at nzb+1 has been set zero further above
639                volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
640                                           v_nzb_p1_for_vfc(i)
[1]641             ENDDO
642          ENDIF
643
644#if defined( __parallel )
645          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
646                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
647          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
648                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
649#else
650          volume_flow_initial = volume_flow_initial_l
651          volume_flow_area    = volume_flow_area_l
652#endif 
653       ENDIF
654
655!
656!--    For the moment, perturbation pressure and vertical velocity are zero
657       p = 0.0; w = 0.0
658
659!
660!--    Initialize array sums (must be defined in first call of pres)
661       sums = 0.0
662
663!
[72]664!--    Treating cloud physics, liquid water content and precipitation amount
665!--    are zero at beginning of the simulation
666       IF ( cloud_physics )  THEN
667          ql = 0.0
668          IF ( precipitation )  precipitation_amount = 0.0
669       ENDIF
[1]670
671!
672!--    Initialize spectra
673       IF ( dt_dosp /= 9999999.9 )  THEN
674          spectrum_x = 0.0
675          spectrum_y = 0.0
676       ENDIF
677
678!
679!--    Impose vortex with vertical axis on the initial velocity profile
680       IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 )  THEN
681          CALL init_rankine
682       ENDIF
683
684!
685!--    Impose temperature anomaly (advection test only)
686       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 )  THEN
687          CALL init_pt_anomaly
688       ENDIF
689
690!
691!--    If required, change the surface temperature at the start of the 3D run
692       IF ( pt_surface_initial_change /= 0.0 )  THEN
693          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
694       ENDIF
695
696!
697!--    If required, change the surface humidity/scalar at the start of the 3D
698!--    run
[75]699       IF ( ( humidity .OR. passive_scalar ) .AND. &
[1]700            q_surface_initial_change /= 0.0 )  THEN
701          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
702       ENDIF
703
704!
705!--    Initialize the random number generator (from numerical recipes)
706       CALL random_function_ini
707
708!
709!--    Impose random perturbation on the horizontal velocity field and then
710!--    remove the divergences from the velocity field
711       IF ( create_disturbances )  THEN
[75]712          CALL disturb_field( nzb_u_inner, tend, u )
713          CALL disturb_field( nzb_v_inner, tend, v )
[1]714          n_sor = nsor_ini
715          CALL pres
716          n_sor = nsor
717       ENDIF
718
719!
720!--    Once again set the perturbation pressure explicitly to zero in order to
721!--    assure that it does not generate any divergences in the first time step.
722!--    At t=0 the velocity field is free of divergence (as constructed above).
723!--    Divergences being created during a time step are not yet known and thus
724!--    cannot be corrected during the time step yet.
725       p = 0.0
726
727!
728!--    Initialize old and new time levels.
729       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
730          e_m = e; pt_m = pt; u_m = u; v_m = v; w_m = w; kh_m = kh; km_m = km
731       ELSE
732          te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
733       ENDIF
734       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
735
[75]736       IF ( humidity  .OR.  passive_scalar )  THEN
[1]737          IF ( ASSOCIATED( q_m ) )  q_m = q
738          IF ( timestep_scheme(1:5) == 'runge' )  tq_m = 0.0
739          q_p = q
[75]740          IF ( humidity  .AND.  ASSOCIATED( vpt_m ) )  vpt_m = vpt
[1]741       ENDIF
742
[73]743!
744!--    Initialize old timelevels needed for radiation boundary conditions
745       IF ( outflow_l )  THEN
[75]746          u_m_l(:,:,:) = u(:,:,-1:1)
747          v_m_l(:,:,:) = v(:,:,-1:1)
748          w_m_l(:,:,:) = w(:,:,-1:1)
[73]749       ENDIF
750       IF ( outflow_r )  THEN
[75]751          u_m_r(:,:,:) = u(:,:,nx-1:nx+1)
752          v_m_r(:,:,:) = v(:,:,nx-1:nx+1)
753          w_m_r(:,:,:) = w(:,:,nx-1:nx+1)
[73]754       ENDIF
755       IF ( outflow_s )  THEN
[75]756          u_m_s(:,:,:) = u(:,-1:1,:)
757          v_m_s(:,:,:) = v(:,-1:1,:)
758          w_m_s(:,:,:) = w(:,-1:1,:)
[73]759       ENDIF
760       IF ( outflow_n )  THEN
[75]761          u_m_n(:,:,:) = u(:,ny-1:ny+1,:)
762          v_m_n(:,:,:) = v(:,ny-1:ny+1,:)
763          w_m_n(:,:,:) = w(:,ny-1:ny+1,:)
[73]764       ENDIF
765
[1]766    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' ) &
767    THEN
768!
769!--    Read binary data from restart file
770       CALL read_3d_binary
771
772!
773!--    Calculate initial temperature field and other constants used in case
774!--    of a sloping surface
775       IF ( sloping_surface )  CALL init_slope
776
777!
778!--    Initialize new time levels (only done in order to set boundary values
779!--    including ghost points)
780       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
[75]781       IF ( humidity  .OR.  passive_scalar )  q_p = q
[1]782
783    ELSE
784!
785!--    Actually this part of the programm should not be reached
786       IF ( myid == 0 )  PRINT*,'+++ init_3d_model: unknown initializing ', &
787                                                    'problem'
788       CALL local_stop
789    ENDIF
790
791!
792!-- If required, initialize dvrp-software
793    IF ( dt_dvrp /= 9999999.9 )  CALL init_dvrp
794
795!
796!-- If required, initialize quantities for handling cloud physics
797!-- This routine must be called before init_particles, because
798!-- otherwise, array pt_d_t, needed in data_output_dvrp (called by
799!-- init_particles) is not defined.
800    CALL init_cloud_physics
801
802!
803!-- If required, initialize particles
[63]804    IF ( particle_advection )  CALL init_particles
[1]805
806!
807!-- Initialize quantities for special advections schemes
808    CALL init_advec
809
810!
811!-- Initialize Rayleigh damping factors
812    rdf = 0.0
813    IF ( rayleigh_damping_factor /= 0.0 )  THEN
814       DO  k = nzb+1, nzt
815          IF ( zu(k) >= rayleigh_damping_height )  THEN
816             rdf(k) = rayleigh_damping_factor * &
817                      ( SIN( pi * 0.5 * ( zu(k) - rayleigh_damping_height )    &
818                                      / ( zu(nzt) - rayleigh_damping_height ) )&
819                      )**2
820          ENDIF
821       ENDDO
822    ENDIF
823
824!
825!-- Initialize diffusivities used within the outflow damping layer in case of
826!-- non-cyclic lateral boundaries. A linear increase is assumed over the first
827!-- half of the width of the damping layer
[73]828    IF ( bc_lr == 'dirichlet/radiation' )  THEN
[1]829
830       DO  i = nxl-1, nxr+1
[73]831          IF ( i >= nx - outflow_damping_width )  THEN
832             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
833                            ( i - ( nx - outflow_damping_width ) ) /   &
834                            REAL( outflow_damping_width/2 )            &
835                                             )
836          ELSE
837             km_damp_x(i) = 0.0
838          ENDIF
839       ENDDO
[1]840
[73]841    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
[1]842
[73]843       DO  i = nxl-1, nxr+1
844          IF ( i <= outflow_damping_width )  THEN
845             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
846                            ( outflow_damping_width - i ) /            &
847                            REAL( outflow_damping_width/2 )            &
848                                             )
849          ELSE
850             km_damp_x(i) = 0.0
851          ENDIF
852       ENDDO
[1]853
[73]854    ENDIF
[1]855
[73]856    IF ( bc_ns == 'dirichlet/radiation' )  THEN
[1]857
[73]858       DO  j = nys-1, nyn+1
859          IF ( j >= ny - outflow_damping_width )  THEN
860             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
861                            ( j - ( ny - outflow_damping_width ) ) /   &
862                            REAL( outflow_damping_width/2 )            &
863                                             )
864          ELSE
865             km_damp_y(j) = 0.0
[1]866          ENDIF
867       ENDDO
868
[73]869    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
[1]870
871       DO  j = nys-1, nyn+1
[73]872          IF ( j <= outflow_damping_width )  THEN
873             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
874                            ( outflow_damping_width - j ) /            &
875                            REAL( outflow_damping_width/2 )            &
876                                             )
877          ELSE
878             km_damp_y(j) = 0.0
[1]879          ENDIF
[73]880       ENDDO
[1]881
882    ENDIF
883
884!
885!-- Initialize local summation arrays for UP flow_statistics. This is necessary
886!-- because they may not yet have been initialized when they are called from
887!-- flow_statistics (or - depending on the chosen model run - are never
888!-- initialized)
889    sums_divnew_l      = 0.0
890    sums_divold_l      = 0.0
891    sums_l_l           = 0.0
892    sums_up_fraction_l = 0.0
893    sums_wsts_bc_l     = 0.0
894
895!
896!-- Pre-set masks for regional statistics. Default is the total model domain.
897    rmask = 1.0
898
899!
[51]900!-- User-defined initializing actions. Check afterwards, if maximum number
901!-- of allowed timeseries is not exceeded
[1]902    CALL user_init
903
[51]904    IF ( dots_num > dots_max )  THEN
905       IF ( myid == 0 )  THEN
906          PRINT*, '+++ user_init: number of time series quantities exceeds', &
907                  ' its maximum of dots_max = ', dots_max
908          PRINT*, '    Please increase dots_max in modules.f90.'
909       ENDIF
910       CALL local_stop
911    ENDIF
912
[1]913!
914!-- Input binary data file is not needed anymore. This line must be placed
915!-- after call of user_init!
916    CALL close_file( 13 )
917
918!
919!-- Compute total sum of active mask grid points
920!-- ngp_2dh: number of grid points of a horizontal cross section through the
921!--          total domain
922!-- ngp_3d:  number of grid points of the total domain
923!-- Note: The lower vertical index nzb_s_outer imposes a small error on the 2D
924!-- ----  averages of staggered variables such as u and v due to the topography
925!--       arrangement on the staggered grid. Maybe revise later.
926    ngp_2dh_outer_l = 0
927    ngp_2dh_outer   = 0
928    ngp_2dh_l       = 0
929    ngp_2dh         = 0
930    ngp_3d_inner_l  = 0
931    ngp_3d_inner    = 0
932    ngp_3d          = 0
933    ngp_sums        = ( nz + 2 ) * var_sum
934
935    DO  sr = 0, statistic_regions
936       DO  i = nxl, nxr
937          DO  j = nys, nyn
938             IF ( rmask(j,i,sr) == 1.0 )  THEN
939!
940!--             All xy-grid points
941                ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
942!
943!--             xy-grid points above topography
944                DO  k = nzb_s_outer(j,i), nz + 1
945                   ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr) + 1
946                ENDDO
947!
948!--             All grid points of the total domain above topography
949                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + &
950                                     ( nz - nzb_s_inner(j,i) + 2 )
951             ENDIF
952          ENDDO
953       ENDDO
954    ENDDO
955
956    sr = statistic_regions + 1
957#if defined( __parallel )
958    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,  &
959                        comm2d, ierr )
960    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr, &
961                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
962    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner(0), sr, MPI_INTEGER, &
963                        MPI_SUM, comm2d, ierr )
964#else
965    ngp_2dh       = ngp_2dh_l
966    ngp_2dh_outer = ngp_2dh_outer_l
967    ngp_3d_inner  = ngp_3d_inner_l
968#endif
969
970    ngp_3d = ngp_2dh * ( nz + 2 )
971
972!
973!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
974!-- buoyancy, etc. A zero value will occur for cases where all grid points of
975!-- the respective subdomain lie below the surface topography
976    ngp_2dh_outer = MAX( 1, ngp_2dh_outer(:,:) ) 
977    ngp_3d_inner  = MAX( 1, ngp_3d_inner(:)    )
978
979    DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l )
980
981
982 END SUBROUTINE init_3d_model
Note: See TracBrowser for help on using the repository browser.