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

Last change on this file since 808 was 808, checked in by maronga, 12 years ago

last commit documented

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