source: palm/trunk/SOURCE/init_masks.f90

Last change on this file was 4896, checked in by raasch, 7 months ago

small re-formatting to follow the coding standard, typo in file appendix removed, more meaningful variable names assigned, redundant code removed

  • Property svn:keywords set to Id
File size: 36.8 KB
Line 
1!> @file init_masks.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2021 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: init_masks.f90 4896 2021-03-03 16:10:18Z banzhafs $
26! small re-formatting to follow the coding standard
27!
28! 4895 2021-03-03 15:39:08Z suehring
29! Allow only mask_k_over_surface >= 1 for terrain-following mask output
30!
31! 4828 2021-01-05 11:21:41Z Giersch
32! Implement snow and graupel (bulk microphysics)
33!
34! 4648 2020-08-25 07:52:08Z raasch
35! file re-formatted to follow the PALM coding standard
36!
37! 4521 2020-05-06 11:39:49Z schwenkel
38! Rename variable
39!
40! 4502 2020-04-17 16:14:16Z schwenkel
41! Implementation of ice microphysics
42!
43! 4444 2020-03-05 15:59:50Z raasch
44! bugfix: cpp-directives for serial mode added
45!
46! 4360 2020-01-07 11:25:50Z suehring
47! Corrected "Former revisions" section
48!
49! 4069 2019-07-01 14:05:51Z Giersch
50! Masked output running index mid has been introduced as a local variable to avoid runtime error
51! (Loop variable has been modified) in time_integration
52!
53! 3766 2019-02-26 16:23:41Z raasch
54! unused variables removed
55!
56! 3687 2019-01-22 10:42:06Z knoop
57! unused variables removed
58!
59! 3655 2019-01-07 16:51:22Z knoop
60! Move the control parameter "salsa" from salsa_mod to control_parameters (M. Kurppa)
61!
62! 410 2009-12-04 17:05:40Z letzel
63! Initial revision
64!
65!
66! Description:
67! ------------
68!> Initialize masked data output
69!--------------------------------------------------------------------------------------------------!
70 SUBROUTINE init_masks
71
72    USE arrays_3d,                                                                                 &
73        ONLY:  zu, zw
74
75    USE bulk_cloud_model_mod,                                                                      &
76        ONLY: bulk_cloud_model, microphysics_ice_phase, microphysics_morrison,                     &
77              microphysics_seifert, snow, graupel
78
79
80    USE control_parameters,                                                                        &
81        ONLY:  constant_diffusion, cloud_droplets, data_output_masks, data_output_masks_user, doav,&
82               doav_n, domask, domask_no, dz, dz_stretch_level_start, humidity, mask, masks,       &
83               mask_scale, mask_i, mask_i_global, mask_j, mask_j_global, mask_k, mask_k_global,    &
84               mask_k_over_surface, mask_loop, mask_size, mask_size_l, mask_start_l, mask_surface, &
85               mask_x, mask_x_loop, mask_xyz_dimension, mask_y, mask_y_loop, mask_z, mask_z_loop,  &
86               max_masks,  message_string, passive_scalar, ocean_mode, varnamelength
87
88    USE grid_variables,                                                                            &
89        ONLY:  dx, dy
90
91    USE indices,                                                                                   &
92        ONLY:  nx, nxl, nxr, ny, nyn, nys, nz, nzb, nzt
93
94    USE kinds
95
96    USE module_interface,                                                                          &
97        ONLY:  module_interface_init_masks
98
99    USE netcdf_interface,                                                                          &
100        ONLY:  domask_unit, netcdf_data_format
101
102    USE particle_attributes,                                                                       &
103        ONLY:  particle_advection
104
105    USE pegrid
106
107    IMPLICIT NONE
108
109    CHARACTER (LEN=varnamelength) ::  var  !< contains variable name
110    CHARACTER (LEN=7)             ::  unit !< contains unit of variable
111
112    CHARACTER (LEN=varnamelength), DIMENSION(max_masks,100) ::  do_mask      !< list of output variables
113    CHARACTER (LEN=varnamelength), DIMENSION(max_masks,100) ::  do_mask_user !< list of user-specified output variables
114
115    INTEGER(iwp) ::  count        !< counting masking indices along a dimension
116    INTEGER(iwp) ::  i            !< loop index
117    INTEGER(iwp) ::  ilen         !< length of string saved in 'do_mask'
118    INTEGER(iwp) ::  ind_array(1) !< array index
119    INTEGER(iwp) ::  j            !< loop index
120    INTEGER(iwp) ::  k            !< loop index
121    INTEGER(iwp) ::  m            !< mask index
122    INTEGER(iwp) ::  mid            !< masked output running index
123#if defined( __parallel )
124    INTEGER(iwp) ::  ind(6)       !< index limits (lower/upper bounds) of output array
125    INTEGER(iwp) ::  n            !< loop index
126    INTEGER(iwp) ::  sender       !< PE id of sending PE
127#endif
128
129    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  tmp_array !< temporary 1D array
130
131    LOGICAL ::  found !< true if variable is found
132
133!
134!-- Initial values are explicitly set here due to a bug in the Cray compiler in case of assignments
135!-- of initial values in declaration statements for arrays with more than 9999 elements
136!-- (appears with -eD only)
137    domask = ' '
138
139!
140!-- Allocation and initialization
141    ALLOCATE( tmp_array( MAX(nx,ny,nz)+2 ) )
142
143    ALLOCATE( mask_i(max_masks,nxr-nxl+2),                                                         &
144              mask_j(max_masks,nyn-nys+2),                                                         &
145              mask_k(max_masks,nzt-nzb+2) )
146!
147!-- internal mask arrays ("mask,dimension,selection")
148    ALLOCATE( mask(max_masks,3,mask_xyz_dimension), mask_loop(max_masks,3,3) )
149
150!
151!-- Parallel mask output not yet supported. In check_parameters data format is restricted and is
152!-- switched back to non-parallel output. Therefore the following error can not occur at the moment.
153    IF ( netcdf_data_format > 4 )  THEN
154       message_string = 'netCDF file formats '//                                                   &
155                        '5 and 6 (with parallel I/O support)'//                                    &
156                        ' are currently not supported.'
157       CALL message( 'init_masks', 'PA0328', 1, 2, 0, 6, 0 )
158    ENDIF
159
160!
161!-- Store data output parameters for masked data output in few shared arrays
162    DO  mid = 1, masks
163
164       do_mask(mid,:)      = data_output_masks(mid,:)
165       do_mask_user(mid,:) = data_output_masks_user(mid,:)
166
167       mask(mid,1,:) = mask_x(mid,:)
168       mask(mid,2,:) = mask_y(mid,:)
169       mask(mid,3,:) = mask_z(mid,:)
170!
171!--    Flag a mask as terrain following
172       IF ( mask_k_over_surface(mid,1) /= -1_iwp )  THEN
173          mask_surface(mid) = .TRUE.
174       ENDIF
175
176       IF ( mask_x_loop(mid,1) == -1.0_wp  .AND.  mask_x_loop(mid,2) == -1.0_wp  .AND.             &
177            mask_x_loop(mid,3) == -1.0_wp )  THEN
178          mask_loop(mid,1,1:2) = -1.0_wp
179          mask_loop(mid,1,3)   =  0.0_wp
180       ELSE
181          mask_loop(mid,1,:) = mask_x_loop(mid,:)
182       ENDIF
183       IF ( mask_y_loop(mid,1) == -1.0_wp  .AND.  mask_y_loop(mid,2) == -1.0_wp  .AND.             &
184            mask_y_loop(mid,3) == -1.0_wp )  THEN
185          mask_loop(mid,2,1:2) = -1.0_wp
186          mask_loop(mid,2,3)   =  0.0_wp
187       ELSE
188          mask_loop(mid,2,:) = mask_y_loop(mid,:)
189       ENDIF
190       IF ( mask_z_loop(mid,1) == -1.0_wp  .AND.  mask_z_loop(mid,2) == -1.0_wp  .AND.             &
191            mask_z_loop(mid,3) == -1.0_wp )  THEN
192          mask_loop(mid,3,1:2) = -1.0_wp
193          mask_loop(mid,3,3)   =  0.0_wp
194       ELSE
195          mask_loop(mid,3,:) = mask_z_loop(mid,:)
196       ENDIF
197
198    ENDDO
199
200    mask_i = -1; mask_j = -1; mask_k = -1
201
202!
203!-- Global arrays are required by define_netcdf_header.
204    IF ( myid == 0  .OR.  netcdf_data_format > 4 )  THEN
205       ALLOCATE( mask_i_global(max_masks,nx+2),                                                    &
206                 mask_j_global(max_masks,ny+2),                                                    &
207                 mask_k_global(max_masks,nz+2) )
208       mask_i_global = -1; mask_j_global = -1; mask_k_global = -1
209    ENDIF
210
211!
212!-- Determine variable names for each mask
213    DO  mid = 1, masks
214!
215!--    Append user-defined data output variables to the standard data output
216       IF ( do_mask_user(mid,1) /= ' ' )  THEN
217          i = 1
218          DO  WHILE ( do_mask(mid,i) /= ' '  .AND.  i <= 100 )
219             i = i + 1
220          ENDDO
221          j = 1
222          DO  WHILE ( do_mask_user(mid,j) /= ' '  .AND.  j <= 100 )
223             IF ( i > 100 )  THEN
224                WRITE ( message_string, * ) 'number of output quantitities ',                      &
225                                            'given by data_output_mask and data_output_mask_user ',&
226                                            'exceeds the limit of 100'
227                CALL message( 'init_masks', 'PA0329', 1, 2, 0, 6, 0 )
228             ENDIF
229             do_mask(mid,i) = do_mask_user(mid,j)
230             i = i + 1
231             j = j + 1
232          ENDDO
233       ENDIF
234
235!
236!--    Check and set steering parameters for mask data output and averaging
237       i   = 1
238       DO WHILE ( do_mask(mid,i) /= ' '  .AND.  i <= 100 )
239!
240!--       Check for data averaging
241          ilen = LEN_TRIM( do_mask(mid,i) )
242          j = 0                                              ! no data averaging
243          IF ( ilen > 3 )  THEN
244             IF ( do_mask(mid,i)(ilen-2:ilen) == '_av' )  THEN
245                j = 1                                           ! data averaging
246                do_mask(mid,i) = do_mask(mid,i)(1:ilen-3)
247             ENDIF
248          ENDIF
249          var = TRIM( do_mask(mid,i) )
250!
251!--       Check for allowed value and set units
252          SELECT CASE ( TRIM( var ) )
253
254             CASE ( 'e' )
255                IF ( constant_diffusion )  THEN
256                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
257                                               '" requires constant_diffusion = .FALSE.'
258                   CALL message( 'init_masks', 'PA0103', 1, 2, 0, 6, 0 )
259                ENDIF
260                unit = 'm2/s2'
261
262             CASE ( 'thetal' )
263                IF ( .NOT. bulk_cloud_model )  THEN
264                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
265                                               '" requires bulk_cloud_model = .TRUE.'
266                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
267                ENDIF
268                unit = 'K'
269
270             CASE ( 'nc' )
271                IF ( .NOT. bulk_cloud_model )  THEN
272                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
273                                               '" requires bulk_cloud_model = .TRUE.'
274                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
275                ELSEIF ( .NOT. microphysics_morrison )  THEN
276                   message_string = 'output of "' // TRIM( var ) // '" ' // 'requires  = morrison'
277                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
278                ENDIF
279                unit = '1/m3'
280
281             CASE ( 'ni' )
282                IF ( .NOT. bulk_cloud_model )  THEN
283                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
284                                               '" requires bulk_cloud_model = .TRUE.'
285                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
286                 ELSEIF ( .NOT. microphysics_ice_phase )  THEN
287                   message_string = 'output of "' // TRIM( var ) // '" ' //                        &
288                                    'requires  microphysics_ice_phase = .TRUE.'
289                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
290                ENDIF
291                unit = '1/m3'
292
293             CASE ( 'ng' )
294                IF ( .NOT. bulk_cloud_model )  THEN
295                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
296                        '" requires bulk_cloud_model = .TRUE.'
297                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
298                 ELSEIF ( .NOT. microphysics_ice_phase  .OR.  .NOT.  graupel )  THEN
299                   message_string = 'output of "' // TRIM( var ) // '" ' //                        &
300                         'requires  microphysics_ice_phase = .TRUE.'
301                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
302                ENDIF
303                unit = '1/m3'
304
305
306             CASE ( 'nr' )
307                IF ( .NOT. bulk_cloud_model )  THEN
308                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
309                                               '" requires bulk_cloud_model = .TRUE.'
310                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
311                ELSEIF ( .NOT. microphysics_seifert )  THEN
312                   message_string = 'output of "' // TRIM( var ) // '"' //                         &
313                                    'requires cloud_scheme = seifert_beheng'
314                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
315                ENDIF
316                unit = '1/m3'
317
318             CASE ( 'ns' )
319                IF ( .NOT. bulk_cloud_model )  THEN
320                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
321                        '" requires bulk_cloud_model = .TRUE.'
322                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
323                 ELSEIF ( .NOT. microphysics_seifert  .OR.  .NOT.  snow)  THEN
324                   message_string = 'output of "' // TRIM( var ) // '"' //                         &
325                         'requires cloud_scheme = seifert_beheng'
326                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
327                ENDIF
328                unit = '1/m3'
329
330
331             CASE ( 'pc', 'pr' )
332                IF ( .NOT. particle_advection )  THEN
333                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
334                                               '" requires a "particles_par"-NAMELIST in the ',    &
335                                               'parameter file (PARIN)'
336                   CALL message( 'init_masks', 'PA0104', 1, 2, 0, 6, 0 )
337                ENDIF
338                IF ( TRIM( var ) == 'pc' )  unit = 'number'
339                IF ( TRIM( var ) == 'pr' )  unit = 'm'
340
341             CASE ( 'q', 'thetav' )
342                IF ( .NOT. humidity )  THEN
343                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
344                                               '" requires humidity = .TRUE.'
345                   CALL message( 'init_masks', 'PA0105', 1, 2, 0, 6, 0 )
346                ENDIF
347                IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
348                IF ( TRIM( var ) == 'thetav' )  unit = 'K'
349
350             CASE ( 'qc' )
351                IF ( .NOT. bulk_cloud_model )  THEN
352                   message_string = 'output of "' // TRIM( var ) // '"' //                         &
353                                    'requires bulk_cloud_model = .TRUE.'
354                   CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
355                ENDIF
356                unit = 'kg/kg'
357
358             CASE ( 'ql' )
359                IF ( .NOT. ( bulk_cloud_model  .OR.  cloud_droplets ) )  THEN
360                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
361                                               '" requires bulk_cloud_model = .TRUE. or ',         &
362                                               'cloud_droplets = .TRUE.'
363                   CALL message( 'init_masks', 'PA0106', 1, 2, 0, 6, 0 )
364                ENDIF
365                unit = 'kg/kg'
366
367             CASE ( 'ql_c', 'ql_v', 'ql_vp' )
368                IF ( .NOT. cloud_droplets )  THEN
369                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
370                                               '" requires cloud_droplets = .TRUE.'
371                   CALL message( 'init_masks', 'PA0107', 1, 2, 0, 6, 0 )
372                ENDIF
373                IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
374                IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
375                IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
376
377             CASE ( 'qv' )
378                IF ( .NOT. bulk_cloud_model )  THEN
379                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
380                        '                      " requires bulk_cloud_model = .TRUE.'
381                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
382                ENDIF
383                unit = 'kg/kg'
384
385             CASE ( 'qg' )
386                IF ( .NOT. bulk_cloud_model )  THEN
387                   message_string = 'output of "' // TRIM( var ) // '" ' //                        &
388                            'requires bulk_cloud_model = .TRUE.'
389                   CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
390                ELSEIF ( .NOT. microphysics_ice_phase  .OR.  .NOT.  graupel )  THEN
391                   message_string = 'output of "' // TRIM( var ) // '" ' //                        &
392                            'requires microphysics_ice_phase = .TRUE.'
393                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
394                ENDIF
395                unit = 'kg/kg'
396
397             CASE ( 'qi' )
398                IF ( .NOT. bulk_cloud_model )  THEN
399                   message_string = 'output of "' // TRIM( var ) // '" ' //                        &
400                                    'requires bulk_cloud_model = .TRUE.'
401                   CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
402                ELSEIF ( .NOT. microphysics_ice_phase ) THEN
403                   message_string = 'output of "' // TRIM( var ) // '" ' //                        &
404                                    'requires microphysics_ice_phase = .TRUE.'
405                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
406                ENDIF
407                unit = 'kg/kg'
408
409             CASE ( 'qr' )
410                IF ( .NOT. bulk_cloud_model )  THEN
411                   message_string = 'output of "' // TRIM( var ) // '" ' //                        &
412                                    'requires bulk_cloud_model = .TRUE.'
413                   CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
414                ELSEIF ( .NOT. microphysics_seifert ) THEN
415                   message_string = 'output of "' // TRIM( var ) // '" ' //                        &
416                                    'requires cloud_scheme = seifert_beheng'
417                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
418                ENDIF
419                unit = 'kg/kg'
420
421             CASE ( 'qs' )
422                IF ( .NOT. bulk_cloud_model )  THEN
423                   message_string = 'output of "' // TRIM( var ) // '" ' //                        &
424                            'requires bulk_cloud_model = .TRUE.'
425                   CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
426                ELSEIF ( .NOT. microphysics_seifert  .OR.  snow )  THEN
427                   message_string = 'output of "' // TRIM( var ) // '" ' //                        &
428                            'requires cloud_scheme = seifert_beheng'
429                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
430                ENDIF
431                unit = 'kg/kg'
432
433
434             CASE ( 'rho_sea_water' )
435                IF ( .NOT. ocean_mode )  THEN
436                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
437                                               '" requires ocean mode'
438                   CALL message( 'init_masks', 'PA0109', 1, 2, 0, 6, 0 )
439                ENDIF
440                unit = 'kg/m3'
441
442             CASE ( 's' )
443                IF ( .NOT. passive_scalar )  THEN
444                   WRITE ( message_string, * ) 'output of "', TRIM( var ),                         &
445                                               '" requires passive_scalar = .TRUE.'
446                   CALL message( 'init_masks', 'PA0110', 1, 2, 0, 6, 0 )
447                ENDIF
448                unit = 'conc'
449
450             CASE ( 'sa' )
451                IF ( .NOT. ocean_mode )  THEN
452                   WRITE ( message_string, * ) 'output of "', TRIM( var ), '" requires ocean mode'
453                   CALL message( 'init_masks', 'PA0109', 1, 2, 0, 6, 0 )
454                ENDIF
455                unit = 'psu'
456
457             CASE ( 'us*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*', 'z0h*' )
458                WRITE ( message_string, * ) 'illegal value for data_', 'output: "', TRIM( var ),   &
459                                            '" is only allowed', 'for horizontal cross section'
460                CALL message( 'init_masks', 'PA0111', 1, 2, 0, 6, 0 )
461
462             CASE ( 'p', 'theta', 'u', 'v', 'w' )
463                IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
464                IF ( TRIM( var ) == 'theta' )  unit = 'K'
465                IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
466                IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
467                IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
468                CONTINUE
469
470             CASE DEFAULT
471!
472!--             Allocate arrays for other modules
473                CALL module_interface_init_masks( var, unit )
474
475                IF ( unit == 'illegal' )  THEN
476                   IF ( do_mask_user(mid,1) /= ' ' )  THEN
477                      WRITE ( message_string, * ) 'illegal value for data_',                       &
478                                                  'output_masks or data_output_masks_user: "',     &
479                                                  TRIM( do_mask(mid,i) ), '"'
480                      CALL message( 'init_masks', 'PA0018', 1, 2, 0, 6, 0 )
481                   ELSE
482                      WRITE ( message_string, * ) 'illegal value for data_',                       &
483                                                  ' output_masks : "', TRIM( do_mask(mid,i) ), '"'
484                      CALL message( 'init_masks', 'PA0330', 1, 2, 0, 6, 0 )
485                   ENDIF
486                ENDIF
487
488          END SELECT
489!
490!--       Set the internal steering parameters appropriately
491          domask_no(mid,j)                    = domask_no(mid,j) + 1
492          domask(mid,j,domask_no(mid,j))      = do_mask(mid,i)
493          domask_unit(mid,j,domask_no(mid,j)) = unit
494
495          IF ( j == 1 )  THEN
496!
497!--          Check, if variable is already subject to averaging
498             found = .FALSE.
499             DO  k = 1, doav_n
500                IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
501             ENDDO
502
503             IF ( .NOT. found )  THEN
504                doav_n = doav_n + 1
505                doav(doav_n) = var
506             ENDIF
507          ENDIF
508
509          i = i + 1
510
511       ENDDO   ! do_mask(mid,i)
512    ENDDO   ! mid
513
514
515!
516!-- Determine mask locations for each mask
517    DO  mid = 1, masks
518!
519!--    Set local masks for each subdomain along all three dimensions
520       CALL set_mask_locations( 1, dx, 'dx', nx, 'nx', nxl, nxr )
521       CALL set_mask_locations( 2, dy, 'dy', ny, 'ny', nys, nyn )
522       IF ( .NOT. mask_surface(mid) )  THEN
523          CALL set_mask_locations( 3, dz(1), 'dz', nz, 'nz', nzb, nzt )
524       ELSE
525!
526!--       Set vertical mask locations and size in case of terrain-following output
527          count = 0
528          DO  WHILE ( mask_k_over_surface(mid, count+1) >= 0 )
529             m = mask_k_over_surface(mid, count+1)
530             IF ( m < 1 )  THEN
531                WRITE ( message_string, '(I3,A,I3,A)' )  m, ' in mask ', mid, ' must be > 0.'
532                CALL message( 'init_masks', 'PA0743', 1, 2, 0, 6, 0 )
533             ENDIF
534             IF ( m > nz+1 )  THEN
535                WRITE ( message_string, '(I3,A,I3,A,I1,3A,I3)' )  m, ' in mask ', mid,             &
536                                                                  ' along dimension ', 3,          &
537                                                                  ' exceeds (nz+1) = ', nz+1
538                CALL message( 'init_masks', 'PA0331', 1, 2, 0, 6, 0 )
539             ENDIF
540             count = count + 1
541             mask_k(mid,count) = mask_k_over_surface(mid, count)
542             IF ( count == mask_xyz_dimension )  EXIT
543          ENDDO
544          mask_start_l(mid,3) = 1
545          mask_size(mid,3)    = count
546          mask_size_l(mid,3)  = count
547       ENDIF
548!
549!--    Set global masks along all three dimensions (required by define_netcdf_header).
550#if defined( __parallel )
551!
552!--    PE0 receives partial arrays from all processors of the respective mask and outputs them. Here
553!--    a barrier has to be set, because otherwise "-MPI- FATAL: Remote protocol queue full" may
554!--    occur.
555
556       CALL MPI_BARRIER( comm2d, ierr )
557
558       IF ( myid == 0 )  THEN
559!
560!--       Local arrays can be relocated directly.
561          mask_i_global(mid,mask_start_l(mid,1): &
562                       mask_start_l(mid,1)+mask_size_l(mid,1)-1) = &
563                       mask_i(mid,:mask_size_l(mid,1))
564          mask_j_global(mid,mask_start_l(mid,2): &
565                       mask_start_l(mid,2)+mask_size_l(mid,2)-1) = &
566                       mask_j(mid,:mask_size_l(mid,2))
567          mask_k_global(mid,mask_start_l(mid,3): &
568                       mask_start_l(mid,3)+mask_size_l(mid,3)-1) = &
569                       mask_k(mid,:mask_size_l(mid,3))
570!
571!--       Receive data from all other PEs.
572          DO  n = 1, numprocs-1
573!
574!--          Receive index limits first, then arrays.
575!--          Index limits are received in arbitrary order from the PEs.
576             CALL MPI_RECV( ind(1), 6, MPI_INTEGER, MPI_ANY_SOURCE, 0, comm2d, status, ierr )
577!
578!--          Not all PEs have data for the mask.
579             IF ( ind(1) /= -9999 )  THEN
580                sender = status(MPI_SOURCE)
581                CALL MPI_RECV( tmp_array(ind(1)), ind(2)-ind(1)+1, MPI_INTEGER, sender, 1, comm2d, &
582                               status, ierr )
583                mask_i_global(mid,ind(1):ind(2)) = tmp_array(ind(1):ind(2))
584                CALL MPI_RECV( tmp_array(ind(3)), ind(4)-ind(3)+1, MPI_INTEGER, sender, 2, comm2d, &
585                               status, ierr )
586                mask_j_global(mid,ind(3):ind(4)) = tmp_array(ind(3):ind(4))
587                CALL MPI_RECV( tmp_array(ind(5)), ind(6)-ind(5)+1, MPI_INTEGER, sender, 3, comm2d, &
588                               status, ierr )
589                mask_k_global(mid,ind(5):ind(6)) = tmp_array(ind(5):ind(6))
590             ENDIF
591          ENDDO
592
593       ELSE
594!
595!--       If at least part of the mask resides on the PE, send the index limits for the target
596!--       array, otherwise send -9999 to PE0.
597          IF ( mask_size_l(mid,1) > 0  .AND.  mask_size_l(mid,2) > 0  .AND.                        &
598               mask_size_l(mid,3) > 0  )  THEN
599             ind(1) = mask_start_l(mid,1)
600             ind(2) = mask_start_l(mid,1) + mask_size_l(mid,1) - 1
601             ind(3) = mask_start_l(mid,2)
602             ind(4) = mask_start_l(mid,2) + mask_size_l(mid,2) - 1
603             ind(5) = mask_start_l(mid,3)
604             ind(6) = mask_start_l(mid,3) + mask_size_l(mid,3) - 1
605          ELSE
606             ind(1) = -9999; ind(2) = -9999
607             ind(3) = -9999; ind(4) = -9999
608             ind(5) = -9999; ind(6) = -9999
609          ENDIF
610          CALL MPI_SEND( ind(1), 6, MPI_INTEGER, 0, 0, comm2d, ierr )
611!
612!--       If applicable, send data to PE0.
613          IF ( ind(1) /= -9999 )  THEN
614             tmp_array(:mask_size_l(mid,1)) = mask_i(mid,:mask_size_l(mid,1))
615             CALL MPI_SEND( tmp_array(1), mask_size_l(mid,1), MPI_INTEGER, 0, 1, comm2d, ierr )
616             tmp_array(:mask_size_l(mid,2)) = mask_j(mid,:mask_size_l(mid,2))
617             CALL MPI_SEND( tmp_array(1), mask_size_l(mid,2), MPI_INTEGER, 0, 2, comm2d, ierr )
618             tmp_array(:mask_size_l(mid,3)) = mask_k(mid,:mask_size_l(mid,3))
619             CALL MPI_SEND( tmp_array(1), mask_size_l(mid,3), MPI_INTEGER, 0, 3, comm2d, ierr )
620          ENDIF
621       ENDIF
622!
623!--    A barrier has to be set, because otherwise some PEs may proceed too fast so that PE0 may
624!--    receive wrong data on tag 0.
625       CALL MPI_BARRIER( comm2d, ierr )
626
627       IF ( netcdf_data_format > 4 )  THEN
628
629          CALL MPI_BCAST( mask_i_global(mid,:), nx+2, MPI_INTEGER, 0, comm2d, ierr )
630          CALL MPI_BCAST( mask_j_global(mid,:), ny+2, MPI_INTEGER, 0, comm2d, ierr )
631          CALL MPI_BCAST( mask_k_global(mid,:), nz+2, MPI_INTEGER, 0, comm2d, ierr )
632
633       ENDIF
634
635#else
636!
637!--    Local arrays can be relocated directly.
638       mask_i_global(mid,:) = mask_i(mid,:)
639       mask_j_global(mid,:) = mask_j(mid,:)
640       mask_k_global(mid,:) = mask_k(mid,:)
641#endif
642    ENDDO   ! mid
643
644    DEALLOCATE( tmp_array )
645!
646!-- Internal mask arrays cannot be deallocated on PE 0 because they are required for header output
647!-- on PE 0.
648    IF ( myid /= 0 )  DEALLOCATE( mask, mask_loop )
649
650 CONTAINS
651
652!--------------------------------------------------------------------------------------------------!
653! Description:
654! ------------
655!> Set local mask for each subdomain along 'dim' direction.
656!--------------------------------------------------------------------------------------------------!
657    SUBROUTINE set_mask_locations( dim, dxyz, dxyz_string, nxyz, nxyz_string, lb, ub )
658
659       IMPLICIT NONE
660
661       CHARACTER (LEN=2) ::  dxyz_string !<
662       CHARACTER (LEN=2) ::  nxyz_string !<
663
664       INTEGER(iwp)  ::  count       !<
665       INTEGER(iwp)  ::  count_l     !<
666       INTEGER(iwp)  ::  dim         !<
667       INTEGER(iwp)  ::  m           !<
668       INTEGER(iwp)  ::  loop_begin  !<
669       INTEGER(iwp)  ::  loop_end    !<
670       INTEGER(iwp)  ::  loop_stride !<
671       INTEGER(iwp)  ::  lb          !<
672       INTEGER(iwp)  ::  nxyz        !<
673       INTEGER(iwp)  ::  ub          !<
674
675       REAL(wp)      ::  dxyz  !<
676       REAL(wp)      ::  ddxyz !<
677       REAL(wp)      ::  tmp1  !<
678       REAL(wp)      ::  tmp2  !<
679
680       count = 0;  count_l = 0
681       ddxyz = 1.0_wp / dxyz
682       tmp1  = 0.0_wp
683       tmp2  = 0.0_wp
684
685       IF ( mask(mid,dim,1) >= 0.0_wp )  THEN
686!
687!--       Use predefined mask_* array
688          DO  WHILE ( mask(mid,dim,count+1) >= 0.0_wp )
689             count = count + 1
690             IF ( dim == 1  .OR.  dim == 2 )  THEN
691                m = NINT( mask(mid,dim,count) * mask_scale(dim) * ddxyz - 0.5_wp )
692                IF ( m < 0 )  m = 0  ! avoid negative values
693             ELSEIF ( dim == 3 )  THEN
694                ind_array =  &
695                     MINLOC( ABS( mask(mid,dim,count) * mask_scale(dim) - zu ) )
696                m = ind_array(1) - 1 + nzb  ! MINLOC uses lower array bound 1
697             ENDIF
698             IF ( m > (nxyz+1) )  THEN
699                WRITE ( message_string, '(I3,A,I3,A,I1,3A,I3)' )  m, ' in mask ', mid,             &
700                                                                  ' along dimension ' ,dim,        &
701                                                                  ' exceeds (' ,nxyz_string,       &
702                                                                  '+1) = ', nxyz+1
703                CALL message( 'init_masks', 'PA0331', 1, 2, 0, 6, 0 )
704             ENDIF
705             IF ( ( m >= lb  .AND.  m <= ub )  .OR.  ( m == (nxyz+1)  .AND.  ub == nxyz )  )  THEN
706                IF ( count_l == 0 )  mask_start_l(mid,dim) = count
707                count_l = count_l + 1
708                IF ( dim == 1 )  THEN
709                   mask_i(mid,count_l) = m
710                ELSEIF ( dim == 2 )  THEN
711                   mask_j(mid,count_l) = m
712                ELSEIF ( dim == 3 )  THEN
713                   mask_k(mid,count_l) = m
714                ENDIF
715             ENDIF
716             IF ( count == mask_xyz_dimension )  EXIT
717          ENDDO
718          mask_size(mid,dim)   = count
719          mask_size_l(mid,dim) = count_l
720
721       ELSE
722!
723!--       Use predefined mask_loop_* array, or use the default (all grid points along this
724!--       direction)
725          IF ( mask_loop(mid,dim,1) < 0.0_wp )  THEN
726             tmp1 = mask_loop(mid,dim,1)
727             mask_loop(mid,dim,1) = zw(nzb)  !   lowest level  (default)
728          ENDIF
729          IF ( dim == 1 .OR. dim == 2 )  THEN
730             IF ( mask_loop(mid,dim,2) < 0.0_wp )  THEN
731                tmp2 = mask_loop(mid,dim,2)
732                mask_loop(mid,dim,2) = nxyz*dxyz / mask_scale(dim)   ! (default)
733             ENDIF
734             IF ( MAXVAL( mask_loop(mid,dim,1:2) )  &
735                  > (nxyz+1) * dxyz / mask_scale(dim) )  THEN
736                WRITE ( message_string, '(2(A,I3,A,I1,A,F9.3),5A,I1,A,F9.3)' )                     &
737                     'mask_loop(', mid, ',', dim, ',1)=', mask_loop(mid,dim,1),                    &
738                     ' and/or mask_loop(', mid, ',', dim, ',2)=', mask_loop(mid,dim,2),            &
739                     ' exceed (', nxyz_string,'+1)*',dxyz_string,'/mask_scale(',dim,')=',          &
740                     (nxyz+1)*dxyz/mask_scale(dim)
741                CALL message( 'init_masks', 'PA0332', 1, 2, 0, 6, 0 )
742             ENDIF
743             loop_begin  = NINT( mask_loop(mid,dim,1) * mask_scale(dim) * ddxyz - 0.5_wp )
744             loop_end    = NINT( mask_loop(mid,dim,2) * mask_scale(dim) * ddxyz - 0.5_wp )
745             loop_stride = NINT( mask_loop(mid,dim,3) * mask_scale(dim) * ddxyz )
746             IF ( loop_begin == -1 )  loop_begin = 0  ! avoid negative values
747          ELSEIF ( dim == 3 )  THEN
748             IF ( mask_loop(mid,dim,2) < 0.0_wp )  THEN
749                tmp2 = mask_loop(mid,dim,2)
750                mask_loop(mid,dim,2) = zu(nz+1) / mask_scale(dim)   ! (default)
751             ENDIF
752             IF ( MAXVAL( mask_loop(mid,dim,1:2) ) > zu(nz+1) / mask_scale(dim) )  THEN
753                WRITE ( message_string, '(2(A,I3,A,I1,A,F9.3),A,I1,A,F9.3)' )                      &
754                      'mask_loop(', mid, ',', dim, ',1)=', mask_loop(mid,dim,1),                   &
755                      ' and/or mask_loop(', mid, ',', dim, ',2)=', mask_loop(mid,dim,2),           &
756                      ' exceed zu(nz+1)/mask_scale(', dim, ')=',zu(nz+1)/mask_scale(dim)
757                CALL message( 'init_masks', 'PA0333', 1, 2, 0, 6, 0 )
758             ENDIF
759             ind_array  = MINLOC( ABS( mask_loop(mid,dim,1) * mask_scale(dim) - zu ) )
760             loop_begin = ind_array(1) - 1 + nzb ! MINLOC uses lower array bound 1
761             ind_array  = MINLOC( ABS( mask_loop(mid,dim,2) * mask_scale(dim) - zu ) )
762             loop_end   = ind_array(1) - 1 + nzb ! MINLOC uses lower array bound 1
763!
764!--          The following line assumes a constant vertical grid spacing within the vertical mask
765!--          range; it fails for vertical grid stretching.
766!--          Maybe revise later. Issue warning but continue execution. ABS(...) within the IF
767!--          statement is necessary because the default value of dz_stretch_level_start is
768!--          -9999999.9_wp.
769             loop_stride = NINT( mask_loop(mid,dim,3) * mask_scale(dim) * ddxyz )
770
771             IF ( mask_loop(mid,dim,2) * mask_scale(dim) > ABS( dz_stretch_level_start(1) ) )  THEN
772                WRITE ( message_string, '(A,I3,A,I1,A,F9.3,A,F8.2,3A)' )                           &
773                      'mask_loop(', mid, ',', dim, ',2)=', mask_loop(mid,dim,2),                   &
774                      ' exceeds dz_stretch_level=', dz_stretch_level_start(1),                     &
775                      '.&Vertical mask locations will not ',                                       &
776                      'match the desired heights within the stretching ', 'region.'
777                CALL message( 'init_masks', 'PA0334', 0, 1, 0, 6, 0 )
778             ENDIF
779
780          ENDIF
781!
782!--       If necessary, reset mask_loop(mid,dim,1) and mask_loop(mid,dim,2).
783          IF ( tmp1 < 0.0_wp )  mask_loop(mid,dim,1) = tmp1
784          IF ( tmp2 < 0.0_wp )  mask_loop(mid,dim,2) = tmp2
785!
786!--       The default stride +/-1 (every grid point) applies if mask_loop(mid,dim,3) is not
787!--       specified (its default is zero).
788          IF ( loop_stride == 0 )  THEN
789             IF ( loop_end >= loop_begin )  THEN
790                loop_stride =  1
791             ELSE
792                loop_stride = -1
793             ENDIF
794          ENDIF
795          DO  m = loop_begin, loop_end, loop_stride
796             count = count + 1
797             IF ( ( m >= lb  .AND.  m <= ub ) .OR.  ( m == (nxyz+1) .AND. ub == nxyz )  )  THEN
798                IF ( count_l == 0 )  mask_start_l(mid,dim) = count
799                count_l = count_l + 1
800                IF ( dim == 1 )  THEN
801                   mask_i(mid,count_l) = m
802                ELSEIF ( dim == 2 )  THEN
803                   mask_j(mid,count_l) = m
804                ELSEIF ( dim == 3 )  THEN
805                   mask_k(mid,count_l) = m
806                ENDIF
807             ENDIF
808          ENDDO
809          mask_size(mid,dim)   = count
810          mask_size_l(mid,dim) = count_l
811
812       ENDIF
813
814    END SUBROUTINE set_mask_locations
815
816 END SUBROUTINE init_masks
Note: See TracBrowser for help on using the repository browser.