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

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

more preliminary uncomplete changes for ocean version

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