source: palm/trunk/SOURCE/init_masks.f90 @ 553

Last change on this file since 553 was 553, checked in by weinreis, 14 years ago

several parameters for masked output are replaced by arrays

  • Property svn:keywords set to Id
File size: 23.0 KB
Line 
1 SUBROUTINE init_masks
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! parameters for masked output are replaced by arrays
7!
8! Former revisions:
9! -----------------
10! $Id: init_masks.f90 553 2010-09-01 14:09:06Z weinreis $
11!
12! 493 2010-03-01 08:30:24Z raasch
13! netcdf_format_mask* and format_parallel_io replaced by netcdf_data_format
14!
15! 410 2009-12-04 17:05:40Z letzel
16! Initial revision
17!
18!
19! Description:
20! ------------
21! Initialize masked data output
22!------------------------------------------------------------------------------!
23
24    USE arrays_3d
25    USE control_parameters
26    USE grid_variables
27    USE indices
28    USE netcdf_control
29    USE particle_attributes
30    USE pegrid
31
32    IMPLICIT NONE
33
34    CHARACTER (LEN=6)   ::  var
35    CHARACTER (LEN=7)   ::  unit
36    CHARACTER (LEN=10), DIMENSION(max_masks,100) ::  do_mask, do_mask_user
37
38    INTEGER :: i, ilen, ind(6), ind_array(1), j, k, n, sender
39    INTEGER, DIMENSION(:), ALLOCATABLE ::  tmp_array
40
41    LOGICAL ::  found
42!
43!-- Allocation and initialization
44    ALLOCATE( tmp_array( MAX(nx,ny,nz)+1 ) )
45
46    ALLOCATE( mask_i(max_masks,nxr-nxl+1), &
47              mask_j(max_masks,nyn-nys+1), &
48              mask_k(max_masks,nzt-nzb+1) )
49!
50!-- internal mask arrays ("mask,dimension,selection")
51    ALLOCATE( mask(max_masks,3,mask_xyz_dimension), &
52              mask_loop(max_masks,3,3) )
53   
54!
55!-- Parallel mask output not yet tested
56    IF ( netcdf_data_format > 2 )  THEN
57       message_string = 'NetCDF file formats '//                         &
58                        '3 (NetCDF 4) and 4 (NetCDF 4 Classic model)'//  &
59                        '&are currently not supported (not yet tested).'
60       CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
61    ENDIF
62
63!
64!-- Store data output parameters for masked data output in few shared arrays
65    DO mid = 1, masks
66       do_mask     (mid,:) = data_output_masks(mid,:)
67       do_mask_user(mid,:) = data_output_masks_user(mid,:)
68       mask      (mid,1,:) = mask_x(mid,:) 
69       mask      (mid,2,:) = mask_y(mid,:)
70       mask      (mid,3,:) = mask_z(mid,:) 
71       mask_loop (mid,1,:) = mask_x_loop(mid,:) 
72       mask_loop (mid,2,:) = mask_y_loop(mid,:) 
73       mask_loop (mid,3,:) = mask_z_loop(mid,:) 
74    ENDDO
75   
76    mask_i = -1; mask_j = -1; mask_k = -1
77   
78!
79!-- Global arrays are required by define_netcdf_header.
80    IF ( myid == 0  .OR.  netcdf_data_format > 2 )  THEN
81       ALLOCATE( mask_i_global(max_masks,nx+1), &
82                 mask_j_global(max_masks,ny+1), &
83                 mask_k_global(max_masks,nz+1) )
84       mask_i_global = -1; mask_j_global = -1; mask_k_global = -1
85    ENDIF
86
87!
88!-- Determine variable names for each mask
89    DO  mid = 1, masks
90!
91!--    Append user-defined data output variables to the standard data output
92       IF ( do_mask_user(mid,1) /= ' ' )  THEN
93          i = 1
94          DO  WHILE ( do_mask(mid,i) /= ' '  .AND.  i <= 100 )
95             i = i + 1
96          ENDDO
97          j = 1
98          DO  WHILE ( do_mask_user(mid,j) /= ' '  .AND.  j <= 100 )
99             IF ( i > 100 )  THEN
100                WRITE ( message_string, * ) 'number of output quantitities ', &
101                     'given by data_output_mask and data_output_mask_user ', &
102                     'exceeds the limit of 100'
103                CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
104             ENDIF
105             do_mask(mid,i) = do_mask_user(mid,j)
106             i = i + 1
107             j = j + 1
108          ENDDO
109       ENDIF
110
111!
112!--    Check and set steering parameters for mask data output and averaging
113       i   = 1
114       DO  WHILE ( do_mask(mid,i) /= ' '  .AND.  i <= 100 )
115!
116!--       Check for data averaging
117          ilen = LEN_TRIM( do_mask(mid,i) )
118          j = 0                                              ! no data averaging
119          IF ( ilen > 3 )  THEN
120             IF ( do_mask(mid,i)(ilen-2:ilen) == '_av' )  THEN
121                j = 1                                           ! data averaging
122                do_mask(mid,i) = do_mask(mid,i)(1:ilen-3)
123             ENDIF
124          ENDIF
125          var = TRIM( do_mask(mid,i) )
126!
127!--       Check for allowed value and set units
128          SELECT CASE ( TRIM( var ) )
129
130             CASE ( 'e' )
131                IF ( constant_diffusion )  THEN
132                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
133                        '" requires constant_diffusion = .FALSE.'
134                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
135                ENDIF
136                unit = 'm2/s2'
137
138             CASE ( 'pc', 'pr' )
139                IF ( .NOT. particle_advection )  THEN
140                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
141                        '" requires a "particles_par"-NAMELIST in the ', &
142                        'parameter file (PARIN)'
143                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
144                ENDIF
145                IF ( TRIM( var ) == 'pc' )  unit = 'number'
146                IF ( TRIM( var ) == 'pr' )  unit = 'm'
147
148             CASE ( 'q', 'vpt' )
149                IF ( .NOT. humidity )  THEN
150                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
151                        '" requires humidity = .TRUE.'
152                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
153                ENDIF
154                IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
155                IF ( TRIM( var ) == 'vpt' )  unit = 'K'
156
157             CASE ( 'ql' )
158                IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
159                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
160                        '" requires cloud_physics = .TRUE. or cloud_droplets', &
161                        ' = .TRUE.'
162                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
163                ENDIF
164                unit = 'kg/kg'
165
166             CASE ( 'ql_c', 'ql_v', 'ql_vp' )
167                IF ( .NOT. cloud_droplets )  THEN
168                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
169                        '" requires cloud_droplets = .TRUE.'
170                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
171                ENDIF
172                IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
173                IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
174                IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
175
176             CASE ( 'qv' )
177                IF ( .NOT. cloud_physics )  THEN
178                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
179                        '" requires cloud_physics = .TRUE.'
180                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
181                ENDIF
182                unit = 'kg/kg'
183
184             CASE ( 'rho' )
185                IF ( .NOT. ocean )  THEN
186                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
187                        '" requires ocean = .TRUE.'
188                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
189                ENDIF
190                unit = 'kg/m3'
191
192             CASE ( 's' )
193                IF ( .NOT. passive_scalar )  THEN
194                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
195                        '" requires passive_scalar = .TRUE.'
196                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
197                ENDIF
198                unit = 'conc'
199
200             CASE ( 'sa' )
201                IF ( .NOT. ocean )  THEN
202                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
203                        '" requires ocean = .TRUE.'
204                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
205                ENDIF
206                unit = 'psu'
207
208             CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*' )
209                WRITE ( message_string, * ) 'illegal value for data_',&
210                     'output: "', TRIM( var ), '" is only allowed', &
211                     'for horizontal cross section'
212                CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
213!                IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
214!                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
215!                         '" requires cloud_physics = .TRUE.'
216!                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
217!                ENDIF
218!                IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
219!                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
220!                         '" requires precipitation = .TRUE.'
221!                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
222!                ENDIF
223!                IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
224!                   WRITE ( message_string, * ) 'temporal averaging of ', &
225!                         ' precipitation amount "', TRIM( var ),         &
226!                         '" not possible'
227!                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
228!                ENDIF
229!                IF ( TRIM( var ) == 'prr*'  .AND.  .NOT. precipitation )  THEN
230!                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
231!                         '" requires precipitation = .TRUE.'
232!                   CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
233!                ENDIF
234!
235!                IF ( TRIM( var ) == 'u*'   )  unit = 'm/s'
236!                IF ( TRIM( var ) == 't*'   )  unit = 'K'
237!                IF ( TRIM( var ) == 'lwp*' )  unit = 'kg/kg*m'
238!                IF ( TRIM( var ) == 'pra*' )  unit = 'mm'
239!                IF ( TRIM( var ) == 'prr*' )  unit = 'mm/s'
240!                IF ( TRIM( var ) == 'z0*'  )  unit = 'm'
241
242             CASE ( 'p', 'pt', 'u', 'v', 'w' )
243                IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
244                IF ( TRIM( var ) == 'pt' )  unit = 'K'
245                IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
246                IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
247                IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
248                CONTINUE
249
250             CASE DEFAULT
251                CALL user_check_data_output( var, unit )
252
253                IF ( unit == 'illegal' )  THEN
254                   IF ( do_mask_user(mid,1) /= ' ' )  THEN
255                      WRITE ( message_string, * ) 'illegal value for data_',&
256                           'output or data_output_user: "', &
257                           TRIM( do_mask(mid,i) ), '"'
258                      CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
259                   ELSE
260                      WRITE ( message_string, * ) 'illegal value for data_',&
261                           'output: "', TRIM( do_mask(mid,i) ), '"'
262                      CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
263                   ENDIF
264                ENDIF
265
266          END SELECT
267!
268!--       Set the internal steering parameters appropriately
269          domask_no(mid,j)                    = domask_no(mid,j) + 1
270          domask(mid,j,domask_no(mid,j))      = do_mask(mid,i)
271          domask_unit(mid,j,domask_no(mid,j)) = unit
272
273          IF ( j == 1 )  THEN
274!
275!--          Check, if variable is already subject to averaging
276             found = .FALSE.
277             DO  k = 1, doav_n
278                IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
279             ENDDO
280
281             IF ( .NOT. found )  THEN
282                doav_n = doav_n + 1
283                doav(doav_n) = var
284             ENDIF
285          ENDIF
286
287          i = i + 1
288
289       ENDDO   ! do_mask(mid,i)
290    ENDDO   ! mid
291
292
293!
294!-- Determine mask locations for each mask
295    DO  mid = 1, masks
296!
297!--    Set local masks for each subdomain along all three dimensions
298       CALL set_mask_locations( 1, dx, 'dx', nx, 'nx', nxl, nxr )
299       CALL set_mask_locations( 2, dy, 'dy', ny, 'ny', nys, nyn )
300       CALL set_mask_locations( 3, dz, 'dz', nz, 'nz', nzb, nzt )
301!
302!--    Set global masks along all three dimensions (required by
303!--    define_netcdf_header).
304#if defined( __parallel )
305!
306!--    PE0 receives partial arrays from all processors of the respective mask
307!--    and outputs them. Here a barrier has to be set, because otherwise
308!--    "-MPI- FATAL: Remote protocol queue full" may occur.
309       CALL MPI_BARRIER( comm2d, ierr )
310
311       IF ( myid == 0 )  THEN
312!
313!--       Local arrays can be relocated directly.
314          mask_i_global(mid,mask_start_l(mid,1): &
315                       mask_start_l(mid,1)+mask_size_l(mid,1)-1) = &
316                       mask_i(mid,:mask_size_l(mid,1))
317          mask_j_global(mid,mask_start_l(mid,2): &
318                       mask_start_l(mid,2)+mask_size_l(mid,2)-1) = &
319                       mask_j(mid,:mask_size_l(mid,2))
320          mask_k_global(mid,mask_start_l(mid,3): &
321                       mask_start_l(mid,3)+mask_size_l(mid,3)-1) = &
322                       mask_k(mid,:mask_size_l(mid,3))
323!
324!--       Receive data from all other PEs.
325          DO  n = 1, numprocs-1
326!
327!--          Receive index limits first, then arrays.
328!--          Index limits are received in arbitrary order from the PEs.
329             CALL MPI_RECV( ind(1), 6, MPI_INTEGER, MPI_ANY_SOURCE, 0,  &
330                  comm2d, status, ierr )
331!
332!--          Not all PEs have data for the mask.
333             IF ( ind(1) /= -9999 )  THEN
334                sender = status(MPI_SOURCE)
335                CALL MPI_RECV( tmp_array(ind(1)), ind(2)-ind(1)+1,  &
336                               MPI_INTEGER, sender, 1, comm2d, status, ierr )
337                mask_i_global(mid,ind(1):ind(2)) = tmp_array(ind(1):ind(2))
338                CALL MPI_RECV( tmp_array(ind(3)), ind(4)-ind(3)+1,  &
339                               MPI_INTEGER, sender, 2, comm2d, status, ierr )
340                mask_j_global(mid,ind(3):ind(4)) = tmp_array(ind(3):ind(4))
341                CALL MPI_RECV( tmp_array(ind(5)), ind(6)-ind(5)+1,  &
342                               MPI_INTEGER, sender, 3, comm2d, status, ierr )
343                mask_k_global(mid,ind(5):ind(6)) = tmp_array(ind(5):ind(6))
344             ENDIF
345          ENDDO
346
347       ELSE
348!
349!--       If at least part of the mask resides on the PE, send the index limits
350!--       for the target array, otherwise send -9999 to PE0.
351          IF ( mask_size_l(mid,1) > 0 .AND.  mask_size_l(mid,2) > 0 .AND. &
352               mask_size_l(mid,3) > 0  ) &
353               THEN
354             ind(1) = mask_start_l(mid,1)
355             ind(2) = mask_start_l(mid,1) + mask_size_l(mid,1) - 1
356             ind(3) = mask_start_l(mid,2)
357             ind(4) = mask_start_l(mid,2) + mask_size_l(mid,2) - 1
358             ind(5) = mask_start_l(mid,3)
359             ind(6) = mask_start_l(mid,3) + mask_size_l(mid,3) - 1
360          ELSE
361             ind(1) = -9999; ind(2) = -9999
362             ind(3) = -9999; ind(4) = -9999
363             ind(5) = -9999; ind(6) = -9999
364          ENDIF
365          CALL MPI_SEND( ind(1), 6, MPI_INTEGER, 0, 0, comm2d, ierr )
366!
367!--       If applicable, send data to PE0.
368          IF ( ind(1) /= -9999 )  THEN
369             tmp_array(:mask_size_l(mid,1)) = mask_i(mid,:mask_size_l(mid,1))
370             CALL MPI_SEND( tmp_array(1), mask_size_l(mid,1),  &
371                            MPI_INTEGER, 0, 1, comm2d, ierr )
372             tmp_array(:mask_size_l(mid,2)) = mask_j(mid,:mask_size_l(mid,2))
373             CALL MPI_SEND( tmp_array(1), mask_size_l(mid,2),  &
374                            MPI_INTEGER, 0, 2, comm2d, ierr )
375             tmp_array(:mask_size_l(mid,3)) = mask_k(mid,:mask_size_l(mid,3))
376             CALL MPI_SEND( tmp_array(1), mask_size_l(mid,3),  &
377                            MPI_INTEGER, 0, 3, comm2d, ierr )
378          ENDIF
379       ENDIF
380!
381!--    A barrier has to be set, because otherwise some PEs may proceed too fast
382!--    so that PE0 may receive wrong data on tag 0.
383       CALL MPI_BARRIER( comm2d, ierr )
384       
385       IF ( netcdf_data_format > 2 )  THEN 
386         
387          CALL MPI_BCAST( mask_i_global(mid,:), nx+1, MPI_INTEGER, 0, comm2d, &
388                          ierr )
389          CALL MPI_BCAST( mask_j_global(mid,:), ny+1, MPI_INTEGER, 0, comm2d, &
390                          ierr )
391          CALL MPI_BCAST( mask_k_global(mid,:), nz+1, MPI_INTEGER, 0, comm2d, &
392                          ierr ) 
393         
394       ENDIF
395
396#else
397!
398!--    Local arrays can be relocated directly.
399       mask_i_global(mid,:) = mask_i(mid,:)
400       mask_j_global(mid,:) = mask_j(mid,:)
401       mask_k_global(mid,:) = mask_k(mid,:)
402#endif
403    ENDDO   ! mid
404
405    DEALLOCATE( tmp_array )
406!
407!-- Internal mask arrays cannot be deallocated on PE 0 because they are
408!-- required for header output on PE 0.
409    IF ( myid /= 0 )  DEALLOCATE( mask, mask_loop )
410
411 CONTAINS
412
413    SUBROUTINE set_mask_locations( dim, dxyz, dxyz_string, nxyz, nxyz_string, &
414                                   lb, ub )
415!------------------------------------------------------------------------------!
416!
417! Description:
418! ------------
419! Set local mask for each subdomain along 'dim' direction.
420!------------------------------------------------------------------------------!
421
422       IMPLICIT NONE
423
424       CHARACTER (LEN=2) :: dxyz_string, nxyz_string
425       INTEGER  ::  count, count_l, dim, m, loop_begin, loop_end, loop_stride, &
426                    lb, nxyz, ub
427       REAL     ::  dxyz, ddxyz, tmp1, tmp2
428
429       count = 0;  count_l = 0;  ddxyz = 1.0 / dxyz;  tmp1 = 0.0;  tmp2 = 0.0
430
431       IF ( mask(mid,dim,1) >= 0.0 )  THEN
432!
433!--       use predefined mask_* array
434          DO  WHILE ( mask(mid,dim,count+1) >= 0.0 )
435             count = count + 1
436             IF ( dim == 1 .OR. dim == 2 )  THEN
437                m = NINT( mask(mid,dim,count) * mask_scale(dim) * ddxyz )
438             ELSEIF ( dim == 3 )  THEN
439                ind_array =  &
440                     MINLOC( ABS( mask(mid,dim,count) * mask_scale(dim) - zw ) )
441                m = ind_array(1) - 1 + nzb  ! MINLOC uses lower array bound 1
442             ENDIF
443             IF ( m > nxyz )  THEN
444                WRITE ( message_string, '(I3,A,I3,A,I1,3A,I3)' )  &
445                     m,' in mask ',mid,' along dimension ',dim,  &
446                     ' exceeds ',nxyz_string,' = ',nxyz
447                CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
448             ENDIF
449             IF ( m >= lb .AND. m <= ub )  THEN
450                IF ( count_l == 0 )  mask_start_l(mid,dim) = count
451                count_l = count_l + 1
452                IF ( dim == 1 )  THEN
453                   mask_i(mid,count_l) = m
454                ELSEIF ( dim == 2 )  THEN
455                   mask_j(mid,count_l) = m
456                ELSEIF ( dim == 3 )  THEN
457                   mask_k(mid,count_l) = m
458                ENDIF
459             ENDIF
460             IF ( count == mask_xyz_dimension )  EXIT
461          ENDDO
462          mask_size(mid,dim)   = count
463          mask_size_l(mid,dim) = count_l
464
465       ELSE
466!
467!--       use predefined mask_loop_* array, or use the default (all grid points
468!--       along this direction)
469          IF ( mask_loop(mid,dim,1) < 0.0 )  THEN
470             tmp1 = mask_loop(mid,dim,1)
471             mask_loop(mid,dim,1) = 0.0   ! (default)
472          ENDIF
473          IF ( dim == 1 .OR. dim == 2 )  THEN
474             IF ( mask_loop(mid,dim,2) < 0.0 )  THEN
475                tmp2 = mask_loop(mid,dim,2)
476                mask_loop(mid,dim,2) = nxyz*dxyz / mask_scale(dim)   ! (default)
477             ENDIF
478             IF ( MAXVAL( mask_loop(mid,dim,1:2) )  &
479                  > nxyz * dxyz / mask_scale(dim) )  THEN
480                WRITE ( message_string, '(2(A,I3,A,I1,A,F9.3),5A,I1,A,F9.3)' ) &
481                     'mask_loop(',mid,',',dim,',1)=',mask_loop(mid,dim,1), &
482                     ' and/or mask_loop(',mid,',',dim,',2)=', &
483                     mask_loop(mid,dim,2),'&exceed ', &
484                     nxyz_string,'*',dxyz_string,'/mask_scale(',dim,')=', &
485                     nxyz*dxyz/mask_scale(dim)
486                CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
487             ENDIF
488             loop_begin  = NINT( mask_loop(mid,dim,1) * mask_scale(dim) * ddxyz )
489             loop_end    = NINT( mask_loop(mid,dim,2) * mask_scale(dim) * ddxyz )
490             loop_stride = NINT( mask_loop(mid,dim,3) * mask_scale(dim) * ddxyz )
491          ELSEIF ( dim == 3 )  THEN
492             IF ( mask_loop(mid,dim,2) < 0.0 )  THEN
493                tmp2 = mask_loop(mid,dim,2)
494                mask_loop(mid,dim,2) = zw(nz) / mask_scale(dim)   ! (default)
495             ENDIF
496             IF ( MAXVAL( mask_loop(mid,dim,1:2) )  &
497                  > zw(nz) / mask_scale(dim) )  THEN
498                WRITE ( message_string, '(2(A,I3,A,I1,A,F9.3),A,I1,A,F9.3)' ) &
499                     'mask_loop(',mid,',',dim,',1)=',mask_loop(mid,dim,1), &
500                     ' and/or mask_loop(',mid,',',dim,',2)=', &
501                     mask_loop(mid,dim,2),'&exceed zw(nz)/mask_scale(',dim, &
502                     ')=',zw(nz)/mask_scale(dim)
503                CALL message( 'init_masks', 'PA9998', 1, 2, 0, 6, 0 )
504             ENDIF
505             ind_array =  &
506                  MINLOC( ABS( mask_loop(mid,dim,1) * mask_scale(dim) - zw ) )
507             loop_begin =  &
508                  ind_array(1) - 1 + nzb ! MINLOC uses lower array bound 1
509             ind_array =  &
510                  MINLOC( ABS( mask_loop(mid,dim,2) * mask_scale(dim) - zw ) )
511             loop_end = ind_array(1) - 1 + nzb ! MINLOC uses lower array bound 1
512!
513!--          The following line assumes a constant vertical grid spacing within
514!--          the vertical mask range; it fails for vertical grid stretching.
515!--          Maybe revise later. Issue warning but continue execution.
516             loop_stride = NINT( mask_loop(mid,dim,3) * mask_scale(dim) * ddxyz )
517
518             IF ( mask_loop(mid,dim,2) * mask_scale(dim) > dz_stretch_level )  &
519                  THEN
520                WRITE ( message_string, '(A,I3,A,I1,A,F9.3,A,F8.2,3A)' )  &
521                     'mask_loop(',mid,dim,',2)=', mask_loop(mid,dim,2),&
522                     ' exceeds dz_stretch_level=',dz_stretch_level, &
523                     '.&Vertical mask locations will not ', &
524                     'match the desired heights&within the stretching ', &
525                     'region. Recommendation: use mask instead of mask_loop.'
526                CALL message( 'init_masks', 'PA9998', 0, 1, 0, 6, 0 )
527             ENDIF
528
529          ENDIF
530!
531!--       If necessary, reset mask_loop(mid,dim,1) and mask_loop(mid,dim,2).
532          IF ( tmp1 < 0.0 )  mask_loop(mid,dim,1) = tmp1
533          IF ( tmp2 < 0.0 )  mask_loop(mid,dim,2) = tmp2
534!
535!--       The default stride +/-1 (every grid point) applies if
536!--       mask_loop(mid,dim,3) is not specified (its default is zero).
537          IF ( loop_stride == 0 )  THEN
538             IF ( loop_end >= loop_begin )  THEN
539                loop_stride =  1
540             ELSE
541                loop_stride = -1
542             ENDIF
543          ENDIF
544          DO  m = loop_begin, loop_end, loop_stride
545             count = count + 1
546             IF ( m >= lb .AND. m <= ub )  THEN
547                IF ( count_l == 0 )  mask_start_l(mid,dim) = count
548                count_l = count_l + 1
549                IF ( dim == 1 )  THEN
550                   mask_i(mid,count_l) = m
551                ELSEIF ( dim == 2 )  THEN
552                   mask_j(mid,count_l) = m
553                ELSEIF ( dim == 3 )  THEN
554                   mask_k(mid,count_l) = m
555                ENDIF
556             ENDIF
557          ENDDO
558          mask_size(mid,dim)   = count
559          mask_size_l(mid,dim) = count_l
560       ENDIF
561
562    END SUBROUTINE set_mask_locations
563
564 END SUBROUTINE init_masks
Note: See TracBrowser for help on using the repository browser.