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

Last change on this file since 108 was 108, checked in by letzel, 14 years ago
  • Improved coupler: evaporation - salinity-flux coupling for humidity = .T.,

avoid MPI hangs when coupled runs terminate, add DOC/app/chapter_3.8;

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