Changeset 16b61a5 in flexpart.git for src/init_domainfill_mpi.f90


Ignore:
Timestamp:
Oct 14, 2016, 3:19:00 PM (8 years ago)
Author:
Espen Sollum ATMOS <eso@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
4c64400
Parents:
861805a
Message:

Reworked the domain-filling option (MPI). Fixed a slow loop which had errors in loop counter (MPI)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/init_domainfill_mpi.f90

    r0f7835d r16b61a5  
    4747!                                                                            *
    4848!*****************************************************************************
     49! MPI version:
     50!
     51!   -Root process allocates temporary arrays holding properties for
     52!    all particles in the simulation
     53!   -An index array is used to assign 1st particle to 1st process, 2nd particle
     54!    to 2nd process and so on so that they are evenly distibuted geographically
     55!   -Inititialization for domain-filling is done as in the serial code
     56!   -Root process distributes particles evenly to other processes
     57!
     58!*****************************************************************************
    4959
    5060  use point_mod
    … …  
    5666  implicit none
    5767
    58   integer :: j,ix,jy,kz,ncolumn,numparttot
     68! ncolumn_mpi,numparttot_mpi        ncolumn,numparttot per process
     69  integer :: j,ix,jy,kz,ncolumn,numparttot,ncolumn_mpi,numparttot_mpi, arr_size
    5970  real :: gridarea(0:nymax-1),pp(nzmax),ylat,ylatp,ylatm,hzone
    6071  real :: cosfactm,cosfactp,deltacol,dz1,dz2,dz,pnew,fractus
    … …  
    6677
    6778  integer :: idummy = -11
     79  integer,allocatable,dimension(:) :: idx ! index array
     80  integer :: stride
     81  integer, parameter :: nullsize=0
    6882  logical :: first_call=.true.
    6983
    70   ! Use different seed for each process
     84! Use different seed for each process ! TODO: not needed anymore
    7185  if (first_call) then
    7286    idummy=idummy+mp_seed
    … …  
    102116
    103117
     118! This section only done by the root process
     119!*******************************************
     120  if (lroot) then
     121! Arrays for particles to be released. Add a small number to npart(1) in case of
     122! round-off errors
     123    arr_size = npart(1) + mp_np
     124    allocate(itra1_tmp(arr_size),npoint_tmp(arr_size),nclass_tmp(arr_size),&
     125         & idt_tmp(arr_size),itramem_tmp(arr_size),itrasplit_tmp(arr_size),&
     126         & xtra1_tmp(arr_size),ytra1_tmp(arr_size),ztra1_tmp(arr_size),&
     127         & xmass1_tmp(arr_size, maxspec))
     128
     129! Index array for particles. This is to avoid having particles
     130! near edges of domain all on one process.
     131!****************************************************************************
     132    allocate(idx(npart(1)))
     133    stride = npart(1)/mp_partgroup_np
     134
     135    jj=0
     136    do j=1, stride
     137      do i=0, mp_partgroup_np-1
     138        jj = jj+1
     139        if (jj.gt.npart(1)) exit
     140        idx(jj) = i*stride+j
     141      end do
     142    end do
     143
     144! Add extra indices if npart(1) not evenly divisible by number of processes
     145    do i=1, mod(npart(1),mp_partgroup_np)
     146      jj = jj+1
     147      if (jj.gt.npart(1)) exit
     148      idx(jj) = jj
     149    end do
     150
     151! Initialize all particles as non-existent   
     152    itra1_tmp(:)=-999999999
     153
    104154! Calculate area of grid cell with formula M=2*pi*R*h*dx/360,
    105155! see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90
    106156!************************************************************
    107157
    108   do jy=ny_sn(1),ny_sn(2)      ! loop about latitudes
    109     ylat=ylat0+real(jy)*dy
    110     ylatp=ylat+0.5*dy
    111     ylatm=ylat-0.5*dy
    112     if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
    113       hzone=1./dyconst
    114     else
     158    do jy=ny_sn(1),ny_sn(2)      ! loop about latitudes
     159      ylat=ylat0+real(jy)*dy
     160      ylatp=ylat+0.5*dy
     161      ylatm=ylat-0.5*dy
     162      if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
     163        hzone=1./dyconst
     164      else
     165        cosfactp=cos(ylatp*pih)*r_earth
     166        cosfactm=cos(ylatm*pih)*r_earth
     167        if (cosfactp.lt.cosfactm) then
     168          hzone=sqrt(r_earth**2-cosfactp**2)- &
     169               sqrt(r_earth**2-cosfactm**2)
     170        else
     171          hzone=sqrt(r_earth**2-cosfactm**2)- &
     172               sqrt(r_earth**2-cosfactp**2)
     173        endif
     174      endif
     175      gridarea(jy)=2.*pi*r_earth*hzone*dx/360.
     176    end do
     177
     178! Do the same for the south pole
     179
     180    if (sglobal) then
     181      ylat=ylat0
     182      ylatp=ylat+0.5*dy
     183      ylatm=ylat
     184      cosfactm=0.
    115185      cosfactp=cos(ylatp*pih)*r_earth
     186      hzone=sqrt(r_earth**2-cosfactm**2)- &
     187           sqrt(r_earth**2-cosfactp**2)
     188      gridarea(0)=2.*pi*r_earth*hzone*dx/360.
     189    endif
     190
     191! Do the same for the north pole
     192
     193    if (nglobal) then
     194      ylat=ylat0+real(nymin1)*dy
     195      ylatp=ylat
     196      ylatm=ylat-0.5*dy
     197      cosfactp=0.
    116198      cosfactm=cos(ylatm*pih)*r_earth
    117       if (cosfactp.lt.cosfactm) then
    118         hzone=sqrt(r_earth**2-cosfactp**2)- &
    119              sqrt(r_earth**2-cosfactm**2)
    120       else
    121         hzone=sqrt(r_earth**2-cosfactm**2)- &
    122              sqrt(r_earth**2-cosfactp**2)
    123       endif
     199      hzone=sqrt(r_earth**2-cosfactp**2)- &
     200           sqrt(r_earth**2-cosfactm**2)
     201      gridarea(nymin1)=2.*pi*r_earth*hzone*dx/360.
    124202    endif
    125     gridarea(jy)=2.*pi*r_earth*hzone*dx/360.
    126   end do
    127 
    128 ! Do the same for the south pole
    129 
    130   if (sglobal) then
    131     ylat=ylat0
    132     ylatp=ylat+0.5*dy
    133     ylatm=ylat
    134     cosfactm=0.
    135     cosfactp=cos(ylatp*pih)*r_earth
    136     hzone=sqrt(r_earth**2-cosfactm**2)- &
    137          sqrt(r_earth**2-cosfactp**2)
    138     gridarea(0)=2.*pi*r_earth*hzone*dx/360.
    139   endif
    140 
    141 ! Do the same for the north pole
    142 
    143   if (nglobal) then
    144     ylat=ylat0+real(nymin1)*dy
    145     ylatp=ylat
    146     ylatm=ylat-0.5*dy
    147     cosfactp=0.
    148     cosfactm=cos(ylatm*pih)*r_earth
    149     hzone=sqrt(r_earth**2-cosfactp**2)- &
    150          sqrt(r_earth**2-cosfactm**2)
    151     gridarea(nymin1)=2.*pi*r_earth*hzone*dx/360.
    152   endif
    153203
    154204
    … …  
    156206!*********************************************************************
    157207
    158   colmasstotal=0.
    159   do jy=ny_sn(1),ny_sn(2)          ! loop about latitudes
    160     do ix=nx_we(1),nx_we(2)      ! loop about longitudes
    161       pp(1)=rho(ix,jy,1,1)*r_air*tt(ix,jy,1,1)
    162       pp(nz)=rho(ix,jy,nz,1)*r_air*tt(ix,jy,nz,1)
    163       ! Each MPI process is assigned an equal share of particles
    164       colmass(ix,jy)=(pp(1)-pp(nz))/ga*gridarea(jy)/mp_partgroup_np
    165       colmasstotal=colmasstotal+colmass(ix,jy)
    166 
    167     end do
    168   end do
    169 
    170   if (lroot) write(*,*) 'Atm. mass: ',colmasstotal
    171 
    172 
    173   if (ipin.eq.0) numpart=0
     208    colmasstotal=0.
     209    do jy=ny_sn(1),ny_sn(2)          ! loop about latitudes
     210      do ix=nx_we(1),nx_we(2)      ! loop about longitudes
     211        pp(1)=rho(ix,jy,1,1)*r_air*tt(ix,jy,1,1)
     212        pp(nz)=rho(ix,jy,nz,1)*r_air*tt(ix,jy,nz,1)
     213        colmass(ix,jy)=(pp(1)-pp(nz))/ga*gridarea(jy)
     214        colmasstotal=colmasstotal+colmass(ix,jy)
     215
     216      end do
     217    end do
     218
     219    write(*,*) 'Atm. mass: ',colmasstotal
     220
     221
     222    if (ipin.eq.0) numpart=0
    174223
    175224! Determine the particle positions
    176225!*********************************
    177226
    178   numparttot=0
    179   numcolumn=0
    180   do jy=ny_sn(1),ny_sn(2)      ! loop about latitudes
    181     ylat=ylat0+real(jy)*dy
    182     do ix=nx_we(1),nx_we(2)      ! loop about longitudes
    183       ncolumn=nint(0.999*real(npart(1)/mp_partgroup_np)*colmass(ix,jy)/ &
    184            colmasstotal)
    185       if (ncolumn.eq.0) goto 30
    186       if (ncolumn.gt.numcolumn) numcolumn=ncolumn
     227    numparttot=0
     228    numcolumn=0
     229    do jy=ny_sn(1),ny_sn(2)      ! loop about latitudes
     230      ylat=ylat0+real(jy)*dy
     231      do ix=nx_we(1),nx_we(2)      ! loop about longitudes
     232        ncolumn=nint(0.999*real(npart(1))*colmass(ix,jy)/ &
     233             colmasstotal)
     234        if (ncolumn.eq.0) goto 30
     235        if (ncolumn.gt.numcolumn) numcolumn=ncolumn
    187236
    188237! Calculate pressure at the altitudes of model surfaces, using the air density
    … …  
    190239!*****************************************************************************
    191240
    192       do kz=1,nz
    193         pp(kz)=rho(ix,jy,kz,1)*r_air*tt(ix,jy,kz,1)
    194       end do
    195 
    196 
    197       deltacol=(pp(1)-pp(nz))/real(ncolumn)
    198       pnew=pp(1)+deltacol/2.
    199       jj=0
    200       do j=1,ncolumn
    201         jj=jj+1
     241        do kz=1,nz
     242          pp(kz)=rho(ix,jy,kz,1)*r_air*tt(ix,jy,kz,1)
     243        end do
     244
     245
     246        deltacol=(pp(1)-pp(nz))/real(ncolumn)
     247        pnew=pp(1)+deltacol/2.
     248        jj=0
     249        do j=1,ncolumn
     250          jj=jj+1
    202251
    203252
    … …  
    208257
    209258
    210         if (ncolumn.gt.20) then
    211           pnew=pnew-deltacol
    212         else
    213           pnew=pp(1)-ran1(idummy)*(pp(1)-pp(nz))
    214         endif
    215 
    216         do kz=1,nz-1
    217           if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew)) then
    218             dz1=pp(kz)-pnew
    219             dz2=pnew-pp(kz+1)
    220             dz=1./(dz1+dz2)
     259          if (ncolumn.gt.20) then
     260            pnew=pnew-deltacol
     261          else
     262            pnew=pp(1)-ran1(idummy)*(pp(1)-pp(nz))
     263          endif
     264
     265          do kz=1,nz-1
     266            if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew)) then
     267              dz1=pp(kz)-pnew
     268              dz2=pnew-pp(kz+1)
     269              dz=1./(dz1+dz2)
    221270
    222271! Assign particle position
    … …  
    224273! Do the following steps only if particles are not read in from previous model run
    225274!*****************************************************************************
    226             if (ipin.eq.0) then
    227               xtra1(numpart+jj)=real(ix)-0.5+ran1(idummy)
    228               if (ix.eq.0) xtra1(numpart+jj)=ran1(idummy)
    229               if (ix.eq.nxmin1) xtra1(numpart+jj)= &
    230                    real(nxmin1)-ran1(idummy)
    231               ytra1(numpart+jj)=real(jy)-0.5+ran1(idummy)
    232               ztra1(numpart+jj)=(height(kz)*dz2+height(kz+1)*dz1)*dz
    233               if (ztra1(numpart+jj).gt.height(nz)-0.5) &
    234                    ztra1(numpart+jj)=height(nz)-0.5
     275              if (ipin.eq.0) then
     276                xtra1_tmp(idx(numpart+jj))=real(ix)-0.5+ran1(idummy)
     277                if (ix.eq.0) xtra1_tmp(idx(numpart+jj))=ran1(idummy)
     278                if (ix.eq.nxmin1) xtra1_tmp(idx(numpart+jj))= &
     279                     real(nxmin1)-ran1(idummy)
     280                ytra1_tmp(idx(numpart+jj))=real(jy)-0.5+ran1(idummy)
     281                ztra1_tmp(idx(numpart+jj))=(height(kz)*dz2+height(kz+1)*dz1)*dz
     282                if (ztra1_tmp(idx(numpart+jj)).gt.height(nz)-0.5) &
     283                     ztra1_tmp(idx(numpart+jj))=height(nz)-0.5
    235284
    236285
    237286! Interpolate PV to the particle position
    238287!****************************************
    239               ixm=int(xtra1(numpart+jj))
    240               jym=int(ytra1(numpart+jj))
    241               ixp=ixm+1
    242               jyp=jym+1
    243               ddx=xtra1(numpart+jj)-real(ixm)
    244               ddy=ytra1(numpart+jj)-real(jym)
    245               rddx=1.-ddx
    246               rddy=1.-ddy
    247               p1=rddx*rddy
    248               p2=ddx*rddy
    249               p3=rddx*ddy
    250               p4=ddx*ddy
    251               do i=2,nz
    252                 if (height(i).gt.ztra1(numpart+jj)) then
    253                   indzm=i-1
    254                   indzp=i
    255                   goto 6
    256                 endif
    257               end do
    258 6             continue
    259               dz1=ztra1(numpart+jj)-height(indzm)
    260               dz2=height(indzp)-ztra1(numpart+jj)
    261               dz=1./(dz1+dz2)
    262               do in=1,2
    263                 indzh=indzm+in-1
    264                 y1(in)=p1*pv(ixm,jym,indzh,1) &
    265                      +p2*pv(ixp,jym,indzh,1) &
    266                      +p3*pv(ixm,jyp,indzh,1) &
    267                      +p4*pv(ixp,jyp,indzh,1)
    268               end do
    269               pvpart=(dz2*y1(1)+dz1*y1(2))*dz
    270               if (ylat.lt.0.) pvpart=-1.*pvpart
     288                ixm=int(xtra1_tmp(idx(numpart+jj)))
     289                jym=int(ytra1_tmp(idx(numpart+jj)))
     290                ixp=ixm+1
     291                jyp=jym+1
     292                ddx=xtra1_tmp(idx(numpart+jj))-real(ixm)
     293                ddy=ytra1_tmp(idx(numpart+jj))-real(jym)
     294                rddx=1.-ddx
     295                rddy=1.-ddy
     296                p1=rddx*rddy
     297                p2=ddx*rddy
     298                p3=rddx*ddy
     299                p4=ddx*ddy
     300                do i=2,nz
     301                  if (height(i).gt.ztra1_tmp(idx(numpart+jj))) then
     302                    indzm=i-1
     303                    indzp=i
     304                    goto 6
     305                  endif
     306                end do
     3076               continue
     308                dz1=ztra1_tmp(idx(numpart+jj))-height(indzm)
     309                dz2=height(indzp)-ztra1_tmp(idx(numpart+jj))
     310                dz=1./(dz1+dz2)
     311                do in=1,2
     312                  indzh=indzm+in-1
     313                  y1(in)=p1*pv(ixm,jym,indzh,1) &
     314                       +p2*pv(ixp,jym,indzh,1) &
     315                       +p3*pv(ixm,jyp,indzh,1) &
     316                       +p4*pv(ixp,jyp,indzh,1)
     317                end do
     318                pvpart=(dz2*y1(1)+dz1*y1(2))*dz
     319                if (ylat.lt.0.) pvpart=-1.*pvpart
    271320
    272321
    … …  
    274323!*****************************************************************************
    275324
    276               if (((ztra1(numpart+jj).gt.3000.).and. &
    277                    (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then
     325                if (((ztra1_tmp(idx(numpart+jj)).gt.3000.).and. &
     326                     (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then
    278327
    279328! Assign certain properties to the particle
    280329!******************************************
    281                 nclass(numpart+jj)=min(int(ran1(idummy)* &
    282                      real(nclassunc))+1,nclassunc)
    283                 numparticlecount=numparticlecount+1
    284                 npoint(numpart+jj)=numparticlecount
    285                 idt(numpart+jj)=mintime
    286                 itra1(numpart+jj)=0
    287                 itramem(numpart+jj)=0
    288                 itrasplit(numpart+jj)=itra1(numpart+jj)+ldirect* &
    289                      itsplit
    290                 xmass1(numpart+jj,1)=colmass(ix,jy)/real(ncolumn)
    291                 if (mdomainfill.eq.2) xmass1(numpart+jj,1)= &
    292                      xmass1(numpart+jj,1)*pvpart*48./29.*ozonescale/10.**9
    293               else
    294                 jj=jj-1
     330                  nclass_tmp(idx(numpart+jj))=min(int(ran1(idummy)* &
     331                       real(nclassunc))+1,nclassunc)
     332                  numparticlecount=numparticlecount+1
     333                  npoint_tmp(idx(numpart+jj))=numparticlecount
     334                  idt_tmp(idx(numpart+jj))=mintime
     335                  itra1_tmp(idx(numpart+jj))=0
     336                  itramem_tmp(idx(numpart+jj))=0
     337                  itrasplit_tmp(idx(numpart+jj))=itra1_tmp(idx(numpart+jj))+ldirect* &
     338                       itsplit
     339                  xmass1_tmp(idx(numpart+jj),1)=colmass(ix,jy)/real(ncolumn)
     340                  if (mdomainfill.eq.2) xmass1_tmp(idx(numpart+jj),1)= &
     341                       xmass1_tmp(idx(numpart+jj),1)*pvpart*48./29.*ozonescale/10.**9
     342                else
     343                  jj=jj-1
     344                endif
    295345              endif
    296346            endif
    297           endif
     347          end do
    298348        end do
     349        numparttot=numparttot+ncolumn
     350        if (ipin.eq.0) numpart=numpart+jj
     35130      continue
    299352      end do
    300       numparttot=numparttot+ncolumn
    301       if (ipin.eq.0) numpart=numpart+jj
    302 30    continue
    303     end do
    304   end do
    305 
    306 
    307 ! Check whether numpart is really smaller than maxpart per process
    308 !*****************************************************************
    309 
    310   if (numpart.gt.maxpart_mpi) then
    311     write(*,*) 'numpart too large: change source in init_atm_mass.f'
    312     write(*,*) 'numpart: ',numpart,' maxpart: ',maxpart_mpi
    313   endif
    314 
    315 
    316   xmassperparticle=colmasstotal/real(numparttot)
     353    end do
     354
     355
     356! Check whether numpart is really smaller than maxpart
     357!*****************************************************
     358
     359! ESO :TODO: this warning need to be moved further up, else out-of-bounds error earlier
     360    if (numpart.gt.maxpart) then
     361      write(*,*) 'numpart too large: change source in init_atm_mass.f'
     362      write(*,*) 'numpart: ',numpart,' maxpart: ',maxpart
     363    endif
     364
     365
     366    xmassperparticle=colmasstotal/real(numparttot)
    317367
    318368
    … …  
    320370!***********************************************
    321371
    322   do j=1,numpart
    323     if ((xtra1(j).lt.0.).or.(xtra1(j).ge.real(nxmin1)).or. &
    324          (ytra1(j).lt.0.).or.(ytra1(j).ge.real(nymin1))) then
    325       itra1(j)=-999999999
    326     endif
    327   end do
     372!    do j=1,numpart
     373    do j=1,npoint(1)
     374      if ((xtra1_tmp(j).lt.0.).or.(xtra1_tmp(j).ge.real(nxmin1)).or. &
     375           (ytra1_tmp(j).lt.0.).or.(ytra1_tmp(j).ge.real(nymin1))) then
     376        itra1_tmp(j)=-999999999
     377      endif
     378    end do
    328379
    329380
    … …  
    340391!****************************************************************************
    341392
    342   fractus=real(numcolumn)/real(nz)
    343   write(*,*) 'Total number of particles at model start: ',numpart*mp_partgroup_np
    344   write(*,*) 'Maximum number of particles per column: ',numcolumn*mp_partgroup_np
    345   write(*,*) 'If ',fractus*mp_partgroup_np,' <1, better use more particles'
    346   fractus=sqrt(max(fractus,1.))/2.
    347 
    348   do jy=ny_sn(1),ny_sn(2)      ! loop about latitudes
    349     do ix=nx_we(1),nx_we(2)      ! loop about longitudes
    350       ncolumn=nint(0.999/fractus*real(npart(1)/mp_partgroup_np)*colmass(ix,jy) &
    351            /colmasstotal)
    352       if (ncolumn.gt.maxcolumn) stop 'maxcolumn too small'
    353       if (ncolumn.eq.0) goto 80
     393    fractus=real(numcolumn)/real(nz)
     394    write(*,*) 'Total number of particles at model start: ',numpart
     395    write(*,*) 'Maximum number of particles per column: ',numcolumn
     396    write(*,*) 'If ',fractus,' <1, better use more particles'
     397    fractus=sqrt(max(fractus,1.))/2.
     398
     399    do jy=ny_sn(1),ny_sn(2)      ! loop about latitudes
     400      do ix=nx_we(1),nx_we(2)      ! loop about longitudes
     401        ncolumn=nint(0.999/fractus*real(npart(1))*colmass(ix,jy) &
     402             /colmasstotal)
     403        if (ncolumn.gt.maxcolumn) stop 'maxcolumn too small'
     404        if (ncolumn.eq.0) goto 80
    354405
    355406
    … …  
    359410!************************************************************************
    360411
    361       if (ix.eq.nx_we(1)) numcolumn_we(1,jy)=ncolumn
    362       if (ix.eq.nx_we(2)) numcolumn_we(2,jy)=ncolumn
    363       if (jy.eq.ny_sn(1)) numcolumn_sn(1,ix)=ncolumn
    364       if (jy.eq.ny_sn(2)) numcolumn_sn(2,ix)=ncolumn
     412        if (ix.eq.nx_we(1)) numcolumn_we(1,jy)=ncolumn
     413        if (ix.eq.nx_we(2)) numcolumn_we(2,jy)=ncolumn
     414        if (jy.eq.ny_sn(1)) numcolumn_sn(1,ix)=ncolumn
     415        if (jy.eq.ny_sn(2)) numcolumn_sn(2,ix)=ncolumn
    365416
    366417! Calculate pressure at the altitudes of model surfaces, using the air density
    … …  
    368419!*****************************************************************************
    369420
    370       do kz=1,nz
    371         pp(kz)=rho(ix,jy,kz,1)*r_air*tt(ix,jy,kz,1)
    372       end do
     421        do kz=1,nz
     422          pp(kz)=rho(ix,jy,kz,1)*r_air*tt(ix,jy,kz,1)
     423        end do
    373424
    374425! Determine the reference starting altitudes
    375426!*******************************************
    376427
    377       deltacol=(pp(1)-pp(nz))/real(ncolumn)
    378       pnew=pp(1)+deltacol/2.
    379       do j=1,ncolumn
    380         pnew=pnew-deltacol
    381         do kz=1,nz-1
    382           if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew)) then
    383             dz1=pp(kz)-pnew
    384             dz2=pnew-pp(kz+1)
    385             dz=1./(dz1+dz2)
    386             zposition=(height(kz)*dz2+height(kz+1)*dz1)*dz
    387             if (zposition.gt.height(nz)-0.5) zposition=height(nz)-0.5
     428        deltacol=(pp(1)-pp(nz))/real(ncolumn)
     429        pnew=pp(1)+deltacol/2.
     430        do j=1,ncolumn
     431          pnew=pnew-deltacol
     432          do kz=1,nz-1
     433            if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew)) then
     434              dz1=pp(kz)-pnew
     435              dz2=pnew-pp(kz+1)
     436              dz=1./(dz1+dz2)
     437              zposition=(height(kz)*dz2+height(kz+1)*dz1)*dz
     438              if (zposition.gt.height(nz)-0.5) zposition=height(nz)-0.5
    388439
    389440! Memorize vertical positions where particles are introduced
    … …  
    391442!***********************************************************
    392443
    393             if (ix.eq.nx_we(1)) zcolumn_we(1,jy,j)=zposition
    394             if (ix.eq.nx_we(2)) zcolumn_we(2,jy,j)=zposition
    395             if (jy.eq.ny_sn(1)) zcolumn_sn(1,ix,j)=zposition
    396             if (jy.eq.ny_sn(2)) zcolumn_sn(2,ix,j)=zposition
     444              if (ix.eq.nx_we(1)) zcolumn_we(1,jy,j)=zposition
     445              if (ix.eq.nx_we(2)) zcolumn_we(2,jy,j)=zposition
     446              if (jy.eq.ny_sn(1)) zcolumn_sn(1,ix,j)=zposition
     447              if (jy.eq.ny_sn(2)) zcolumn_sn(2,ix,j)=zposition
    397448
    398449! Initialize mass that has accumulated at boundary to zero
    399450!*********************************************************
    400451
    401             acc_mass_we(1,jy,j)=0.
    402             acc_mass_we(2,jy,j)=0.
    403             acc_mass_sn(1,jy,j)=0.
    404             acc_mass_sn(2,jy,j)=0.
    405           endif
     452              acc_mass_we(1,jy,j)=0.
     453              acc_mass_we(2,jy,j)=0.
     454              acc_mass_sn(1,jy,j)=0.
     455              acc_mass_sn(2,jy,j)=0.
     456            endif
     457          end do
    406458        end do
     45980      continue
    407460      end do
    408 80    continue
    409     end do
    410   end do
     461    end do
    411462
    412463! If particles shall be read in to continue an existing run,
    … …  
    415466!***************************************************************************
    416467
    417 ! :TODO: eso: parallelize
    418   if (ipin.eq.1) then
    419     open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
    420          form='unformatted')
    421     read(unitboundcond) numcolumn_we,numcolumn_sn, &
    422          zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn
    423     close(unitboundcond)
    424   endif
    425 
    426 
     468! eso TODO: only needed for root process
     469    if (ipin.eq.1) then
     470      open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
     471           form='unformatted')
     472      read(unitboundcond) numcolumn_we,numcolumn_sn, &
     473           zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn
     474      close(unitboundcond)
     475    endif
     476
     477    numpart = npart(1)/mp_partgroup_np
     478    if (mod(numpart,mp_partgroup_np).ne.0) numpart=numpart+1
     479
     480  else ! Allocate dummy arrays for receiving processes
     481    allocate(itra1_tmp(nullsize),npoint_tmp(nullsize),nclass_tmp(nullsize),&
     482         & idt_tmp(nullsize),itramem_tmp(nullsize),itrasplit_tmp(nullsize),&
     483         & xtra1_tmp(nullsize),ytra1_tmp(nullsize),ztra1_tmp(nullsize),&
     484         & xmass1_tmp(nullsize, nullsize))
     485   
     486  end if ! end if(lroot) 
     487
     488
     489! Distribute particles to other processes (numpart is 'per-process', not total)
     490    call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
     491! eso TODO: xmassperparticle: not necessary to send
     492    call MPI_Bcast(xmassperparticle, 1, mp_sp, id_root, mp_comm_used, mp_ierr)
     493    call mpif_send_part_properties(numpart)
     494
     495! Deallocate the temporary arrays used for all particles
     496    deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,&
     497         & itrasplit_tmp,xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp)
    427498
    428499
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG