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

Last change on this file since 1551 was 1551, checked in by maronga, 8 years ago

land surface model released

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