!*********************************************************************** !* Copyright 2012,2013 * !* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, * !* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,* !* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, * !* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, * !* * !* This file is part of FLEXPART WRF * !* * !* FLEXPART is free software: you can redistribute it and/or modify * !* it under the terms of the GNU General Public License as published by* !* the Free Software Foundation, either version 3 of the License, or * !* (at your option) any later version. * !* * !* FLEXPART is distributed in the hope that it will be useful, * !* but WITHOUT ANY WARRANTY; without even the implied warranty of * !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * !* GNU General Public License for more details. * !* * !* You should have received a copy of the GNU General Public License * !* along with FLEXPART. If not, see . * !*********************************************************************** subroutine timemanager_mpi(mts) !******************************************************************************* ! * ! Handles the computation of trajectories, i.e. determines which * ! trajectories have to be computed at what time. * ! Manages dry+wet deposition routines, radioactive decay and the computation * ! of concentrations. * ! * ! Author: A. Stohl * ! * ! 20 May 1996 * ! * ! Dec 2005, J. Fast - Only call conccalc & concoutput when (iout.ge.1) * ! Aug 2007, W. Wang - call KFeta convection scheme (lconvection=2or3) ! Note, backward is unavailabe for lconvection=2 ! Mar 2012, J. Brioude: modifications to handle openmp and mpi * !******************************************************************************* ! Changes, Bernd C. Krueger, Feb. 2001: * ! Call of convmix when new windfield is read * !------------------------------------ * ! Changes Petra Seibert, Sept 2002 * ! fix wet scavenging problem * ! Code may not be correct for decay of deposition! * ! Changes Petra Seibert, Nov 2002 * ! call convection BEFORE new fields are read in BWD mode * ! Changes Caroline Forster, Feb 2005 ! new interface between flexpart and convection scheme ! Emanuel's latest subroutine convect43c.f is used !******************************************************************************* ! * ! Variables: * ! DEP .true. if either wet or dry deposition is switched on * ! decay(maxspec) [1/s] decay constant for radioactive decay * ! DRYDEP .true. if dry deposition is switched on * ! ideltas [s] modelling period * ! jtime [s] actual temporal position of calculation * ! ldeltat [s] time since computation of radioact. decay of depositions * ! loutaver [s] averaging period for concentration calculations * ! loutend [s] end of averaging for concentration calculations * ! loutnext [s] next time at which output fields shall be centered * ! loutsample [s] sampling interval for averaging of concentrations * ! loutstart [s] start of averaging for concentration calculations * ! loutstep [s] time interval for which concentrations shall be calculated* ! npoint(maxpart) index, which starting point the trajectory has * ! starting positions of trajectories * ! nstop serves as indicator for fate of particles * ! in the particle loop * ! nstop1 serves as indicator for wind fields (see getfields) * ! outnum number of samples for each concentration calculation * ! outnum number of samples for each concentration calculation * ! prob probability of absorption at ground due to dry deposition * ! WETDEP .true. if wet deposition is switched on * ! weight weight for each concentration sample (1/2 or 1) * ! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to turbulence * ! us(maxpart),vs(maxpart),ws(maxpart) = random velocities due to interpolation * ! xtra1(maxpart), ytra1(maxpart), ztra1(maxpart) = * ! spatial positions of trajectories * ! * ! Constants: * ! maxpart maximum number of trajectories * ! * !******************************************************************************* ! include 'includepar' ! include 'includecom' use unc_mod use point_mod ! use xmass_mod use flux_mod use outg_mod use oh_mod use par_mod use com_mod use mpi_mod use mt_stream ! use ran_mod ! use interpol_mod implicit none include 'mpif.h' integer :: ix,jy,j,ks,kp,l,n,jtime,nstop,nstop1 ! integer :: MPI_COMM_WORLD ! integer :: ksp integer :: loutnext,loutstart,loutend,jj,chunksize !,chunksize2 integer :: chunksize3,omp_get_num_threads integer :: ldeltat,itage,nage,th_itra1,i real :: outnum,weight,prob(maxspec),nrand,decfact ! real :: uap(maxpart),ucp(maxpart),uzp(maxpart) ! real :: us(maxpart),vs(maxpart),ws(maxpart) ! integer(kind=2) :: cbt(maxpart) ! real,allocatable, dimension (:) :: uap,ucp,uzp ! real,allocatable, dimension (:) :: us,vs,ws ! integer(kind=2),allocatable, dimension (:) :: cbt real :: drydeposit(maxspec),gridtotalunc,wetgridtotalunc real :: drygridtotalunc,xold,yold,zold,xmassfract ! integer j,k,l,n,jtime,nstop,nstop1 ! integer loutnext,loutstart,loutend ! integer ix,jy,ldeltat,itage,nage ! real outnum,weight,prob(maxspec) ! real uap(maxpart),ucp(maxpart),uzp(maxpart),decfact ! real us(maxpart),vs(maxpart),ws(maxpart),cbt(maxpart) ! real drydeposit(maxspec),gridtotalunc,wetgridtotalunc ! real drygridtotalunc,xold,yold,zold ! real xm,xm1 integer :: th_npoint,th_idt,th_itramem,jdeb,jfin,stat,th_nclass ! integer,save :: cpt(24)=0 integer,save :: cpt(maxomp)=0 real(kind=dp) :: th_xtra1,th_ytra1 real :: th_ztra1,th_uap,th_ucp,th_uzp real :: th_us,th_vs,th_ws,ran3 integer(kind=2) :: th_cbt integer :: OMP_GET_THREAD_NUM,from real :: p1,p2,p3,p4,ddx,ddy,rddx,rddy,dtt,dt1,dt2 integer :: ixp,jyp,ngrid,indz,indzp,nbp,jj2,ii,offset logical :: depoindicator(maxspec) logical,save :: indzindicator(nzmax) real :: ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw real :: sigw,dsigwdz,dsigw2dz,th_xmass1(maxspec) real :: start, finish real :: uprof(nzmax),vprof(nzmax),wprof(nzmax) real :: usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax) real :: rhoprof(nzmax),rhogradprof(nzmax) real :: tkeprof(nzmax),pttprof(nzmax) real :: u,v,w,usig,vsig,wsig,pvi integer*4 :: now(3) integer :: ttime,cpttra integer, dimension(MPI_STATUS_SIZE) :: status integer :: myid,ntasks,ierr,islave,tag2,ompid,n_threads,tag3,i_omp type (mt_state) :: mts (0: MAX_STREAM) !************************ !JB call MPI_COMM_RANK ( MPI_COMM_WORLD, myid, ierr ) call MPI_COMM_SIZE ( MPI_COMM_WORLD, ntasks, ierr ) ! myid gives the info on the node id loutnext=loutstep/2 outnum=0. loutstart=loutnext-loutaver/2 loutend=loutnext+loutaver/2 ! if (myid.eq.0) then allocate(uap(maxpart) ,stat=stat) allocate(ucp(maxpart) ,stat=stat) allocate(uzp(maxpart) ,stat=stat) allocate(us(maxpart) ,stat=stat) allocate(vs(maxpart) ,stat=stat) allocate(ws(maxpart) ,stat=stat) allocate(cbt(maxpart) ,stat=stat) ! endif ! if (chunksize2.eq.0) chunksize2=1 chunksize2=int(maxpart/ntasks)+1 !if sent homogeneously ! print*,'chunk',myid,chunksize2,numpart if (ntasks.gt.1) then allocate(mpi_npoint(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not 1' allocate(mpi_idt(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not 2' allocate(mpi_itra1(chunksize2) ,stat=stat) allocate(mpi_itramem(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not 3' allocate(mpi_uap(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not 4' allocate(mpi_ucp(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not 5' allocate(mpi_uzp(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not 6' allocate(mpi_us(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not 7' allocate(mpi_vs(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not 8' allocate(mpi_ws(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not 82' allocate(mpi_xtra1(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not 9' allocate(mpi_ytra1(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not10' allocate(mpi_ztra1(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not11' allocate(mpi_cbt(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not12' allocate(mpi_xmass1(chunksize2,nspec) ,stat=stat) ! allocate(mpi_drydep1(chunksize2,nspec) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not13' allocate(mpi_nclass(chunksize2) ,stat=stat) allocate(dummyi2(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not14' allocate(dummyr2(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not15' allocate(dummyi22(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not16' allocate(dummyr22(chunksize2) ,stat=stat) if (stat.ne.0) write(*,*)'ERROR: could not17' chunksize2=chunksize endif !********************************************************************** ! Loop over the whole modelling period in time steps of mintime seconds !********************************************************************** ! print*,'time',myid,ideltas,lsynctime do jtime=0,ideltas,lsynctime ! print*,'jtime',jtime ! Computation of wet deposition every lsynctime seconds ! maybe wet depo frequency can be relaxed later but better be on safe side ! wetdepo must be called BEFORE new fields are read in but should not ! be called in the very beginning before any fields are loaded, or ! before particles are in the system ! Code may not be correct for decay of deposition ! changed by Petra Seibert 9/02 !******************************************************************** if (WETDEP .and. jtime .ne. 0 .and. numpart .gt. 0) & call wetdepo(jtime,lsynctime,loutnext) if (OHREA .and. jtime .ne. 0 .and. numpart .gt. 0) & call ohreaction(jtime,lsynctime,loutnext) ! compute convection for backward runs !************************************* ! if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(jtime.lt.0)) ! & call convmix(jtime) if ((ldirect.eq.-1).and.(jtime.lt.0)) then if (lconvection .eq. 1) call convmix(jtime) if (lconvection .eq. 2 .or. lconvection .eq. 3) & call convmix_kfeta(jtime) endif ! Get necessary wind fields if not available !******************************************* ! call itime(now) ! ttime=now(1)*3600+now(2)*60+now(3) call cpu_time(start) call getfields(jtime,nstop1) if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE' ! call itime(now) ! ttime=now(1)*3600+now(2)*60+now(3)-ttime call cpu_time(finish) ! print*,'read wind time',ttime if (option_verbose.eq.1) then print*,'read wind time',finish-start endif ! Release particles !****************** !JB if (myid.eq.0) then ! I let only the master thread releasing the particles and calculate the output ! call itime(now) call cpu_time(start) if (mdomainfill.ge.1) then if (jtime.eq.0) then call init_domainfill() else call boundcond_domainfill(jtime,loutend) endif else if (numpoint_option.eq.0) then call releaseparticles_irreg(jtime) elseif (numpoint_option.eq.1) then ! print*,'avant release' call releaseparticles_reg(jtime) endif endif ! do i=1,numpart ! print*,'ipart 2',myid,i,ztra1(i) ! enddo ! print*,'test rel',npoint(1),npoint(2),npoint(3) ! print*,'test rel1',npoint(5139),npoint(6002),npoint(100003) ! Compute convective mixing for forward runs ! for backward runs it is done before next windfield is read in !************************************************************** ! if ((ldirect.eq.1).and.(lconvection.eq.1)) & ! call convmix(jtime) if (ldirect.eq.1) then if (lconvection.eq.1)call convmix(jtime) if (lconvection.eq.2 .or. lconvection .eq. 3) & call convmix_kfeta(jtime) endif ! print*,'intermediaire' ! If middle of averaging period of output fields is reached, accumulated ! deposited mass radioactively decays !*********************************************************************** if (DEP.and.(jtime.eq.loutnext).and.(ldirect.gt.0)) then do ks=1,nspec do kp=1,maxpointspec_act if (decay(ks).gt.0.) then do nage=1,nageclass do l=1,nclassunc ! Mother output grid do jy=0,numygrid-1 do ix=0,numxgrid-1 wetgridunc(ix,jy,ks,kp,l,nage)= & wetgridunc(ix,jy,ks,kp,l,nage)* & exp(-1.*outstep*decay(ks)) drygridunc(ix,jy,ks,kp,l,nage)= & drygridunc(ix,jy,ks,kp,l,nage)* & exp(-1.*outstep*decay(ks)) end do end do ! Nested output grid if (nested_output.eq.1) then do jy=0,numygridn-1 do ix=0,numxgridn-1 wetgriduncn(ix,jy,ks,kp,l,nage)= & wetgriduncn(ix,jy,ks,kp,l,nage)* & exp(-1.*outstep*decay(ks)) drygriduncn(ix,jy,ks,kp,l,nage)= & drygriduncn(ix,jy,ks,kp,l,nage)* & exp(-1.*outstep*decay(ks)) end do end do endif end do end do endif end do end do endif !!! CHANGE: These lines may be switched on to check the conservation !!! of mass within FLEXPART ! if (mod(jtime,loutsample).eq.0) then ! xm=0. ! xm1=0. ! do 247 j=1,numpart !47 if (itra1(j).eq.jtime) xm1=xm1+xmass1(j,1) ! xm=xm1 ! do 248 nage=1,nageclass ! do 248 ix=0,numxgrid-1 ! do 248 jy=0,numygrid-1 ! do 248 l=1,nclassunc !48 xm=xm+wetgridunc(ix,jy,1,l,nage)+drygridunc(ix,jy,1,l,nage) ! write(*,'(i6,4f10.3)') jtime,xm,xm1 ! endif !!! CHANGE ! Check whether concentrations are to be calculated !************************************************** if ((ldirect*jtime.ge.ldirect*loutstart).and. & (ldirect*jtime.le.ldirect*loutend)) then ! add to grid if (mod(jtime-loutstart,loutsample).eq.0) then ! If we are exactly at the start or end of the concentration averaging interval, ! give only half the weight to this sample !******************************************************************************* if ((jtime.eq.loutstart).or.(jtime.eq.loutend)) then weight=0.5 else weight=1.0 endif ! print*,'avant conccalc' outnum=outnum+weight if(iout.ge.1) then if (outgrid_option.eq.0) then call conccalc_irreg(jtime,weight) elseif (outgrid_option.eq.1) then call conccalc_reg(jtime,weight) endif endif endif ! print*,'apres conccalc' ! if ((mquasilag.eq.1).and.(jtime.eq.(loutstart+loutend)/2)) & ! call partoutput_short(jtime) ! dump particle positions in extremely compressed format ! Output and reinitialization of grid ! If necessary, first sample of new grid is also taken !***************************************************** if ((jtime.eq.loutend).and.(outnum.gt.0.)) then ! print*,'iout',iout,ipout,outgrid_option if ((iout.le.3.).or.(iout.eq.5)) then if(iout.ge.1) then if (outgrid_option.eq.0) then call concoutput_irreg(jtime,outnum,gridtotalunc, & wetgridtotalunc,drygridtotalunc) if (nested_output.eq.1) call concoutput_nest_irreg(jtime,outnum) elseif (outgrid_option.eq.1) then call concoutput_reg(jtime,outnum,gridtotalunc, & wetgridtotalunc,drygridtotalunc) if (nested_output.eq.1) call concoutput_nest_reg(jtime,outnum) endif endif ! print*,'apres concout' ! if (nested_output.eq.1.and.iout.ge.1) ! + call concoutput_nest(jtime,outnum) outnum=0. endif if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(jtime) if (iflux.eq.1) call fluxoutput(jtime) write(*,45) jtime,numpart,gridtotalunc,wetgridtotalunc, & drygridtotalunc 45 format(i9,' SECONDS SIMULATED: ',i9, & ' PARTICLES: Uncertainty: ',3f7.3) if (ipout.ge.1) call partoutput(jtime) ! dump particle positions loutnext=loutnext+loutstep loutstart=loutnext-loutaver/2 loutend=loutnext+loutaver/2 if (jtime.eq.loutstart) then weight=0.5 outnum=outnum+weight if(iout.ge.1) then if (outgrid_option.eq.0) then call conccalc_irreg(jtime,weight) elseif (outgrid_option.eq.1) then call conccalc_reg(jtime,weight) endif endif endif ! Check, whether particles are to be split: ! If so, create new particles and attribute all information from the old ! particles also to the new ones; old and new particles both get half the ! mass of the old ones !************************************************************************ if (ldirect*jtime.ge.ldirect*itsplit) then n=numpart do j=1,numpart if (ldirect*jtime.ge.ldirect*itrasplit(j)) then if (n.lt.maxpart) then n=n+1 itrasplit(j)=2*(itrasplit(j)-itramem(j))+itramem(j) itrasplit(n)=itrasplit(j) itramem(n)=itramem(j) itra1(n)=itra1(j) idt(n)=idt(j) npoint(n)=npoint(j) nclass(n)=nclass(j) xtra1(n)=xtra1(j) ytra1(n)=ytra1(j) ztra1(n)=ztra1(j) uap(n)=uap(j) ucp(n)=ucp(j) uzp(n)=uzp(j) us(n)=us(j) vs(n)=vs(j) ws(n)=ws(j) cbt(n)=cbt(j) do ks=1,nspec xmass1(j,ks)=xmass1(j,ks)/2. xmass1(n,ks)=xmass1(j,ks) end do endif endif end do numpart=n endif endif endif ! Loop over all particles !************************ ! chunksize=int(numpart/ntasks)+1 !if sent homogeneously chunksize=int(numpart/ntasks) !if sent homogeneously ! call itime(now) ! ttime=now(1)*3600+now(2)*60+now(3)-ttime call cpu_time(finish) ! print*,'processing time',ttime if (option_verbose.eq.1) then print*,'processing time',finish-start endif endif !over myid !JB ! at this stage, I assume that each node has the same shared memory because they run getfields. ! now we need to split the trajectories into pieces for each node ! if (myid.eq.0) then if (jtime.eq.ideltas) exit ! Compute interval since radioactive decay of deposited mass was computed !************************************************************************ if (jtime.lt.loutnext) then ldeltat=jtime-(loutnext-loutstep) else ! first half of next interval ldeltat=jtime-loutnext endif if (myid.eq.0) then ! call itime(now) ! ttime=now(1)*3600+now(2)*60+now(3) call cpu_time(start) do ii=1,ntasks-1 call MPI_SEND(chunksize,1, MPI_INTEGER, ii,3001, MPI_COMM_WORLD, ierr) call MPI_SEND(numpart,1, MPI_INTEGER, ii,3002, MPI_COMM_WORLD, ierr) enddo else call MPI_RECV(chunksize,1, MPI_INTEGER, 0,3001, MPI_COMM_WORLD,status, ierr) call MPI_RECV(numpart,1, MPI_INTEGER, 0,3002, MPI_COMM_WORLD,status, ierr) endif ! print*,'numpart',numpart ! chunksize2=chunksize ! print*,'chunksize 0',chunksize2,ntasks !JB ! here I am going to send the infos to each slave nodes. if (numpart.gt.0 .and. ntasks.gt.1 ) then call MPI_BARRIER(MPI_COMM_WORLD,ierr) call sendint_mpi(1,numpart,chunksize,0) ! print*,'after npoint',myid,numpart,chunksize call sendint_mpi(2,numpart,chunksize,0) ! print*,'after idt',myid,numpart,chunksize call sendint_mpi(3,numpart,chunksize,0) ! print*,'after itra1',myid,numpart,chunksize call sendreal_mpi(4,numpart,chunksize,0) ! print*,'after uap',myid,numpart,chunksize call sendreal_mpi(5,numpart,chunksize,0) ! print*,'after ucp',myid,numpart,chunksize call sendreal_mpi(6,numpart,chunksize,0) ! print*,'after uzp',myid,numpart,chunksize call sendreal_mpi(7,numpart,chunksize,0) ! print*,'after us',myid,numpart,chunksize call sendreal_mpi(8,numpart,chunksize,0) ! print*,'after vs',myid,numpart,chunksize call sendreal_mpi(9,numpart,chunksize,0) ! print*,'after ws',myid,numpart,chunksize call sendreal_mpi(10,numpart,chunksize,0) ! print*,'after ztra1',myid,numpart,chunksize call senddouble_mpi(11,numpart,chunksize,0) ! print*,'after xtra1',myid,numpart,chunksize call senddouble_mpi(12,numpart,chunksize,0) ! print*,'after ytra1',myid,numpart,chunksize call sendint2_mpi(13,numpart,chunksize,0) ! print*,'after cbt',myid,numpart,chunksize call sendint_mpi(14,numpart,chunksize,0) ! print*,'after itramem',myid,numpart,chunksize call sendint_mpi(15,numpart,chunksize,0) ! print*,'after nclass',myid,numpart,chunksize call sendreal2d_mpi(20,numpart,nspec,chunksize,0) ! print*,'after xmass1',myid,numpart,chunksize if (myid.eq.0) then ! call itime(now) ! ttime=now(1)*3600+now(2)*60+now(3)-ttime call cpu_time(finish) ! print*,'sending time',ttime if (option_verbose.eq.1) then print*,'sending time',finish-start endif else chunksize2=chunksize endif endif !if statement on numpart et ntasks ! print*,'debut chunksize',chunksize,chunksize2,myid ! sigw,dsigwdz,dsigw2dz,cpt(nbp),ompid) ! initialize the temporary drydeposition grid if (DRYDEP.and.ldirect.gt.0) then do ks=1,nspec do kp=1,maxpointspec_act do nage=1,nageclass do jy=0,numygrid-1 do ix=0,numxgrid-1 do l=1,nclassunc drygridunc2(ix,jy,ks,kp,l,nage)=0. end do end do end do if (nested_output.eq.1) then do jy=0,numygridn-1 do ix=0,numxgridn-1 do l=1,nclassunc drygriduncn2(ix,jy,ks,kp,l,nage)=0. end do end do end do endif end do end do end do endif ! !JB ! now we are entering the openmp zone. ! if (myid.eq.0) then ! print*,'itra11',numpart,itra1(numpart-2:numpart+2) ! print*,'itra12',chunksize2,mpi_itra1(chunksize2-2:chunksize2+1) ! else ! print*,'itra13',mpi_itra1(chunksize-2:chunksize+1) ! endif ! print*,'continue',myid,chunksize2 !!!$OMP PARALLEL NUM_THREADS(10) DEFAULT(SHARED) & !$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(jj, th_npoint, th_idt, th_uap, th_ucp, & !$OMP th_uzp, th_us, th_vs, th_ws, th_xtra1, th_ytra1, th_ztra1,decfact, & !$OMP th_cbt, xold, yold, zold, kp, itage, prob, nstop, xmassfract, & !$OMP th_nclass,chunksize3,start,finish,ngrid,ompid,depoindicator,nbp, & !$OMP indzindicator,cpttra,th_xmass1,th_itra1,th_itramem,drydeposit,n_threads) & !$OMP SHARED(height,rho,tt,vsetaver,dquer,xtra1,ytra1,ztra1, & !$OMP density,cunningham,itra1,ioutputforeachrelease,cbt,iflux, & !$OMP uun,vvn,wwn,ustar,wstar,oli,uupol,vvpol,uu,vv,ww,drhodz,ptt,tke, & !$OMP rhon,drhodzn,pttn,tken,vdep,vdepn,itramem,nageclass,lage, & !$OMP jtime,ldirect,memind,nglobal,switchnorthg,m_x,m_y,m_xn,m_yn, & !$OMP switchsouthg,numbnests,xln,xrn,yln,yrn,memtime,xresoln, & !$OMP yresoln,hmix,hmixn,tropopause, & !$OMP tropopausen,lsynctime,dxconst,dyconst,mdomainfill, & !$OMP turb_option,turbswitch,ifine,chunksize,chunksize2, & !!!maxrand, & !$OMP xmass,xmass1,DRYDEP,DRYDEPSPEC,nspec,rannumb,uniform_rannumb,cpt, & !$OMP lwindinterv,npart,npoint,idt,uap,ucp,uzp,us,vs,ws, & !$OMP linit_cond,decay,ldeltat,nclass,nested_output,numpart, & !$OMP mpi_npoint,mpi_idt, mpi_uap, mpi_ucp, mpi_uzp, mpi_us, mpi_vs, & !$OMP mpi_ws, mpi_xtra1, mpi_ytra1, mpi_ztra1, & !$OMP mpi_cbt,drygridunc2,drygriduncn2, & !$OMP mpi_xmass1, mpi_itra1,myid,mpi_itramem,mpi_nclass, & !$OMP mts) ! call itime(now) ! ttime=now(1)*3600+now(2)*60+now(3) call cpu_time(start) ! chunksize3=int(chunksize2/omp_get_num_threads())+1 n_threads=omp_get_num_threads() ! ompid=omp_get_num_threads() ompid=OMP_GET_THREAD_NUM() chunksize3=int(real(chunksize2)/real(n_threads)/20.)+1 !more efficient if (ompid+1.gt.maxomp) then print*,'problem with maxomp. modify par_mod.f90',maxomp,ompid+1 stop endif cpttra=0 ! print*,'chunksi',chunksize2,chunksize3,myid if (chunksize2.gt.0 .and. numpart.gt.0) then ! print*,'test rel2',npoint(5139),npoint(6002),npoint(100003) !!!$OMP DO SCHEDULE(STATIC,10) !$OMP DO SCHEDULE(STATIC,chunksize3) !!!$OMP DO SCHEDULE(GUIDED,1) !!!$OMP DO ! do jj=1,numpart ! do jj=numpart,1,-1 ! print*,jj do jj=1,chunksize2 ! If integration step is due, do it !********************************** !$OMP CRITICAL if (ntasks.gt.1) then th_itra1=mpi_itra1(jj) th_itramem=mpi_itramem(jj) th_npoint=mpi_npoint(jj) else th_itra1=itra1(jj) th_itramem=itramem(jj) th_npoint=npoint(jj) endif !$OMP END CRITICAL ! if (th_itra1(jj).eq.jtime) then ! if (mod(jj,30000).eq.1) print*,'middle',jj,th_itra1,jtime ! if (th_itra1.lt.-9999) print*,'middle',jj,th_itra1,jtime ! if (th_itra1.eq.jtime) then ! if (th_itra1.eq.jtime .and. th_npoint.gt.0 .and. th_npoint.le.numpoint) then if (th_itra1.eq.jtime ) then cpttra=cpttra+1 ! print*,'avant init',j ! Initialize newly released particle !*********************************** !$OMP CRITICAL if (ntasks.eq.1) then ! th_npoint=npoint(jj) th_idt=idt(jj) th_uap=uap(jj) th_ucp=ucp(jj) th_uzp=uzp(jj) th_us=us(jj) th_vs=vs(jj) th_ws=ws(jj) th_xtra1=xtra1(jj) th_ytra1=ytra1(jj) th_ztra1=ztra1(jj) th_nclass=nclass(jj) th_cbt=cbt(jj) do ks=1,nspec th_xmass1(ks)=xmass1(jj,ks) enddo if (ioutputforeachrelease.eq.1) then kp=npoint(jj) else kp=1 endif else ! th_npoint=mpi_npoint(jj) th_idt=mpi_idt(jj) th_uap=mpi_uap(jj) th_ucp=mpi_ucp(jj) th_uzp=mpi_uzp(jj) th_us=mpi_us(jj) th_vs=mpi_vs(jj) th_ws=mpi_ws(jj) th_xtra1=mpi_xtra1(jj) th_ytra1=mpi_ytra1(jj) th_ztra1=mpi_ztra1(jj) th_nclass=mpi_nclass(jj) th_cbt=mpi_cbt(jj) do ks=1,nspec th_xmass1(ks)=mpi_xmass1(jj,ks) enddo if (ioutputforeachrelease.eq.1) then kp=mpi_npoint(jj) else kp=1 endif endif !$OMP END CRITICAL ! Determine age class of the particle ! itage=abs(itra1(jj)-itramem(jj)) itage=abs(th_itra1-th_itramem) do nage=1,nageclass if (itage.lt.lage(nage)) exit enddo ! if (jj.lt.5) print*,'xmass1',th_xmass1(1) ! ompid=OMP_GET_THREAD_NUM() nbp=ompid+1 ! print*,th_npoint,jj,npoint(jj) ! print*,'befo',th_xtra1,th_ytra1,th_ztra1 ! iff=0 ! if ((itramem(jj).eq.jtime).or.(jtime.eq.0)) & if ((th_itramem.eq.jtime).or.(jtime.eq.0)) then ! call initialize(jtime,idt(j),uap(j),ucp(j),uzp(j), & ! us(j),vs(j),ws(j),xtra1(j),ytra1(j),ztra1(j),cbt(j)) call initialize(jtime,th_idt,th_uap,th_ucp,th_uzp, & th_us,th_vs,th_ws,th_xtra1,th_ytra1,th_ztra1,th_cbt, & ngrid,depoindicator,indzindicator,cpt(nbp),ompid,myid,n_threads,mts ) endif ! print*,'after',th_xtra1,th_ytra1,th_ztra1 ! Memorize particle positions !**************************** ! if (mod(jj,100000).eq.1) print*,'middle',jj,myid,ompid ! xold=xtra1(j) ! yold=ytra1(j) ! zold=ztra1(j) xold=th_xtra1 yold=th_ytra1 zold=th_ztra1 ! Integrate Lagevin equation for lsynctime seconds !************************************************* ! write(*,*)'numpart,jtime, particle #=',numpart,jtime,j ! call advance(jtime,npoint(j),idt(j),uap(j),ucp(j),uzp(j),us(j), & ! vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob,cbt(j)) ! if ( abs(xold).gt.1000. .or. xold.ne.xold .or. th_xtra1.ne.th_xtra1 ) & ! print*,'pb avant 0',xold,yold,th_itramem, & ! jtime,jj,chunksize2 ! print*,'npoint',mpi_npoint(jj),th_npoint,jj call advance(jtime,th_npoint,th_idt,th_uap,th_ucp,th_uzp, & th_us,th_vs,th_ws,nstop,th_xtra1,& th_ytra1,th_ztra1,prob,th_cbt, & ngrid,depoindicator,indzindicator,cpt(nbp),ompid,myid,n_threads,mts ) ! if ( abs(xold).gt.1000. .or. xold.ne.xold .or. th_xtra1.ne.th_xtra1 ) & ! print*,'pb avant 1',xold,yold,th_itramem, & ! jtime,jj,chunksize2 ! Calculate the gross fluxes across layer interfaces !*************************************************** if (iflux.eq.1) call calcfluxes(nage,jj,xold,yold,zold) ! if (jj.lt.5) print*,'coord after',myid,th_itra1,th_xmass1(1),DRYDEPSPEC(ks) ! Determine, when next time step is due ! If trajectory is terminated, mark it !************************************** !!!$OMP CRITICAL do ks=1,nspec drydeposit(ks)=0. enddo if (nstop.gt.1) then if (linit_cond.ge.1) call initial_cond_calc(jtime,jj) ! itra1(jj)=-999999999 th_itra1=-999999999 if (option_verbose.gt.1) print*,'out of domain',th_xtra1,th_ytra1,th_ztra1 else ! itra1(jj)=jtime+lsynctime th_itra1=jtime+lsynctime ! if (jj.lt.5) print*,'coord after2',myid,th_itra1,th_xmass1(1),DRYDEPSPEC(ks) ! Dry deposition and radioactive decay for each species !****************************************************** xmassfract=0. do ks=1,nspec if (decay(ks).gt.0.) then ! radioactive decay decfact=exp(-real(abs(lsynctime))*decay(ks)) else decfact=1. endif if (DRYDEPSPEC(ks)) then ! dry deposition ! drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact drydeposit(ks)=th_xmass1(ks)*prob(ks)*decfact ! xmass1(j,ks)=xmass1(j,ks)*(1.-prob(ks))*decfact th_xmass1(ks)=th_xmass1(ks)*(1.-prob(ks))*decfact if (decay(ks).gt.0.) then ! correct for decay (see wetdepo) drydeposit(ks)=drydeposit(ks)* & exp(real(abs(ldeltat))*decay(ks)) endif else ! no dry deposition ! xmass1(j,ks)=xmass1(j,ks)*decfact th_xmass1(ks)=th_xmass1(ks)*decfact endif ! if (jj.lt.5) print*,'coord after3',myid,th_itra1,th_xmass1(1),DRYDEPSPEC(ks),xmass(th_npoint,1) if (mdomainfill.eq.0) then if (xmass(th_npoint,ks).gt.0.) & ! xmassfract=max(xmassfract,real(npart(npoint(jj)))* & xmassfract=max(xmassfract,real(npart(th_npoint))* & ! xmass1(j,ks)/xmass(npoint(j),ks)) th_xmass1(ks)/xmass(th_npoint,ks)) else xmassfract=1. endif end do if (ipin.eq.0 .and. xmassfract.lt.0.000001) then ! terminate all particles carrying less mass ! itra1(jj)=-999999999 th_itra1=-999999999 endif ! Sabine Eckhardt, June 2008 ! don't create depofield for backward runs if (DRYDEP.AND.(ldirect.eq.1)) then ! call drydepokernel(nclass(jj),drydeposit,real(xtra1(jj)), & call drydepokernel(th_nclass,drydeposit,real(th_xtra1), & ! real(ytra1(jj)),nage,kp) real(th_ytra1),itage,nage,kp) if (nested_output.eq.1) call drydepokernel_nest( & ! nclass(jj),drydeposit,real(xtra1(jj)),real(ytra1(jj)), & th_nclass,drydeposit,real(th_xtra1),real(th_ytra1), & itage,nage,kp) endif ! Terminate trajectories that are older than maximum allowed age !*************************************************************** ! if (abs(itra1(jj)-itramem(jj)).ge.lage(nageclass)) then if (abs(th_itra1-th_itramem).ge.lage(nageclass)) then if (linit_cond.ge.1) & call initial_cond_calc(jtime+lsynctime,jj) ! itra1(jj)=-999999999 th_itra1=-999999999 endif endif !! print*,xtra1(j),th_xtra1,OMP_GET_THREAD_NUM() !$OMP CRITICAL if (ntasks.eq.1) then npoint(jj)=th_npoint idt(jj)=th_idt uap(jj)=th_uap ucp(jj)=th_ucp uzp(jj)=th_uzp us(jj)=th_us vs(jj)=th_vs ws(jj)=th_ws xtra1(jj)=th_xtra1 ytra1(jj)=th_ytra1 ztra1(jj)=th_ztra1 cbt(jj)=th_cbt do ks=1,nspec xmass1(jj,ks)=th_xmass1(ks) ! drydep1(jj,ks)=drydeposit(ks) enddo ! itramem(jj)=th_itramem itra1(jj)=th_itra1 else mpi_npoint(jj)=th_npoint mpi_idt(jj)=th_idt mpi_uap(jj)=th_uap mpi_ucp(jj)=th_ucp mpi_uzp(jj)=th_uzp mpi_us(jj)=th_us mpi_vs(jj)=th_vs mpi_ws(jj)=th_ws mpi_xtra1(jj)=th_xtra1 mpi_ytra1(jj)=th_ytra1 mpi_ztra1(jj)=th_ztra1 ! mpi_nclass(jj)=th_nclass mpi_cbt(jj)=th_cbt do ks=1,nspec mpi_xmass1(jj,ks)=th_xmass1(ks) ! mpi_drydep1(jj,ks)=drydeposit(ks) enddo ! mpi_itramem(jj)=th_itramem mpi_itra1(jj)=th_itra1 endif ! if (jj.lt.5) print*,'coord cont',th_itra1 !$OMP END CRITICAL endif end do !loop over particles !$OMP END DO endif !!!$OMP !!!$OMP FLUSH(npoint,idt,uap,ucp,uzp,us,vs,ws,xtra1,ytra1,ztra1,cbt,xmass1,itra1) !!!$OMP FLUSH ! call itime(now) ! ttime=now(1)*3600+now(2)*60+now(3)-ttime call cpu_time(finish) ! print*,'time',ttime,cpttra,myid,OMP_GET_THREAD_NUM() if (option_verbose.eq.1) then print*,'time',finish-start,cpttra,myid,ompid endif !$OMP END PARALLEL !JB ! the openmp is done. the output above gives how long it takes to finish the loop over the particles. good benchmark ! here we use mpi to use the mpi_* arrays to update the big ones in the master ! thread ! call MPI_REDUCE (mypi ,pi ,1, MPI_DOUBLE_PRECISION , MPI_SUM ,0, & ! MPI_COMM_WORLD , ierr ) ! print*,'after loop',myid,chunksize if (chunksize.gt.0 .and. ntasks.gt.1) then !JB ! I need a barrier so each node is a the same place ! I am going to send the new results to the master thread now. call MPI_BARRIER(MPI_COMM_WORLD,ierr) if (myid.eq.0) then ! call itime(now) ! ttime=now(1)*3600+now(2)*60+now(3) call cpu_time(start) endif call sendint_mpi(1,numpart,chunksize,1) call sendint_mpi(2,numpart,chunksize,1) call sendint_mpi(3,numpart,chunksize,1) call sendreal_mpi(4,numpart,chunksize,1) call sendreal_mpi(5,numpart,chunksize,1) call sendreal_mpi(6,numpart,chunksize,1) call sendreal_mpi(7,numpart,chunksize,1) call sendreal_mpi(8,numpart,chunksize,1) call sendreal_mpi(9,numpart,chunksize,1) call sendreal_mpi(10,numpart,chunksize,1) call senddouble_mpi(11,numpart,chunksize,1) call senddouble_mpi(12,numpart,chunksize,1) call sendint2_mpi(13,numpart,chunksize,1) call sendint_mpi(14,numpart,chunksize,1) ! call sendint_mpi(15,nclass(1:numpart),numpart,chunksize2,1) call sendreal2d_mpi(99,numpart,nspec,chunksize,1) ! if (DRYDEP) call sendreal2d_mpi(21,numpart,nspec,chunksize,1) if (DRYDEP.and.(ldirect.eq.1)) call senddrydep_mpi(numxgrid*numygrid) if (DRYDEP.and.(ldirect.eq.1).and.nested_output.eq.1) & call senddrydep_nest_mpi(numxgridn*numygridn) ! print*,'size of vector',chunksize,chunksize2,myid,numpart ! ! if (myid.eq.0) then ! print*,'itra1',numpart,itra1(numpart-2:numpart+2) ! print*,'itra2',chunksize2,mpi_itra1(chunksize2-2:chunksize2+1) ! else ! print*,'itra3',mpi_itra1(chunksize-2:chunksize+1) ! endif if (myid.eq.0) then ! call itime(now) ! ttime=now(1)*3600+now(2)*60+now(3)-ttime call cpu_time(finish) ! print*,'receiving time',ttime if (option_verbose.eq.1) then print*,'receiving time',finish-start endif endif endif !finish transfering between nodes ! update the drydepo if (DRYDEP.and.ldirect.gt.0) then do ks=1,nspec do kp=1,maxpointspec_act do nage=1,nageclass do jy=0,numygrid-1 do ix=0,numxgrid-1 do l=1,nclassunc drygridunc(ix,jy,ks,kp,l,nage)=drygridunc(ix,jy,ks,kp,l,nage) & +drygridunc2(ix,jy,ks,kp,l,nage) end do end do end do if (nested_output.eq.1) then do jy=0,numygridn-1 do ix=0,numxgridn-1 do l=1,nclassunc drygriduncn(ix,jy,ks,kp,l,nage)=drygriduncn(ix,jy,ks,kp,l,nage) & +drygriduncn2(ix,jy,ks,kp,l,nage) end do end do end do endif end do end do end do endif ! if (DRYDEP.AND.(ldirect.eq.1)) then !! call drydepokernel(nclass(jj),drydeposit,real(xtra1(jj)), & ! call drydepokernel(th_nclass,drydeposit,real(th_xtra1), & !! real(ytra1(jj)),nage,kp) ! real(th_ytra1),itage,nage,kp) ! if (nested_output.eq.1) call drydepokernel_nest( & !! nclass(jj),drydeposit,real(xtra1(jj)),real(ytra1(jj)), & ! th_nclass,drydeposit,real(th_xtra1),real(th_ytra1), & ! itage,nage,kp) ! endif end do !loop over time ! Complete the calculation of initial conditions for particles not yet terminated !***************************************************************************** call MPI_BARRIER(MPI_COMM_WORLD,ierr) do j=1,numpart if (linit_cond.ge.1) call initial_cond_calc(jtime,j) end do if (ipout.eq.2) call partoutput(jtime) ! dump particle positions if (linit_cond.ge.1) call initial_cond_output(jtime) ! dump initial cond. field close(104) ! De-allocate memory and end !*************************** if (iflux.eq.1) then deallocate(flux) endif if (ntasks.gt.1) then ! print*,'before dealloc',myid deallocate(mpi_npoint,mpi_idt,mpi_itra1) ! print*,'after dealloc, part1',myid deallocate(mpi_uap,mpi_ucp,mpi_uzp) ! print*,'after dealloc, part12',myid deallocate(mpi_us,mpi_vs,mpi_ws,mpi_ztra1) ! print*,'after dealloc, part2',myid deallocate(mpi_xtra1,mpi_ytra1) ! print*,'after dealloc, part3',myid deallocate(mpi_cbt) ! print*,'after dealloc, part4',myid ! deallocate(mpi_xmass1,mpi_drydep1,mpi_nclass) deallocate(mpi_xmass1,mpi_nclass) ! print*,'after dealloc',myid deallocate(dummyi2) deallocate(dummyr2) deallocate(dummyi22) deallocate(dummyr22) endif if (OHREA.eqv..TRUE.) then deallocate(OH_field,OH_field_height) endif deallocate(gridunc) deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass) deallocate(ireleasestart,ireleaseend,npart,kindz) ! deallocate(xmasssave) if (myid.eq.0) then if (nested_output.eq.1) then deallocate(orooutn, arean, volumen) if (ldirect.gt.0) then deallocate(griduncn,drygriduncn,wetgriduncn,drygriduncn2) endif endif if (ldirect.gt.0) then if (allocated(drygridunc)) deallocate(drygridunc) if (allocated(wetgridunc)) deallocate(wetgridunc) if (allocated(drygridunc2)) deallocate(drygridunc2) if (allocated(drygriduncn2)) deallocate(drygriduncn2) endif deallocate(outheight,outheighthalf) deallocate(oroout, area, volume) endif end subroutine timemanager_mpi