Changeset 4600 for palm/trunk/SOURCE


Ignore:
Timestamp:
Jul 13, 2020 6:50:12 PM (4 years ago)
Author:
suehring
Message:

adjustmens for mpi-io - surface data is transformed to a 2D-based surface array before writing; bugfix in counting of surface elements; bugfix in data-output of averaged surface data in case of restarts

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/surface_data_output_mod.f90

    r4577 r4600  
    2020! Current revisions:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! - Change: adjustmens for mpi-io - surface data is transformed to a 2D-based surface array
     28!   before writing.
     29! - Bugfix in counting of surface elements
     30! - Bugfix in data-output of averaged surface data in case of restarts
     31!
     32! 4577 2020-06-25 09:53:58Z raasch
    2733! File re-formatted to follow the PALM coding standard
    2834!
     
    178184
    179185   USE restart_data_mpi_io_mod,                                                                    &
    180        ONLY:  rd_mpi_io_check_array,                                                               &
    181               rrd_mpi_io,                                                                          &
    182               wrd_mpi_io
     186       ONLY:  rrd_mpi_io,                                                                          &
     187              rd_mpi_io_check_array,                                                               &
     188              rrd_mpi_io_surface,                                                                  &
     189              rd_mpi_io_surface_filetypes,                                                         &
     190              wrd_mpi_io,                                                                          &
     191              wrd_mpi_io_surface
    183192
    184193   USE surface_mod,                                                                                &
     
    213222      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zs        !< z-coordinate for NetCDF output
    214223      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  zenith    !< zenith orientation coordinate for NetCDF output
    215       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_out   !< output variables
    216       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_av    !< variables used for averaging
     224      REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_out   !< output variable
     225      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_av    !< variable used for averaging
    217226      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  points    !< points  / vertices of a surface element
    218227      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  polygons  !< polygon data of a surface element
     
    398407
    399408#if defined( __netcdf4_parallel )
    400     CHARACTER (LEN=100)  ::  filename            !< name of output file
    401     CHARACTER (LEN=80)   ::   time_average_text  !< string written to file attribute time_avg
    402     CHARACTER (LEN=4000) ::   var_list           !< list of variables written to NetCDF file
     409    CHARACTER (LEN=100)  ::  filename           !< name of output file
     410    CHARACTER (LEN=80)   ::  time_average_text  !< string written to file attribute time_avg
     411    CHARACTER (LEN=4000) ::  var_list           !< list of variables written to NetCDF file
    403412
    404413    INTEGER(iwp) ::  av                 !< flag for averaged (=1) and non-averaged (=0) data
     
    436445       surfaces%npoints = 0
    437446       DO  l = 0, 1
    438           DO  m = 1, surf_def_h(0)%ns
     447          DO  m = 1, surf_def_h(l)%ns
    439448!
    440449!--          Determine the indices of the respective grid cell inside the topography
    441              i = surf_def_h(0)%i(m) + surf_def_h(0)%ioff
    442              j = surf_def_h(0)%j(m) + surf_def_h(0)%joff
    443              k = surf_def_h(0)%k(m) + surf_def_h(0)%koff
     450             i = surf_def_h(l)%i(m) + surf_def_h(l)%ioff
     451             j = surf_def_h(l)%j(m) + surf_def_h(l)%joff
     452             k = surf_def_h(l)%k(m) + surf_def_h(l)%koff
    444453!
    445454!--          Check if the vertices that define the surface element are already defined, if not,
     
    45794588       CASE ( 'average_count_surf' )
    45804589          READ ( 13 )  average_count_surf
     4590       CASE ( 'time_dosurf_av' )
     4591          READ ( 13 )  time_dosurf_av
    45814592
    45824593       CASE DEFAULT
     
    45984609
    45994610    CALL rrd_mpi_io( 'average_count_surf', average_count_surf )
     4611    CALL rrd_mpi_io( 'time_dosurf_av', time_dosurf_av )
    46004612
    46014613 END SUBROUTINE surface_data_output_rrd_global_mpi
     
    46454657    IMPLICIT NONE
    46464658
    4647     LOGICAL ::  array_found  !<
    4648 
    4649 
    4650     CALL rd_mpi_io_check_array( 'surfaces%var_av' , found = array_found )
    4651 
    4652 !> does not work this way: surface%var_av has non-standard dimensions
    4653 !    IF ( array_found )  THEN
    4654 !       IF ( .NOT. ALLOCATED( surfaces%var_av ) )  ALLOCATE( ....... )
    4655 !       CALL rrd_mpi_io( 'surfaces%var_av', surfaces%var_av )
    4656 !    ENDIF
     4659    CHARACTER(LEN=3) ::  dum !< dummy string to create output-variable name
     4660
     4661    INTEGER(iwp) ::  i  !< grid index in x-direction
     4662    INTEGER(iwp) ::  j  !< grid index in y-direction
     4663    INTEGER(iwp) ::  l  !< running index surface orientation
     4664    INTEGER(iwp) ::  m  !< running index surface elements
     4665    INTEGER(iwp) ::  n  !< counting variable
     4666    INTEGER(iwp) ::  nv !< running index over number of variables
     4667
     4668    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  end_index           !< end index of surface data at (j,i)
     4669    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  global_start_index  !< index array for surface data (MPI-IO)
     4670    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  num_surf            !< number of surface data at (j,i)
     4671    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  start_index         !< start index of surface data at (j,i)
     4672
     4673    LOGICAL ::  array_found  !< flag indicating whether variable is in restart data or not
     4674    LOGICAL ::  ldum         !< dummy variable
     4675
     4676    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surf_in !< input array in expected restart format
     4677
     4678!
     4679!-- Note, surface data which is written to file is organized in a different way than
     4680!-- the output surface data. The output surface data is a concatenated array of the
     4681!-- different surface types and orientations, while the mpi-io expects surface data that
     4682!-- is consecutive in terms of start- and end-index, i.e. organized along the (j,i)
     4683!-- grid index. Hence, data need to be tranformed back to the output surface data.
     4684    ALLOCATE( end_index(nys:nyn,nxl:nxr)          )
     4685    ALLOCATE( num_surf(nys:nyn,nxl:nxr)           )
     4686    ALLOCATE( start_index(nys:nyn,nxl:nxr)        )
     4687    ALLOCATE( global_start_index(nys:nyn,nxl:nxr) )
     4688
     4689    ALLOCATE( surf_in(1:surfaces%ns) )
     4690
     4691    CALL rd_mpi_io_check_array( 'surfaces%start_index', found = array_found )
     4692    IF ( array_found )  CALL rrd_mpi_io( 'surfaces%start_index', start_index )
     4693
     4694    CALL rd_mpi_io_check_array( 'surfaces%end_index', found = array_found )
     4695    IF ( array_found )  CALL rrd_mpi_io( 'surfaces%end_index', end_index )
     4696
     4697    CALL rd_mpi_io_check_array( 'surfaces%global_start_index', found = array_found )
     4698    IF ( array_found )  CALL rrd_mpi_io( 'surfaces%global_start_index', global_start_index )
     4699
     4700    CALL rd_mpi_io_surface_filetypes( start_index, end_index, ldum, global_start_index )
     4701
     4702    DO  nv = 1, dosurf_no(1)
     4703       IF ( nv < 10                     )  WRITE( dum, '(I1)') nv
     4704       IF ( nv < 100   .AND.  nv >= 10  )  WRITE( dum, '(I2)') nv
     4705       IF ( nv < 1000  .AND.  nv >= 100 )  WRITE( dum, '(I3)') nv
     4706
     4707       CALL rd_mpi_io_check_array( 'surfaces%var_av' // TRIM( dum ), found = array_found )
     4708
     4709       IF ( array_found )  THEN
     4710
     4711          CALL rrd_mpi_io_surface( 'surfaces%var_av' // TRIM(dum), surf_in )
     4712!
     4713!--       Write temporary input variable back to surface-output data array.
     4714          n = 0
     4715          num_surf = 0
     4716          DO  l = 0, 1
     4717             DO  m = 1, surf_def_h(l)%ns
     4718                i = surf_def_h(l)%i(m)
     4719                j = surf_def_h(l)%j(m)
     4720                n                     = n + 1
     4721                surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
     4722                num_surf(j,i)         = num_surf(j,i) + 1
     4723             ENDDO
     4724          ENDDO
     4725          DO  m = 1, surf_lsm_h%ns
     4726             i = surf_lsm_h%i(m)
     4727             j = surf_lsm_h%j(m)
     4728             n                     = n + 1
     4729             surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
     4730             num_surf(j,i)         = num_surf(j,i) + 1
     4731          ENDDO
     4732          DO  m = 1, surf_usm_h%ns
     4733             i = surf_usm_h%i(m)
     4734             j = surf_usm_h%j(m)
     4735             n                     = n + 1
     4736             surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
     4737             num_surf(j,i)         = num_surf(j,i) + 1
     4738          ENDDO
     4739
     4740          DO  l = 0, 3
     4741             DO  m = 1, surf_def_v(l)%ns
     4742                i = surf_def_v(l)%i(m)
     4743                j = surf_def_v(l)%j(m)
     4744                n                     = n + 1
     4745                surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
     4746                num_surf(j,i)         = num_surf(j,i) + 1
     4747             ENDDO
     4748             DO  m = 1, surf_lsm_v(l)%ns
     4749                i = surf_lsm_v(l)%i(m)
     4750                j = surf_lsm_v(l)%j(m)
     4751                n                     = n + 1
     4752                surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
     4753                num_surf(j,i)         = num_surf(j,i) + 1
     4754             ENDDO
     4755             DO  m = 1, surf_usm_v(l)%ns
     4756                i = surf_usm_v(l)%i(m)
     4757                j = surf_usm_v(l)%j(m)
     4758                n                     = n + 1
     4759                surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
     4760                num_surf(j,i)         = num_surf(j,i) + 1
     4761             ENDDO
     4762          ENDDO
     4763       ENDIF
     4764    ENDDO
     4765
    46574766
    46584767 END SUBROUTINE surface_data_output_rrd_local_mpi
     
    46734782       WRITE ( 14 )  average_count_surf
    46744783
     4784       CALL wrd_write_string( 'time_dosurf_av' )
     4785       WRITE ( 14 )  time_dosurf_av
     4786
    46754787    ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
    46764788
    46774789      CALL wrd_mpi_io( 'average_count_surf', average_count_surf )
     4790      CALL wrd_mpi_io( 'time_dosurf_av', time_dosurf_av )
    46784791
    46794792    ENDIF
     
    46914804    IMPLICIT NONE
    46924805
     4806    CHARACTER(LEN=3) ::  dum !< dummy string to create output-variable name
     4807
     4808    INTEGER(iwp) ::  i                      !< grid index in x-direction
     4809    INTEGER(iwp) ::  j                      !< grid index in y-direction
     4810    INTEGER(iwp) ::  l                      !< running index surface orientation
     4811    INTEGER(iwp) ::  m                      !< running index surface elements
     4812    INTEGER(iwp) ::  n                      !< counting variable
     4813    INTEGER(iwp) ::  nv                     !< running index over number of variables
     4814    INTEGER(iwp) ::  start_index_aggregated !< sum of start-index at (j,i) over all surface types
     4815
     4816    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  end_index           !< end index of surface data at (j,i)
     4817    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  global_start_index  !< index array for surface data (MPI-IO)
     4818    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  num_surf            !< number of surface data at (j,i)
     4819    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  start_index         !< start index of surface data at (j,i)
     4820
     4821    LOGICAL ::  surface_data_to_write  !< switch for MPI-I/O if PE has surface data to write
     4822
     4823    REAL(wp), DIMENSION(:), ALLOCATABLE ::  surf_out  !< surface data in expected restart format
     4824
     4825
    46934826    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
    46944827
     
    46994832
    47004833    ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
    4701 
    4702        IF ( ALLOCATED( surfaces%var_av ) )  CALL wrd_mpi_io( 'surfaces%var_av', surfaces%var_av )
     4834!
     4835!--    Note, surface data which is written to file is organized in a different way than
     4836!--    the output surface data. The output surface data is a concatenated array of the
     4837!--    different surface types and orientations, while the mpi-io expects surface data that
     4838!--    is consecutive in terms of start- and end-index, i.e. organized along the (j,i)
     4839!--    grid index. Hence, data need to be tranformed before it can be written to file.
     4840       IF ( ALLOCATED( surfaces%var_av ) )  THEN
     4841          ALLOCATE( end_index(nys:nyn,nxl:nxr)          )
     4842          ALLOCATE( num_surf(nys:nyn,nxl:nxr)           )
     4843          ALLOCATE( start_index(nys:nyn,nxl:nxr)        )
     4844          ALLOCATE( global_start_index(nys:nyn,nxl:nxr) )
     4845          ALLOCATE( surf_out(1:surfaces%ns)             )
     4846!
     4847!--       Determine the start and end index at each (j,i)-pair and resort the surface data
     4848          start_index            = 1
     4849          end_index              = 0
     4850          start_index_aggregated = 1
     4851          num_surf = 0
     4852          DO  l = 0, 1
     4853             DO  m = 1, surf_def_h(l)%ns
     4854                i = surf_def_h(l)%i(m)
     4855                j = surf_def_h(l)%j(m)
     4856                num_surf(j,i) = num_surf(j,i) + 1
     4857             ENDDO
     4858          ENDDO
     4859          DO  m = 1, surf_lsm_h%ns
     4860             i = surf_lsm_h%i(m)
     4861             j = surf_lsm_h%j(m)
     4862             num_surf(j,i) = num_surf(j,i) + 1
     4863          ENDDO
     4864          DO  m = 1, surf_usm_h%ns
     4865             i = surf_usm_h%i(m)
     4866             j = surf_usm_h%j(m)
     4867             num_surf(j,i) = num_surf(j,i) + 1
     4868          ENDDO
     4869
     4870          DO  l = 0, 3
     4871             DO  m = 1, surf_def_v(l)%ns
     4872                i = surf_def_v(l)%i(m)
     4873                j = surf_def_v(l)%j(m)
     4874                num_surf(j,i) = num_surf(j,i) + 1
     4875             ENDDO
     4876             DO  m = 1, surf_lsm_v(l)%ns
     4877                i = surf_lsm_v(l)%i(m)
     4878                j = surf_lsm_v(l)%j(m)
     4879                num_surf(j,i) = num_surf(j,i) + 1
     4880             ENDDO
     4881             DO  m = 1, surf_usm_v(l)%ns
     4882                i = surf_usm_v(l)%i(m)
     4883                j = surf_usm_v(l)%j(m)
     4884                num_surf(j,i) = num_surf(j,i) + 1
     4885             ENDDO
     4886          ENDDO
     4887
     4888          start_index = 0
     4889          end_index   = 0
     4890          start_index_aggregated = 1
     4891          DO  i = nxl, nxr
     4892             DO  j = nys, nyn
     4893                start_index(j,i)       = start_index_aggregated
     4894                end_index(j,i)         = start_index(j,i) + num_surf(j,i) - 1
     4895                start_index_aggregated = start_index_aggregated + num_surf(j,i)
     4896             ENDDO
     4897          ENDDO
     4898
     4899          CALL rd_mpi_io_surface_filetypes( start_index, end_index, surface_data_to_write,         &
     4900                                            global_start_index )
     4901          CALL wrd_mpi_io( 'surfaces%start_index',        start_index        )
     4902          CALL wrd_mpi_io( 'surfaces%end_index',          end_index          )
     4903          CALL wrd_mpi_io( 'surfaces%global_start_index', global_start_index )
     4904
     4905          DO  nv = 1, dosurf_no(1)
     4906             n = 0
     4907             num_surf = 0
     4908             DO  l = 0, 1
     4909                DO  m = 1, surf_def_h(l)%ns
     4910                   i = surf_def_h(l)%i(m)
     4911                   j = surf_def_h(l)%j(m)
     4912                   n                                          = n + 1
     4913                   surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
     4914                   num_surf(j,i)                              = num_surf(j,i) + 1
     4915                ENDDO
     4916             ENDDO
     4917             DO  m = 1, surf_lsm_h%ns
     4918                i = surf_lsm_h%i(m)
     4919                j = surf_lsm_h%j(m)
     4920                n                                          = n + 1
     4921                surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
     4922                num_surf(j,i)                              = num_surf(j,i) + 1
     4923             ENDDO
     4924             DO  m = 1, surf_usm_h%ns
     4925                i = surf_usm_h%i(m)
     4926                j = surf_usm_h%j(m)
     4927                n                                          = n + 1
     4928                surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
     4929                num_surf(j,i)                              = num_surf(j,i) + 1
     4930             ENDDO
     4931
     4932             DO  l = 0, 3
     4933                DO  m = 1, surf_def_v(l)%ns
     4934                   i = surf_def_v(l)%i(m)
     4935                   j = surf_def_v(l)%j(m)
     4936                   n                                          = n + 1
     4937                   surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
     4938                   num_surf(j,i)                              = num_surf(j,i) + 1
     4939                ENDDO
     4940                DO  m = 1, surf_lsm_v(l)%ns
     4941                   i = surf_lsm_v(l)%i(m)
     4942                   j = surf_lsm_v(l)%j(m)
     4943                   n                                          = n + 1
     4944                   surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
     4945                   num_surf(j,i)                              = num_surf(j,i) + 1
     4946                ENDDO
     4947                DO  m = 1, surf_usm_v(l)%ns
     4948                   i = surf_usm_v(l)%i(m)
     4949                   j = surf_usm_v(l)%j(m)
     4950                   n                                          = n + 1
     4951                   surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
     4952                   num_surf(j,i)                              = num_surf(j,i) + 1
     4953                ENDDO
     4954             ENDDO
     4955
     4956             IF ( nv < 10                     )  WRITE( dum, '(I1)') nv
     4957             IF ( nv < 100   .AND.  nv >= 10  )  WRITE( dum, '(I2)') nv
     4958             IF ( nv < 1000  .AND.  nv >= 100 )  WRITE( dum, '(I3)') nv
     4959
     4960             CALL wrd_mpi_io_surface( 'surfaces%var_av' // TRIM( dum ), surf_out )
     4961          ENDDO
     4962
     4963       ENDIF
    47034964
    47044965    ENDIF
Note: See TracChangeset for help on using the changeset viewer.