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

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

updating comments and rc-file

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