module my_mpi ! ! Module containing MPI routines needed for decomposition of ! geographic and magnetic grids, and for message passing required ! by the edynamo code. Buffer allocations are done once per run. ! use shr_kind_mod ,only: r8 => shr_kind_r8 ! 8-byte reals use params ,only: iulog use geogrid ,only: nlon,nlat,nlev use maggrid ,only: nmlon,nmlonp1,nmlat,nmlath,nmlev ! note nmlev is not a parameter use namelist ,only: model_name use timing implicit none #include public ! ! Number of MPI tasks and current task id (geo or mag): ! integer :: & ntask, & ! number of mpi tasks mytid ! my task id ! ! Geographic subdomains for current task: ! integer :: & ntaski, & ! number of tasks in lon dimension ntaskj, & ! number of tasks in lat dimension mytidi, & ! i coord for current task in task table mytidj, & ! j coord for current task in task table lat0,lat1, & ! first and last lats for each task lon0,lon1, & ! first and last lons for each task lev0,lev1, & ! first and last levs for each task (not distributed) mxlon,mxlat ! max number of subdomain lon,lat points among all tasks ! ! Magnetic subdomains for current task: ! integer :: & nmagtaski, & ! number of tasks in mag lon dimension nmagtaskj, & ! number of tasks in mag lat dimension magtidi, & ! i coord for current task in task table magtidj, & ! j coord for current task in task table mlat0,mlat1, & ! first and last mag lats for each task mlon0,mlon1, & ! first and last mag lons for each task mlev0,mlev1, & ! first and last mag levs (not distributed) mxmaglon, & ! max number of mag subdomain lon points among all tasks mxmaglat ! max number of mag subdomain lat points among all tasks integer,allocatable,save :: & itask_table_geo(:,:), & ! 2d table of tasks on geographic grid (i,j) itask_table_mag(:,:) ! 2d table of tasks on mag grid (i,j) integer :: cols_comm ! communicators for each task column integer :: rows_comm ! communicators for each task row ! ! Task type: subdomain information for all tasks, known by all tasks: ! type task integer :: mytid ! task id ! ! Geographic subdomains in task structure: integer :: mytidi ! task coord in longitude dimension of task table integer :: mytidj ! task coord in latitude dimension of task table integer :: nlats ! number of latitudes calculated by this task integer :: nlons ! number of longitudes calculated by this task integer :: lat0,lat1 ! first and last latitude indices integer :: lon0,lon1 ! first and last longitude indices ! ! Magnetic subdomains in task structure: integer :: magtidi ! task coord in mag longitude dimension of task table integer :: magtidj ! task coord in mag latitude dimension of task table integer :: nmaglats ! number of mag latitudes calculated by this task integer :: nmaglons ! number of mag longitudes calculated by this task integer :: mlat0,mlat1 ! first and last latitude indices integer :: mlon0,mlon1 ! first and last longitude indices end type task ! ! type(task) :: tasks(ntask) will be made available to all tasks ! (so each task has information about all tasks) ! type(task),allocatable,save :: tasks(:) ! ! Conjugate points in mag subdomains, for mp_mag_foldhem ! integer,allocatable,dimension(:),save :: & ! (ntask) nsend_south, & ! number of south lats to send to north (each task) nrecv_north ! number of north lats to send to south (each task) integer,allocatable,dimension(:,:),save :: & ! (mxlats,ntask) send_south_coords, & ! south j lats to send to north recv_north_coords ! north j lats to recv from south ! contains !----------------------------------------------------------------------- subroutine mp_init ! ! Initialize MPI, and allocate task table. ! integer :: ier call mpi_init(ier) call mpi_comm_size(MPI_COMM_WORLD,ntask,ier) call mpi_comm_rank(MPI_COMM_WORLD,mytid,ier) ! ! Allocate array of task structures: ! allocate(tasks(0:ntask-1),stat=ier) if (ier /= 0) then write(iulog,"('>>> mp_init: error allocating tasks(',i3,')')") ntask call shutdown('tasks') endif end subroutine mp_init !----------------------------------------------------------------------- subroutine mp_decomp(model) ! ! Decompose geographic and geomagnetic domains: ! ! Args: character(len=*),intent(in) :: model ! ! Local: integer :: ier if (model=='TIMEGCM') then call mkntask(ntask,nlon+4,nlat,ntaskj,ntaski,ier) call mp_distribute_geo(nlon+4,nlat) else call mkntask(ntask,nlon,nlat,ntaskj,ntaski,ier) call mp_distribute_geo(nlon,nlat) endif call mkntask(ntask,nmlonp1,nmlat,nmagtaskj,nmagtaski,ier) call mp_distribute_mag call mp_exchange_tasks(1) end subroutine mp_decomp !----------------------------------------------------------------------- subroutine mp_close integer :: ier ! call mpi_finalize(ier) if (ier /= 0) then write(iulog,"(/,'>>> WARNING: error from mp_finalize: ier=',i3)")ier call shutdown('mp_finalize error') endif end subroutine mp_close !----------------------------------------------------------------------- subroutine mp_distribute_geo(nlon,nlat) ! ! Args: integer,intent(in) :: nlon,nlat ! ! Local: integer :: i,j,n,irank,ier,tidrow ! ! Vertical dimension is not distributed, so every mpi task gets the full column. ! For TGCM, k-index lev0==1 is the bottom boundary, and k-index lev1=nlev is ! the top boundary, but this is reversed in CAM/WACCM. ! ! Note that nlev as read from TIMEGCM history file (the lev dimension) ! (module geogrid) is the same as nlevp1 in the timegcm model, so ! nlev=nlevp1=lev=49 (nlevp1 is not used in edynamo) ! lev0 = 1 lev1 = nlev ! ! Allocate and set 2d table of tasks: ! allocate(itask_table_geo(-1:ntaski,-1:ntaskj),stat=ier) if (ier /= 0) then write(iulog,"('>>> Error allocating itable: ntaski,j=',2i4)") ntaski,ntaskj call shutdown('itask_table_geo') endif itask_table_geo(:,:) = MPI_PROC_NULL irank = 0 do j = 0,ntaskj-1 do i = 0,ntaski-1 itask_table_geo(i,j) = irank if (mytid == irank) then mytidi = i mytidj = j endif irank = irank+1 enddo enddo ! ! Calculate starting and ending indices in lon,lat dimensions for each task: ! call distribute_1d(1,nlon,ntaski,mytidi,lon0,lon1) call distribute_1d(1,nlat,ntaskj,mytidj,lat0,lat1) ! ! Define all task structures with current task values ! (redundant for alltoall): ! do n=0,ntask-1 tasks(n)%mytid = mytid tasks(n)%mytidi = mytidi tasks(n)%mytidj = mytidj tasks(n)%nlats = lat1-lat0+1 tasks(n)%nlons = lon1-lon0+1 tasks(n)%lat0 = lat0 tasks(n)%lat1 = lat1 tasks(n)%lon0 = lon0 tasks(n)%lon1 = lon1 enddo ! ! All tasks must have at least 4 longitudes: ! do n=0,ntask-1 if (tasks(n)%nlons < 4) then write(iulog,"(/,i2,i2)") '>>> mp_distribute_geo: each task must carry ',& 'at least 4 longitudes. task=',n,' nlons=',tasks(n)%nlons call shutdown('nlons per task') endif enddo ! ! Create sub-communicators for each task row (used by mp_geopole_3d): ! call mpi_comm_split(MPI_COMM_WORLD,mytidj,mytid,rows_comm,ier) call MPI_Comm_rank(rows_comm,tidrow,ier) end subroutine mp_distribute_geo !----------------------------------------------------------------------- subroutine mp_distribute_mag ! ! Local: integer :: i,j,n,irank,ier,tidcol ! ! Vertical dimension is not distributed: mlev0 = 1 mlev1 = nmlev ! ! Allocate and set 2d table of tasks: allocate(itask_table_mag(-1:nmagtaski,-1:nmagtaskj),stat=ier) if (ier /= 0) then write(iulog,"('>>> Error allocating itable: nmagtaski,j=',2i3)") & nmagtaski,nmagtaskj call shutdown('itask_table_mag') endif itask_table_mag(:,:) = MPI_PROC_NULL irank = 0 do j = 0,nmagtaskj-1 do i = 0,nmagtaski-1 itask_table_mag(i,j) = irank if (mytid == irank) then magtidi = i magtidj = j endif irank = irank+1 enddo ! ! Tasks are periodic in longitude: ! itask_table_mag(-1,j) = itask_table_mag(nmagtaski-1,j) itask_table_mag(nmagtaski,j) = itask_table_mag(0,j) enddo ! ! Calculate start and end indices in mag lon,lat dimensions for each task: ! call distribute_1d(1,nmlonp1,nmagtaski,magtidi,mlon0,mlon1) call distribute_1d(1,nmlat ,nmagtaskj,magtidj,mlat0,mlat1) ! ! Define all task structures with current task values ! (redundant for alltoall): ! do n=0,ntask-1 tasks(n)%magtidi = magtidi tasks(n)%magtidj = magtidj tasks(n)%nmaglats = mlat1-mlat0+1 tasks(n)%nmaglons = mlon1-mlon0+1 tasks(n)%mlat0 = mlat0 tasks(n)%mlat1 = mlat1 tasks(n)%mlon0 = mlon0 tasks(n)%mlon1 = mlon1 enddo ! ! All tasks must have at least 4 longitudes: do n=0,ntask-1 if (tasks(n)%nmaglons < 4) then write(iulog,"(/,i2,i2)") '>>> mp_distribute_mag: each task must carry ',& 'at least 4 longitudes. task=',n,' nmaglons=',tasks(n)%nmaglons call shutdown('nmaglons per task') endif enddo ! ! Create subgroup communicators for each task column: ! These communicators will be used by sub mp_mag_jslot (mpi.F). ! call mpi_comm_split(MPI_COMM_WORLD,mod(mytid,nmagtaski),mytid,& cols_comm,ier) call MPI_Comm_rank(cols_comm,tidcol,ier) end subroutine mp_distribute_mag !----------------------------------------------------------------------- subroutine distribute_1d(n1,n2,nprocs,myrank,istart,iend) ! ! Distribute work across a 1d vector(n1->n2) to nprocs. ! Return start and end indices for proc myrank. ! ! Args: integer,intent(in) :: n1,n2,nprocs,myrank integer,intent(out) :: istart,iend ! ! Local: integer :: lenproc,iremain,n ! n = n2-n1+1 lenproc = n/nprocs iremain = mod(n,nprocs) istart = n1 + myrank*lenproc + min(myrank,iremain) iend = istart+lenproc-1 if (iremain > myrank) iend = iend+1 end subroutine distribute_1d !----------------------------------------------------------------------- subroutine mkntask(ntask,nlon,nlat,ntaskj,ntaski,ier) ! ! Args: integer,intent(in) :: ntask,nlon,nlat integer,intent(out) :: ntaskj,ntaski,ier ! ! Local: integer :: i,j,ntlat(nlon),ntlon(nlon),nchoice,ngap,& ngap_prev,ii ntaskj = 0 ntaski = 0 nchoice=0 do i=1,nlon do j=1,nlat if (i*j==ntask) then nchoice = nchoice+1 ntlat(nchoice)=j ntlon(nchoice)=i endif enddo enddo ! ! Choose combinations of ntlat*ntlon==ntask in which ntlon==ntlat, ! or if that does not exist, use combination with smallest delta ! ntlon-ntlat. ! do i=1,nchoice if (ntlon(i) == ntlat(i)) then ntaskj = ntlat(i) ntaski = ntlon(i) endif enddo if (ntaskj==0.or.ntaski==0) then ngap_prev = nlon*nlat do i=1,nchoice ngap = ntlon(i)-ntlat(i) if (abs(ngap) < ngap_prev) then ngap = ntlon(i)-ntlat(i) ngap_prev = abs(ngap) ii = i endif enddo ntaskj = ntlat(ii) ntaski = ntlon(ii) ! ! If they are not equal, ntlon should be > ntlat because nlon > nlat. if (ntaski < ntaskj) then i = ntaskj ntaskj = ntaski ntaski = i endif endif ier =0 if (ntaskj==0.or.ntaski==0.or.ntaskj*ntaski /= ntask) & ier = 1 if (ntaskj*ntaski /= ntask) then write(iulog,"(i4,i4,i4)") '>>> mkntask: ntaskj*ntaski /= ntask: ',& ' ntaskj=',ntaskj,' ntaski=',ntaski,' ntask=',ntask stop 'mkntask' endif end subroutine mkntask !----------------------------------------------------------------------- subroutine mp_exchange_tasks(iprint) ! ! Args: integer,intent(in) :: iprint ! ! Local: ! itasks_send(len_task_type,ntask) will be used to send tasks(:) info ! to all tasks (directly passing mpi derived data types is reportedly ! not stable, or not available until MPI 2.x). ! integer :: n,ier,i integer,parameter :: len_task_type = 17 ! see type task above integer,allocatable,save :: & itasks_send(:,:), & ! send buffer itasks_recv(:,:) ! send buffer ! ! Pack tasks(mytid) into itasks_send: allocate(itasks_send(len_task_type,0:ntask-1),stat=ier) if (ier /= 0) then write(iulog,"(i4,i4)") '>>> Error allocating itasks_send: len_task_type=',& len_task_type,' ntask=',ntask endif allocate(itasks_recv(len_task_type,0:ntask-1),stat=ier) if (ier /= 0) then write(iulog,"(i4,i4)") '>>> Error allocating itasks_recv: len_task_type=',& len_task_type,' ntask=',ntask endif do n=0,ntask-1 itasks_send(1,n) = tasks(mytid)%mytid itasks_send(2,n) = tasks(mytid)%mytidi itasks_send(3,n) = tasks(mytid)%mytidj itasks_send(4,n) = tasks(mytid)%nlats itasks_send(5,n) = tasks(mytid)%nlons itasks_send(6,n) = tasks(mytid)%lat0 itasks_send(7,n) = tasks(mytid)%lat1 itasks_send(8,n) = tasks(mytid)%lon0 itasks_send(9,n) = tasks(mytid)%lon1 itasks_send(10,n) = tasks(mytid)%magtidi itasks_send(11,n) = tasks(mytid)%magtidj itasks_send(12,n) = tasks(mytid)%nmaglats itasks_send(13,n) = tasks(mytid)%nmaglons itasks_send(14,n) = tasks(mytid)%mlat0 itasks_send(15,n) = tasks(mytid)%mlat1 itasks_send(16,n) = tasks(mytid)%mlon0 itasks_send(17,n) = tasks(mytid)%mlon1 enddo ! ! Send itasks_send and receive itasks_recv: call mpi_alltoall(itasks_send,len_task_type,MPI_INTEGER,& itasks_recv,len_task_type,MPI_INTEGER,& MPI_COMM_WORLD,ier) if (ier /= 0) & call handle_mpi_err(ier,'mpi_alltoall to send/recv itasks') ! ! Unpack itasks_recv into tasks(n) ! do n=0,ntask-1 tasks(n)%mytid = itasks_recv(1,n) tasks(n)%mytidi = itasks_recv(2,n) tasks(n)%mytidj = itasks_recv(3,n) tasks(n)%nlats = itasks_recv(4,n) tasks(n)%nlons = itasks_recv(5,n) tasks(n)%lat0 = itasks_recv(6,n) tasks(n)%lat1 = itasks_recv(7,n) tasks(n)%lon0 = itasks_recv(8,n) tasks(n)%lon1 = itasks_recv(9,n) tasks(n)%magtidi = itasks_recv(10,n) tasks(n)%magtidj = itasks_recv(11,n) tasks(n)%nmaglats = itasks_recv(12,n) tasks(n)%nmaglons = itasks_recv(13,n) tasks(n)%mlat0 = itasks_recv(14,n) tasks(n)%mlat1 = itasks_recv(15,n) tasks(n)%mlon0 = itasks_recv(16,n) tasks(n)%mlon1 = itasks_recv(17,n) ! ! Report to stdout: ! if (n==mytid.and.iprint > 0) then write(iulog,"(/,'Task ',i3,':')") n write(iulog,"(/,'Subdomain on geographic grid:')") write(iulog,"('tasks(',i3,')%mytid =',i3)") n,tasks(n)%mytid write(iulog,"('tasks(',i3,')%mytidi=',i3)") n,tasks(n)%mytidi write(iulog,"('tasks(',i3,')%mytidj=',i3)") n,tasks(n)%mytidj write(iulog,"('tasks(',i3,')%nlats =',i3)") n,tasks(n)%nlats write(iulog,"('tasks(',i3,')%nlons =',i3)") n,tasks(n)%nlons write(iulog,"('tasks(',i3,')%lat0 =',i3)") n,tasks(n)%lat0 write(iulog,"('tasks(',i3,')%lat1 =',i3)") n,tasks(n)%lat1 write(iulog,"('tasks(',i3,')%lon0 =',i3)") n,tasks(n)%lon0 write(iulog,"('tasks(',i3,')%lon1 =',i3)") n,tasks(n)%lon1 write(iulog,"('Number of geo subdomain grid points = ',i6)") & tasks(n)%nlons * tasks(n)%nlats write(iulog,"(/,'Subdomain on geomagnetic grid:')") write(iulog,"('tasks(',i3,')%magtidi=',i3)") n,tasks(n)%magtidi write(iulog,"('tasks(',i3,')%magtidj=',i3)") n,tasks(n)%magtidj write(iulog,"('tasks(',i3,')%nmaglats =',i3)") n,tasks(n)%nmaglats write(iulog,"('tasks(',i3,')%nmaglons =',i3)") n,tasks(n)%nmaglons write(iulog,"('tasks(',i3,')%mlat0 =',i3)") n,tasks(n)%mlat0 write(iulog,"('tasks(',i3,')%mlat1 =',i3)") n,tasks(n)%mlat1 write(iulog,"('tasks(',i3,')%mlon0 =',i3)") n,tasks(n)%mlon0 write(iulog,"('tasks(',i3,')%mlon1 =',i3)") n,tasks(n)%mlon1 write(iulog,"('Number of mag subdomain grid points = ',i6)") & tasks(n)%nmaglons * tasks(n)%nmaglats endif enddo ! ! Release locally allocated space: ! debug: temporarilly remove these deallocs while using getoverstepped() deallocate(itasks_send) deallocate(itasks_recv) ! ! mxlon,mxlat are maximum number of lons,lats owned by all tasks: mxlon = -9999 do n=0,ntask-1 if (tasks(n)%nlons > mxlon) mxlon = tasks(n)%nlons enddo mxlat = -9999 do n=0,ntask-1 if (tasks(n)%nlats > mxlat) mxlat = tasks(n)%nlats enddo ! ! mxmaglon,mxmaglat are maximum number of mag lons,lats owned by all tasks: mxmaglon = -9999 do n=0,ntask-1 if (tasks(n)%nmaglons > mxmaglon) mxmaglon = tasks(n)%nmaglons enddo mxmaglat = -9999 do n=0,ntask-1 if (tasks(n)%nmaglats > mxmaglat) mxmaglat = tasks(n)%nmaglats enddo ! ! Find conjugate points for folding hemispheres: call conjugate_points end subroutine mp_exchange_tasks !----------------------------------------------------------------------- subroutine mp_geopole_3d(f,ilon0,ilon1,ilat0,ilat1,nlevs,& nglblon,jspole,jnpole,fpole_jpm2,nf,fnames) ! ! Args: integer,intent(in) :: ilon0,ilon1,ilat0,ilat1,nlevs,nglblon,& jspole,jnpole,nf real(r8),intent(in) :: f(ilon0:ilon1,ilat0:ilat1,nlevs,nf) real(r8),intent(out) :: fpole_jpm2(nglblon,4,nlevs,nf) character(len=*),intent(in) :: fnames(nf) ! ! Local: integer :: ier,tidrow,len,j,n,k,itask,i0,i1,j0,j1,wid real(r8),allocatable,save :: sndbuf(:,:,:,:,:),rcvbuf(:,:,:,:,:) integer,parameter :: mxnf=7 real(r8) :: starttime,endtime logical :: found1,found2 if (mpi_timing) starttime = mpi_wtime() ! ! Max number fields allowed is mxnf (local parameter): if (nf > mxnf) then write(iulog,"(a,i4,a,i4)") '>>> mp_geopole_3d: nf=',nf,& ' but cannot be called with greater than mxnf=',mxnf call shutdown('mp_geopole_3d') endif ! ! Need only northern-most and southern-most processor rows, ! all other tasks return: ! if (mytidj /= 0.and.mytidj /= ntaskj-1) return ! ! Confirm that lats 1,2 and nlat-1,nlat are available on the ! remaining tasks: ! found1=.false. ; found2=.false. do j=ilat0,ilat1 if (j==jspole) then ! south pole if (ilat0 > jspole+2) & ! need lats 1,2 for south pole call shutdown('mp_geopole_3d: ilat1 < 2') found1=.true. elseif (j==jnpole) then ! north pole if (ilat1 < jnpole-2) & ! need top 2 lats for north pole call shutdown('mp_geopole_3d: ilat1 > jnpole-2') found2=.true. endif enddo if ((mytidj==0.and..not.found1).or.(mytidj==ntaskj-1.and..not.found2)) then write(iulog,"('>>> Could not find poles: ilat0,1=',2i4,' jspole,jnpole=',2i4,' jspole+2=',i4,' jnpole-2=',i4)") & ilat0,ilat1,jspole,jnpole,jspole+2,jnpole-2 call shutdown('mp_geopole_3d: need poles') endif ! ! Using rows_comm communicators, so only need ntaski tasks: ! Note mxnf dimension, so this can be called w/ different number ! of fields, but only allocate once. ! if (.not.allocated(sndbuf)) then allocate(sndbuf(nglblon,4,nlevs,mxnf,0:ntaski-1),stat=ier) if (ier /= 0) call shutdown('mp_geopole_3d: error allocating sndbuf') endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(nglblon,4,nlevs,mxnf,0:ntaski-1),stat=ier) if (ier /= 0) call shutdown('mp_geopole_3d: error allocating rcvbuf') endif sndbuf = 0._r8 rcvbuf = 0._r8 ! ! Messages are built for only nf fields, not mxnf: len = nglblon*nlevs*4*nf ! ! Task id of my row communicator: call MPI_Comm_rank(rows_comm,tidrow,ier) wid = itask_table_geo(tidrow,mytidj) ! world-comm id for this task ! ! Init and load send buffers: fpole_jpm2=0._r8 do itask=0,ntaski-1 do n=1,nf do j=ilat0,ilat1 if (j==jspole+1) then ! south pole +1 sndbuf(ilon0:ilon1,1,:,n,itask) = f(ilon0:ilon1,j,:,n) elseif (j==jspole+2) then ! south pole +2 sndbuf(ilon0:ilon1,2,:,n,itask) = f(ilon0:ilon1,j,:,n) elseif (j==jnpole-1) then ! north pole -1 sndbuf(ilon0:ilon1,3,:,n,itask) = f(ilon0:ilon1,j,:,n) elseif (j==jnpole-2) then ! north pole -2 sndbuf(ilon0:ilon1,4,:,n,itask) = f(ilon0:ilon1,j,:,n) endif enddo ! j=ilat0,ilat1 enddo ! n=1,nf enddo ! itask=0,ntaski-1 ! ! Do the exchange within my row: call mpi_alltoall(sndbuf(:,:,:,1:nf,:),len,MPI_REAL8,& rcvbuf(:,:,:,1:nf,:),len,MPI_REAL8,rows_comm,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_geopole_3d call mpi_alltoall') ! ! Unpack rcvbuf into output fpole_jpm2(nglblon,4,nlevs,nf) do itask=0,ntaski-1 i0 = tasks(itask)%lon0 i1 = tasks(itask)%lon1 if (trim(model_name)=='TIMEGCM') then if (i0 /= 1) i0 = i0-2 if (i1 /= nlon+4) then i1 = i1-2 else i1 = i1-4 endif endif j0 = tasks(wid)%lat0 j1 = tasks(wid)%lat1 do n=1,nf do j=j0,j1 if (j==jspole+1) then ! south pole +1 fpole_jpm2(i0:i1,1,:,n) = rcvbuf(i0:i1,1,:,n,itask) elseif (j==jspole+2) then ! south pole +2 fpole_jpm2(i0:i1,2,:,n) = rcvbuf(i0:i1,2,:,n,itask) elseif (j==jnpole-1) then ! north pole -1 fpole_jpm2(i0:i1,3,:,n) = rcvbuf(i0:i1,3,:,n,itask) elseif (j==jnpole-2) then ! north pole -2 fpole_jpm2(i0:i1,4,:,n) = rcvbuf(i0:i1,4,:,n,itask) endif enddo ! do j=j0,j1 enddo ! do n=1,nf enddo ! itask=0,ntaski-1 if (mpi_timing) then endtime = mpi_wtime() time_geopole_3d=time_geopole_3d+(endtime-starttime) endif end subroutine mp_geopole_3d !----------------------------------------------------------------------- subroutine mp_geopole_2d(f,ilon0,ilon1,ilat0,ilat1,& nglblon,jspole,jnpole,fpole_jpm2,nf) ! ! Return fpole_jpm2(nglblon,1->4,nf) as: ! 1: j = jspole+1 (spole+1) ! 2: j = jspole+2 (spole+2) ! 3: j = jnpole-1 (npole-1) ! 4: j = jnpole-2 (npole-2) ! ! Args: integer,intent(in) :: ilon0,ilon1,ilat0,ilat1,nglblon,jspole,jnpole,nf real(r8),intent(in) :: f(ilon0:ilon1,ilat0:ilat1,nf) real(r8),intent(out) :: fpole_jpm2(nglblon,4,nf) ! ! Local: integer :: j,k,ier,len,itask,i0,i1,j0,j1 real(r8),allocatable,save :: sndbuf(:,:,:,:),rcvbuf(:,:,:,:) real(r8) :: starttime,endtime if (mpi_timing) starttime = mpi_wtime() ! if (.not.allocated(sndbuf)) then allocate(sndbuf(nglblon,4,nf,0:ntask-1),stat=ier) if (ier /= 0) call shutdown('mp_geopole_2d: error allocating sndbuf') endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(nglblon,4,nf,0:ntask-1),stat=ier) if (ier /= 0) call shutdown('mp_geopole_2d: error allocating rcvbuf') endif sndbuf = 0._r8 rcvbuf = 0._r8 len = nglblon*4*nf ! ! Load send buffer with values at poles +/- 2 for current subdomain ! (redundant to all tasks for alltoall) ! do itask=0,ntask-1 do j=ilat0,ilat1 if (j==jspole+1) then ! south pole +1 sndbuf(ilon0:ilon1,1,:,itask) = f(ilon0:ilon1,j,:) elseif (j==jspole+2) then ! south pole +2 sndbuf(ilon0:ilon1,2,:,itask) = f(ilon0:ilon1,j,:) elseif (j==jnpole-1) then ! north pole -1 sndbuf(ilon0:ilon1,3,:,itask) = f(ilon0:ilon1,j,:) elseif (j==jnpole-2) then ! north pole -2 sndbuf(ilon0:ilon1,4,:,itask) = f(ilon0:ilon1,j,:) endif enddo enddo ! ! Do the exchange: call mpi_alltoall(sndbuf,len,MPI_REAL8, & rcvbuf,len,MPI_REAL8,MPI_COMM_WORLD,ier) if (ier /= 0) & call handle_mpi_err(ier,'mp_geopole_2d call mpi_alltoall') ! ! Unpack rcvbuf into output fpole_jpm2(nglblon,4,nf) do itask=0,ntask-1 i0 = tasks(itask)%lon0 i1 = tasks(itask)%lon1 if (trim(model_name)=='TIMEGCM') then if (i0 /= 1) i0 = i0-2 if (i1 /= nlon+4) then i1 = i1-2 else i1 = i1-4 endif endif j0 = tasks(itask)%lat0 j1 = tasks(itask)%lat1 do j=j0,j1 if (j==jspole+1) then ! south pole +1 fpole_jpm2(i0:i1,1,:) = rcvbuf(i0:i1,1,:,itask) elseif (j==jspole+2) then ! south pole +2 fpole_jpm2(i0:i1,2,:) = rcvbuf(i0:i1,2,:,itask) elseif (j==jnpole-1) then ! north pole -1 fpole_jpm2(i0:i1,3,:) = rcvbuf(i0:i1,3,:,itask) elseif (j==jnpole-2) then ! north pole -2 fpole_jpm2(i0:i1,4,:) = rcvbuf(i0:i1,4,:,itask) endif enddo enddo if (mpi_timing) then endtime = mpi_wtime() time_geopole_2d=time_geopole_2d+(endtime-starttime) endif end subroutine mp_geopole_2d !----------------------------------------------------------------------- subroutine mp_periodic_f3d(f,lev0,lev1,lon0,lon1,lat0,lat1,nf) ! ! Define periodic points for field f(nlevels,lon0:lon1,lat0:lat1): ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 ! ! Args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1,nf real(r8),intent(inout) :: f(lev0:lev1,lon0:lon1,lat0:lat1,nf) ! ! Local: integer :: k,j,idest,isrc,isend,irecv,ier,len,nlevs,nlons,nlats,n real(r8),allocatable,save :: sndbuf(:,:,:,:),rcvbuf(:,:,:,:) ! (nlevs,2,mxlat,nf) integer :: nfsave=0 real(r8) :: starttime,endtime integer :: irstat(MPI_STATUS_SIZE) ! mpi receive status ! if (mpi_timing) starttime = mpi_wtime() ! ! Serial run of an MPI job (1 task): if (ntask==1) then do n=1,nf do j=lat0,lat1 f(:,lon0:lon0+1,j,n) = f(:,nlon+1:nlon+2,j,n) f(:,lon1-1:lon1,j,n) = f(:,lon0+2:lon0+3,j,n) enddo ! j=lat0,lat1 enddo return endif ! ! The restriction that every task must carry at least 4 longitudes ! means that tasks with mytidi==0 will have lons 1->4 and tasks ! with mytidi==ntaski-1 will have lons nlonp4-3->nlonp4. ! if (mytidi /= 0 .and. mytidi /= ntaski-1) return ! ! Allocate send and receive buffers: nlevs = lev1-lev0+1 nlons = lon1-lon0+1 nlats = lat1-lat0+1 if (nf > nfsave) then if (allocated(sndbuf)) deallocate(sndbuf) if (allocated(rcvbuf)) deallocate(rcvbuf) nfsave = nf endif if (.not.allocated(sndbuf)) then allocate(sndbuf(nlevs,2,mxlat,nf),stat=ier) if (ier /= 0) call shutdown('mp_periodic_f3d: error allocating sndbuf') endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(nlevs,2,mxlat,nf),stat=ier) if (ier /= 0) call shutdown('mp_periodic_f3d: error allocating rcvbuf') endif sndbuf = 0._r8 ; rcvbuf = 0._r8 len = nlevs*2*mxlat*nf ! ! The restriction that every task must carry at least 4 longitudes ! means that tasks with mytidi==0 will have lons 1->4 and tasks ! with mytidi==ntaski-1 will have lons nlonp4-3->nlonp4. ! if (mytidi /= 0 .and. mytidi /= ntaski-1) return ! ! Allocate send and receive buffers: nlevs = lev1-lev0+1 nlons = lon1-lon0+1 nlats = lat1-lat0+1 if (nf > nfsave) then if (allocated(sndbuf)) deallocate(sndbuf) if (allocated(rcvbuf)) deallocate(rcvbuf) nfsave = nf endif if (.not.allocated(sndbuf)) then allocate(sndbuf(nlevs,2,mxlat,nf),stat=ier) if (ier /= 0) call shutdown('mp_periodic_f3d: error allocating sndbuf') endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(nlevs,2,mxlat,nf),stat=ier) if (ier /= 0) call shutdown('mp_periodic_f3d: error allocating rcvbuf') endif sndbuf = 0._r8 ; rcvbuf = 0._r8 len = nlevs*2*mxlat*nf ! ! Send lons 3,4 to lons nlonp4-1,nlonp4 at task with mytidi==ntaski-1, ! and receive lons 1,2 from the same task. if (mytidi==0) then idest = itask_table_geo(ntaski-1,mytidj) do n=1,nf do j=lat0,lat1 do k=lev0,lev1 sndbuf(k-lev0+1,1:2,j-lat0+1,n) = f(k,3:4,j,n) enddo enddo ! j=lat0,lat1 enddo call mpi_isend(sndbuf(:,:,:,1:nf),len,MPI_REAL8,idest,1,& MPI_COMM_WORLD,isend,ier) if (ier /= 0) & call handle_mpi_err(ier,'mp_periodic_f3d send to idest') isrc = idest call mpi_irecv(rcvbuf(:,:,:,1:nf),len,MPI_REAL8,isrc,1,& MPI_COMM_WORLD,irecv,ier) if (ier /= 0) & call handle_mpi_err(ier,'mp_periodic_f3d recv fm isrc') ! ! Send lons nlonp4-3,nlonp4-2 to lons 1,2 at task with mytidi==0, ! and receive lons nlonp4-1,nlonp4 from same task: elseif (mytidi==ntaski-1) then idest = itask_table_geo(0,mytidj) do n=1,nf do j=lat0,lat1 do k=lev0,lev1 sndbuf(k-lev0+1,1:2,j-lat0+1,n) = & f(k,nlon+1:nlon+2,j,n) enddo enddo ! j=lat0,lat1 enddo call mpi_isend(sndbuf(:,:,:,1:nf),len,MPI_REAL8,idest,1, & MPI_COMM_WORLD,isend,ier) if (ier /= 0) & call handle_mpi_err(ier,'mp_periodic_f3d send to idest') isrc = idest call mpi_irecv(rcvbuf(:,:,:,1:nf),len,MPI_REAL8,isrc,1, & MPI_COMM_WORLD,irecv,ier) if (ier /= 0) & call handle_mpi_err(ier,'mp_periodic_f3d recv fm isrc') endif ! ! Wait for completions: call mpi_wait(isend,irstat,ier) if (ier /= 0) & call handle_mpi_err(ier,'mp_periodic_f3d wait for send') call mpi_wait(irecv,irstat,ier) if (ier /= 0) & call handle_mpi_err(ier,'mp_periodic_f3d wait for recv') ! ! Copy from receive buffer: if (mytidi==0) then ! I am on west end do n=1,nf do j=lat0,lat1 do k=lev0,lev1 f(k,1:2,j,n) = rcvbuf(k-lev0+1,1:2,j-lat0+1,n) enddo ! k=lev0,lev1 enddo ! j=lat0,lat1 enddo elseif (mytidi==ntaski-1) then ! I am on east end do n=1,nf do j=lat0,lat1 do k=lev0,lev1 f(k,nlon+3:nlon+4,j,n) = rcvbuf(k-lev0+1,1:2,j-lat0+1,n) enddo ! k=lev0,lev1 enddo ! j=lat0,lat1 enddo endif if (mpi_timing) then endtime = mpi_wtime() time_periodic_f3d=time_periodic_f3d+(endtime-starttime) endif end subroutine mp_periodic_f3d !----------------------------------------------------------------------- subroutine mp_periodic_f2d(f,lon0,lon1,lat0,lat1,nf) ! ! Define periodic points for field f(lon0:lon1,lat0:lat1): ! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4 ! ! Args: integer,intent(in) :: lon0,lon1,lat0,lat1 real(r8),intent(inout) :: f(lon0:lon1,lat0:lat1,nf) ! ! Local: integer :: j,idest,isrc,isend,irecv,ier,len,nlons,nlats,nf,n real(r8),allocatable,save :: sndbuf(:,:,:),rcvbuf(:,:,:) ! (2,nlats,nf) integer :: nfsave=0 real(r8) :: starttime,endtime integer :: irstat(MPI_STATUS_SIZE) ! mpi receive status if (mpi_timing) starttime = mpi_wtime() ! ! The restriction that every task must carry at least 4 longitudes ! means that tasks with mytidi==0 will have lons 1->4 and tasks ! with mytidi==ntaski-1 will have lons nlonp4-3->nlonp4. ! if (mytidi /= 0 .and. mytidi /= ntaski-1) return ! ! Allocate send and receive buffers: nlons = lon1-lon0+1 nlats = lat1-lat0+1 if (nf > nfsave) then if (allocated(sndbuf)) deallocate(sndbuf) if (allocated(rcvbuf)) deallocate(rcvbuf) nfsave = nf endif if (.not.allocated(sndbuf)) then allocate(sndbuf(2,mxlat,nf),stat=ier) if (ier /= 0) & call shutdown('mp_periodic_f2d: error allocating sndbuf') endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(2,mxlat,nf),stat=ier) if (ier /= 0) & call shutdown('mp_periodic_f2d: error allocating rcvbuf') endif sndbuf = 0._r8 ; rcvbuf = 0._r8 len = 2*mxlat*nf ! ! Send lons 3,4 to lons nlonp4-1,nlonp4 at task with mytidi==ntaski-1, ! and receive lons 1,2 from the same task. if (mytidi==0) then idest = itask_table_geo(ntaski-1,mytidj) do n=1,nf do j=lat0,lat1 sndbuf(1:2,j-lat0+1,n) = f(3:4,j,n) enddo ! j=lat0,lat1 enddo call mpi_isend(sndbuf(:,:,1:nf),len,MPI_REAL8,idest,1, & MPI_COMM_WORLD,isend,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_periodic_f2d send to idest') isrc = idest call mpi_irecv(rcvbuf(:,:,1:nf),len,MPI_REAL8,isrc,1, & MPI_COMM_WORLD,irecv,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_periodic_f2d recv fm isrc') ! ! Send lons nlonp4-3,nlonp4-2 to lons 1,2 at task with mytidi==0, ! and receive lons nlonp4-1,nlonp4 from same task: elseif (mytidi==ntaski-1) then idest = itask_table_geo(0,mytidj) do n=1,nf do j=lat0,lat1 sndbuf(1:2,j-lat0+1,n) = f(nlon+1:nlon+2,j,n) enddo ! j=lat0,lat1 enddo call mpi_isend(sndbuf(:,:,1:nf),len,MPI_REAL8,idest,1, & MPI_COMM_WORLD,isend,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_periodic_f2d send to idest') isrc = idest call mpi_irecv(rcvbuf(:,:,1:nf),len,MPI_REAL8,isrc,1, & MPI_COMM_WORLD,irecv,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_periodic_f2d recv fm isrc') endif ! ! Wait for completions: call mpi_wait(isend,irstat,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_periodic_f2d wait for send') call mpi_wait(irecv,irstat,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_periodic_f2d wait for recv') ! ! Copy from receive buffer: if (mytidi==0) then do n=1,nf do j=lat0,lat1 f(1:2,j,n) = rcvbuf(1:2,j-lat0+1,n) enddo ! j=lat0,lat1 enddo elseif (mytidi==ntaski-1) then do n=1,nf do j=lat0,lat1 f(nlon+3:nlon+4,j,n) = rcvbuf(1:2,j-lat0+1,n) enddo ! j=lat0,lat1 enddo endif if (mpi_timing) then endtime = mpi_wtime() time_periodic_f2d=time_periodic_f2d+(endtime-starttime) endif end subroutine mp_periodic_f2d !----------------------------------------------------------------------- subroutine mp_mageq(fin,fout,nf,mlon0,mlon1,mlat0,mlat1,nmlev) ! ! Each task needs values of conductivities and adotv1,2 fields at the ! at the mag equator for its longitude subdomain (and all levels), for ! the fieldline integrations. ! ! On input, fin is ped_mag, hal_mag, adotv1_mag, adotv2_mag ! on (i,j,k) magnetic subdomain. ! On output, fout(mlon0:mlon1,nmlev,nf) is ped_meq, hal_meq, adotv1_meq, ! adotv2_meq at mag equator at longitude subdomain and all levels. ! ! Args: integer :: mlon0,mlon1,mlat0,mlat1,nmlev,nf real(r8),intent(in) :: fin (mlon0:mlon1,mlat0:mlat1,nmlev,nf) real(r8),intent(out) :: fout(mlon0:mlon1,nmlev,nf) ! ! Local: real(r8),allocatable,save :: & ! mpi buffers sndbuf(:,:,:), & ! mxmaglon,nmlev,nf rcvbuf(:,:,:) ! mxmaglon,nmlev,nf integer :: i,j,n,itask,ier,len,jlateq,ireqsend,ireqrecv integer :: irstat(MPI_STATUS_SIZE) ! mpi receive status real(r8) :: starttime,endtime logical :: have_eq if (mpi_timing) starttime = mpi_wtime() ! ! Allocate buffers once per run: if (.not.allocated(sndbuf)) then allocate(sndbuf(mxmaglon,nmlev,nf),stat=ier) if (ier /= 0) then write(iulog,"('>>> mp_mageq: error allocating sndbuf')") call shutdown('mp_mageq') endif endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(mxmaglon,nmlev,nf),stat=ier) if (ier /= 0) then write(iulog,"('>>> mp_mageq: error allocating rcvbuf')") call shutdown('mp_mageq') endif endif sndbuf = 0._r8 rcvbuf = 0._r8 len = mxmaglon*nmlev*nf ! ! If mag equator is in current subdomain, load it into sndbuf ! and send to other tasks in my task column (mytidi) ! jlateq = (nmlat+1)/2 ! lat index of mag equator (49) have_eq = .false. do j=mlat0,mlat1 if (j == jlateq) then ! load send buffer w/ data at equator have_eq = .true. do i=mlon0,mlon1 sndbuf(i-mlon0+1,:,:) = fin(i,j,:,:) enddo ! ! Send mag equator data to other tasks in my task column (mytidi): do itask=0,ntask-1 if (itask /= mytid.and.tasks(itask)%mytidi==mytidi) then call mpi_isend(sndbuf,len,MPI_REAL8,itask,1, & MPI_COMM_WORLD,ireqsend,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mageq isend') call mpi_wait(ireqsend,irstat,ier) endif ! another task in mytidi enddo ! itask=0,ntask-1 endif ! j==jlateq enddo ! j=mlat0,mlat1 ! ! Receive by other tasks in the sending task's column: fout = 0._r8 if (.not.have_eq) then ! find task to receive from do itask=0,ntask-1 do j=tasks(itask)%mlat0,tasks(itask)%mlat1 if (j == jlateq.and.tasks(itask)%mytidi==mytidi) then call mpi_irecv(rcvbuf,len,MPI_REAL8,itask,1, & MPI_COMM_WORLD,ireqrecv,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mageq irecv') call mpi_wait(ireqrecv,irstat,ier) do n=1,nf do i=mlon0,mlon1 fout(i,:,n) = rcvbuf(i-mlon0+1,:,n) enddo enddo endif ! itask has mag eq and is in my column (sending task) enddo ! scan itask latitudes enddo ! task table search ! ! If I am the sending task, set fout to equator values of input array: else do n=1,nf do i=mlon0,mlon1 fout(i,:,n) = fin(i,jlateq,:,n) enddo enddo endif ! I am receiving or sending task if (mpi_timing) then endtime = mpi_wtime() time_mageq=time_mageq+(endtime-starttime) endif end subroutine mp_mageq !----------------------------------------------------------------------- subroutine mp_mageq_jpm1(f,mlon0,mlon1,mlat0,mlat1,nmlonp1,feq_jpm1,nf) ! ! All tasks need data at mag latitudes equator-1, equator+1 at global ! longitudes. ! On input: f is 6 fields on mag subdomains: zigm11,zigm22,zigmc,zigm2,rim1,rim2 ! On output: feq_jpm1(nmlonp1,2,nf) ! ! Args: integer,intent(in) :: mlon0,mlon1,mlat0,mlat1,nmlonp1,nf real(r8),intent(in) :: f(mlon0:mlon1,mlat0:mlat1,nf) real(r8),intent(out) :: feq_jpm1(nmlonp1,2,nf) ! eq-1,eq+1 ! ! Local: integer :: ifld,j,ier,len,jlateq,itask,i0,i1 real(r8),allocatable,save :: sndbuf(:,:,:,:), rcvbuf(:,:,:,:) real(r8) :: starttime,endtime if (mpi_timing) starttime = mpi_wtime() ! ! Allocate buffers once per run: if (.not.allocated(sndbuf)) then allocate(sndbuf(nmlonp1,2,nf,0:ntask-1),stat=ier) if (ier /= 0) call shutdown('mp_mageq_jpm1 error allocating sndbuf') endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(nmlonp1,2,nf,0:ntask-1),stat=ier) if (ier /= 0) call shutdown('mp_mageq_jpm1 error allocating rcvbuf') endif sndbuf = 0._r8 rcvbuf = 0._r8 len = nmlonp1*2*nf ! ! Load send buffer w/ eq +/- 1 for current subdomain ! (redundant to all tasks for alltoall) ! jlateq = (nmlat+1)/2 do itask=0,ntask-1 do j=mlat0,mlat1 if (j == jlateq+1) then ! equator+1 sndbuf(mlon0:mlon1,1,:,itask) = f(mlon0:mlon1,j,:) elseif (j == jlateq-1) then ! equator-1 sndbuf(mlon0:mlon1,2,:,itask) = f(mlon0:mlon1,j,:) endif ! j==jlateq enddo ! j=mlat0,mlat1 enddo ! itask ! ! Do the exchange: call mpi_alltoall(sndbuf,len,MPI_REAL8, & rcvbuf,len,MPI_REAL8,MPI_COMM_WORLD,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mageq_jpm1 call mpi_alltoall') ! ! Unpack receive buffer to output feq_jpm1(nmlonp1,2,nf) ! rcvbuf(nmlonp1,2,nf,0:ntask-1) do itask=0,ntask-1 i0 = tasks(itask)%mlon0 ; i1 = tasks(itask)%mlon1 do j=tasks(itask)%mlat0,tasks(itask)%mlat1 if (j == jlateq+1) then ! equator+1 feq_jpm1(i0:i1,1,:) = rcvbuf(i0:i1,1,:,itask) elseif (j == jlateq-1) then ! equator-1 feq_jpm1(i0:i1,2,:) = rcvbuf(i0:i1,2,:,itask) endif enddo ! j=mlat0,mlat1 enddo ! itask ! ! Periodic point: feq_jpm1(nmlonp1,:,:) = feq_jpm1(1,:,:) if (mpi_timing) then endtime = mpi_wtime() time_mageq_jpm1=time_mageq_jpm1+(endtime-starttime) endif end subroutine mp_mageq_jpm1 !----------------------------------------------------------------------- subroutine mp_mageq_jpm3(f,mlon0,mlon1,mlat0,mlat1,nmlonp1,feq_jpm3,nf) ! ! All tasks need global longitudes at mag latitudes equator, ! and equator +/- 1,2,3 ! On input: f is nf fields on mag subdomains ! On output: feq_jpm3(nmlonp1,-3:3,nf) has global lons at eq, eq +/- 1,2,3 ! 2nd dimension of feq_jpm3 (and send/recv buffers) is as follows: ! +3: eq+3 ! +2: eq+2 ! +1: eq+1 ! 0: eq ! -1: eq-1 ! -2: eq-2 ! -3: eq-3 ! ! Args: integer,intent(in) :: mlon0,mlon1,mlat0,mlat1,nmlonp1,nf real(r8),intent(in) :: f(mlon0:mlon1,mlat0:mlat1,nf) real(r8),intent(out) :: feq_jpm3(nmlonp1,-3:3,nf) ! ! Local: integer :: ifld,j,ier,len,jlateq,itask,i0,i1 real(r8),allocatable,save :: sndbuf(:,:,:,:), rcvbuf(:,:,:,:) real(r8) :: starttime,endtime integer,parameter :: mxnf=6 if (mpi_timing) starttime = mpi_wtime() if (nf > mxnf) then write(iulog,"('>>> mp_mageq_jpm3: nf=',i4,' but cannot be called with greater than mxnf=',i4)") & nf,mxnf call shutdown('mp_mageq_jpm3') endif ! ! Allocate buffers once per run (max number of fields): if (.not.allocated(sndbuf)) then allocate(sndbuf(nmlonp1,-3:3,mxnf,0:ntask-1),stat=ier) if (ier /= 0) & call shutdown('mp_mageq_jpm3 error allocating sndbuf') endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(nmlonp1,-3:3,mxnf,0:ntask-1),stat=ier) if (ier /= 0) & call shutdown('mp_mageq_jpm3 error allocating rcvbuf') endif sndbuf = 0._r8 rcvbuf = 0._r8 len = nmlonp1*7*nf ! ! Load send buffer w/ eq +/- 3 for current subdomain ! (redundant to all tasks for alltoall) ! jlateq = (nmlat+1)/2 do itask=0,ntask-1 do j=mlat0,mlat1 if (j == jlateq-3) then ! equator-3 sndbuf(mlon0:mlon1,-3,1:nf,itask) = f(mlon0:mlon1,j,:) elseif (j == jlateq-2) then ! equator-2 sndbuf(mlon0:mlon1,-2,1:nf,itask) = f(mlon0:mlon1,j,:) elseif (j == jlateq-1) then ! equator-1 sndbuf(mlon0:mlon1,-1,1:nf,itask) = f(mlon0:mlon1,j,:) elseif (j == jlateq) then ! equator sndbuf(mlon0:mlon1,0,1:nf,itask) = f(mlon0:mlon1,j,:) elseif (j == jlateq+1) then ! equator+1 sndbuf(mlon0:mlon1,1,1:nf,itask) = f(mlon0:mlon1,j,:) elseif (j == jlateq+2) then ! equator+2 sndbuf(mlon0:mlon1,2,1:nf,itask) = f(mlon0:mlon1,j,:) elseif (j == jlateq+3) then ! equator+3 sndbuf(mlon0:mlon1,3,1:nf,itask) = f(mlon0:mlon1,j,:) endif ! j==jlateq enddo ! j=mlat0,mlat1 enddo ! itask ! ! Do the exchange: call mpi_alltoall(sndbuf(:,:,1:nf,:),len,MPI_REAL8, & rcvbuf(:,:,1:nf,:),len,MPI_REAL8,MPI_COMM_WORLD,ier) if (ier /= 0) & call handle_mpi_err(ier,'mp_mageq_jpm3 call mpi_alltoall') ! ! Unpack receive buffer to output feq_jpm3(nmlonp1,-3:3,nf) ! rcvbuf(nmlonp1,-3:3,nf,0:ntask-1) do itask=0,ntask-1 i0 = tasks(itask)%mlon0 ; i1 = tasks(itask)%mlon1 do j=tasks(itask)%mlat0,tasks(itask)%mlat1 if (j == jlateq-3) then ! equator-3 feq_jpm3(i0:i1,-3,:) = rcvbuf(i0:i1,-3,1:nf,itask) elseif (j == jlateq-2) then ! equator-2 feq_jpm3(i0:i1,-2,:) = rcvbuf(i0:i1,-2,1:nf,itask) elseif (j == jlateq-1) then ! equator-1 feq_jpm3(i0:i1,-1,:) = rcvbuf(i0:i1,-1,1:nf,itask) elseif (j == jlateq) then ! equator feq_jpm3(i0:i1,0,:) = rcvbuf(i0:i1,0,1:nf,itask) elseif (j == jlateq+1) then ! equator+1 feq_jpm3(i0:i1,1,:) = rcvbuf(i0:i1,1,1:nf,itask) elseif (j == jlateq+2) then ! equator+2 feq_jpm3(i0:i1,2,:) = rcvbuf(i0:i1,2,1:nf,itask) elseif (j == jlateq+3) then ! equator+3 feq_jpm3(i0:i1,3,:) = rcvbuf(i0:i1,3,1:nf,itask) endif enddo ! j=mlat0,mlat1 enddo ! itask ! ! Periodic point: feq_jpm3(nmlonp1,:,:) = feq_jpm3(1,:,:) if (mpi_timing) then endtime = mpi_wtime() time_mageq_jpm3=time_mageq_jpm3+(endtime-starttime) endif end subroutine mp_mageq_jpm3 !----------------------------------------------------------------------- subroutine mp_magpole_2d(f,ilon0,ilon1,ilat0,ilat1, & nglblon,jspole,jnpole,fpole_jpm2,nf) ! ! Similiar to mp_geopole_2d, but for mag pole ! Return fpole_jpm2(nglblon,1->4,nf) as: ! 1: j = jspole+1 (spole+1) ! 2: j = jspole+2 (spole+2) ! 3: j = jnpole-1 (npole-1) ! 4: j = jnpole-2 (npole-2) ! This can be called with different number of fields nf, but cannot ! be called w/ > mxnf fields. Buffers are allocated once per run for ! mxnf fields, then passed and processed for nf fields. ! ! Args: integer,intent(in) :: ilon0,ilon1,ilat0,ilat1,nglblon,jspole,jnpole,nf real(r8),intent(in) :: f(ilon0:ilon1,ilat0:ilat1,nf) real(r8),intent(out) :: fpole_jpm2(nglblon,4,nf) ! ! Local: integer :: j,k,ier,len,itask,i0,i1,j0,j1 real(r8),allocatable,save :: sndbuf(:,:,:,:),rcvbuf(:,:,:,:) real(r8) :: starttime,endtime integer,parameter :: mxnf=6 if (mpi_timing) starttime = mpi_wtime() if (nf > mxnf) then write(iulog,"('>>> mp_magpole_2d: nf=',i4,' but cannot be called with greater than mxnf=',i4)") & nf,mxnf call shutdown('mp_magpole_2d') endif ! ! Allocate the max number of fields, only once per run: if (.not.allocated(sndbuf)) then allocate(sndbuf(nglblon,4,mxnf,0:ntask-1),stat=ier) if (ier /= 0) call shutdown('mp_magpole_2d: error allocating sndbuf') endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(nglblon,4,mxnf,0:ntask-1),stat=ier) if (ier /= 0) call shutdown('mp_magpole_2d: error allocating rcvbuf') endif sndbuf = 0._r8 rcvbuf = 0._r8 len = nglblon*4*nf ! ! Load send buffer with values at poles +/- 2 for current subdomain ! (redundant to all tasks for alltoall) ! do itask=0,ntask-1 do j=ilat0,ilat1 if (j==jspole+1) then ! south pole +1 sndbuf(ilon0:ilon1,1,1:nf,itask) = f(ilon0:ilon1,j,:) elseif (j==jspole+2) then ! south pole +2 sndbuf(ilon0:ilon1,2,1:nf,itask) = f(ilon0:ilon1,j,:) elseif (j==jnpole-1) then ! north pole -1 sndbuf(ilon0:ilon1,3,1:nf,itask) = f(ilon0:ilon1,j,:) elseif (j==jnpole-2) then ! north pole -2 sndbuf(ilon0:ilon1,4,1:nf,itask) = f(ilon0:ilon1,j,:) endif enddo enddo ! ! Do the exchange: call mpi_alltoall(sndbuf(:,:,1:nf,:),len,MPI_REAL8, & rcvbuf(:,:,1:nf,:),len,MPI_REAL8,MPI_COMM_WORLD,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_magpole_2d call mpi_alltoall') ! ! Unpack rcvbuf into output fpole_jpm2(nglblon,4,nf) do itask=0,ntask-1 i0 = tasks(itask)%mlon0 i1 = tasks(itask)%mlon1 j0 = tasks(itask)%mlat0 j1 = tasks(itask)%mlat1 do j=j0,j1 if (j==jspole+1) then ! south pole +1 fpole_jpm2(i0:i1,1,:) = rcvbuf(i0:i1,1,1:nf,itask) elseif (j==jspole+2) then ! south pole +2 fpole_jpm2(i0:i1,2,:) = rcvbuf(i0:i1,2,1:nf,itask) elseif (j==jnpole-1) then ! north pole -1 fpole_jpm2(i0:i1,3,:) = rcvbuf(i0:i1,3,1:nf,itask) elseif (j==jnpole-2) then ! north pole -2 fpole_jpm2(i0:i1,4,:) = rcvbuf(i0:i1,4,1:nf,itask) endif enddo enddo if (mpi_timing) then endtime = mpi_wtime() time_magpole_2d=time_magpole_2d+(endtime-starttime) endif end subroutine mp_magpole_2d !----------------------------------------------------------------------- subroutine mp_magpole_3d(f,ilon0,ilon1,ilat0,ilat1,nlev,& nglblon,jspole,jnpole,fpole_jpm2,nf) ! ! Similiar to mp_geopole_2d, but for mag pole ! Return fpole_jpm2(nglblon,1->4,nlev,nf) as: ! 1: j = jspole+1 (spole+1) ! 2: j = jspole+2 (spole+2) ! 3: j = jnpole-1 (npole-1) ! 4: j = jnpole-2 (npole-2) ! This can be called with different number of fields nf, but cannot ! be called w/ > mxnf fields. Buffers are allocated once per run for ! mxnf fields, then passed and processed for nf fields. ! ! Args: integer,intent(in) :: ilon0,ilon1,ilat0,ilat1,nglblon,& jspole,jnpole,nf,nlev real(r8),intent(in) :: f(ilon0:ilon1,ilat0:ilat1,nlev,nf) real(r8),intent(out) :: fpole_jpm2(nglblon,4,nlev,nf) ! ! Local: integer :: j,k,ier,len,itask,i0,i1,j0,j1 real(r8),allocatable,save :: sndbuf(:,:,:,:,:),rcvbuf(:,:,:,:,:) real(r8) :: starttime,endtime integer,parameter :: mxnf=6 if (mpi_timing) starttime = mpi_wtime() if (nf > mxnf) then write(iulog,"('>>> mp_magpole_3d: nf=',i4,' but cannot be called with greater than mxnf=',i4)") & nf,mxnf call shutdown('mp_magpole_3d') endif ! ! Allocate the max number of fields, only once per run: if (.not.allocated(sndbuf)) then allocate(sndbuf(nglblon,4,nlev,mxnf,0:ntask-1),stat=ier) if (ier /= 0) call shutdown('mp_magpole_3d: error allocating sndbuf') endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(nglblon,4,nlev,mxnf,0:ntask-1),stat=ier) if (ier /= 0) call shutdown('mp_magpole_3d: error allocating rcvbuf') endif sndbuf = 0._r8 rcvbuf = 0._r8 len = nglblon*4*nlev*nf ! ! Load send buffer with values at poles +/- 2 for current subdomain ! (redundant to all tasks for alltoall) ! do itask=0,ntask-1 do j=ilat0,ilat1 do k=1,nlev if (j==jspole+1) then ! south pole +1 sndbuf(ilon0:ilon1,1,k,1:nf,itask) = f(ilon0:ilon1,j,k,:) elseif (j==jspole+2) then ! south pole +2 sndbuf(ilon0:ilon1,2,k,1:nf,itask) = f(ilon0:ilon1,j,k,:) elseif (j==jnpole-1) then ! north pole -1 sndbuf(ilon0:ilon1,3,k,1:nf,itask) = f(ilon0:ilon1,j,k,:) elseif (j==jnpole-2) then ! north pole -2 sndbuf(ilon0:ilon1,4,k,1:nf,itask) = f(ilon0:ilon1,j,k,:) endif enddo enddo enddo ! ! Do the exchange: call mpi_alltoall(sndbuf(:,:,:,1:nf,:),len,MPI_REAL8,& rcvbuf(:,:,:,1:nf,:),len,MPI_REAL8,MPI_COMM_WORLD,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_magpole_3d call mpi_alltoall') ! ! Unpack rcvbuf into output fpole_jpm2(nglblon,4,nf) do itask=0,ntask-1 i0 = tasks(itask)%mlon0 i1 = tasks(itask)%mlon1 j0 = tasks(itask)%mlat0 j1 = tasks(itask)%mlat1 do j=j0,j1 if (j==jspole+1) then ! south pole +1 fpole_jpm2(i0:i1,1,:,:) = rcvbuf(i0:i1,1,:,1:nf,itask) elseif (j==jspole+2) then ! south pole +2 fpole_jpm2(i0:i1,2,:,:) = rcvbuf(i0:i1,2,:,1:nf,itask) elseif (j==jnpole-1) then ! north pole -1 fpole_jpm2(i0:i1,3,:,:) = rcvbuf(i0:i1,3,:,1:nf,itask) elseif (j==jnpole-2) then ! north pole -2 fpole_jpm2(i0:i1,4,:,:) = rcvbuf(i0:i1,4,:,1:nf,itask) endif enddo enddo if (mpi_timing) then endtime = mpi_wtime() time_magpole_3d=time_magpole_3d+(endtime-starttime) endif end subroutine mp_magpole_3d !----------------------------------------------------------------------- subroutine mp_magpoles(f,ilon0,ilon1,ilat0,ilat1,nglblon,& jspole,jnpole,fpoles,nf) ! ! Similiar to mp_magpole_2d, but returns global longitudes for ! j==1 and j==nmlat (not for poles +/- 2) ! Return fpoles(nglblon,2,nf) as: ! 1: j = jspole (spole) ! 2: j = jnpole (npole) ! This can be called with different number of fields nf, but cannot ! be called w/ > mxnf fields. Buffers are allocated once per run for ! mxnf fields, then passed and processed for nf fields. ! ! Args: integer,intent(in) :: ilon0,ilon1,ilat0,ilat1,nglblon, & jspole,jnpole,nf real(r8),intent(in) :: f(ilon0:ilon1,ilat0:ilat1,nf) real(r8),intent(out) :: fpoles(nglblon,2,nf) ! ! Local: integer :: j,k,ier,len,itask,i0,i1,j0,j1 real(r8),allocatable,save :: sndbuf(:,:,:,:),rcvbuf(:,:,:,:) real(r8) :: starttime,endtime integer,save :: mxnf if (mpi_timing) starttime = mpi_wtime() mxnf = nmlev if (nf > mxnf) then write(iulog,"('>>> mp_magpoles: nf=',i4,' but cannot be called with greater than mxnf=',i4)") & nf,mxnf call shutdown('mp_magpoles') endif ! ! Allocate the max number of fields, only once per run: if (.not.allocated(sndbuf)) then allocate(sndbuf(nglblon,2,mxnf,0:ntask-1),stat=ier) if (ier /= 0) & call shutdown('mp_magpoles: error allocating sndbuf') endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(nglblon,2,mxnf,0:ntask-1),stat=ier) if (ier /= 0) & call shutdown('mp_magpoles: error allocating rcvbuf') endif sndbuf = 0._r8 rcvbuf = 0._r8 len = nglblon*2*nf ! ! Load send buffer with values at poles +/- 2 for current subdomain ! (redundant to all tasks for alltoall) ! do itask=0,ntask-1 do j=ilat0,ilat1 if (j==jspole) then ! south pole sndbuf(ilon0:ilon1,1,1:nf,itask) = f(ilon0:ilon1,j,:) elseif (j==jnpole) then ! npole pole sndbuf(ilon0:ilon1,2,1:nf,itask) = f(ilon0:ilon1,j,:) endif enddo enddo ! ! Do the exchange: call mpi_alltoall(sndbuf(:,:,1:nf,:),len,MPI_REAL8, & rcvbuf(:,:,1:nf,:),len,MPI_REAL8,MPI_COMM_WORLD,ier) if (ier /= 0) & call handle_mpi_err(ier,'mp_magpoles call mpi_alltoall') ! ! Unpack rcvbuf into output fpoles(nglblon,2,nf) do itask=0,ntask-1 i0 = tasks(itask)%mlon0 i1 = tasks(itask)%mlon1 j0 = tasks(itask)%mlat0 j1 = tasks(itask)%mlat1 do j=j0,j1 if (j==jspole) then ! south pole fpoles(i0:i1,1,:) = rcvbuf(i0:i1,1,1:nf,itask) elseif (j==jnpole) then ! north pole fpoles(i0:i1,2,:) = rcvbuf(i0:i1,2,1:nf,itask) endif enddo enddo if (mpi_timing) then endtime = mpi_wtime() time_magpoles=time_magpoles+(endtime-starttime) endif end subroutine mp_magpoles !----------------------------------------------------------------------- subroutine conjugate_points use maggrid,only: gmlat ! ! Local: integer :: ier,i,j,js,jn,itask,ii,jj real :: starttime,endtime if (mpi_timing) starttime = mpi_wtime() ! ! nsend_south(ntask): number of lats in south to send north ! nrecv_north(ntask): number of lats in north to recv from south ! allocate(nsend_south(0:ntask-1),stat=ier) allocate(nrecv_north(0:ntask-1),stat=ier) ! ! send_south_coords: south j lats to send north ! recv_north_coords: north j lats to recv from south ! allocate(send_south_coords(mxmaglat,0:ntask-1),stat=ier) allocate(recv_north_coords(mxmaglat,0:ntask-1),stat=ier) nsend_south(:) = 0 nrecv_north(:) = 0 send_south_coords(:,:) = 0 recv_north_coords(:,:) = 0 magloop: do j=mlat0,mlat1 ! ! In north hem: find tasks w/ conjugate points in south to recv: ! (nmlath is in params module) if (gmlat(j) > 0._r8) then ! in north hem of current task js = nmlath-(j-nmlath) ! j index to south conjugate point (should be -j) do itask=0,ntask-1 do jj = tasks(itask)%mlat0,tasks(itask)%mlat1 ! ! Receive these north coords from the south: if (jj==js.and.mlon0==tasks(itask)%mlon0.and. & mlon1==tasks(itask)%mlon1) then nrecv_north(itask) = nrecv_north(itask)+1 recv_north_coords(nrecv_north(itask),itask) = j endif enddo ! jj of remote task enddo ! itask=0,ntask-1 if (all(nrecv_north==0)) & write(iulog,"(2a,i4,a,f8.2)") '>>> WARNING: could not find north conjugate',& ' points corresponding to south latitude js=',js,' gmlat(js)=',gmlat(js) ! ! In south hem: find tasks w/ conjugate points in north to send: elseif (gmlat(j) < 0._r8.and.j /= nmlath) then ! in south hem jn = nmlath+(nmlath-j) ! j index of north conjugate point do itask=0,ntask-1 do jj = tasks(itask)%mlat0,tasks(itask)%mlat1 if (jj==jn.and.mlon0==tasks(itask)%mlon0.and. & mlon1==tasks(itask)%mlon1) then nsend_south(itask) = nsend_south(itask)+1 ! Send these south coords to the north: send_south_coords(nsend_south(itask),itask) = j endif enddo ! jj of remote task enddo ! itask=0,ntask-1 if (all(nsend_south==0)) & write(iulog,"(2a,i4,a,f8.2)") '>>> WARNING: could not find south conjugate',& ' points corresponding to north latitude jn=',jn,' gmlat(jn)=',gmlat(jn) endif ! in north or south hem enddo magloop ! j=mlat0,mlat1 ! do itask=0,ntask-1 ! if (nrecv_north(itask) > 0) then ! write(iulog,"('I will receive ',i4,' north conjugate lats from task ',i4,':',/,'recv_north_coords =',/,(8i5))") & ! nrecv_north(itask),itask,recv_north_coords(1:nrecv_north(itask),itask) ! endif ! if (nsend_south(itask) > 0) then ! write(iulog,"('I will send ',i4,' south conjugate lats to task ',i4,':',/,'send_south_coords =',/,(8i5))") & ! nsend_south(itask),itask,send_south_coords(1:nsend_south(itask),itask) ! endif ! enddo if (mpi_timing) then endtime = mpi_wtime() time_conjugate_points=time_conjugate_points+(endtime-starttime) endif end subroutine conjugate_points !----------------------------------------------------------------------- subroutine mp_mag_foldhem(f,mlon0,mlon1,mlat0,mlat1,nf) ! ! For each point in northern hemisphere (if any) of the current task ! subdomain, receive data from conjugate point in the south (from the ! south task that owns it), and sum it to the north point data. ! Do this for nf fields. Conjugate point indices to send/recv to/from ! each task were determined by sub conjugate_points (this module). ! nsend_south, ! number of south lats to send to north (each task) ! nrecv_north ! number of north lats to send to south (each task) ! ! This routine is called from edynamo at every timestep. ! Sub conjugate_points is called once per run, from mp_distribute. ! ! Args: integer,intent(in) :: mlon0,mlon1,mlat0,mlat1,nf real(r8),intent(inout) :: f(mlon0:mlon1,mlat0:mlat1,nf) ! ! Local: integer :: i,j,n,len,itask,ifld,ier,nmlons real(r8) :: sndbuf(mxmaglon,mxmaglat,nf,0:ntask-1) real(r8) :: rcvbuf(mxmaglon,mxmaglat,nf,0:ntask-1) integer :: jsend(0:ntask-1),jrecv(0:ntask-1) real(r8) :: starttime,endtime integer :: irstat(MPI_STATUS_SIZE) ! mpi receive status if (mpi_timing) starttime = mpi_wtime() ! sndbuf = 0._r8 ; rcvbuf = 0._r8 jsend = 0 ; jrecv = 0 len = mxmaglon*mxmaglat*nf nmlons = mlon1-mlon0+1 ! ! Send south data to north itask: ! (To avoid deadlock, do not send if north task is also myself. This will ! happen when there is an odd number of tasks in the latitude dimension, ! e.g., ntask == 12, 30, etc) ! do itask=0,ntask-1 ! Attempt to fetch from allocatable variable NSEND_SOUTH when it is not allocated if (nsend_south(itask) > 0 .and. itask /= mytid) then do ifld = 1,nf do n=1,nsend_south(itask) sndbuf(1:nmlons,n,ifld,itask) = & f(:,send_south_coords(n,itask),ifld) enddo enddo ! ifld=1,nf call mpi_isend(sndbuf(1,1,1,itask),len,MPI_REAL8, & itask,1,MPI_COMM_WORLD,jsend(itask),ier) call mpi_wait(jsend(itask),irstat,ier) endif ! nsend_south(itask) > 0 enddo ! itask=0,ntask-1 ! ! Receive north data from south itask and add to north, ! i.e., north = north+south. (do not receive if south task is ! also myself, but do add south data to my north points, see below) ! do itask=0,ntask-1 if (nrecv_north(itask) > 0 .and. itask /= mytid) then call mpi_irecv(rcvbuf(1,1,1,itask),len,MPI_REAL8, & itask,1,MPI_COMM_WORLD,jrecv(itask),ier) call mpi_wait(jrecv(itask),irstat,ier) do ifld=1,nf do n=1,nrecv_north(itask) ! ! Receive lats in reverse order: f(mlon0:mlon1, & recv_north_coords(nrecv_north(itask)-n+1,itask),ifld) = & f(mlon0:mlon1, & recv_north_coords(nrecv_north(itask)-n+1,itask),ifld) + & rcvbuf(1:nmlons,n,ifld,itask) enddo ! n=1,nrecv_north(itask) enddo ! ifld=1,nf ! ! If I am send *and* receive task, simply add my south data to my north points: elseif (nrecv_north(itask) > 0 .and. itask == mytid) then do ifld=1,nf do n=1,nrecv_north(itask) f(mlon0:mlon1, & recv_north_coords(nrecv_north(itask)-n+1,itask),ifld) = & f(mlon0:mlon1, & recv_north_coords(nrecv_north(itask)-n+1,itask),ifld) + & f(mlon0:mlon1,send_south_coords(n,itask),ifld) enddo ! n=1,nrecv_north(itask) enddo ! ifld=1,nf endif ! nrecv_north(itask) > 0 enddo ! itask=0,ntask-1 ! ! Mag equator is also "folded", but not included in conjugate points, ! so double it here: do j=mlat0,mlat1 if (j==nmlath) then do ifld=1,nf f(:,j,ifld) = f(:,j,ifld)+f(:,j,ifld) enddo endif enddo if (mpi_timing) then endtime = mpi_wtime() time_mag_foldhem=time_mag_foldhem+(endtime-starttime) endif end subroutine mp_mag_foldhem !----------------------------------------------------------------------- subroutine mp_mag_periodic_f2d(f,mlon0,mlon1,mlat0,mlat1,nf) ! ! Args: integer,intent(in) :: mlon0,mlon1,mlat0,mlat1,nf real(r8),intent(inout) :: f(mlon0:mlon1,mlat0:mlat1,nf) ! ! Local: integer :: j,ier,idest,isrc,len,ireqsend,ireqrecv,ifld,msgtag real(r8) :: sndbuf(mxmaglat,nf),rcvbuf(mxmaglat,nf) real(r8) :: starttime,endtime integer :: irstat(MPI_STATUS_SIZE) ! mpi receive status if (mpi_timing) starttime = mpi_wtime() len = mxmaglat*nf ! ! I am a western-most task. Send lon 1 to eastern-most tasks: if (mytidi==0) then idest = itask_table_mag(ntaski-1,mytidj) do j=mlat0,mlat1 sndbuf(j-mlat0+1,:) = f(1,j,:) enddo msgtag = mytid call mpi_isend(sndbuf,len,MPI_REAL8,idest,msgtag,MPI_COMM_WORLD,& ireqsend,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_periodic_f2d send to idest') call mpi_wait(ireqsend,irstat,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_periodic_f2d wait for send') ! ! I am eastern-most task. Receive lon 1 from western-most tasks, ! and assign to nmlonp1: elseif (mytidi==ntaski-1) then isrc = itask_table_mag(0,mytidj) msgtag = isrc call mpi_irecv(rcvbuf,len,MPI_REAL8,isrc,msgtag,MPI_COMM_WORLD,& ireqrecv,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_periodic_f2d recv from isrc') call mpi_wait(ireqrecv,irstat,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_periodic_f2d wait for recv') do j=mlat0,mlat1 f(nmlonp1,j,:) = rcvbuf(j-mlat0+1,:) enddo endif ! mytidi == 0 or ntaski-1 if (mpi_timing) then endtime = mpi_wtime() time_mag_periodic_f2d=time_mag_periodic_f2d+(endtime-starttime) endif end subroutine mp_mag_periodic_f2d !----------------------------------------------------------------------- subroutine mp_mag_periodic_f3d(f,mlon0,mlon1,mlat0,mlat1,nlev,nf) ! ! Args: integer,intent(in) :: mlon0,mlon1,mlat0,mlat1,nlev,nf real(r8),intent(inout) :: f(mlon0:mlon1,mlat0:mlat1,nlev,nf) ! ! Local: integer :: j,ier,idest,isrc,len,ireqsend,ireqrecv,ifld,msgtag real(r8) :: sndbuf(mxmaglat,nlev,nf),rcvbuf(mxmaglat,nlev,nf) real(r8) :: starttime,endtime integer :: irstat(MPI_STATUS_SIZE) ! mpi receive status if (mpi_timing) starttime = mpi_wtime() len = mxmaglat*nlev*nf ! ! I am a western-most task. Send lon 1 to eastern-most tasks: if (mytidi==0) then idest = itask_table_mag(ntaski-1,mytidj) do j=mlat0,mlat1 sndbuf(j-mlat0+1,:,:) = f(1,j,:,:) enddo msgtag = mytid call mpi_isend(sndbuf,len,MPI_REAL8,idest,msgtag,MPI_COMM_WORLD,& ireqsend,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_periodic_f3d send to idest') call mpi_wait(ireqsend,irstat,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_periodic_f3d wait for send') ! ! I am eastern-most task. Receive lon 1 from western-most tasks, ! and assign to nmlonp1: elseif (mytidi==ntaski-1) then isrc = itask_table_mag(0,mytidj) msgtag = isrc call mpi_irecv(rcvbuf,len,MPI_REAL8,isrc,msgtag,MPI_COMM_WORLD,& ireqrecv,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_periodic_f2d recv from isrc') call mpi_wait(ireqrecv,irstat,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_periodic_f2d wait for recv') do j=mlat0,mlat1 f(nmlonp1,j,:,:) = rcvbuf(j-mlat0+1,:,:) enddo endif ! mytidi == 0 or ntaski-1 if (mpi_timing) then endtime = mpi_wtime() time_mag_periodic_f3d=time_mag_periodic_f3d+(endtime-starttime) endif end subroutine mp_mag_periodic_f3d !----------------------------------------------------------------------- subroutine mp_gather2root(fsub,ilon0,ilon1,ilat0,ilat1,fglb,nglblon,nglblat,nlev,geo) ! ! Args: integer,intent(in) :: ilon0,ilon1,ilat0,ilat1,nglblon,nglblat,nlev real(r8),intent(inout) :: fsub(ilon0:ilon1,ilat0:ilat1,nlev) real(r8),intent(inout) :: fglb(nglblon,nglblat,nlev) logical,intent(in) :: geo ! ! Local: integer :: len,i,j,ier,i0,i1,j0,j1,itask real(r8),dimension(nglblon,nglblat,nlev) :: sndbuf real(r8),dimension(nglblon,nglblat,nlev,0:ntask-1) :: rcvbuf sndbuf = 0._r8 rcvbuf = 0._r8 do j=ilat0,ilat1 sndbuf(ilon0:ilon1,j,:) = fsub(ilon0:ilon1,j,:) enddo ! write(iulog,"('mp_gather2root: min,max sndbuf=',2e12.4)") & ! minval(sndbuf(ilon0:ilon1,ilat0:ilat1,1:nlev)), & ! maxval(sndbuf(ilon0:ilon1,ilat0:ilat1,1:nlev)) ! ! Gather to root task: len = nglblon*nglblat*nlev call mpi_gather(sndbuf,len,MPI_REAL8,rcvbuf,len,MPI_REAL8,0, & MPI_COMM_WORLD,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_gather2root') ! ! Root task unpack from receive buffer: if (mytid==0) then do itask=0,ntask-1 if (geo) then i0 = tasks(itask)%lon0 ; i1 = tasks(itask)%lon1 j0 = tasks(itask)%lat0 ; j1 = tasks(itask)%lat1 else i0 = tasks(itask)%mlon0 ; i1 = tasks(itask)%mlon1 j0 = tasks(itask)%mlat0 ; j1 = tasks(itask)%mlat1 endif do j=j0,j1 fglb(i0:i1,j,:) = rcvbuf(i0:i1,j,:,itask) enddo enddo endif end subroutine mp_gather2root !----------------------------------------------------------------------- subroutine mp_mag_halos(fmsub,mlon0,mlon1,mlat0,mlat1,nf) ! ! Exchange halo/ghost points between magnetic grid subdomains for nf fields. ! Only a single halo point is required in both lon and lat dimensions. ! Note that all tasks in any row of the task matrix have the same ! mlat0,mlat1, and that all tasks in any column of the task matrix ! have the same mlon0,mlon1. ! Longitude halos are done first, exchanging mlat0:mlat1, then latitude ! halos are done, exchanging mlon0-1:mlon1+1 (i.e., including the ! longitude halos that were defined first). ! ! Args: integer,intent(in) :: mlon0,mlon1,mlat0,mlat1,nf real(r8),intent(inout) :: fmsub(mlon0-1:mlon1+1,mlat0-1:mlat1+1,nf) ! ! Local: integer :: i,j,ifld,west,east,north,south,len,isend0,isend1, & irecv0,irecv1,ier,nmlats,istat(MPI_STATUS_SIZE,4),ireq(4),nmlons real(r8),dimension(mlat1-mlat0+1,nf)::sndlon0,sndlon1,rcvlon0,rcvlon1 real(r8),dimension((mlon1+1)-(mlon0-1)+1,nf) :: & sndlat0,sndlat1,rcvlat0,rcvlat1 real(r8) :: starttime,endtime if (mpi_timing) starttime = mpi_wtime() ! ! Init send/recv buffers for lon halos: sndlon0 = 0._r8 ; rcvlon0 = 0._r8 sndlon1 = 0._r8 ; rcvlon1 = 0._r8 ! ! Identify east and west neightbors: west = itask_table_mag(mytidi-1,mytidj) east = itask_table_mag(mytidi+1,mytidj) ! ! Exchange mlat0:mlat1 (lat halos are not yet defined): nmlats = mlat1-mlat0+1 len = nmlats*nf ! ! Send mlon0 to the west neighbor, and mlon1 to the east. ! However, tasks are periodic in longitude (see itask_table_mag), ! and far west tasks send mlon0+1, and far east tasks send mlon1-1 ! do ifld=1,nf ! Far west tasks send mlon0+1 to far east (periodic) tasks: if (mytidi==0) then sndlon0(:,ifld) = fmsub(mlon0+1,mlat0:mlat1,ifld) ! Interior tasks send mlon0 to west neighbor: else sndlon0(:,ifld) = fmsub(mlon0,mlat0:mlat1,ifld) endif ! Far east tasks send mlon1-1 to far west (periodic) tasks: if (mytidi==nmagtaski-1) then sndlon1(:,ifld) = fmsub(mlon1-1,mlat0:mlat1,ifld) ! Interior tasks send mlon1 to east neighbor: else sndlon1(:,ifld) = fmsub(mlon1,mlat0:mlat1,ifld) endif enddo ! ifld=1,nf ! ! Send mlon0 to the west: call mpi_isend(sndlon0,len,MPI_REAL8,west,1,MPI_COMM_WORLD,isend0,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_halos send mlon0 to west') ! ! Send mlon1 to the east: call mpi_isend(sndlon1,len,MPI_REAL8,east,1,MPI_COMM_WORLD,isend1,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_halos send mlon1 to east') ! ! Recv mlon0-1 from west: call mpi_irecv(rcvlon0,len,MPI_REAL8,west,1,MPI_COMM_WORLD,irecv0,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_halos recv mlon0 from west') ! ! Recv mlon1+1 from east: call mpi_irecv(rcvlon1,len,MPI_REAL8,east,1,MPI_COMM_WORLD,irecv1,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_halos recv mlon1 from east') ! ! Wait for completions: ireq = (/isend0,isend1,irecv0,irecv1/) istat = 0 call mpi_waitall(4,ireq,istat,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_halos waitall for lons') ! ! Copy mlon0-1 from rcvlon0, and mlon1+1 from rcvlon1: do ifld=1,nf fmsub(mlon0-1,mlat0:mlat1,ifld) = rcvlon0(:,ifld) fmsub(mlon1+1,mlat0:mlat1,ifld) = rcvlon1(:,ifld) ! ! Fix special case of 2 tasks in longitude dimension: if (east == west) then fmsub(mlon0-1,mlat0:mlat1,ifld) = rcvlon1(:,ifld) fmsub(mlon1+1,mlat0:mlat1,ifld) = rcvlon0(:,ifld) endif enddo ! ifld=1,nf ! ! Now exchange latitudes: sndlat0 = 0._r8 ; rcvlat0 = 0._r8 sndlat1 = 0._r8 ; rcvlat1 = 0._r8 south = itask_table_mag(mytidi,mytidj-1) ! neighbor to south north = itask_table_mag(mytidi,mytidj+1) ! neighbor to north ! ! Include halo longitudes that were defined by the exchanges above: nmlons = (mlon1+1)-(mlon0-1)+1 len = nmlons*nf ! ! Send mlat0 to south neighbor, and mlat1 to north: do ifld=1,nf sndlat0(:,ifld) = fmsub(:,mlat0,ifld) sndlat1(:,ifld) = fmsub(:,mlat1,ifld) enddo ! ! Send mlat0 to south: call mpi_isend(sndlat0,len,MPI_REAL8,south,1,MPI_COMM_WORLD,isend0,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_halos send mlat0 to south') ! ! Send mlat1 to north: call mpi_isend(sndlat1,len,MPI_REAL8,north,1,MPI_COMM_WORLD,isend1,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_halos send mlat1 to north') ! ! Recv mlat0-1 from south: call mpi_irecv(rcvlat0,len,MPI_REAL8,south,1,MPI_COMM_WORLD,irecv0,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_halos recv mlat0-1 from south') ! ! Recv mlat1+1 from north: call mpi_irecv(rcvlat1,len,MPI_REAL8,north,1,MPI_COMM_WORLD,irecv1,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_halos recv mlat1+1 from north') ! ! Wait for completions: ireq = (/isend0,isend1,irecv0,irecv1/) istat = 0 call mpi_waitall(4,ireq,istat,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_mag_halos waitall for lats') ! ! Copy mlat0-1 from rcvlat0, and mlat1+1 from rcvlat1: do ifld=1,nf fmsub(:,mlat0-1,ifld) = rcvlat0(:,ifld) fmsub(:,mlat1+1,ifld) = rcvlat1(:,ifld) enddo ! ifld=1,nf if (mpi_timing) then endtime = mpi_wtime() time_mag_halos=time_mag_halos+(endtime-starttime) endif end subroutine mp_mag_halos !----------------------------------------------------------------------- subroutine mp_gather_edyn(fmsub,mlon0,mlon1,mlat0,mlat1,fmglb,nmlonp1,nmlat,nf) ! ! Gather fields on mag subdomains to root task, so root task can ! complete non-parallel portion of dynamo (starting after rhspde) ! ! Args: integer,intent(in) :: mlon0,mlon1,mlat0,mlat1,nmlonp1,nmlat,nf real(r8),intent(in) :: fmsub(mlon0:mlon1,mlat0:mlat1,nf) real(r8),intent(out) :: fmglb(nmlonp1,nmlat,nf) ! ! Local: integer :: len,i,j,ifld,ier,itask,i0,i1,j0,j1 real(r8),dimension(nmlonp1,nmlat,nf) :: sndbuf real(r8),dimension(nmlonp1,nmlat,nf,0:ntask-1) :: rcvbuf real(r8) :: starttime,endtime if (mpi_timing) starttime = mpi_wtime() sndbuf = 0._r8 rcvbuf = 0._r8 len = nmlonp1*nmlat*nf ! ! Load send buffer with my subdomain: do ifld=1,nf do j=mlat0,mlat1 sndbuf(mlon0:mlon1,j,ifld) = fmsub(:,j,ifld) enddo enddo ! ! Gather to root: call mpi_gather(sndbuf,len,MPI_REAL8,rcvbuf,len,MPI_REAL8,0, & MPI_COMM_WORLD,ier) if (ier /= 0) call handle_mpi_err(ier,'mp_gather_edyn: mpi_gather to root') ! ! Root task unpack from receive buffer: if (mytid==0) then do itask=0,ntask-1 do ifld=1,nf i0 = tasks(itask)%mlon0 ; i1 = tasks(itask)%mlon1 j0 = tasks(itask)%mlat0 ; j1 = tasks(itask)%mlat1 do j=j0,j1 fmglb(i0:i1,j,ifld) = rcvbuf(i0:i1,j,ifld,itask) enddo ! j=mlat0,mlat1 enddo ! ifld=1,nf enddo ! itask=0,ntask-1 endif ! root task if (mpi_timing) then endtime = mpi_wtime() time_gather_edyn=time_gather_edyn+(endtime-starttime) endif end subroutine mp_gather_edyn !----------------------------------------------------------------------- subroutine mp_scatter_phim(phim_glb,phim) real(r8),intent(in) :: phim_glb(nmlonp1,nmlat) real(r8),intent(out) :: phim(mlon0:mlon1,mlat0:mlat1) ! ! Local: integer :: ier,len,i,j real(r8) :: starttime,endtime if (mpi_timing) starttime = mpi_wtime() ! ! Broadcast global phim (from pdynamo phim(nmlonp1,nmlat)): len = nmlat*nmlonp1 call mpi_bcast(phim_glb,len,MPI_REAL8,0,MPI_COMM_WORLD,ier) if (ier /= 0) & call handle_mpi_err(ier,'mp_scatter_phim: bcast global phim') ! ! Define subdomains: do j=mlat0,mlat1 do i=mlon0,mlon1 phim(i,j) = phim_glb(i,j) enddo enddo if (mpi_timing) then endtime = mpi_wtime() time_scatter_phim=time_scatter_phim+(endtime-starttime) endif end subroutine mp_scatter_phim !----------------------------------------------------------------------- subroutine mp_mag_jslot(fin,mlon00,mlon11,mlat00,mlat11, & fout,jneed,mxneed,nf) ! ! Current task needs to receive (from other tasks) field f at (non-zero) ! latitude indices in jneed, at all longitudes in the current subdomain. ! Note subdomains include halo points mlon0-1 and mlat1+1. Data in f also ! includes halo points (will need the lat data at halo-longitudes) ! ! Args: integer,intent(in) :: mlon00,mlon11,mlat00,mlat11 ! subdomains w/ halos integer,intent(in) :: nf ! number of fields integer,intent(in) :: mxneed ! max number of needed lats (nmlat+2) integer,intent(in) :: jneed(mxneed) ! j-indices of needed lats (where /= -1) real(r8),intent(in) :: fin(mlon00:mlon11,mlat00:mlat11,nf) ! data at current subdomain real(r8),intent(out) :: fout(mlon00:mlon11,mxneed,nf) ! returned data at needed lats ! ! Local: integer :: ier,njneed,i,j,n,nj,idest,ireq(2), & isend(ntask),irecv(ntask),icount, & len,istat,nlons,isrc,msgid,ifld integer :: tij ! rank in cols_comm (0 to nmagtaskj-1) integer :: jhave(mxneed),njhave,wid integer :: peersneed(mxneed,0:nmagtaskj-1) integer :: jneedall (mxneed,0:nmagtaskj-1) real(r8),allocatable :: sndbuf(:,:,:) real(r8),allocatable :: rcvbuf(:,:,:) integer :: irstat(MPI_STATUS_SIZE) ! mpi receive status real(r8) :: starttime,endtime ! ! External: integer,external :: ixfind if (mpi_timing) starttime = mpi_wtime() ! if (.not.allocated(sndbuf)) then allocate(sndbuf(mxmaglon+2,mxneed,nf),stat=ier) ! +2 for single halos if (ier /= 0) & call shutdown('mp_mag_jslot: error allocating sndbuf') endif if (.not.allocated(rcvbuf)) then allocate(rcvbuf(mxmaglon+2,mxneed,nf),stat=ier) if (ier /= 0) & call shutdown('mp_mag_jslot: error allocating rcvbuf') endif sndbuf = 0._r8 rcvbuf = 0._r8 njneed = 0 do j=1,mxneed if (jneed(j) /= -1) njneed=njneed+1 enddo if (any(jneed(1:njneed)==-1)) call shutdown('mp_mag_jslot jneed') ! call MPI_Comm_rank(cols_comm,tij,ier) ! ! Send needed lat indices to all tasks in my column: ! (redundant for alltoall) do n=0,nmagtaskj-1 jneedall(:,n) = jneed(:) enddo call mpi_alltoall(jneedall,mxneed,MPI_INTEGER, & peersneed,mxneed,MPI_INTEGER,cols_comm,ier) if (ier /= 0) & call handle_mpi_err(ier,'mp_mag_jslot call mpi_alltoall') ! ! Check if I have any needed lats, and who to send to: do n=0,nmagtaskj-1 if (n==tij) cycle njhave = 0 do j=1,mxneed if (peersneed(j,n) >= mlat00.and.peersneed(j,n) <= mlat11)then njhave = njhave+1 jhave(njhave) = peersneed(j,n) idest = n wid = itask_table_geo(mytidi,idest) endif enddo if (njhave > 0) then ! ! Load send buffer: nlons = mlon11-mlon00+1 do ifld=1,nf do j=1,njhave do i=mlon00,mlon11 sndbuf(i-mlon00+1,j,ifld) = fin(i,jhave(j),ifld) enddo enddo enddo len = nlons*njhave*nf msgid = mytid+wid call mpi_send(sndbuf(1:nlons,1:njhave,:),len,MPI_REAL8, & idest,1,cols_comm,ier) endif enddo ! n=0,nmagtaskj-1 ! ! Determine which tasks to receive which lats from. Task to ! receive from must be in same task column magtidi as I am. if (njneed > 0) then njhave = 0 jhave(:) = -1 do n=0,ntask-1 njhave = 0 do j=1,njneed if (jneed(j) >= tasks(n)%mlat0-1 .and. & jneed(j) <= tasks(n)%mlat1+1) then njhave = njhave+1 jhave(njhave) = jneed(j) endif enddo if (njhave > 0 .and. tasks(n)%magtidi==magtidi) then isrc = tasks(n)%magtidj ! task id in cols_comm to recv from nlons = mlon11-mlon00+1 len = nlons*njhave*nf msgid = mytid+n call mpi_recv(rcvbuf(1:nlons,1:njhave,:),len,MPI_REAL8, & isrc,1,cols_comm,irstat,ier) ! ! Get data from receive buffer: ! real,intent(out) :: fout(mlon00:mlon11,mxneed) ! returned data at needed lats do ifld=1,nf do j=1,njhave nj = ixfind(jneed,mxneed,jhave(j),icount) if (nj==0) call shutdown('jhave(j) not in jneed') do i=mlon00,mlon11 fout(i,nj,ifld) = rcvbuf(i-mlon00+1,j,ifld) enddo enddo ! j=1,njhave enddo ! ifld=1,nf endif ! jhave > 0 enddo ! n=0,ntask-1 endif ! njneed > 0 ! if (mpi_timing) then endtime = mpi_wtime() time_mag_jslot=time_mag_jslot+(endtime-starttime) endif end subroutine mp_mag_jslot !----------------------------------------------------------------------- subroutine handle_mpi_err(ierrcode,string) ! ! Args: integer,intent(in) :: ierrcode character(len=*) :: string ! ! Local: character(len=80) :: errstring integer :: len_errstring ! call mpi_error_string(ierrcode,errstring,len_errstring) write(iulog,"(/,'>>> mpi error: ',a)") trim(string) write(iulog,"(' ierrcode=',i3,': ',a)") trim(errstring) end subroutine handle_mpi_err !----------------------------------------------------------------------- end module my_mpi