00001 module seq_rearr_mod
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 use shr_sys_mod
00014 use mct_mod
00015 use seq_flds_mod
00016 use seq_cdata_mod
00017 use seq_comm_mct
00018 use seq_diag_mct
00019 use shr_kind_mod, only: R8 => SHR_KIND_R8
00020
00021 implicit none
00022 save
00023 private
00024
00025
00026
00027
00028
00029 public :: seq_rearr_init
00030 public :: seq_rearr_gsmapIdentical
00031 private :: seq_rearr_gsmapExtend
00032 private :: seq_rearr_gsmapCreate
00033 private :: seq_rearr_avExtend
00034 private :: seq_rearr_avCreate
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #ifdef CPP_VECTOR
00045 logical :: usevector = .true.
00046 #else
00047 logical :: usevector = .false.
00048 #endif
00049
00050 #ifdef SYSUNICOS
00051 logical :: usealltoall = .true.
00052 #else
00053 logical :: usealltoall = .false.
00054 #endif
00055
00056 character(*),parameter :: subName = '(seq_rearr_mct)'
00057 real(r8),parameter :: c1 = 1.0_r8
00058
00059 #include <mpif.h>
00060
00061
00062 contains
00063
00064
00065 subroutine seq_rearr_init( cdata_old, AV1_old, AV2_old, ID_old, &
00066 cdata_new, AV1_new, AV2_new, ID_new, ID_join, &
00067 Re_old2new, Re_new2old)
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089 implicit none
00090
00091
00092
00093
00094
00095 type(seq_cdata),intent(in) :: cdata_old
00096 type(mct_aVect),intent(inout),optional :: AV1_old
00097 type(mct_aVect),intent(inout),optional :: AV2_old
00098 integer ,intent(in) :: ID_old
00099 type(seq_cdata),intent(inout) :: cdata_new
00100 type(mct_aVect),intent(inout),optional :: AV1_new
00101 type(mct_aVect),intent(inout),optional :: AV2_new
00102 integer ,intent(in) :: ID_new
00103 integer ,intent(in) :: ID_join
00104 type(mct_rearr),intent(inout) :: Re_old2new
00105 type(mct_rearr),intent(inout) :: Re_new2old
00106
00107
00108
00109 character(len=*),parameter :: subname = "(seq_rearr_init) "
00110 type(mct_gsMap), pointer :: gsmap_new
00111 type(mct_gsMap), pointer :: gsmap_old
00112 type(mct_gGrid), pointer :: dom_new
00113 type(mct_gGrid), pointer :: dom_old
00114
00115 integer :: lsize
00116 integer :: mpicom_old
00117 integer :: mpicom_new
00118 integer :: mpicom_join
00119 logical :: master_join
00120 integer :: mpigrp_old, mpigrp_join
00121 integer :: ierr
00122 integer :: ID
00123
00124 type(mct_gsMap) :: gsmap_new_join
00125 type(mct_gsMap) :: gsmap_old_join
00126 type(mct_gGrid) :: dom_join
00127 type(mct_gGrid) :: dom_empty
00128
00129
00130
00131
00132 call seq_cdata_setptrs(cdata_old, gsMap=gsmap_old, dom=dom_old, ID=ID)
00133 call seq_cdata_setptrs(cdata_new, gsMap=gsmap_new, dom=dom_new)
00134
00135 if (ID /= ID_old) then
00136 write(logunit,*) subname,' ERROR: inconsistent ID',ID,ID_old
00137 call shr_sys_abort()
00138 endif
00139
00140 call seq_comm_setptrs(ID_old ,mpicom=mpicom_old ,mpigrp=mpigrp_old )
00141 call seq_comm_setptrs(ID_new ,mpicom=mpicom_new )
00142 call seq_comm_setptrs(ID_join,mpicom=mpicom_join,mpigrp=mpigrp_join,iamroot=master_join)
00143
00144
00145
00146
00147
00148
00149 call seq_rearr_gsmapExtend(gsmap_old , mpicom_old , gsmap_old_join, mpicom_join, ID_join)
00150 call seq_rearr_gsmapCreate(gsmap_old_join, mpicom_join , gsmap_new , mpicom_new , ID_new )
00151 call seq_rearr_gsmapExtend(gsmap_new , mpicom_new , gsmap_new_join, mpicom_join, ID_join)
00152
00153
00154
00155 call mct_rearr_init(gsmap_old_join, gsmap_new_join, mpicom_join, Re_old2new)
00156 call mct_rearr_init(gsmap_new_join, gsmap_old_join, mpicom_join, Re_new2old)
00157
00158
00159
00160 lsize = mct_gsMap_lsize(gsmap_new_join, mpicom_join)
00161 call mct_gGrid_init( GGrid=dom_new, CoordChars=seq_flds_dom_coord, OtherChars=seq_flds_dom_other, lsize=lsize )
00162 call mct_avect_zero(dom_new%data)
00163 call seq_rearr_avExtend(dom_old%data, ID_old, ID_join)
00164 call mct_rearr_rearrange(dom_old%data, dom_new%data, Re_old2new, VECTOR=usevector, ALLTOALL=usealltoall)
00165
00166
00167
00168 call mct_gsMap_clean(gsmap_old_join)
00169 call mct_gsMap_clean(gsmap_new_join)
00170
00171
00172
00173 lsize = 0
00174 if (seq_comm_iamin(ID_new)) then
00175 lsize = mct_gsMap_lsize(gsMap_new,mpicom_new)
00176 endif
00177 if (present(AV1_old)) then
00178 call seq_rearr_avExtend(AV1_old, ID_old, ID_join)
00179 if (present(AV1_new)) then
00180 call seq_rearr_avCreate(AV1_old, ID_old, AV1_new, ID_join, lsize)
00181 endif
00182 endif
00183 if (present(AV2_old)) then
00184 call seq_rearr_avExtend(AV2_old, ID_old, ID_join)
00185 if (present(AV2_new)) then
00186 call seq_rearr_avCreate(AV2_old, ID_old, AV2_new, ID_join, lsize)
00187 endif
00188 endif
00189
00190 end subroutine seq_rearr_init
00191
00192
00193
00194 subroutine seq_rearr_gsmapExtend(gsmapi, mpicomi, gsmapo, mpicomo, compido)
00195
00196
00197
00198
00199
00200
00201 implicit none
00202 type(mct_gsMap), intent(IN) :: gsmapi
00203 integer , intent(IN) :: mpicomi
00204 type(mct_gsMap), intent(OUT):: gsmapo
00205 integer , intent(IN) :: mpicomo
00206 integer , intent(IN) :: compido
00207
00208 character(len=*),parameter :: subname = "(seq_rearr_gsmapExtend) "
00209 integer :: n
00210 integer :: ngseg
00211 integer :: gsize
00212 integer :: msizei,msizeo
00213 integer :: mrank,mranko,mrankog
00214 integer :: mpigrpi,mpigrpo
00215 integer :: ierr
00216 integer, pointer :: pei(:),peo(:)
00217 integer, pointer :: start(:),length(:),peloc(:)
00218
00219 mranko = -1
00220
00221
00222
00223 if (mpicomi /= MPI_COMM_NULL) then
00224 call mpi_comm_rank(mpicomi,mrank,ierr)
00225 call shr_mpi_chkerr(ierr,subname//' gsm_cop mpi_comm_rank i')
00226 if (mrank == 0) then
00227 call mpi_comm_group(mpicomi,mpigrpi,ierr)
00228 call shr_mpi_chkerr(ierr,subname//' gsm_cop mpi_comm_group i')
00229 call mpi_comm_group(mpicomo,mpigrpo,ierr)
00230 call shr_mpi_chkerr(ierr,subname//' gsm_cop mpi_comm_group o')
00231 call mpi_comm_size(mpicomi,msizei,ierr)
00232 call shr_mpi_chkerr(ierr,subname//' gsm_cop mpi_comm_size i')
00233 call mpi_comm_size(mpicomo,msizeo,ierr)
00234 call shr_mpi_chkerr(ierr,subname//' gsm_cop mpi_comm_size o')
00235
00236
00237
00238
00239 allocate(pei(0:msizei-1),peo(0:msizei-1))
00240 do n = 0,msizei-1
00241 pei(n) = n
00242 enddo
00243
00244 peo = -1
00245 call mpi_group_translate_ranks(mpigrpi,msizei,pei,mpigrpo,peo,ierr)
00246 call shr_mpi_chkerr(ierr,subname//' gsm_cop mpi_group_translate_ranks')
00247
00248 do n = 0,msizei-1
00249 if (peo(n) < 0 .or. peo(n) > msizeo-1) then
00250 write(logunit,*) subname,' peo out of bounds ',peo(n),msizeo
00251 call shr_sys_abort()
00252 endif
00253 enddo
00254
00255 mranko = peo(0)
00256
00257
00258
00259
00260 ngseg = gsmapi%ngseg
00261 gsize = gsmapi%gsize
00262 allocate(start(ngseg),length(ngseg),peloc(ngseg))
00263 do n = 1,ngseg
00264 start(n) = gsmapi%start(n)
00265 length(n) = gsmapi%length(n)
00266 peloc(n) = peo(gsmapi%pe_loc(n))
00267 enddo
00268
00269
00270
00271 call mct_gsmap_init(gsmapo,compido,ngseg,gsize,start,length,peloc)
00272
00273 deallocate(pei,peo,start,length,peloc)
00274 endif
00275 endif
00276
00277
00278
00279
00280 call mpi_allreduce(mranko,mrankog,1,MPI_INTEGER,MPI_MAX,mpicomo,ierr)
00281 call shr_mpi_chkerr(ierr,subname//' gsm_cop mpi_allreduce max')
00282
00283
00284
00285 call mct_gsmap_bcast(gsmapo, mrankog, mpicomo)
00286
00287
00288 #if (1 == 0)
00289 write(logunit,*) trim(subname),'tcxa ',mpicomi,mpicomo
00290 call shr_sys_flush(logunit)
00291 call mpi_barrier(mpicomo,ierr)
00292
00293 if (mpicomi /= MPI_COMM_NULL) then
00294 call mpi_comm_rank(mpicomi,mrank,ierr)
00295 write(logunit,*) 'tcxbi ',mrank
00296 if (mrank == 0) then
00297 write(logunit,*) 'tcxci ',gsmapi%ngseg,size(gsmapi%start),gsmapi%gsize,gsmapi%comp_id
00298 do n = 1,gsmapi%ngseg
00299 write(logunit,*) 'tcx gsmti ',n,gsmapi%start(n),gsmapi%length(n),gsmapi%pe_loc(n)
00300 enddo
00301 call shr_sys_flush(logunit)
00302 endif
00303 endif
00304
00305 if (mpicomo /= MPI_COMM_NULL) then
00306 call mpi_comm_rank(mpicomo,mrank,ierr)
00307 write(logunit,*) 'tcxbo ',mrank
00308 if (mrank == 0) then
00309 write(logunit,*) 'tcxco ',gsmapo%ngseg,size(gsmapo%start),gsmapo%gsize,gsmapo%comp_id
00310 do n = 1,gsmapo%ngseg
00311 write(logunit,*) 'tcx gsmto ',n,gsmapo%start(n),gsmapo%length(n),gsmapo%pe_loc(n)
00312 enddo
00313 call shr_sys_flush(logunit)
00314 endif
00315 endif
00316
00317 call shr_sys_flush(logunit)
00318 call mpi_barrier(mpicomo,ierr)
00319 #endif
00320
00321
00322 end subroutine seq_rearr_gsmapExtend
00323
00324
00325 subroutine seq_rearr_gsmapCreate(gsmapi, mpicomi, gsmapo, mpicomo, compido)
00326
00327
00328
00329
00330
00331 implicit none
00332 type(mct_gsMap), intent(IN) :: gsmapi
00333 integer , intent(IN) :: mpicomi
00334 type(mct_gsMap), intent(OUT):: gsmapo
00335 integer , intent(IN) :: mpicomo
00336 integer , intent(IN) :: compido
00337
00338 character(len=*),parameter :: subname = "(seq_rearr_gsmapCreate) "
00339 integer :: n,m,k
00340 integer :: ktot
00341 integer :: apesi, apeso
00342 integer :: lsizeo
00343 integer :: ngsegi,ngsego
00344 integer :: gsizei,gsizeo
00345 integer :: msizei,msizeo
00346 integer :: mranki,mranko
00347 integer :: ierr
00348 integer :: decomp_type
00349 integer, pointer :: start(:),length(:),peloc(:),perm(:),gindex(:),lindex(:)
00350 real(r8):: rpeloc
00351
00352
00353
00354
00355
00356 if (mpicomo /= MPI_COMM_NULL) then
00357 call mpi_comm_rank(mpicomi,mranki,ierr)
00358 call shr_mpi_chkerr(ierr,subname//' mpi_comm_rank i')
00359 call mpi_comm_size(mpicomi,msizei,ierr)
00360 call shr_mpi_chkerr(ierr,subname//' mpi_comm_size i')
00361 call mpi_comm_rank(mpicomo,mranko,ierr)
00362 call shr_mpi_chkerr(ierr,subname//' mpi_comm_rank o')
00363 call mpi_comm_size(mpicomo,msizeo,ierr)
00364 call shr_mpi_chkerr(ierr,subname//' mpi_comm_size o')
00365
00366 ngsegi = gsmapi%ngseg
00367 gsizei = gsmapi%gsize
00368 gsizeo = gsizei
00369 call mct_gsMap_activepes(gsmapi,apesi)
00370
00371 decomp_type = 0
00372
00373 if (msizeo == apesi) then
00374
00375 decomp_type = 2
00376 elseif (ngsegi >= msizeo) then
00377 decomp_type = 2
00378 else
00379 decomp_type = 3
00380 endif
00381
00382
00383
00384
00385 select case (decomp_type)
00386
00387 case(1)
00388
00389
00390 call mct_gsMap_copy(gsmapi,gsmapo)
00391 ngsego = ngsegi
00392 do n = 1,ngsego
00393 gsmapo%pe_loc(n) = mod(gsmapo%pe_loc(n),msizeo)
00394 enddo
00395
00396 case(2)
00397
00398
00399 ngsego = ngsegi
00400 allocate(start(ngsego),length(ngsego),peloc(ngsego),perm(ngsego))
00401 do n = 1,ngsego
00402 start(n) = gsmapi%start(n)
00403 length(n) = gsmapi%length(n)
00404 enddo
00405
00406 call mct_indexset(perm)
00407 call mct_indexsort(ngsego,perm,start)
00408 call mct_permute(start,perm,ngsego)
00409 call mct_permute(length,perm,ngsego)
00410
00411 do n = 1,ngsego
00412 rpeloc = (((msizeo*c1)*((n-1)*c1))/(ngsego*c1))
00413 peloc(n) = int(rpeloc)
00414 enddo
00415 call mct_gsmap_init(gsmapo,ngsego,start,length,peloc,0,mpicomo,compido,gsizeo)
00416 deallocate(start,length,peloc,perm)
00417
00418 case(3)
00419
00420
00421
00422 k = 0
00423 do n = 1,ngsegi
00424 do m = 1,gsmapi%length(n)
00425 k = k + 1
00426 if (k > gsizei) then
00427 write(logunit,*) trim(subname),' ERROR in gindex ',k,gsizei
00428 call shr_sys_abort()
00429 endif
00430 enddo
00431 enddo
00432 ktot = k
00433
00434 allocate(gindex(ktot),perm(ktot))
00435
00436 k = 0
00437 do n = 1,ngsegi
00438 do m = 1,gsmapi%length(n)
00439 k = k + 1
00440 gindex(k) = gsmapi%start(n) + m - 1
00441 enddo
00442 enddo
00443 call mct_indexset(perm)
00444 call mct_indexsort(ktot,perm,gindex)
00445 call mct_permute(gindex,perm,ktot)
00446
00447 k = 0
00448 do m = 0,msizeo-1
00449 lsizeo = ktot/msizeo
00450 if (m < (ktot - lsizeo*msizeo)) lsizeo = lsizeo + 1
00451 if (mranko == m) then
00452 allocate(lindex(lsizeo))
00453 if (k+lsizeo > ktot) then
00454 write(logunit,*) trim(subname),' ERROR: decomp out of bounds ',mranko,k,lsizeo,ktot
00455 call shr_sys_abort()
00456 endif
00457 lindex(1:lsizeo) = gindex(k+1:k+lsizeo)
00458
00459 endif
00460 k = k + lsizeo
00461 enddo
00462 if (k /= ktot) then
00463 write(logunit,*) trim(subname),' ERROR: decomp incomplete ',k,ktot
00464 call shr_sys_abort()
00465 endif
00466
00467 call mct_gsmap_init(gsmapo,lindex,mpicomo,compido,size(lindex),gsizeo)
00468 deallocate(gindex,perm,lindex)
00469
00470 end select
00471
00472 if (mranko == 0) then
00473 write(logunit,102) trim(subname),' created new gsmap decomp_type =',decomp_type
00474 write(logunit,102) trim(subname),' ngseg/gsize = ', &
00475 mct_gsmap_ngseg(gsmapo),mct_gsmap_gsize(gsmapo)
00476 call mct_gsmap_activepes(gsmapo,apeso)
00477 write(logunit,102) trim(subname),' mpisize/active_pes = ', &
00478 msizeo,apeso
00479 write(logunit,102) trim(subname),' avg seg per pe/ape = ', &
00480 mct_gsmap_ngseg(gsmapo)/msizeo,mct_gsmap_ngseg(gsmapo)/apeso
00481 write(logunit,102) trim(subname),' nlseg/maxnlsegs = ', &
00482 mct_gsmap_nlseg(gsmapo,0),mct_gsmap_maxnlseg(gsmapo)
00483 102 format(2A,2I8)
00484 endif
00485
00486
00487
00488
00489
00490
00491 endif
00492
00493 end subroutine seq_rearr_gsmapCreate
00494
00495
00496
00497 subroutine seq_rearr_avExtend(AVin,IDin,ID)
00498
00499
00500
00501
00502
00503
00504 implicit none
00505 type(mct_aVect), intent(INOUT):: AVin
00506 integer ,intent(IN) :: IDin
00507 integer , intent(IN) :: ID
00508
00509
00510
00511 character(len=*),parameter :: subname = "(seq_rearr_avExtend) "
00512 integer :: mpicom
00513 integer :: rank,rank2
00514 integer :: lsizei, lsizen
00515 integer :: srank,srankg
00516 integer :: ierr
00517 integer :: nints
00518 character(len=1024) :: iList,rList
00519
00520
00521 call seq_comm_setptrs(ID,mpicom=mpicom,iam=rank)
00522
00523
00524
00525
00526 lsizei = -1
00527 if (seq_comm_iamin(IDin)) lsizei = mct_aVect_lsize(AVin)
00528 lsizen = 0
00529
00530
00531
00532
00533 srank = -1
00534 srankg = -1
00535 if (lsizei > 0) srank = rank
00536
00537 call mpi_allreduce(srank,srankg,1,MPI_INTEGER,MPI_MAX,mpicom,ierr)
00538 call shr_mpi_chkerr(ierr,subname//' mpi_allreduce max')
00539
00540 if (srankg < 0) then
00541 write(logunit,*) subname,' WARNING AVin empty '
00542 return
00543 endif
00544
00545
00546
00547
00548 iList = " "
00549 rList = " "
00550 if (rank == srankg) then
00551 if (mct_aVect_nIAttr(AVin) /= 0) iList = mct_aVect_ExportIList2c(AVin)
00552 if (mct_aVect_nRattr(AVin) /= 0) rList = mct_aVect_ExportRList2c(AVin)
00553 endif
00554
00555 call mpi_bcast(iList,len(iList),MPI_CHARACTER,srankg,mpicom,ierr)
00556 call mpi_bcast(rList,len(rList),MPI_CHARACTER,srankg,mpicom,ierr)
00557
00558
00559
00560
00561
00562 if (lsizei <= 0) then
00563 if(len_trim(iList) > 0 .and. len_trim(rList) > 0) then
00564 call mct_aVect_init(AVin,iList=iList,rList=rList,lsize=lsizen)
00565 elseif (len_trim(iList) > 0 .and. len_trim(rList) == 0) then
00566 call mct_aVect_init(AVin,iList=iList,lsize=lsizen)
00567 elseif (len_trim(iList) == 0 .and. len_trim(rList) > 0) then
00568 call mct_aVect_init(AVin,rList=rList,lsize=lsizen)
00569 endif
00570 endif
00571
00572 end subroutine seq_rearr_avExtend
00573
00574
00575
00576 subroutine seq_rearr_avCreate(AVin,IDin,AVout,ID,lsize)
00577
00578
00579
00580
00581
00582
00583 implicit none
00584 type(mct_aVect), intent(INOUT):: AVin
00585 integer ,intent(IN) :: IDin
00586 type(mct_aVect), intent(INOUT):: AVout
00587 integer , intent(IN) :: ID
00588 integer , intent(IN) :: lsize
00589
00590
00591
00592 character(len=*),parameter :: subname = "(seq_rearr_avCreate) "
00593 integer :: mpicom
00594 integer :: rank,rank2
00595 integer :: lsizei, lsizen
00596 integer :: srank,srankg
00597 integer :: ierr
00598 integer :: nints
00599 character(len=1024) :: iList,rList
00600
00601 call seq_comm_setptrs(ID,mpicom=mpicom,iam=rank)
00602
00603
00604
00605 lsizei = -1
00606 if (seq_comm_iamin(IDin)) lsizei = mct_aVect_lsize(AVin)
00607 lsizen = lsize
00608
00609
00610
00611
00612 srank = -1
00613 srankg = -1
00614 if (lsizei > 0) srank = rank
00615
00616 call mpi_allreduce(srank,srankg,1,MPI_INTEGER,MPI_MAX,mpicom,ierr)
00617 call shr_mpi_chkerr(ierr,subname//' mpi_allreduce max')
00618
00619 if (srankg < 0) then
00620 write(logunit,*) subname,' ERROR AVin not initialized '
00621 call shr_sys_abort()
00622 endif
00623
00624
00625
00626
00627 iList = " "
00628 rList = " "
00629 if (rank == srankg) then
00630 if (mct_aVect_nIAttr(AVin) /= 0) iList = mct_aVect_ExportIList2c(AVin)
00631 if (mct_aVect_nRattr(AVin) /= 0) rList = mct_aVect_ExportRList2c(AVin)
00632 endif
00633
00634 call mpi_bcast(iList,len(iList),MPI_CHARACTER,srankg,mpicom,ierr)
00635 call mpi_bcast(rList,len(rList),MPI_CHARACTER,srankg,mpicom,ierr)
00636
00637
00638
00639
00640 if(len_trim(iList) > 0 .and. len_trim(rList) > 0) then
00641 call mct_aVect_init(AVout,iList=iList,rList=rList,lsize=lsizen)
00642 elseif (len_trim(iList) > 0 .and. len_trim(rList) == 0) then
00643 call mct_aVect_init(AVout,iList=iList,lsize=lsizen)
00644 elseif (len_trim(iList) == 0 .and. len_trim(rList) > 0) then
00645 call mct_aVect_init(AVout,rList=rList,lsize=lsizen)
00646 endif
00647
00648 end subroutine seq_rearr_avCreate
00649
00650
00651 logical function seq_rearr_gsmapIdentical(gsmap1,gsmap2)
00652
00653 implicit none
00654 type(mct_gsMap), intent(IN):: gsmap1
00655 type(mct_gsMap), intent(IN):: gsmap2
00656
00657
00658
00659 character(len=*),parameter :: subname = "(seq_rearr_gsmapIdentical) "
00660 integer :: n
00661 logical :: identical
00662
00663
00664
00665 identical = .true.
00666
00667
00668 if (identical) then
00669 if (mct_gsMap_gsize(gsmap1) /= mct_gsMap_gsize(gsmap2)) identical = .false.
00670 if (mct_gsMap_ngseg(gsmap1) /= mct_gsMap_ngseg(gsmap2)) identical = .false.
00671 endif
00672
00673
00674 if (identical) then
00675 do n = 1,mct_gsMap_ngseg(gsmap1)
00676 if (gsmap1%start(n) /= gsmap2%start(n) ) identical = .false.
00677 if (gsmap1%length(n) /= gsmap2%length(n)) identical = .false.
00678 if (gsmap1%pe_loc(n) /= gsmap2%pe_loc(n)) identical = .false.
00679 enddo
00680 endif
00681
00682 seq_rearr_gsmapIdentical = identical
00683
00684 end function seq_rearr_gsmapIdentical
00685
00686
00687
00688 end module seq_rearr_mod