SUBROUTINE data_output_2d( mode, av ) !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! simulated_time in NetCDF output replaced by time_since_reference_point. ! Output of NetCDF messages with aid of message handling routine. ! Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) ! Output of messages replaced by message handling routine. ! Output of user defined 2D (XY) arrays at z=nzb+1 is now possible ! ! ! Former revisions: ! ----------------- ! $Id: data_output_2d.f90 343 2009-06-24 12:59:09Z raasch $ ! ! 215 2008-11-18 09:54:31Z raasch ! Bugfix: no output of particle concentration and radius unless particles ! have been started ! ! 96 2007-06-04 08:07:41Z raasch ! Output of density and salinity ! ! 75 2007-03-22 09:54:05Z raasch ! Output of precipitation amount/rate and roughness length, ! 2nd+3rd argument removed from exchange horiz ! ! RCS Log replace by Id keyword, revision history cleaned up ! ! Revision 1.5 2006/08/22 13:50:29 raasch ! xz and yz cross sections now up to nzt+1 ! ! Revision 1.2 2006/02/23 10:19:22 raasch ! Output of time-averaged data, output of averages along x, y, or z, ! output of user-defined quantities, ! section data are copied from local_pf to local_2d before they are output, ! output of particle concentration and mean radius, ! Former subroutine plot_2d renamed data_output_2d, pl2d.. renamed do2d.., ! anz renamed ngp, ebene renamed section, pl2d_.._anz renamed do2d_.._n ! ! Revision 1.1 1997/08/11 06:24:09 raasch ! Initial revision ! ! ! Description: ! ------------ ! Data output of horizontal cross-sections in NetCDF format or binary format ! compatible to old graphic software iso2d. ! Attention: The position of the sectional planes is still not always computed ! --------- correctly. (zu is used always)! !------------------------------------------------------------------------------! USE arrays_3d USE averaging USE cloud_parameters USE control_parameters USE cpulog USE grid_variables USE indices USE interfaces USE netcdf_control USE particle_attributes USE pegrid IMPLICIT NONE CHARACTER (LEN=2) :: do2d_mode, mode CHARACTER (LEN=4) :: grid CHARACTER (LEN=25) :: section_chr CHARACTER (LEN=50) :: rtext INTEGER :: av, ngp, file_id, i, if, is, j, k, l, layer_xy, n, psi, s, & sender, & ind(4) LOGICAL :: found, resorted, two_d REAL :: mean_r, s_r3, s_r4 REAL, DIMENSION(:), ALLOCATABLE :: level_z REAL, DIMENSION(:,:), ALLOCATABLE :: local_2d, local_2d_l REAL, DIMENSION(:,:,:), ALLOCATABLE :: local_pf #if defined( __parallel ) REAL, DIMENSION(:,:), ALLOCATABLE :: total_2d #endif REAL, DIMENSION(:,:,:), POINTER :: to_be_resorted NAMELIST /LOCAL/ rtext CALL cpu_log (log_point(3),'data_output_2d','start') ! !-- Immediate return, if no output is requested (no respective sections !-- found in parameter data_output) IF ( mode == 'xy' .AND. .NOT. data_output_xy(av) ) RETURN IF ( mode == 'xz' .AND. .NOT. data_output_xz(av) ) RETURN IF ( mode == 'yz' .AND. .NOT. data_output_yz(av) ) RETURN two_d = .FALSE. ! local variable to distinguish between output of pure 2D ! arrays and cross-sections of 3D arrays. ! !-- Depending on the orientation of the cross-section, the respective output !-- files have to be opened. SELECT CASE ( mode ) CASE ( 'xy' ) s = 1 ALLOCATE( level_z(0:nzt+1), local_2d(nxl-1:nxr+1,nys-1:nyn+1) ) #if defined( __netcdf ) IF ( myid == 0 .AND. netcdf_output ) CALL check_open( 101+av*10 ) #endif IF ( data_output_2d_on_each_pe ) THEN CALL check_open( 21 ) ELSE IF ( myid == 0 ) THEN IF ( iso2d_output ) CALL check_open( 21 ) #if defined( __parallel ) ALLOCATE( total_2d(-1:nx+1,-1:ny+1) ) #endif ENDIF ENDIF CASE ( 'xz' ) s = 2 ALLOCATE( local_2d(nxl-1:nxr+1,nzb:nzt+1) ) #if defined( __netcdf ) IF ( myid == 0 .AND. netcdf_output ) CALL check_open( 102+av*10 ) #endif IF ( data_output_2d_on_each_pe ) THEN CALL check_open( 22 ) ELSE IF ( myid == 0 ) THEN IF ( iso2d_output ) CALL check_open( 22 ) #if defined( __parallel ) ALLOCATE( total_2d(-1:nx+1,nzb:nzt+1) ) #endif ENDIF ENDIF CASE ( 'yz' ) s = 3 ALLOCATE( local_2d(nys-1:nyn+1,nzb:nzt+1) ) #if defined( __netcdf ) IF ( myid == 0 .AND. netcdf_output ) CALL check_open( 103+av*10 ) #endif IF ( data_output_2d_on_each_pe ) THEN CALL check_open( 23 ) ELSE IF ( myid == 0 ) THEN IF ( iso2d_output ) CALL check_open( 23 ) #if defined( __parallel ) ALLOCATE( total_2d(-1:ny+1,nzb:nzt+1) ) #endif ENDIF ENDIF CASE DEFAULT message_string = 'unknown cross-section: ' // TRIM( mode ) CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 ) END SELECT ! !-- Allocate a temporary array for resorting (kji -> ijk). ALLOCATE( local_pf(nxl-1:nxr+1,nys-1:nyn+1,nzb:nzt+1) ) ! !-- Loop of all variables to be written. !-- Output dimensions chosen if = 1 l = MAX( 2, LEN_TRIM( do2d(av,if) ) ) do2d_mode = do2d(av,if)(l-1:l) DO WHILE ( do2d(av,if)(1:1) /= ' ' ) IF ( do2d_mode == mode ) THEN ! !-- Store the array chosen on the temporary array. resorted = .FALSE. SELECT CASE ( TRIM( do2d(av,if) ) ) CASE ( 'e_xy', 'e_xz', 'e_yz' ) IF ( av == 0 ) THEN to_be_resorted => e ELSE to_be_resorted => e_av ENDIF IF ( mode == 'xy' ) level_z = zu CASE ( 'lwp*_xy' ) ! 2d-array IF ( av == 0 ) THEN DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) * & dzw(1:nzt+1) ) ENDDO ENDDO ELSE DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 local_pf(i,j,nzb+1) = lwp_av(j,i) ENDDO ENDDO ENDIF resorted = .TRUE. two_d = .TRUE. level_z(nzb+1) = zu(nzb+1) CASE ( 'p_xy', 'p_xz', 'p_yz' ) IF ( av == 0 ) THEN to_be_resorted => p ELSE to_be_resorted => p_av ENDIF IF ( mode == 'xy' ) level_z = zu CASE ( 'pc_xy', 'pc_xz', 'pc_yz' ) ! particle concentration IF ( av == 0 ) THEN IF ( simulated_time >= particle_advection_start ) THEN tend = prt_count CALL exchange_horiz( tend ) ELSE tend = 0.0 ENDIF DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 DO k = nzb, nzt+1 local_pf(i,j,k) = tend(k,j,i) ENDDO ENDDO ENDDO resorted = .TRUE. ELSE CALL exchange_horiz( pc_av ) to_be_resorted => pc_av ENDIF CASE ( 'pr_xy', 'pr_xz', 'pr_yz' ) ! mean particle radius IF ( av == 0 ) THEN IF ( simulated_time >= particle_advection_start ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 psi = prt_start_index(k,j,i) s_r3 = 0.0 s_r4 = 0.0 DO n = psi, psi+prt_count(k,j,i)-1 s_r3 = s_r3 + particles(n)%radius**3 s_r4 = s_r4 + particles(n)%radius**4 ENDDO IF ( s_r3 /= 0.0 ) THEN mean_r = s_r4 / s_r3 ELSE mean_r = 0.0 ENDIF tend(k,j,i) = mean_r ENDDO ENDDO ENDDO CALL exchange_horiz( tend ) ELSE tend = 0.0 ENDIF DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 DO k = nzb, nzt+1 local_pf(i,j,k) = tend(k,j,i) ENDDO ENDDO ENDDO resorted = .TRUE. ELSE CALL exchange_horiz( pr_av ) to_be_resorted => pr_av ENDIF CASE ( 'pra*_xy' ) ! 2d-array / integral quantity => no av CALL exchange_horiz_2d( precipitation_amount ) DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 local_pf(i,j,nzb+1) = precipitation_amount(j,i) ENDDO ENDDO precipitation_amount = 0.0 ! reset for next integ. interval resorted = .TRUE. two_d = .TRUE. level_z(nzb+1) = zu(nzb+1) CASE ( 'prr*_xy' ) ! 2d-array IF ( av == 0 ) THEN CALL exchange_horiz_2d( precipitation_rate ) DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 local_pf(i,j,nzb+1) = precipitation_rate(j,i) ENDDO ENDDO ELSE CALL exchange_horiz_2d( precipitation_rate_av ) DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 local_pf(i,j,nzb+1) = precipitation_rate_av(j,i) ENDDO ENDDO ENDIF resorted = .TRUE. two_d = .TRUE. level_z(nzb+1) = zu(nzb+1) CASE ( 'pt_xy', 'pt_xz', 'pt_yz' ) IF ( av == 0 ) THEN IF ( .NOT. cloud_physics ) THEN to_be_resorted => pt ELSE DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 DO k = nzb, nzt+1 local_pf(i,j,k) = pt(k,j,i) + l_d_cp * & pt_d_t(k) * & ql(k,j,i) ENDDO ENDDO ENDDO resorted = .TRUE. ENDIF ELSE to_be_resorted => pt_av ENDIF IF ( mode == 'xy' ) level_z = zu CASE ( 'q_xy', 'q_xz', 'q_yz' ) IF ( av == 0 ) THEN to_be_resorted => q ELSE to_be_resorted => q_av ENDIF IF ( mode == 'xy' ) level_z = zu CASE ( 'ql_xy', 'ql_xz', 'ql_yz' ) IF ( av == 0 ) THEN to_be_resorted => ql ELSE to_be_resorted => ql_av ENDIF IF ( mode == 'xy' ) level_z = zu CASE ( 'ql_c_xy', 'ql_c_xz', 'ql_c_yz' ) IF ( av == 0 ) THEN to_be_resorted => ql_c ELSE to_be_resorted => ql_c_av ENDIF IF ( mode == 'xy' ) level_z = zu CASE ( 'ql_v_xy', 'ql_v_xz', 'ql_v_yz' ) IF ( av == 0 ) THEN to_be_resorted => ql_v ELSE to_be_resorted => ql_v_av ENDIF IF ( mode == 'xy' ) level_z = zu CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' ) IF ( av == 0 ) THEN to_be_resorted => ql_vp ELSE to_be_resorted => ql_vp_av ENDIF IF ( mode == 'xy' ) level_z = zu CASE ( 'qv_xy', 'qv_xz', 'qv_yz' ) IF ( av == 0 ) THEN DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 DO k = nzb, nzt+1 local_pf(i,j,k) = q(k,j,i) - ql(k,j,i) ENDDO ENDDO ENDDO resorted = .TRUE. ELSE to_be_resorted => qv_av ENDIF IF ( mode == 'xy' ) level_z = zu CASE ( 'rho_xy', 'rho_xz', 'rho_yz' ) IF ( av == 0 ) THEN to_be_resorted => rho ELSE to_be_resorted => rho_av ENDIF CASE ( 's_xy', 's_xz', 's_yz' ) IF ( av == 0 ) THEN to_be_resorted => q ELSE to_be_resorted => q_av ENDIF CASE ( 'sa_xy', 'sa_xz', 'sa_yz' ) IF ( av == 0 ) THEN to_be_resorted => sa ELSE to_be_resorted => sa_av ENDIF CASE ( 't*_xy' ) ! 2d-array IF ( av == 0 ) THEN DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 local_pf(i,j,nzb+1) = ts(j,i) ENDDO ENDDO ELSE DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 local_pf(i,j,nzb+1) = ts_av(j,i) ENDDO ENDDO ENDIF resorted = .TRUE. two_d = .TRUE. level_z(nzb+1) = zu(nzb+1) CASE ( 'u_xy', 'u_xz', 'u_yz' ) IF ( av == 0 ) THEN to_be_resorted => u ELSE to_be_resorted => u_av ENDIF IF ( mode == 'xy' ) level_z = zu ! !-- Substitute the values generated by "mirror" boundary condition !-- at the bottom boundary by the real surface values. IF ( do2d(av,if) == 'u_xz' .OR. do2d(av,if) == 'u_yz' ) THEN IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0 ENDIF CASE ( 'u*_xy' ) ! 2d-array IF ( av == 0 ) THEN DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 local_pf(i,j,nzb+1) = us(j,i) ENDDO ENDDO ELSE DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 local_pf(i,j,nzb+1) = us_av(j,i) ENDDO ENDDO ENDIF resorted = .TRUE. two_d = .TRUE. level_z(nzb+1) = zu(nzb+1) CASE ( 'v_xy', 'v_xz', 'v_yz' ) IF ( av == 0 ) THEN to_be_resorted => v ELSE to_be_resorted => v_av ENDIF IF ( mode == 'xy' ) level_z = zu ! !-- Substitute the values generated by "mirror" boundary condition !-- at the bottom boundary by the real surface values. IF ( do2d(av,if) == 'v_xz' .OR. do2d(av,if) == 'v_yz' ) THEN IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0 ENDIF CASE ( 'vpt_xy', 'vpt_xz', 'vpt_yz' ) IF ( av == 0 ) THEN to_be_resorted => vpt ELSE to_be_resorted => vpt_av ENDIF IF ( mode == 'xy' ) level_z = zu CASE ( 'w_xy', 'w_xz', 'w_yz' ) IF ( av == 0 ) THEN to_be_resorted => w ELSE to_be_resorted => w_av ENDIF IF ( mode == 'xy' ) level_z = zw CASE ( 'z0*_xy' ) ! 2d-array IF ( av == 0 ) THEN DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 local_pf(i,j,nzb+1) = z0(j,i) ENDDO ENDDO ELSE DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 local_pf(i,j,nzb+1) = z0_av(j,i) ENDDO ENDDO ENDIF resorted = .TRUE. two_d = .TRUE. level_z(nzb+1) = zu(nzb+1) CASE DEFAULT ! !-- User defined quantity CALL user_data_output_2d( av, do2d(av,if), found, grid, & local_pf, two_d ) resorted = .TRUE. IF ( grid == 'zu' ) THEN IF ( mode == 'xy' ) level_z = zu ELSEIF ( grid == 'zw' ) THEN IF ( mode == 'xy' ) level_z = zw ELSEIF ( grid == 'zu1' ) THEN IF ( mode == 'xy' ) level_z(nzb+1) = zu(nzb+1) ENDIF IF ( .NOT. found ) THEN message_string = 'no output provided for: ' // & TRIM( do2d(av,if) ) CALL message( 'data_output_2d', 'PA0181', 0, 0, 0, 6, 0 ) ENDIF END SELECT ! !-- Resort the array to be output, if not done above IF ( .NOT. resorted ) THEN DO i = nxl-1, nxr+1 DO j = nys-1, nyn+1 DO k = nzb, nzt+1 local_pf(i,j,k) = to_be_resorted(k,j,i) ENDDO ENDDO ENDDO ENDIF ! !-- Output of the individual cross-sections, depending on the cross- !-- section mode chosen. is = 1 loop1: DO WHILE ( section(is,s) /= -9999 .OR. two_d ) SELECT CASE ( mode ) CASE ( 'xy' ) ! !-- Determine the cross section index IF ( two_d ) THEN layer_xy = nzb+1 ELSE layer_xy = section(is,s) ENDIF ! !-- Update the NetCDF xy cross section time axis IF ( myid == 0 ) THEN IF ( simulated_time /= do2d_xy_last_time(av) ) THEN do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1 do2d_xy_last_time(av) = simulated_time IF ( .NOT. data_output_2d_on_each_pe .AND. & netcdf_output ) THEN #if defined( __netcdf ) nc_stat = NF90_PUT_VAR( id_set_xy(av), & id_var_time_xy(av), & (/ time_since_reference_point /), & start = (/ do2d_xy_time_count(av) /), & count = (/ 1 /) ) CALL handle_netcdf_error( 'data_output_2d', 53 ) #endif ENDIF ENDIF ENDIF ! !-- If required, carry out averaging along z IF ( section(is,s) == -1 .AND. .NOT. two_d ) THEN local_2d = 0.0 ! !-- Carry out the averaging (all data are on the PE) DO k = nzb, nzt+1 DO j = nys-1, nyn+1 DO i = nxl-1, nxr+1 local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k) ENDDO ENDDO ENDDO local_2d = local_2d / ( nzt -nzb + 2.0 ) ELSE ! !-- Just store the respective section on the local array local_2d = local_pf(:,:,layer_xy) ENDIF #if defined( __parallel ) IF ( data_output_2d_on_each_pe ) THEN ! !-- Output of partial arrays on each PE #if defined( __netcdf ) IF ( netcdf_output .AND. myid == 0 ) THEN WRITE ( 21 ) simulated_time, do2d_xy_time_count(av), & av ENDIF #endif WRITE ( 21 ) nxl-1, nxr+1, nys-1, nyn+1 WRITE ( 21 ) local_2d ELSE ! !-- PE0 receives partial arrays from all processors and then !-- outputs them. Here a barrier has to be set, because !-- otherwise "-MPI- FATAL: Remote protocol queue full" may !-- occur. CALL MPI_BARRIER( comm2d, ierr ) ngp = ( nxr-nxl+3 ) * ( nyn-nys+3 ) IF ( myid == 0 ) THEN ! !-- Local array can be relocated directly. total_2d(nxl-1:nxr+1,nys-1:nyn+1) = local_2d ! !-- Receive data from all other PEs. DO n = 1, numprocs-1 ! !-- Receive index limits first, then array. !-- Index limits are received in arbitrary order from !-- the PEs. CALL MPI_RECV( ind(1), 4, MPI_INTEGER, & MPI_ANY_SOURCE, 0, comm2d, status, & ierr ) sender = status(MPI_SOURCE) DEALLOCATE( local_2d ) ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) ) CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, & MPI_REAL, sender, 1, comm2d, & status, ierr ) total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d ENDDO ! !-- Output of the total cross-section. IF ( iso2d_output ) WRITE (21) total_2d(0:nx+1,0:ny+1) ! !-- Relocate the local array for the next loop increment DEALLOCATE( local_2d ) ALLOCATE( local_2d(nxl-1:nxr+1,nys-1:nyn+1) ) #if defined( __netcdf ) IF ( netcdf_output ) THEN IF ( two_d ) THEN nc_stat = NF90_PUT_VAR( id_set_xy(av), & id_var_do2d(av,if), & total_2d(0:nx+1,0:ny+1), & start = (/ 1, 1, 1, do2d_xy_time_count(av) /), & count = (/ nx+2, ny+2, 1, 1 /) ) ELSE nc_stat = NF90_PUT_VAR( id_set_xy(av), & id_var_do2d(av,if), & total_2d(0:nx+1,0:ny+1), & start = (/ 1, 1, is, do2d_xy_time_count(av) /), & count = (/ nx+2, ny+2, 1, 1 /) ) ENDIF CALL handle_netcdf_error( 'data_output_2d', 54 ) ENDIF #endif ELSE ! !-- First send the local index limits to PE0 ind(1) = nxl-1; ind(2) = nxr+1 ind(3) = nys-1; ind(4) = nyn+1 CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, & ierr ) ! !-- Send data to PE0 CALL MPI_SEND( local_2d(nxl-1,nys-1), ngp, MPI_REAL, & 0, 1, comm2d, ierr ) ENDIF ! !-- A barrier has to be set, because otherwise some PEs may !-- proceed too fast so that PE0 may receive wrong data on !-- tag 0 CALL MPI_BARRIER( comm2d, ierr ) ENDIF #else IF ( iso2d_output ) THEN WRITE (21) local_2d(nxl:nxr+1,nys:nyn+1) ENDIF #if defined( __netcdf ) IF ( netcdf_output ) THEN IF ( two_d ) THEN nc_stat = NF90_PUT_VAR( id_set_xy(av), & id_var_do2d(av,if), & local_2d(nxl:nxr+1,nys:nyn+1), & start = (/ 1, 1, 1, do2d_xy_time_count(av) /), & count = (/ nx+2, ny+2, 1, 1 /) ) ELSE nc_stat = NF90_PUT_VAR( id_set_xy(av), & id_var_do2d(av,if), & local_2d(nxl:nxr+1,nys:nyn+1), & start = (/ 1, 1, is, do2d_xy_time_count(av) /), & count = (/ nx+2, ny+2, 1, 1 /) ) ENDIF CALL handle_netcdf_error( 'data_output_2d', 55 ) ENDIF #endif #endif do2d_xy_n = do2d_xy_n + 1 ! !-- Write LOCAL parameter set for ISO2D. IF ( myid == 0 .AND. iso2d_output ) THEN IF ( section(is,s) /= -1 ) THEN WRITE ( section_chr, '(''z = '',F7.2,'' m (GP '',I3, & &'')'')' & ) level_z(layer_xy), layer_xy ELSE section_chr = 'averaged along z' ENDIF IF ( av == 0 ) THEN rtext = TRIM( do2d(av,if) ) // ' t = ' // & TRIM( simulated_time_chr ) // ' ' // & TRIM( section_chr ) ELSE rtext = TRIM( do2d(av,if) ) // ' averaged t = ' // & TRIM( simulated_time_chr ) // ' ' // & TRIM( section_chr ) ENDIF WRITE (27,LOCAL) ENDIF ! !-- For 2D-arrays (e.g. u*) only one cross-section is available. !-- Hence exit loop of output levels. IF ( two_d ) THEN two_d = .FALSE. EXIT loop1 ENDIF CASE ( 'xz' ) ! !-- Update the NetCDF xz cross section time axis IF ( myid == 0 ) THEN IF ( simulated_time /= do2d_xz_last_time(av) ) THEN do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1 do2d_xz_last_time(av) = simulated_time IF ( .NOT. data_output_2d_on_each_pe .AND. & netcdf_output ) THEN #if defined( __netcdf ) nc_stat = NF90_PUT_VAR( id_set_xz(av), & id_var_time_xz(av), & (/ time_since_reference_point /), & start = (/ do2d_xz_time_count(av) /), & count = (/ 1 /) ) CALL handle_netcdf_error( 'data_output_2d', 56 ) #endif ENDIF ENDIF ENDIF ! !-- If required, carry out averaging along y IF ( section(is,s) == -1 ) THEN ALLOCATE( local_2d_l(nxl-1:nxr+1,nzb:nzt+1) ) local_2d_l = 0.0 ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 ) ! !-- First local averaging on the PE DO k = nzb, nzt+1 DO j = nys, nyn DO i = nxl-1, nxr+1 local_2d_l(i,k) = local_2d_l(i,k) + & local_pf(i,j,k) ENDDO ENDDO ENDDO #if defined( __parallel ) ! !-- Now do the averaging over all PEs along y CALL MPI_ALLREDUCE( local_2d_l(nxl-1,nzb), & local_2d(nxl-1,nzb), ngp, MPI_REAL, & MPI_SUM, comm1dy, ierr ) #else local_2d = local_2d_l #endif local_2d = local_2d / ( ny + 1.0 ) DEALLOCATE( local_2d_l ) ELSE ! !-- Just store the respective section on the local array !-- (but only if it is available on this PE!) IF ( section(is,s) >= nys .AND. section(is,s) <= nyn ) & THEN local_2d = local_pf(:,section(is,s),nzb:nzt+1) ENDIF ENDIF #if defined( __parallel ) IF ( data_output_2d_on_each_pe ) THEN ! !-- Output of partial arrays on each PE. If the cross section !-- does not reside on the PE, output special index values. #if defined( __netcdf ) IF ( netcdf_output .AND. myid == 0 ) THEN WRITE ( 22 ) simulated_time, do2d_xz_time_count(av), & av ENDIF #endif IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn ) .OR.& ( section(is,s) == -1 .AND. nys-1 == -1 ) ) & THEN WRITE (22) nxl-1, nxr+1, nzb, nzt+1 WRITE (22) local_2d ELSE WRITE (22) -1, -1, -1, -1 ENDIF ELSE ! !-- PE0 receives partial arrays from all processors of the !-- respective cross section and outputs them. Here a !-- barrier has to be set, because otherwise !-- "-MPI- FATAL: Remote protocol queue full" may occur. CALL MPI_BARRIER( comm2d, ierr ) ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 ) IF ( myid == 0 ) THEN ! !-- Local array can be relocated directly. IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn ) & .OR. ( section(is,s) == -1 .AND. nys-1 == -1 ) ) & THEN total_2d(nxl-1:nxr+1,nzb:nzt+1) = local_2d ENDIF ! !-- Receive data from all other PEs. DO n = 1, numprocs-1 ! !-- Receive index limits first, then array. !-- Index limits are received in arbitrary order from !-- the PEs. CALL MPI_RECV( ind(1), 4, MPI_INTEGER, & MPI_ANY_SOURCE, 0, comm2d, status, & ierr ) ! !-- Not all PEs have data for XZ-cross-section. IF ( ind(1) /= -9999 ) THEN sender = status(MPI_SOURCE) DEALLOCATE( local_2d ) ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) ) CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, & MPI_REAL, sender, 1, comm2d, & status, ierr ) total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d ENDIF ENDDO ! !-- Output of the total cross-section. IF ( iso2d_output ) THEN WRITE (22) total_2d(0:nx+1,nzb:nzt+1) ENDIF ! !-- Relocate the local array for the next loop increment DEALLOCATE( local_2d ) ALLOCATE( local_2d(nxl-1:nxr+1,nzb:nzt+1) ) #if defined( __netcdf ) IF ( netcdf_output ) THEN nc_stat = NF90_PUT_VAR( id_set_xz(av), & id_var_do2d(av,if), & total_2d(0:nx+1,nzb:nzt+1),& start = (/ 1, is, 1, do2d_xz_time_count(av) /), & count = (/ nx+2, 1, nz+2, 1 /) ) CALL handle_netcdf_error( 'data_output_2d', 57 ) ENDIF #endif ELSE ! !-- If the cross section resides on the PE, send the !-- local index limits, otherwise send -9999 to PE0. IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn ) & .OR. ( section(is,s) == -1 .AND. nys-1 == -1 ) ) & THEN ind(1) = nxl-1; ind(2) = nxr+1 ind(3) = nzb; ind(4) = nzt+1 ELSE ind(1) = -9999; ind(2) = -9999 ind(3) = -9999; ind(4) = -9999 ENDIF CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, & ierr ) ! !-- If applicable, send data to PE0. IF ( ind(1) /= -9999 ) THEN CALL MPI_SEND( local_2d(nxl-1,nzb), ngp, MPI_REAL, & 0, 1, comm2d, ierr ) ENDIF ENDIF ! !-- A barrier has to be set, because otherwise some PEs may !-- proceed too fast so that PE0 may receive wrong data on !-- tag 0 CALL MPI_BARRIER( comm2d, ierr ) ENDIF #else IF ( iso2d_output ) THEN WRITE (22) local_2d(nxl:nxr+1,nzb:nzt+1) ENDIF #if defined( __netcdf ) IF ( netcdf_output ) THEN nc_stat = NF90_PUT_VAR( id_set_xz(av), & id_var_do2d(av,if), & local_2d(nxl:nxr+1,nzb:nzt+1), & start = (/ 1, is, 1, do2d_xz_time_count(av) /), & count = (/ nx+2, 1, nz+2, 1 /) ) CALL handle_netcdf_error( 'data_output_2d', 58 ) ENDIF #endif #endif do2d_xz_n = do2d_xz_n + 1 ! !-- Write LOCAL-parameter set for ISO2D. IF ( myid == 0 .AND. iso2d_output ) THEN IF ( section(is,s) /= -1 ) THEN WRITE ( section_chr, '(''y = '',F8.2,'' m (GP '',I3, & &'')'')' & ) section(is,s)*dy, section(is,s) ELSE section_chr = 'averaged along y' ENDIF IF ( av == 0 ) THEN rtext = TRIM( do2d(av,if) ) // ' t = ' // & TRIM( simulated_time_chr ) // ' ' // & TRIM( section_chr ) ELSE rtext = TRIM( do2d(av,if) ) // ' averaged t = ' // & TRIM( simulated_time_chr ) // ' ' // & TRIM( section_chr ) ENDIF WRITE (28,LOCAL) ENDIF CASE ( 'yz' ) ! !-- Update the NetCDF xy cross section time axis IF ( myid == 0 ) THEN IF ( simulated_time /= do2d_yz_last_time(av) ) THEN do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1 do2d_yz_last_time(av) = simulated_time IF ( .NOT. data_output_2d_on_each_pe .AND. & netcdf_output ) THEN #if defined( __netcdf ) nc_stat = NF90_PUT_VAR( id_set_yz(av), & id_var_time_yz(av), & (/ time_since_reference_point /), & start = (/ do2d_yz_time_count(av) /), & count = (/ 1 /) ) CALL handle_netcdf_error( 'data_output_2d', 59 ) #endif ENDIF ENDIF ENDIF ! !-- If required, carry out averaging along x IF ( section(is,s) == -1 ) THEN ALLOCATE( local_2d_l(nys-1:nyn+1,nzb:nzt+1) ) local_2d_l = 0.0 ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 ) ! !-- First local averaging on the PE DO k = nzb, nzt+1 DO j = nys-1, nyn+1 DO i = nxl, nxr local_2d_l(j,k) = local_2d_l(j,k) + & local_pf(i,j,k) ENDDO ENDDO ENDDO #if defined( __parallel ) ! !-- Now do the averaging over all PEs along x CALL MPI_ALLREDUCE( local_2d_l(nys-1,nzb), & local_2d(nys-1,nzb), ngp, MPI_REAL, & MPI_SUM, comm1dx, ierr ) #else local_2d = local_2d_l #endif local_2d = local_2d / ( nx + 1.0 ) DEALLOCATE( local_2d_l ) ELSE ! !-- Just store the respective section on the local array !-- (but only if it is available on this PE!) IF ( section(is,s) >= nxl .AND. section(is,s) <= nxr ) & THEN local_2d = local_pf(section(is,s),:,nzb:nzt+1) ENDIF ENDIF #if defined( __parallel ) IF ( data_output_2d_on_each_pe ) THEN ! !-- Output of partial arrays on each PE. If the cross section !-- does not reside on the PE, output special index values. #if defined( __netcdf ) IF ( netcdf_output .AND. myid == 0 ) THEN WRITE ( 23 ) simulated_time, do2d_yz_time_count(av), & av ENDIF #endif IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr ) .OR.& ( section(is,s) == -1 .AND. nxl-1 == -1 ) ) & THEN WRITE (23) nys-1, nyn+1, nzb, nzt+1 WRITE (23) local_2d ELSE WRITE (23) -1, -1, -1, -1 ENDIF ELSE ! !-- PE0 receives partial arrays from all processors of the !-- respective cross section and outputs them. Here a !-- barrier has to be set, because otherwise !-- "-MPI- FATAL: Remote protocol queue full" may occur. CALL MPI_BARRIER( comm2d, ierr ) ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 ) IF ( myid == 0 ) THEN ! !-- Local array can be relocated directly. IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr ) & .OR. ( section(is,s) == -1 .AND. nxl-1 == -1 ) ) & THEN total_2d(nys-1:nyn+1,nzb:nzt+1) = local_2d ENDIF ! !-- Receive data from all other PEs. DO n = 1, numprocs-1 ! !-- Receive index limits first, then array. !-- Index limits are received in arbitrary order from !-- the PEs. CALL MPI_RECV( ind(1), 4, MPI_INTEGER, & MPI_ANY_SOURCE, 0, comm2d, status, & ierr ) ! !-- Not all PEs have data for YZ-cross-section. IF ( ind(1) /= -9999 ) THEN sender = status(MPI_SOURCE) DEALLOCATE( local_2d ) ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) ) CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, & MPI_REAL, sender, 1, comm2d, & status, ierr ) total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d ENDIF ENDDO ! !-- Output of the total cross-section. IF ( iso2d_output ) THEN WRITE (23) total_2d(0:ny+1,nzb:nzt+1) ENDIF ! !-- Relocate the local array for the next loop increment DEALLOCATE( local_2d ) ALLOCATE( local_2d(nys-1:nyn+1,nzb:nzt+1) ) #if defined( __netcdf ) IF ( netcdf_output ) THEN nc_stat = NF90_PUT_VAR( id_set_yz(av), & id_var_do2d(av,if), & total_2d(0:ny+1,nzb:nzt+1),& start = (/ is, 1, 1, do2d_yz_time_count(av) /), & count = (/ 1, ny+2, nz+2, 1 /) ) CALL handle_netcdf_error( 'data_output_2d', 60 ) ENDIF #endif ELSE ! !-- If the cross section resides on the PE, send the !-- local index limits, otherwise send -9999 to PE0. IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr ) & .OR. ( section(is,s) == -1 .AND. nxl-1 == -1 ) ) & THEN ind(1) = nys-1; ind(2) = nyn+1 ind(3) = nzb; ind(4) = nzt+1 ELSE ind(1) = -9999; ind(2) = -9999 ind(3) = -9999; ind(4) = -9999 ENDIF CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, & ierr ) ! !-- If applicable, send data to PE0. IF ( ind(1) /= -9999 ) THEN CALL MPI_SEND( local_2d(nys-1,nzb), ngp, MPI_REAL, & 0, 1, comm2d, ierr ) ENDIF ENDIF ! !-- A barrier has to be set, because otherwise some PEs may !-- proceed too fast so that PE0 may receive wrong data on !-- tag 0 CALL MPI_BARRIER( comm2d, ierr ) ENDIF #else IF ( iso2d_output ) THEN WRITE (23) local_2d(nys:nyn+1,nzb:nzt+1) ENDIF #if defined( __netcdf ) IF ( netcdf_output ) THEN nc_stat = NF90_PUT_VAR( id_set_yz(av), & id_var_do2d(av,if), & local_2d(nys:nyn+1,nzb:nzt+1), & start = (/ is, 1, 1, do2d_xz_time_count(av) /), & count = (/ 1, ny+2, nz+2, 1 /) ) CALL handle_netcdf_error( 'data_output_2d', 61 ) ENDIF #endif #endif do2d_yz_n = do2d_yz_n + 1 ! !-- Write LOCAL-parameter set for ISO2D. IF ( myid == 0 .AND. iso2d_output ) THEN IF ( section(is,s) /= -1 ) THEN WRITE ( section_chr, '(''x = '',F8.2,'' m (GP '',I3, & &'')'')' & ) section(is,s)*dx, section(is,s) ELSE section_chr = 'averaged along x' ENDIF IF ( av == 0 ) THEN rtext = TRIM( do2d(av,if) ) // ' t = ' // & TRIM( simulated_time_chr ) // ' ' // & TRIM( section_chr ) ELSE rtext = TRIM( do2d(av,if) ) // ' averaged t = ' // & TRIM( simulated_time_chr ) // ' ' // & TRIM( section_chr ) ENDIF WRITE (29,LOCAL) ENDIF END SELECT is = is + 1 ENDDO loop1 ENDIF if = if + 1 l = MAX( 2, LEN_TRIM( do2d(av,if) ) ) do2d_mode = do2d(av,if)(l-1:l) ENDDO ! !-- Deallocate temporary arrays. IF ( ALLOCATED( level_z ) ) DEALLOCATE( level_z ) DEALLOCATE( local_pf, local_2d ) #if defined( __parallel ) IF ( .NOT. data_output_2d_on_each_pe .AND. myid == 0 ) THEN DEALLOCATE( total_2d ) ENDIF #endif ! !-- Close plot output file. file_id = 20 + s IF ( data_output_2d_on_each_pe ) THEN CALL close_file( file_id ) ELSE IF ( myid == 0 ) CALL close_file( file_id ) ENDIF CALL cpu_log (log_point(3),'data_output_2d','stop','nobarrier') END SUBROUTINE data_output_2d