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

Last change on this file since 4 was 4, checked in by raasch, 15 years ago

Id keyword set as property for all *.f90 files

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