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

Last change on this file since 19 was 19, checked in by raasch, 17 years ago

preliminary version of modified boundary conditions at top

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