source: palm/trunk/SOURCE/pmc_child_mod.f90 @ 4100

Last change on this file since 4100 was 3964, checked in by suehring, 6 years ago

In a nested child domain, distinguish between soil moisture and temperature initialization in case the parent domain is initialized via the dynamic input file; in the offline nesting, add a safety factor for the calculation of bulk Richardson number in order to avoid division by zero which can potentially happen if 3D buildings are located directly at the lateral model boundaries

  • Property svn:keywords set to Id
File size: 30.6 KB
RevLine 
[3962]1MODULE pmc_child
2
3!------------------------------------------------------------------------------!
4! This file is part of the PALM model system.
5!
6! PALM is free software: you can redistribute it and/or modify it under the
7! terms of the GNU General Public License as published by the Free Software
8! Foundation, either version 3 of the License, or (at your option) any later
9! version.
10!
11! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
14!
15! You should have received a copy of the GNU General Public License along with
16! PALM. If not, see <http://www.gnu.org/licenses/>.
17!
18! Copyright 1997-2019 Leibniz Universitaet Hannover
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! ------------------
[3964]23!
24!
[3962]25! Former revisions:
26! -----------------
27! $Id: pmc_child_mod.f90 3964 2019-05-09 09:48:32Z forkel $
[3964]28! Remove unused variable
29!
30! 3963 2019-05-08 20:09:11Z suehring
[3962]31! Bugfixes in initial settings of child and parent communication patterns.
32!
[3963]33! 3945 2019-05-02 11:29:27Z raasch
[3962]34!
35! 3932 2019-04-24 17:31:34Z suehring
36! typo removed
37!
38! 2019-02-25 15:31:42Z raasch
39! statement added to avoid compiler warning
40!
41! 3655 2019-01-07 16:51:22Z knoop
42! explicit kind settings
43!
44! 3241 2018-09-12 15:02:00Z raasch
45! unused variables removed
46!
47! 2841 2018-02-27 15:02:57Z knoop
48! Bugfix: wrong placement of include 'mpif.h' corrected
49!
50! 2809 2018-02-15 09:55:58Z schwenkel
51! Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE
52!
53! 2801 2018-02-14 16:01:55Z thiele
54! Introduce particle transfer in nested models.
55!
56! 2718 2018-01-02 08:49:38Z maronga
57! Corrected "Former revisions" section
58!
59! 2696 2017-12-14 17:12:51Z kanani
60! Change in file header (GPL part)
61!
62! 2599 2017-11-01 13:18:45Z hellstea
63! Some cleanup and commenting improvements only.
64!
65! 2101 2017-01-05 16:42:31Z suehring
66!
67! 2000 2016-08-20 18:09:15Z knoop
68! Forced header and separation lines into 80 columns
69!
70! 1897 2016-05-03 08:10:23Z raasch
71! Module renamed. Code clean up. The words server/client changed to parent/child.
72!
73! 1896 2016-05-03 08:06:41Z raasch
74! re-formatted to match PALM style
75!
76! 1850 2016-04-08 13:29:27Z maronga
77! Module renamed
78!
79!
80! 1833 2016-04-07 14:23:03Z raasch
81! gfortran requires pointer attributes for some array declarations,
82! long line wrapped
83!
84! 1808 2016-04-05 19:44:00Z raasch
85! MPI module used by default on all machines
86!
87! 1797 2016-03-21 16:50:28Z raasch
88! introduction of different datatransfer modes
89!
90! 1791 2016-03-11 10:41:25Z raasch
91! Debug write-statement commented out
92!
93! 1786 2016-03-08 05:49:27Z raasch
94! change in child-parent data transfer: parent now gets data from child
95! instead of that child puts it to the parent
96!
97! 1783 2016-03-06 18:36:17Z raasch
98! Bugfix: wrong data-type in MPI_WIN_CREATE replaced
99!
100! 1779 2016-03-03 08:01:28Z raasch
101! kind=dp replaced by wp, dim_order removed
102! array management changed from linked list to sequential loop
103!
104! 1764 2016-02-28 12:45:19Z raasch
105! cpp-statement added (nesting can only be used in parallel mode),
106! all kinds given in PALM style
107!
108! 1762 2016-02-25 12:31:13Z hellstea
109! Initial revision by K. Ketelsen
110!
111! Description:
112! ------------
113!
114! Child part of Palm Model Coupler
115!------------------------------------------------------------------------------!
116
117#if defined( __parallel )
118
119    USE, INTRINSIC ::  iso_c_binding
120
121#if !defined( __mpifh )
122    USE MPI
123#endif
124
125    USE kinds
126
127    USE pmc_general,                                                           &
128        ONLY:  arraydef, childdef, da_desclen, da_namedef, da_namelen, pedef,  &
129               pmc_da_name_err,  pmc_g_setname, pmc_max_array, pmc_status_ok
130
131    USE pmc_handle_communicator,                                               &
132        ONLY:  m_model_comm, m_model_npes, m_model_rank, m_to_parent_comm
133
134    USE pmc_mpi_wrapper,                                                       &
135        ONLY:  pmc_alloc_mem, pmc_bcast, pmc_inter_bcast, pmc_time
136
137    IMPLICIT NONE
138
139#if defined( __mpifh )
140    INCLUDE "mpif.h"
141#endif
142
143    PRIVATE
144    SAVE
145
146    TYPE(childdef), PUBLIC ::  me   !<
147
148    INTEGER(iwp) ::  myindex = 0         !< counter and unique number for data arrays
149    INTEGER(iwp) ::  next_array_in_list = 0   !<
150
151
152    INTERFACE pmc_childinit
153        MODULE PROCEDURE pmc_childinit
154    END INTERFACE pmc_childinit
155
156    INTERFACE pmc_c_clear_next_array_list
157        MODULE PROCEDURE pmc_c_clear_next_array_list
158    END INTERFACE pmc_c_clear_next_array_list
159
160    INTERFACE pmc_c_getbuffer
161        MODULE PROCEDURE pmc_c_getbuffer
162    END INTERFACE pmc_c_getbuffer
163
164    INTERFACE pmc_c_getnextarray
165        MODULE PROCEDURE pmc_c_getnextarray
166    END INTERFACE pmc_c_getnextarray
167
168    INTERFACE pmc_c_get_2d_index_list
169        MODULE PROCEDURE pmc_c_get_2d_index_list
170    END INTERFACE pmc_c_get_2d_index_list
171
172    INTERFACE pmc_c_putbuffer
173        MODULE PROCEDURE pmc_c_putbuffer
174    END INTERFACE pmc_c_putbuffer
175
176    INTERFACE pmc_c_setind_and_allocmem
177        MODULE PROCEDURE pmc_c_setind_and_allocmem
178    END INTERFACE pmc_c_setind_and_allocmem
179
180    INTERFACE pmc_c_set_dataarray
181        MODULE PROCEDURE pmc_c_set_dataarray_2d
182        MODULE PROCEDURE pmc_c_set_dataarray_3d
183        MODULE PROCEDURE pmc_c_set_dataarray_ip2d
184    END INTERFACE pmc_c_set_dataarray
185
186    INTERFACE pmc_set_dataarray_name
187        MODULE PROCEDURE pmc_set_dataarray_name
188        MODULE PROCEDURE pmc_set_dataarray_name_lastentry
189    END INTERFACE pmc_set_dataarray_name
190
191
192    PUBLIC pmc_childinit, pmc_c_clear_next_array_list, pmc_c_getbuffer,        &
193           pmc_c_getnextarray, pmc_c_putbuffer, pmc_c_setind_and_allocmem,     &
194           pmc_c_set_dataarray, pmc_set_dataarray_name, pmc_c_get_2d_index_list
195
196 CONTAINS
197
198
199
200 SUBROUTINE pmc_childinit
201
202     IMPLICIT NONE
203
204     INTEGER(iwp) ::  i        !<
205     INTEGER(iwp) ::  istat    !<
206
207!
208!--  Get / define the MPI environment
209     me%model_comm = m_model_comm
210     me%inter_comm = m_to_parent_comm
211
212     CALL MPI_COMM_RANK( me%model_comm, me%model_rank, istat )
213     CALL MPI_COMM_SIZE( me%model_comm, me%model_npes, istat )
214     CALL MPI_COMM_REMOTE_SIZE( me%inter_comm, me%inter_npes, istat )
215!
216!--  Intra-communicator is used for MPI_GET
217     CALL MPI_INTERCOMM_MERGE( me%inter_comm, .TRUE., me%intra_comm, istat )
218     CALL MPI_COMM_RANK( me%intra_comm, me%intra_rank, istat )
219
220     ALLOCATE( me%pes(me%inter_npes) )
221!
222!--  Allocate an array of type arraydef for all parent processes to store
223!--  information of then transfer array
224     DO  i = 1, me%inter_npes
225        ALLOCATE( me%pes(i)%array_list(pmc_max_array) )
226     ENDDO
227
228 END SUBROUTINE pmc_childinit
229
230
231
232 SUBROUTINE pmc_set_dataarray_name( parentarraydesc, parentarrayname,          &
233                                    childarraydesc, childarrayname, istat )
234
235    IMPLICIT NONE
236
237    CHARACTER(LEN=*), INTENT(IN) ::  parentarrayname  !<
238    CHARACTER(LEN=*), INTENT(IN) ::  parentarraydesc  !<
239    CHARACTER(LEN=*), INTENT(IN) ::  childarrayname   !<
240    CHARACTER(LEN=*), INTENT(IN) ::  childarraydesc   !<
241
242    INTEGER(iwp), INTENT(OUT) ::  istat  !<
243!
244!-- Local variables
245    TYPE(da_namedef) ::  myname  !<
246
247    INTEGER(iwp) ::  mype  !<
248
249
250    istat = pmc_status_ok
251!
252!-- Check length of array names
253    IF ( LEN( TRIM( parentarrayname) ) > da_namelen  .OR.                      &
254         LEN( TRIM( childarrayname) ) > da_namelen )  THEN
255       istat = pmc_da_name_err
256    ENDIF
257
258    IF ( m_model_rank == 0 )  THEN
259       myindex = myindex + 1
260       myname%couple_index = myindex
261       myname%parentdesc   = TRIM( parentarraydesc )
262       myname%nameonparent = TRIM( parentarrayname )
263       myname%childdesc    = TRIM( childarraydesc )
264       myname%nameonchild  = TRIM( childarrayname )
265    ENDIF
266
267!
268!-- Broadcast to all child processes
269!
270!-- The complete description of an transfer names array is broadcasted
271
272    CALL pmc_bcast( myname%couple_index, 0, comm=m_model_comm )
273    CALL pmc_bcast( myname%parentdesc,   0, comm=m_model_comm )
274    CALL pmc_bcast( myname%nameonparent, 0, comm=m_model_comm )
275    CALL pmc_bcast( myname%childdesc,    0, comm=m_model_comm )
276    CALL pmc_bcast( myname%nameonchild,  0, comm=m_model_comm )
277!
278!-- Broadcast to all parent processes
279!-- The complete description of an transfer array names is broadcasted als to all parent processe
280!   Only the root PE of the broadcasts to parent using intra communicator
281
282    IF ( m_model_rank == 0 )  THEN
283        mype = MPI_ROOT
284    ELSE
285        mype = MPI_PROC_NULL
286    ENDIF
287
288    CALL pmc_bcast( myname%couple_index, mype, comm=m_to_parent_comm )
289    CALL pmc_bcast( myname%parentdesc,   mype, comm=m_to_parent_comm )
290    CALL pmc_bcast( myname%nameonparent, mype, comm=m_to_parent_comm )
291    CALL pmc_bcast( myname%childdesc,    mype, comm=m_to_parent_comm )
292    CALL pmc_bcast( myname%nameonchild,  mype, comm=m_to_parent_comm )
293
294    CALL pmc_g_setname( me, myname%couple_index, myname%nameonchild )
295
296 END SUBROUTINE pmc_set_dataarray_name
297
298
299
300 SUBROUTINE pmc_set_dataarray_name_lastentry( lastentry )
301
302    IMPLICIT NONE
303
304    LOGICAL, INTENT(IN), OPTIONAL ::  lastentry  !<
305!
306!-- Local variables
307    INTEGER ::  idum  !<
308    INTEGER ::  mype  !<
309    TYPE(dA_namedef) ::  myname  !<
310
311    myname%couple_index = -1
312
313    IF ( m_model_rank == 0 )  THEN
314       mype = MPI_ROOT
315    ELSE
316       mype = MPI_PROC_NULL
317    ENDIF
318
319    CALL pmc_bcast( myname%couple_index, mype, comm=m_to_parent_comm )
320
321!
322!-- Next statement is just to avoid compiler warnings about unused variables
323    IF ( PRESENT( lastentry ) )  idum = 1
324
325 END SUBROUTINE pmc_set_dataarray_name_lastentry
326
327
328
329 SUBROUTINE pmc_c_get_2d_index_list
330
331    IMPLICIT NONE
332
333    INTEGER(iwp) :: dummy               !<
334    INTEGER(iwp) :: i, ierr, i2, j, nr  !<
335    INTEGER(iwp) :: indwin              !< MPI window object
336    INTEGER(iwp) :: indwin2             !< MPI window object
337
338    INTEGER(KIND=MPI_ADDRESS_KIND) :: win_size !< Size of MPI window 1 (in bytes)
339    INTEGER(KIND=MPI_ADDRESS_KIND) :: disp     !< Displacement unit (Integer = 4, floating poit = 8
340    INTEGER(KIND=MPI_ADDRESS_KIND) :: winsize  !< Size of MPI window 2 (in bytes)
341
342    INTEGER, DIMENSION(me%inter_npes*2) :: nrele  !< Number of Elements of a
343                                                  !< horizontal slice
344    INTEGER, DIMENSION(:), POINTER ::  myind  !<
345
346    TYPE(pedef), POINTER ::  ape  !> Pointer to pedef structure
347
348
349    win_size = STORAGE_SIZE( dummy )/8
350    CALL MPI_WIN_CREATE( dummy, win_size, iwp, MPI_INFO_NULL, me%intra_comm,   &
351                         indwin, ierr )
352!
353!-- Close window on child side and open on parent side
354    CALL MPI_WIN_FENCE( 0, indwin, ierr )
355
356!   Between the two MPI_WIN_FENCE calls, the parent can fill the RMA window
357
358!-- Close window on parent side and open on child side
359
360    CALL MPI_WIN_FENCE( 0, indwin, ierr )
361
362    DO  i = 1, me%inter_npes
363       disp = me%model_rank * 2
364       CALL MPI_GET( nrele((i-1)*2+1), 2, MPI_INTEGER, i-1, disp, 2,           &
365                     MPI_INTEGER, indwin, ierr )
366    ENDDO
367!
368!-- MPI_GET is non-blocking -> data in nrele is not available until MPI_FENCE is
369!-- called
370    CALL MPI_WIN_FENCE( 0, indwin, ierr )
371!
372!-- Allocate memory for index array
373    winsize = 0
374    DO  i = 1, me%inter_npes
375       ape => me%pes(i)
376       i2 = ( i-1 ) * 2 + 1
377       nr = nrele(i2+1)
378       IF ( nr > 0 )  THEN
379          ALLOCATE( ape%locind(nr) )
380       ELSE
381          NULLIFY( ape%locind )
382       ENDIF
383       winsize = MAX( INT( nr, MPI_ADDRESS_KIND ), winsize )
384    ENDDO
385
386    ALLOCATE( myind(2*winsize) )
387    winsize = 1
388!
389!-- Local buffer used in MPI_GET can but must not be inside the MPI Window.
390!-- Here, we use a dummy for the MPI window because the parent processes do
391!-- not access the RMA window via MPI_GET or MPI_PUT
392    CALL MPI_WIN_CREATE( dummy, winsize, iwp, MPI_INFO_NULL, me%intra_comm,    &
393                         indwin2, ierr )
394!
395!-- MPI_GET is non-blocking -> data in nrele is not available until MPI_FENCE is
396!-- called
397
398    CALL MPI_WIN_FENCE( 0, indwin2, ierr )
399
400!   Between the two MPI_WIN_FENCE calls, the parent can fill the RMA window
401
402    CALL MPI_WIN_FENCE( 0, indwin2, ierr )
403
404    DO  i = 1, me%inter_npes
405       ape => me%pes(i)
406       nr = nrele(i*2)
407       IF ( nr > 0 )  THEN
408          disp = nrele(2*(i-1)+1)
409          CALL MPI_WIN_LOCK( MPI_LOCK_SHARED , i-1, 0, indwin2, ierr )
410          CALL MPI_GET( myind, 2*nr, MPI_INTEGER, i-1, disp, 2*nr,             &
411                        MPI_INTEGER, indwin2, ierr )
412          CALL MPI_WIN_UNLOCK( i-1, indwin2, ierr )
413          DO  j = 1, nr
414             ape%locind(j)%i = myind(2*j-1)
415             ape%locind(j)%j = myind(2*j)
416          ENDDO
417          ape%nrele = nr
418       ELSE
419          ape%nrele = -1
420       ENDIF
421    ENDDO
422!
423!-- Don't know why, but this barrier is necessary before we can free the windows
424    CALL MPI_BARRIER( me%intra_comm, ierr )
425
426    CALL MPI_WIN_FREE( indWin,  ierr )
427    CALL MPI_WIN_FREE( indwin2, ierr )
428    DEALLOCATE( myind )
429
430 END SUBROUTINE pmc_c_get_2d_index_list
431
432
433
434 SUBROUTINE pmc_c_clear_next_array_list
435
436    IMPLICIT NONE
437
438    next_array_in_list = 0
439
440 END SUBROUTINE pmc_c_clear_next_array_list
441
442
443
444 LOGICAL FUNCTION pmc_c_getnextarray( myname )
445
446!
447!--  List handling is still required to get minimal interaction with
448!--  pmc_interface
449     CHARACTER(LEN=*), INTENT(OUT) ::  myname  !<
450!
451!-- Local variables
452    TYPE(pedef), POINTER    :: ape
453    TYPE(arraydef), POINTER :: ar
454
455
456    next_array_in_list = next_array_in_list + 1
457!
458!-- Array names are the same on all child PEs, so take first process to
459!-- get the name   
460    ape => me%pes(1)
461!
462!-- Check if all arrays have been processed
463    IF ( next_array_in_list > ape%nr_arrays )  THEN
464       pmc_c_getnextarray = .FALSE.
465       RETURN
466    ENDIF
467
468    ar => ape%array_list( next_array_in_list )
469
470    myname = ar%name
471!
472!-- Return true if annother array
473!-- If all array have been processed, the RETURN statement a couple of lines above is active
474
475    pmc_c_getnextarray = .TRUE.
476
477 END FUNCTION pmc_c_getnextarray
478
479
480
481 SUBROUTINE pmc_c_set_dataarray_2d( array )
482
483    IMPLICIT NONE
484
485    REAL(wp), INTENT(IN) , DIMENSION(:,:), POINTER ::  array  !<
486
487    INTEGER(iwp)               ::  i       !<
488    INTEGER(iwp)               ::  nrdims  !<
489    INTEGER(iwp), DIMENSION(4) ::  dims    !<
490
491    TYPE(C_PTR)             ::  array_adr
492    TYPE(arraydef), POINTER ::  ar
493    TYPE(pedef), POINTER    ::  ape
494
495
496    dims    = 1
497    nrdims  = 2
498    dims(1) = SIZE( array, 1 )
499    dims(2) = SIZE( array, 2 )
500
501    array_adr = C_LOC( array )
502
503    DO  i = 1, me%inter_npes
504       ape => me%pes(i)
505       ar  => ape%array_list(next_array_in_list)
506       ar%nrdims = nrdims
507       ar%dimkey = nrdims
508       ar%a_dim  = dims
509       ar%data   = array_adr
510    ENDDO
511
512 END SUBROUTINE pmc_c_set_dataarray_2d
513
514 SUBROUTINE pmc_c_set_dataarray_ip2d( array )
515
516    IMPLICIT NONE
517
518    INTEGER(idp), INTENT(IN) , DIMENSION(:,:), POINTER ::  array  !<
519
520    INTEGER(iwp)               ::  i       !<
521    INTEGER(iwp)               ::  nrdims  !<
522    INTEGER(iwp), DIMENSION(4) ::  dims    !<
523
524    TYPE(C_PTR)             ::  array_adr
525    TYPE(arraydef), POINTER ::  ar
526    TYPE(pedef), POINTER    ::  ape
527
528    dims    = 1
529    nrdims  = 2
530    dims(1) = SIZE( array, 1 )
531    dims(2) = SIZE( array, 2 )
532
533    array_adr = C_LOC( array )
534
535    DO  i = 1, me%inter_npes
536       ape => me%pes(i)
537       ar  => ape%array_list(next_array_in_list)
538       ar%nrdims = nrdims
539       ar%dimkey = 22
540       ar%a_dim  = dims
541       ar%data   = array_adr
542    ENDDO
543
544 END SUBROUTINE pmc_c_set_dataarray_ip2d
545
546 SUBROUTINE pmc_c_set_dataarray_3d (array)
547
548    IMPLICIT NONE
549
550    REAL(wp), INTENT(IN), DIMENSION(:,:,:), POINTER ::  array  !<
551
552    INTEGER(iwp)                ::  i
553    INTEGER(iwp)                ::  nrdims
554    INTEGER(iwp), DIMENSION (4) ::  dims
555   
556    TYPE(C_PTR)             ::  array_adr
557    TYPE(pedef), POINTER    ::  ape
558    TYPE(arraydef), POINTER ::  ar
559
560
561    dims    = 1
562    nrdims  = 3
563    dims(1) = SIZE( array, 1 )
564    dims(2) = SIZE( array, 2 )
565    dims(3) = SIZE( array, 3 )
566
567    array_adr = C_LOC( array )
568
569    DO  i = 1, me%inter_npes
570       ape => me%pes(i)
571       ar  => ape%array_list(next_array_in_list)
572       ar%nrdims = nrdims
573       ar%dimkey = nrdims
574       ar%a_dim  = dims
575       ar%data   = array_adr
576    ENDDO
577
578 END SUBROUTINE pmc_c_set_dataarray_3d
579
580
581
582 SUBROUTINE pmc_c_setind_and_allocmem
583
584    IMPLICIT NONE
585
586!
587!-- Naming convention for appendices:  _pc  -> parent to child transfer
588!--                                    _cp  -> child to parent transfer
589!--                                    recv -> parent to child transfer
590!--                                    send -> child to parent transfer
591    INTEGER(iwp) ::  arlen        !<
592    INTEGER(iwp) ::  myindex      !<
593    INTEGER(iwp) ::  i            !<
594    INTEGER(iwp) ::  ierr         !<
595    INTEGER(iwp) ::  istat        !<
596    INTEGER(iwp) ::  j            !<
597    INTEGER(iwp) ::  lo_nr_arrays !<
598    INTEGER(iwp) ::  rcount       !<
599    INTEGER(iwp) ::  tag          !<
600    INTEGER(iwp) ::  total_npes   !<
601
602    INTEGER(iwp), PARAMETER ::  noindex = -1  !<
603
604    INTEGER(idp)                   ::  bufsize  !< size of MPI data window
605    INTEGER(KIND=MPI_ADDRESS_KIND) ::  winsize  !<
606   
607    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  myindex_s
608    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  myindex_r
609
610    REAL(wp), DIMENSION(:), POINTER, SAVE ::  base_array_pc  !< base array
611    REAL(wp), DIMENSION(:), POINTER, SAVE ::  base_array_cp  !< base array
612
613    TYPE(pedef), POINTER    ::  ape       !<
614    TYPE(arraydef), POINTER ::  ar        !<
615    Type(C_PTR)             ::  base_ptr  !<
616
617 
618    CALL MPI_COMM_SIZE (me%intra_comm, total_npes, ierr)
619
620    lo_nr_arrays = me%pes(1)%nr_arrays
621
622    ALLOCATE(myindex_s(lo_nr_arrays,0:total_npes-1))
623    ALLOCATE(myindex_r(lo_nr_arrays,0:total_npes-1))
624
625    myindex_s = 0
626
627!
628!-- Receive indices from child
629    CALL MPI_ALLTOALL( myindex_s, lo_nr_arrays, MPI_INTEGER,                   &
630                       myindex_r, lo_nr_arrays, MPI_INTEGER,                   &
631                       me%intra_comm, ierr )
632
633    myindex = 0
634    bufsize = 8
635!
636!-- Parent to child direction.
637!-- First stride: compute size and set index
638    DO  i = 1, me%inter_npes
639       ape => me%pes(i)
640       DO  j = 1, ape%nr_arrays
641          ar => ape%array_list(j)
642          ar%recvindex = myindex_r(j,i-1)
643!
644!--       Determine max, because child buffer is allocated only once
645!--       All 2D and 3d arrays use the same buffer
646
647          IF ( ar%nrdims == 3 )  THEN
648             bufsize = MAX( bufsize,                                           &
649                            INT( ar%a_dim(1)*ar%a_dim(2)*ar%a_dim(3),          &
650                                 MPI_ADDRESS_KIND ) )
651          ELSE
652             bufsize = MAX( bufsize,                                           &
653                            INT( ar%a_dim(1)*ar%a_dim(2), MPI_ADDRESS_KIND ) )
654          ENDIF
655       ENDDO
656    ENDDO
657
658!
659!-- Create RMA (one sided communication) data buffer.
660!-- The buffer for MPI_GET can be PE local, i.e. it can but must not be part of
661!-- the MPI RMA window
662    CALL pmc_alloc_mem( base_array_pc, bufsize, base_ptr )
663    me%totalbuffersize = bufsize*wp  ! total buffer size in byte
664!
665!-- Second stride: set buffer pointer
666    DO  i = 1, me%inter_npes
667       ape => me%pes(i)
668       DO  j = 1, ape%nr_arrays
669          ar => ape%array_list(j)
670          ar%recvbuf = base_ptr
671       ENDDO
672    ENDDO
673!
674!-- Child to parent direction
675    myindex = 1
676    rcount  = 0
677    bufsize = 8
678
679    myindex_s = 0
680    myindex_r = 0
681
682    DO  i = 1, me%inter_npes
683       ape => me%pes(i)
684       tag = 300
685       DO  j = 1, ape%nr_arrays
686          ar => ape%array_list(j)
687          IF ( ar%nrdims == 2 )  THEN
688             arlen = ape%nrele
689          ELSEIF( ar%nrdims == 3 )  THEN
690             arlen = ape%nrele*ar%a_dim(1)
691          ENDIF
692
693          IF ( ape%nrele > 0 )  THEN
694             ar%sendindex = myindex
695          ELSE
696             ar%sendindex = noindex
697          ENDIF
698
699          myindex_s(j,i-1) = ar%sendindex
700
701          IF ( ape%nrele > 0 )  THEN
702             ar%sendsize = arlen
703             myindex     = myindex + arlen
704             bufsize     = bufsize + arlen
705          ENDIF
706
707       ENDDO
708
709    ENDDO
710!
711!-- Send indices to parent
712
713    CALL MPI_ALLTOALL( myindex_s, lo_nr_arrays, MPI_INTEGER,                   &
714                       myindex_r, lo_nr_arrays, MPI_INTEGER,                   &
715                       me%intra_comm, ierr)
716
717    DEALLOCATE( myindex_s )
718    DEALLOCATE( myindex_r )
719
720!
721!-- Create RMA (one sided communication) window for data buffer child to parent
722!-- transfer.
723!-- The buffer of MPI_GET (counter part of transfer) can be PE-local, i.e. it
724!-- can but must not be part of the MPI RMA window. Only one RMA window is
725!-- required to prepare the data
726!--        for parent -> child transfer on the parent side
727!-- and
728!--        for child -> parent transfer on the child side
729    CALL pmc_alloc_mem( base_array_cp, bufsize )
730    me%totalbuffersize = bufsize * wp  ! total buffer size in byte
731
732    winSize = me%totalbuffersize
733
734    CALL MPI_WIN_CREATE( base_array_cp, winsize, wp, MPI_INFO_NULL,            &
735                         me%intra_comm, me%win_parent_child, ierr )
736    CALL MPI_WIN_FENCE( 0, me%win_parent_child, ierr )
737    CALL MPI_BARRIER( me%intra_comm, ierr )
738!
739!-- Second stride: set buffer pointer
740    DO  i = 1, me%inter_npes
741       ape => me%pes(i)
742       DO  j = 1, ape%nr_arrays
743          ar => ape%array_list(j)         
744          IF ( ape%nrele > 0 )  THEN
745             ar%sendbuf = C_LOC( base_array_cp(ar%sendindex) )
746!
747!--          TODO: if this is an error to be really expected, replace the
748!--                following message by a meaningful standard PALM message using
749!--                the message-routine
750             IF ( ar%sendindex+ar%sendsize > bufsize )  THEN
751                WRITE( 0,'(a,i4,4i7,1x,a)') 'Child buffer too small ', i,      &
752                          ar%sendindex, ar%sendsize, ar%sendindex+ar%sendsize, &
753                          bufsize, TRIM( ar%name )
754                CALL MPI_ABORT( MPI_COMM_WORLD, istat, ierr )
755             ENDIF
756          ENDIF
757       ENDDO
758    ENDDO
759
760 END SUBROUTINE pmc_c_setind_and_allocmem
761
762
763
764 SUBROUTINE pmc_c_getbuffer( waittime, particle_transfer )
765
766    IMPLICIT NONE
767
768    REAL(wp), INTENT(OUT), OPTIONAL ::  waittime  !<
769    LOGICAL, INTENT(IN), OPTIONAL   ::  particle_transfer  !<
770
771    LOGICAL                        ::  lo_ptrans!<
772   
773    INTEGER(iwp)                        ::  ierr    !<
774    INTEGER(iwp)                        ::  ij      !<
775    INTEGER(iwp)                        ::  ip      !<
776    INTEGER(iwp)                        ::  j       !<
777    INTEGER(iwp)                        ::  myindex !<
778    INTEGER(iwp)                        ::  nr      !< number of elements to get
779                                                    !< from parent
780    INTEGER(KIND=MPI_ADDRESS_KIND) ::  target_disp
781    INTEGER,DIMENSION(1)           ::  buf_shape
782
783    REAL(wp)                            ::  t1
784    REAL(wp)                            ::  t2
785
786    REAL(wp), POINTER, DIMENSION(:)     ::  buf
787    REAL(wp), POINTER, DIMENSION(:,:)   ::  data_2d
788    REAL(wp), POINTER, DIMENSION(:,:,:) ::  data_3d
789    TYPE(pedef), POINTER                ::  ape
790    TYPE(arraydef), POINTER             ::  ar
791    INTEGER(idp), POINTER, DIMENSION(:)     ::  ibuf      !<
792    INTEGER(idp), POINTER, DIMENSION(:,:)   ::  idata_2d  !<
793
794!
795!-- Synchronization of the model is done in pmci_synchronize.
796!-- Therefore the RMA window can be filled without
797!-- sychronization at this point and a barrier is not necessary.
798
799!-- In case waittime is present, the following barrier is necessary to
800!-- insure the same number of barrier calls on parent and child
801!-- This means, that here on child side two barriers are call successively
802!-- The parent is filling its buffer between the two barrier calls
803
804!-- Please note that waittime has to be set in pmc_s_fillbuffer AND
805!-- pmc_c_getbuffer
806    IF ( PRESENT( waittime ) )  THEN
807       t1 = pmc_time()
808       CALL MPI_BARRIER( me%intra_comm, ierr )
809       t2 = pmc_time()
810       waittime = t2 - t1
811    ENDIF
812
813    lo_ptrans = .FALSE.
814    IF ( PRESENT( particle_transfer))    lo_ptrans = particle_transfer
815
816!
817!-- Wait for buffer is filled.
818!
819!-- The parent side (in pmc_s_fillbuffer) is filling the buffer in the MPI RMA window
820!-- When the filling is complet, a MPI_BARRIER is called.
821!-- The child is not allowd to access the parent-buffer before it is completely filled
822!-- therefore the following barrier is required.
823
824    CALL MPI_BARRIER( me%intra_comm, ierr )
825
826    DO  ip = 1, me%inter_npes
827       ape => me%pes(ip)
828       DO  j = 1, ape%nr_arrays
829          ar => ape%array_list(j)
830
831          IF ( ar%dimkey == 2 .AND. .NOT.lo_ptrans)  THEN
832             nr = ape%nrele
833          ELSEIF ( ar%dimkey == 3 .AND. .NOT.lo_ptrans)  THEN
834             nr = ape%nrele * ar%a_dim(1)
835          ELSE IF ( ar%dimkey == 22 .AND. lo_ptrans)  THEN
836             nr = ape%nrele
837          ELSE
838             CYCLE                    ! Particle array ar not transferd here
839          ENDIF
840          buf_shape(1) = nr
841          IF ( lo_ptrans )   THEN
842             CALL C_F_POINTER( ar%recvbuf, ibuf, buf_shape )
843          ELSE
844             CALL C_F_POINTER( ar%recvbuf, buf, buf_shape )
845          ENDIF
846!
847!--       MPI passive target RMA
848!--       One data array is fetcht from MPI RMA window on parent
849
850          IF ( nr > 0 )  THEN
851             target_disp = ar%recvindex - 1
852             CALL MPI_WIN_LOCK( MPI_LOCK_SHARED , ip-1, 0,                     &
853                                me%win_parent_child, ierr )
854             IF ( lo_ptrans )   THEN
855                CALL MPI_GET( ibuf, nr*8, MPI_BYTE, ip-1, target_disp, nr*8, MPI_BYTE,  &               !There is no MPI_INTEGER8 datatype
856                                   me%win_parent_child, ierr )
857             ELSE
858                CALL MPI_GET( buf, nr, MPI_REAL, ip-1, target_disp, nr,        &
859                              MPI_REAL, me%win_parent_child, ierr )
860             ENDIF
861             CALL MPI_WIN_UNLOCK( ip-1, me%win_parent_child, ierr )
862          ENDIF
863          myindex = 1
864          IF ( ar%dimkey == 2 .AND. .NOT.lo_ptrans)  THEN
865
866             CALL C_F_POINTER( ar%data, data_2d, ar%a_dim(1:2) )
867             DO  ij = 1, ape%nrele
868                data_2d(ape%locind(ij)%j,ape%locind(ij)%i) = buf(myindex)
869                myindex = myindex + 1
870             ENDDO
871
872          ELSEIF ( ar%dimkey == 3 .AND. .NOT.lo_ptrans)  THEN
873
874             CALL C_F_POINTER( ar%data, data_3d, ar%a_dim(1:3) )
875             DO  ij = 1, ape%nrele
876                data_3d(:,ape%locind(ij)%j,ape%locind(ij)%i) =                 &
877                                              buf(myindex:myindex+ar%a_dim(1)-1)
878                myindex = myindex+ar%a_dim(1)
879             ENDDO
880
881          ELSEIF ( ar%dimkey == 22 .AND. lo_ptrans)  THEN
882             CALL C_F_POINTER( ar%data, idata_2d, ar%a_dim(1:2) )
883
884             DO  ij = 1, ape%nrele
885                idata_2d(ape%locind(ij)%j,ape%locind(ij)%i) = ibuf(myindex)
886                myindex = myindex + 1
887             ENDDO
888
889          ENDIF
890       ENDDO
891    ENDDO
892
893 END SUBROUTINE pmc_c_getbuffer
894
895
896
897 SUBROUTINE pmc_c_putbuffer( waittime , particle_transfer )
898
899    IMPLICIT NONE
900
901    REAL(wp), INTENT(OUT), OPTIONAL ::  waittime  !<
902    LOGICAL, INTENT(IN), OPTIONAL   ::  particle_transfer  !<
903
904    LOGICAL ::  lo_ptrans!<
905   
906    INTEGER(iwp) ::  ierr         !<
907    INTEGER(iwp) ::  ij           !<
908    INTEGER(iwp) ::  ip           !<
909    INTEGER(iwp) ::  j            !<
910    INTEGER(iwp) ::  myindex      !<
911
912    INTEGER(iwp), DIMENSION(1) ::  buf_shape    !<
913
914    REAL(wp) ::  t1  !<
915    REAL(wp) ::  t2  !<
916
917    REAL(wp), POINTER, DIMENSION(:)         ::  buf      !<
918    REAL(wp), POINTER, DIMENSION(:,:)       ::  data_2d  !<
919    REAL(wp), POINTER, DIMENSION(:,:,:)     ::  data_3d  !<
920   
921    INTEGER(idp), POINTER, DIMENSION(:)     ::  ibuf      !<
922    INTEGER(idp), POINTER, DIMENSION(:,:)   ::  idata_2d  !<
923
924    TYPE(pedef), POINTER                    ::  ape  !<
925    TYPE(arraydef), POINTER                 ::  ar   !<
926
927!
928!-- Wait for empty buffer
929!-- Switch RMA epoche
930
931    t1 = pmc_time()
932    CALL MPI_BARRIER( me%intra_comm, ierr )
933    t2 = pmc_time()
934    IF ( PRESENT( waittime ) )  waittime = t2 - t1
935
936    lo_ptrans = .FALSE.
937    IF ( PRESENT( particle_transfer))    lo_ptrans = particle_transfer
938
939    DO  ip = 1, me%inter_npes
940       ape => me%pes(ip)
941       DO  j = 1, ape%nr_arrays
942          ar => aPE%array_list(j)
943          myindex = 1
944
945          IF ( ar%dimkey == 2 .AND. .NOT.lo_ptrans )  THEN
946
947             buf_shape(1) = ape%nrele
948             CALL C_F_POINTER( ar%sendbuf, buf,     buf_shape     )
949             CALL C_F_POINTER( ar%data,    data_2d, ar%a_dim(1:2) )
950             DO  ij = 1, ape%nrele
951                buf(myindex) = data_2d(ape%locind(ij)%j,ape%locind(ij)%i)
952                myindex = myindex + 1
953             ENDDO
954
955          ELSEIF ( ar%dimkey == 3 .AND. .NOT.lo_ptrans )  THEN
956
957             buf_shape(1) = ape%nrele*ar%a_dim(1)
958             CALL C_F_POINTER( ar%sendbuf, buf,     buf_shape     )
959             CALL C_F_POINTER( ar%data,    data_3d, ar%a_dim(1:3) )
960             DO  ij = 1, ape%nrele
961                buf(myindex:myindex+ar%a_dim(1)-1) =                            &
962                                    data_3d(:,ape%locind(ij)%j,ape%locind(ij)%i)
963                myindex = myindex + ar%a_dim(1)
964             ENDDO
965
966          ELSE IF ( ar%dimkey == 22 .AND. lo_ptrans)  THEN
967
968             buf_shape(1) = ape%nrele
969             CALL C_F_POINTER( ar%sendbuf, ibuf,     buf_shape     )
970             CALL C_F_POINTER( ar%data,    idata_2d, ar%a_dim(1:2) )
971
972             DO  ij = 1, ape%nrele
973                ibuf(myindex) = idata_2d(ape%locind(ij)%j,ape%locind(ij)%i)
974                myindex = myindex + 1
975             ENDDO
976
977          ENDIF
978       ENDDO
979    ENDDO
980!
981!-- Buffer is filled
982!-- Switch RMA epoche
983
984    CALL MPI_Barrier(me%intra_comm, ierr)
985
986 END SUBROUTINE pmc_c_putbuffer
987
988#endif
989 END MODULE pmc_child
Note: See TracBrowser for help on using the repository browser.