source: palm/trunk/SOURCE/data_output_3d.f90 @ 1560

Last change on this file since 1560 was 1552, checked in by maronga, 10 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 21.8 KB
Line 
1 SUBROUTINE data_output_3d( av )
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: data_output_3d.f90 1552 2015-03-03 14:27:15Z keck $
27!
28! 1551 2015-03-03 14:18:16Z maronga
29! Added suppport for land surface model and radiation model output. In the course
30! of this action, the limits for vertical loops have been changed (from nzb and
31! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
32! Moreover, a new vertical grid zs was introduced.
33!
34! 1359 2014-04-11 17:15:14Z hoffmann
35! New particle structure integrated.
36!
37! 1353 2014-04-08 15:21:23Z heinze
38! REAL constants provided with KIND-attribute
39!
40! 1327 2014-03-21 11:00:16Z raasch
41! parts concerning avs output removed,
42! -netcdf output queries
43!
44! 1320 2014-03-20 08:40:49Z raasch
45! ONLY-attribute added to USE-statements,
46! kind-parameters added to all INTEGER and REAL declaration statements,
47! kinds are defined in new module kinds,
48! old module precision_kind is removed,
49! revision history before 2012 removed,
50! comment fields (!:) to be used for variable explanations added to
51! all variable declaration statements
52!
53! 1318 2014-03-17 13:35:16Z raasch
54! barrier argument removed from cpu_log,
55! module interfaces removed
56!
57! 1308 2014-03-13 14:58:42Z fricke
58! Check, if the limit of the time dimension is exceeded for parallel output
59! To increase the performance for parallel output, the following is done:
60! - Update of time axis is only done by PE0
61!
62! 1244 2013-10-31 08:16:56Z raasch
63! Bugfix for index bounds in case of 3d-parallel output
64!
65! 1115 2013-03-26 18:16:16Z hoffmann
66! ql is calculated by calc_liquid_water_content
67!
68! 1106 2013-03-04 05:31:38Z raasch
69! array_kind renamed precision_kind
70!
71! 1076 2012-12-05 08:30:18Z hoffmann
72! Bugfix in output of ql
73!
74! 1053 2012-11-13 17:11:03Z hoffmann
75! +nr, qr, prr, qc and averaged quantities
76!
77! 1036 2012-10-22 13:43:42Z raasch
78! code put under GPL (PALM 3.9)
79!
80! 1031 2012-10-19 14:35:30Z raasch
81! netCDF4 without parallel file support implemented
82!
83! 1007 2012-09-19 14:30:36Z franke
84! Bugfix: missing calculation of ql_vp added
85!
86! Revision 1.1  1997/09/03 06:29:36  raasch
87! Initial revision
88!
89!
90! Description:
91! ------------
92! Output of the 3D-arrays in netCDF and/or AVS format.
93!------------------------------------------------------------------------------!
94
95    USE arrays_3d,                                                             &
96        ONLY:  e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho, sa, tend, u, v,   &
97               vpt, w
98       
99    USE averaging
100       
101    USE cloud_parameters,                                                      &
102        ONLY:  l_d_cp, prr, pt_d_t
103       
104    USE control_parameters,                                                    &
105        ONLY:  avs_data_file, cloud_physics, do3d, do3d_avs_n,                 &
106               do3d_no, do3d_time_count, io_blocks, io_group,                  &
107               message_string, netcdf_data_format, ntdim_3d,                   &
108               nz_do3d, plot_3d_precision, psolver, simulated_time,            &
109               simulated_time_chr, skip_do_avs, time_since_reference_point
110       
111    USE cpulog,                                                                &
112        ONLY:  log_point, cpu_log
113       
114    USE indices,                                                               &
115        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzt,  &
116               nzb
117       
118    USE kinds
119   
120    USE land_surface_model_mod,                                                &
121        ONLY: m_soil, m_soil_av, nzb_soil, nzt_soil, t_soil, t_soil_av
122
123    USE netcdf_control
124       
125    USE particle_attributes,                                                   &
126        ONLY:  grid_particles, number_of_particles, particles,                 &
127               particle_advection_start, prt_count
128       
129    USE pegrid
130
131    IMPLICIT NONE
132
133    CHARACTER (LEN=9) ::  simulated_time_mod  !:
134
135    INTEGER(iwp) ::  av        !:
136    INTEGER(iwp) ::  i         !:
137    INTEGER(iwp) ::  if        !:
138    INTEGER(iwp) ::  j         !:
139    INTEGER(iwp) ::  k         !:
140    INTEGER(iwp) ::  n         !:
141    INTEGER(iwp) ::  nzb_do    !: vertical lower limit for data output
142    INTEGER(iwp) ::  nzt_do    !: vertical upper limit for data output
143    INTEGER(iwp) ::  pos       !:
144    INTEGER(iwp) ::  prec      !:
145    INTEGER(iwp) ::  psi       !:
146
147    LOGICAL      ::  found     !:
148    LOGICAL      ::  resorted  !:
149
150    REAL(wp)     ::  mean_r    !:
151    REAL(wp)     ::  s_r2      !:
152    REAL(wp)     ::  s_r3      !:
153
154    REAL(sp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf  !:
155
156    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !:
157
158!
159!-- Return, if nothing to output
160    IF ( do3d_no(av) == 0 )  RETURN
161!
162!-- For parallel netcdf output the time axis must be limited. Return, if this
163!-- limit is exceeded. This could be the case, if the simulated time exceeds
164!-- the given end time by the length of the given output interval.
165    IF ( netcdf_data_format > 4 )  THEN
166       IF ( do3d_time_count(av) + 1 > ntdim_3d(av) )  THEN
167          WRITE ( message_string, * ) 'Output of 3d data is not given at t=',  &
168                                      simulated_time, '&because the maximum ', & 
169                                      'number of output time levels is ',      &
170                                      'exceeded.'
171          CALL message( 'data_output_3d', 'PA0387', 0, 1, 0, 6, 0 )         
172          RETURN
173       ENDIF
174    ENDIF
175
176    CALL cpu_log (log_point(14),'data_output_3d','start')
177
178!
179!-- Open output file.
180!-- Also creates coordinate and fld-file for AVS.
181!-- For classic or 64bit netCDF output or output of other (old) data formats,
182!-- for a run on more than one PE, each PE opens its own file and
183!-- writes the data of its subdomain in binary format (regardless of the format
184!-- the user has requested). After the run, these files are combined to one
185!-- file by combine_plot_fields in the format requested by the user (netcdf
186!-- and/or avs).
187!-- For netCDF4/HDF5 output, data is written in parallel into one file.
188    IF ( netcdf_data_format < 5 )  THEN
189       CALL check_open( 30 )
190       IF ( myid == 0 )  CALL check_open( 106+av*10 )
191    ELSE
192       CALL check_open( 106+av*10 )
193    ENDIF
194
195!
196!-- Update the netCDF time axis
197!-- In case of parallel output, this is only done by PE0 to increase the
198!-- performance.
199#if defined( __netcdf )
200    do3d_time_count(av) = do3d_time_count(av) + 1
201    IF ( myid == 0 )  THEN
202       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_time_3d(av),           &
203                               (/ time_since_reference_point /),            &
204                               start = (/ do3d_time_count(av) /),           &
205                               count = (/ 1 /) )
206       CALL handle_netcdf_error( 'data_output_3d', 376 )
207    ENDIF
208#endif
209
210!
211!-- Loop over all variables to be written.
212    if = 1
213
214    DO  WHILE ( do3d(av,if)(1:1) /= ' ' )
215!
216!--    Store the array chosen on the temporary array.
217       resorted = .FALSE.
218       nzb_do = nzb
219       nzt_do = nz_do3d
220
221!
222!--    Allocate a temporary array with the desired output dimensions.
223       ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
224
225       SELECT CASE ( TRIM( do3d(av,if) ) )
226
227          CASE ( 'e' )
228             IF ( av == 0 )  THEN
229                to_be_resorted => e
230             ELSE
231                to_be_resorted => e_av
232             ENDIF
233
234          CASE ( 'lpt' )
235             IF ( av == 0 )  THEN
236                to_be_resorted => pt
237             ELSE
238                to_be_resorted => lpt_av
239             ENDIF
240
241          CASE ( 'm_soil' )
242             nzb_do = nzb_soil
243             nzt_do = nzt_soil
244!
245!--          For soil model quantities, it is required to re-allocate local_pf
246             DEALLOCATE ( local_pf )
247             ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
248
249             IF ( av == 0 )  THEN
250                to_be_resorted => m_soil
251             ELSE
252                to_be_resorted => m_soil_av
253             ENDIF
254
255          CASE ( 'nr' )
256             IF ( av == 0 )  THEN
257                to_be_resorted => nr
258             ELSE
259                to_be_resorted => nr_av
260             ENDIF
261
262          CASE ( 'p' )
263             IF ( av == 0 )  THEN
264                IF ( psolver /= 'sor' )  CALL exchange_horiz( p, nbgp )
265                to_be_resorted => p
266             ELSE
267                IF ( psolver /= 'sor' )  CALL exchange_horiz( p_av, nbgp )
268                to_be_resorted => p_av
269             ENDIF
270
271          CASE ( 'pc' )  ! particle concentration (requires ghostpoint exchange)
272             IF ( av == 0 )  THEN
273                IF ( simulated_time >= particle_advection_start )  THEN
274                   tend = prt_count
275                   CALL exchange_horiz( tend, nbgp )
276                ELSE
277                   tend = 0.0_wp
278                ENDIF
279                DO  i = nxlg, nxrg
280                   DO  j = nysg, nyng
281                      DO  k = nzb_do, nzt_do
282                         local_pf(i,j,k) = tend(k,j,i)
283                      ENDDO
284                   ENDDO
285                ENDDO
286                resorted = .TRUE.
287             ELSE
288                CALL exchange_horiz( pc_av, nbgp )
289                to_be_resorted => pc_av
290             ENDIF
291
292          CASE ( 'pr' )  ! mean particle radius (effective radius)
293             IF ( av == 0 )  THEN
294                IF ( simulated_time >= particle_advection_start )  THEN
295                   DO  i = nxl, nxr
296                      DO  j = nys, nyn
297                         DO  k = nzb_do, nzt_do
298                            number_of_particles = prt_count(k,j,i)
299                            IF (number_of_particles <= 0)  CYCLE
300                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
301                            s_r2 = 0.0_wp
302                            s_r3 = 0.0_wp
303                            DO  n = 1, number_of_particles
304                               IF ( particles(n)%particle_mask )  THEN
305                                  s_r2 = s_r2 + particles(n)%radius**2 * &
306                                         particles(n)%weight_factor
307                                  s_r3 = s_r3 + particles(n)%radius**3 * &
308                                         particles(n)%weight_factor
309                               ENDIF
310                            ENDDO
311                            IF ( s_r2 > 0.0_wp )  THEN
312                               mean_r = s_r3 / s_r2
313                            ELSE
314                               mean_r = 0.0_wp
315                            ENDIF
316                            tend(k,j,i) = mean_r
317                         ENDDO
318                      ENDDO
319                   ENDDO
320                   CALL exchange_horiz( tend, nbgp )
321                ELSE
322                   tend = 0.0_wp
323                ENDIF
324                DO  i = nxlg, nxrg
325                   DO  j = nysg, nyng
326                      DO  k = nzb_do, nzt_do
327                         local_pf(i,j,k) = tend(k,j,i)
328                      ENDDO
329                   ENDDO
330                ENDDO
331                resorted = .TRUE.
332             ELSE
333                CALL exchange_horiz( pr_av, nbgp )
334                to_be_resorted => pr_av
335             ENDIF
336
337          CASE ( 'prr' )
338             IF ( av == 0 )  THEN
339                CALL exchange_horiz( prr, nbgp )
340                DO  i = nxlg, nxrg
341                   DO  j = nysg, nyng
342                      DO  k = nzb, nzt+1
343                         local_pf(i,j,k) = prr(k,j,i)
344                      ENDDO
345                   ENDDO
346                ENDDO
347             ELSE
348                CALL exchange_horiz( prr_av, nbgp )
349                DO  i = nxlg, nxrg
350                   DO  j = nysg, nyng
351                      DO  k = nzb, nzt+1
352                         local_pf(i,j,k) = prr_av(k,j,i)
353                      ENDDO
354                   ENDDO
355                ENDDO
356             ENDIF
357             resorted = .TRUE.
358
359          CASE ( 'pt' )
360             IF ( av == 0 )  THEN
361                IF ( .NOT. cloud_physics ) THEN
362                   to_be_resorted => pt
363                ELSE
364                   DO  i = nxlg, nxrg
365                      DO  j = nysg, nyng
366                         DO  k = nzb_do, nzt_do
367                            local_pf(i,j,k) = pt(k,j,i) + l_d_cp *             &
368                                                          pt_d_t(k) *          &
369                                                          ql(k,j,i)
370                         ENDDO
371                      ENDDO
372                   ENDDO
373                   resorted = .TRUE.
374                ENDIF
375             ELSE
376                to_be_resorted => pt_av
377             ENDIF
378
379          CASE ( 'q' )
380             IF ( av == 0 )  THEN
381                to_be_resorted => q
382             ELSE
383                to_be_resorted => q_av
384             ENDIF
385
386          CASE ( 'qc' )
387             IF ( av == 0 )  THEN
388                to_be_resorted => qc
389             ELSE
390                to_be_resorted => qc_av
391             ENDIF
392
393          CASE ( 'ql' )
394             IF ( av == 0 )  THEN
395                to_be_resorted => ql
396             ELSE
397                to_be_resorted => ql_av
398             ENDIF
399
400          CASE ( 'ql_c' )
401             IF ( av == 0 )  THEN
402                to_be_resorted => ql_c
403             ELSE
404                to_be_resorted => ql_c_av
405             ENDIF
406
407          CASE ( 'ql_v' )
408             IF ( av == 0 )  THEN
409                to_be_resorted => ql_v
410             ELSE
411                to_be_resorted => ql_v_av
412             ENDIF
413
414          CASE ( 'ql_vp' )
415             IF ( av == 0 )  THEN
416                IF ( simulated_time >= particle_advection_start )  THEN
417                   DO  i = nxl, nxr
418                      DO  j = nys, nyn
419                         DO  k = nzb_do, nzt_do
420                            number_of_particles = prt_count(k,j,i)
421                            IF (number_of_particles <= 0)  CYCLE
422                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
423                            DO  n = 1, number_of_particles
424                               IF ( particles(n)%particle_mask )  THEN
425                                  tend(k,j,i) =  tend(k,j,i) +                 &
426                                                 particles(n)%weight_factor /  &
427                                                 prt_count(k,j,i)
428                               ENDIF
429                            ENDDO
430                         ENDDO
431                      ENDDO
432                   ENDDO
433                   CALL exchange_horiz( tend, nbgp )
434                ELSE
435                   tend = 0.0_wp
436                ENDIF
437                DO  i = nxlg, nxrg
438                   DO  j = nysg, nyng
439                      DO  k = nzb_do, nzt_do
440                         local_pf(i,j,k) = tend(k,j,i)
441                      ENDDO
442                   ENDDO
443                ENDDO
444                resorted = .TRUE.
445             ELSE
446                CALL exchange_horiz( ql_vp_av, nbgp )
447                to_be_resorted => ql_vp_av
448             ENDIF
449
450          CASE ( 'qr' )
451             IF ( av == 0 )  THEN
452                to_be_resorted => qr
453             ELSE
454                to_be_resorted => qr_av
455             ENDIF
456
457          CASE ( 'qv' )
458             IF ( av == 0 )  THEN
459                DO  i = nxlg, nxrg
460                   DO  j = nysg, nyng
461                      DO  k = nzb_do, nzt_do
462                         local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
463                      ENDDO
464                   ENDDO
465                ENDDO
466                resorted = .TRUE.
467             ELSE
468                to_be_resorted => qv_av
469             ENDIF
470
471          CASE ( 'rho' )
472             IF ( av == 0 )  THEN
473                to_be_resorted => rho
474             ELSE
475                to_be_resorted => rho_av
476             ENDIF
477
478          CASE ( 's' )
479             IF ( av == 0 )  THEN
480                to_be_resorted => q
481             ELSE
482                to_be_resorted => s_av
483             ENDIF
484
485          CASE ( 'sa' )
486             IF ( av == 0 )  THEN
487                to_be_resorted => sa
488             ELSE
489                to_be_resorted => sa_av
490             ENDIF
491
492          CASE ( 't_soil' )
493             nzb_do = nzb_soil
494             nzt_do = nzt_soil
495!
496!--          For soil model quantities, it is required to re-allocate local_pf
497             DEALLOCATE ( local_pf )
498             ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )
499
500             IF ( av == 0 )  THEN
501                to_be_resorted => t_soil
502             ELSE
503                to_be_resorted => t_soil_av
504             ENDIF
505
506          CASE ( 'u' )
507             IF ( av == 0 )  THEN
508                to_be_resorted => u
509             ELSE
510                to_be_resorted => u_av
511             ENDIF
512
513          CASE ( 'v' )
514             IF ( av == 0 )  THEN
515                to_be_resorted => v
516             ELSE
517                to_be_resorted => v_av
518             ENDIF
519
520          CASE ( 'vpt' )
521             IF ( av == 0 )  THEN
522                to_be_resorted => vpt
523             ELSE
524                to_be_resorted => vpt_av
525             ENDIF
526
527          CASE ( 'w' )
528             IF ( av == 0 )  THEN
529                to_be_resorted => w
530             ELSE
531                to_be_resorted => w_av
532             ENDIF
533
534          CASE DEFAULT
535!
536!--          User defined quantity
537             CALL user_data_output_3d( av, do3d(av,if), found, local_pf,       &
538                                       nzb_do, nzt_do )
539             resorted = .TRUE.
540
541             IF ( .NOT. found )  THEN
542                message_string =  'no output available for: ' //               &
543                                  TRIM( do3d(av,if) )
544                CALL message( 'data_output_3d', 'PA0182', 0, 0, 0, 6, 0 )
545             ENDIF
546
547       END SELECT
548
549!
550!--    Resort the array to be output, if not done above
551       IF ( .NOT. resorted )  THEN
552          DO  i = nxlg, nxrg
553             DO  j = nysg, nyng
554                DO  k = nzb_do, nzt_do
555                   local_pf(i,j,k) = to_be_resorted(k,j,i)
556                ENDDO
557             ENDDO
558          ENDDO
559       ENDIF
560
561!
562!--    Output of the 3D-array
563#if defined( __parallel )
564       IF ( netcdf_data_format < 5 )  THEN
565!
566!--       Non-parallel netCDF output. Data is output in parallel in
567!--       FORTRAN binary format here, and later collected into one file by
568!--       combine_plot_fields
569          IF ( myid == 0 )  THEN
570             WRITE ( 30 )  time_since_reference_point,                   &
571                           do3d_time_count(av), av
572          ENDIF
573          DO  i = 0, io_blocks-1
574             IF ( i == io_group )  THEN
575                WRITE ( 30 )  nxlg, nxrg, nysg, nyng, nzb_do, nzt_do
576                WRITE ( 30 )  local_pf(:,:,nzb_do:nzt_do)
577             ENDIF
578#if defined( __parallel )
579             CALL MPI_BARRIER( comm2d, ierr )
580#endif
581          ENDDO
582
583       ELSE
584#if defined( __netcdf )
585!
586!--       Parallel output in netCDF4/HDF5 format.
587!--       Do not output redundant ghost point data except for the
588!--       boundaries of the total domain.
589          IF ( nxr == nx  .AND.  nyn /= ny )  THEN
590             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
591                               local_pf(nxl:nxr+1,nys:nyn,nzb_do:nzt_do),    &
592                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
593                count = (/ nxr-nxl+2, nyn-nys+1, nzt_do-nzb_do+1, 1 /) )
594          ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
595             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
596                               local_pf(nxl:nxr,nys:nyn+1,nzb_do:nzt_do),    &
597                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
598                count = (/ nxr-nxl+1, nyn-nys+2, nzt_do-nzb_do+1, 1 /) )
599          ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
600             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
601                             local_pf(nxl:nxr+1,nys:nyn+1,nzb_do:nzt_do  ),  &
602                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
603                count = (/ nxr-nxl+2, nyn-nys+2, nzt_do-nzb_do+1, 1 /) )
604          ELSE
605             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
606                                 local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do),    &
607                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
608                count = (/ nxr-nxl+1, nyn-nys+1, nzt_do-nzb_do+1, 1 /) )
609          ENDIF
610          CALL handle_netcdf_error( 'data_output_3d', 386 )
611#endif
612       ENDIF
613#else
614#if defined( __netcdf )
615       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),        &
616                         local_pf(nxl:nxr+1,nys:nyn+1,nzb_do:nzt_do),        &
617                         start = (/ 1, 1, 1, do3d_time_count(av) /),     &
618                         count = (/ nx+2, ny+2, nzt_do-nzb_do+1, 1 /) )
619       CALL handle_netcdf_error( 'data_output_3d', 446 )
620#endif
621#endif
622
623       if = if + 1
624
625!
626!--    Deallocate temporary array
627       DEALLOCATE ( local_pf )
628
629    ENDDO
630
631    CALL cpu_log( log_point(14), 'data_output_3d', 'stop' )
632
633!
634!-- Formats.
6353300 FORMAT ('variable ',I4,'  file=',A,'  filetype=unformatted  skip=',I12/   &
636             'label = ',A,A)
637
638 END SUBROUTINE data_output_3d
Note: See TracBrowser for help on using the repository browser.