Changeset 20963b1 in flexpart.git


Ignore:
Timestamp:
Apr 13, 2018, 2:33:40 PM (11 months ago)
Author:
Espen Sollum ATMOS <eso@…>
Branches:
dev, release-10, univie
Children:
0ca4976, f5b967d, 0eda008
Parents:
3f149cc
Message:

Backwards deposition for the MPI version

Location:
src
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • src/FLEXPART_MPI.f90

    ra9cf4b1 r20963b1  
    6767  integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
    6868  integer :: detectformat
    69 
     69  integer(selected_int_kind(16)), dimension(maxspec) :: tot_b=0, &
     70       & tot_i=0
    7071
    7172
     
    206207
    207208  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    208     print *,'ECMWF metdata detected'
     209    if (lroot) print *,'ECMWF metdata detected'
    209210  elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
    210     print *,'NCEP metdata detected'
     211    if (lroot) print *,'NCEP metdata detected'
    211212  else
    212     print *,'Unknown metdata format'
     213    if (lroot) print *,'Unknown metdata format'
    213214    stop
    214215  endif
     
    483484
    484485! NIK 16.02.2005
    485   if (lroot) then
    486     call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
     486  if (mp_partgroup_pid.ge.0) then ! Skip for readwind process
     487    call MPI_Reduce(tot_blc_count, tot_b, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
    487488         & mp_comm_used, mp_ierr)
    488     call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
     489    call MPI_Reduce(tot_inc_count, tot_i, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
    489490         & mp_comm_used, mp_ierr)
    490   else
    491     if (mp_partgroup_pid.ge.0) then ! Skip for readwind process
    492       call MPI_Reduce(tot_blc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
    493            & mp_comm_used, mp_ierr)
    494       call MPI_Reduce(tot_inc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
    495            & mp_comm_used, mp_ierr)
    496     end if
    497491  end if
     492 
    498493
    499494  if (lroot) then
     
    502497      write(*,*) 'Scavenging statistics for species ', species(i), ':'
    503498      write(*,*) 'Total number of occurences of below-cloud scavenging', &
    504            & tot_blc_count(i)
     499           & tot_b(i)
     500!           & tot_blc_count(i)
    505501      write(*,*) 'Total number of occurences of in-cloud    scavenging', &
    506            & tot_inc_count(i)
     502           & tot_i(i)
     503!           & tot_inc_count(i)
    507504      write(*,*) '**********************************************'
    508505    end do
  • src/concoutput.f90

    r08a38b5 r20963b1  
    140140  close(unitdates) 
    141141
     142  ! Overwrite existing dates file on first call, later append to it
     143  ! This fixes a bug where the dates file kept growing across multiple runs
    142144  IF (init) THEN
    143145    file_stat='OLD'
     
    313315                   drygridsigma(ix,jy)* &
    314316                   sqrt(real(nclassunc))
    315 125           drygridsigmatotal=drygridsigmatotal+ &
     317              drygridsigmatotal=drygridsigmatotal+ &
    316318                   drygridsigma(ix,jy)
    317319            endif
  • src/concoutput_mpi.f90

    rb1e0742 r20963b1  
    104104  real,parameter :: weightair=28.97
    105105  logical :: sp_zer
    106   LOGICAL,save :: init=.true.
     106  logical,save :: init=.true.
    107107  character :: adate*8,atime*6
    108108  character(len=3) :: anspec
     
    128128! This fixes a bug where the dates file kept growing across multiple runs
    129129
    130 ! If 'dates' file exists, make a backup
     130! If 'dates' file exists in output directory, make a backup
    131131  inquire(file=path(2)(1:length(2))//'dates', exist=ldates_file)
    132132  if (ldates_file.and.init) then
     
    257257
    258258    write(anspec,'(i3.3)') ks
    259     if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
    260       if (ldirect.eq.1) then
    261         open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
     259
     260    if (DRYBKDEP.or.WETBKDEP) then !scavdep output
     261      if (DRYBKDEP) &
     262      open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// &
     263           atime//'_'//anspec,form='unformatted')
     264      if (WETBKDEP) &
     265      open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_'//adate// &
     266           atime//'_'//anspec,form='unformatted')
     267      write(unitoutgrid) itime
     268    else
     269      if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
     270        if (ldirect.eq.1) then
     271          open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
     272               atime//'_'//anspec,form='unformatted')
     273        else
     274          open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
     275               atime//'_'//anspec,form='unformatted')
     276        endif
     277        write(unitoutgrid) itime
     278      endif
     279
     280      if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
     281        open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
    262282             atime//'_'//anspec,form='unformatted')
    263       else
    264         open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
    265              atime//'_'//anspec,form='unformatted')
     283
     284        write(unitoutgridppt) itime
    266285      endif
    267       write(unitoutgrid) itime
    268     endif
    269 
    270     if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
    271       open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
    272            atime//'_'//anspec,form='unformatted')
    273 
    274       write(unitoutgridppt) itime
    275     endif
     286    endif ! if deposition output
    276287
    277288    do kp=1,maxpointspec_act
     
    353364! Concentration output
    354365!*********************
    355         if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
     366        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
    356367
    357368! Wet deposition
     
    448459                  endif
    449460                  sp_count_r=sp_count_r+1
     461                  if (lparticlecountoutput) then
     462                    sparse_dump_r(sp_count_r)= &
     463                         sp_fact* &
     464                         grid(ix,jy,kz)
     465                  else
    450466                  sparse_dump_r(sp_count_r)= &
    451467                       sp_fact* &
    452468                       grid(ix,jy,kz)* &
    453469                       factor3d(ix,jy,kz)/tot_mu(ks,kp)
     470                  end if
     471
     472
    454473!                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
    455474!    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
     
    637656!*************************
    638657
    639   do ks=1,nspec
    640     do kp=1,maxpointspec_act
    641       do i=1,numreceptor
    642         creceptor(i,ks)=0.
    643       end do
    644       do jy=0,numygrid-1
    645         do ix=0,numxgrid-1
    646           do l=1,nclassunc
    647             do nage=1,nageclass
    648               do kz=1,numzgrid
    649                 gridunc(ix,jy,kz,ks,kp,l,nage)=0.
    650               end do
    651             end do
    652           end do
    653         end do
    654       end do
    655     end do
    656   end do
     658  creceptor(:,:)=0.
     659  gridunc(:,:,:,:,:,:,:)=0.
    657660
    658661  if (mp_measure_time) call mpif_mtime('rootonly',1)
  • src/getfields_mpi.f90

    r6ecb30a r20963b1  
    23323340  indmin=indj
    234234
     235   if (WETBKDEP.and.(lmpreader.or.(.not.lmp_use_reader.and.lroot))) then
     236        call writeprecip(itime,memind(1))
     237      endif
     238
    235239  else
    236240
     
    28729160  indmin=indj
    288292
    289     mind3=memstat
     293!   if (WETBKDEP.and.lroot) then
     294    if (WETBKDEP.and.(lmpreader.or.(.not.lmp_use_reader.and.lroot))) then
     295      call writeprecip(itime,memind(1))
     296    endif
    290297
    291298  endif
  • src/readcommand.f90

    ra9cf4b1 r20963b1  
    332332     case (3)  ! 3 .. wet deposition in outputfield
    333333        ind_rel = 3
    334          write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE    #### '
    335          write(*,*) ' #### Releaseheight is forced to 0 - 20km      #### '
    336          write(*,*) ' #### Release is performed above ground lev    #### '
     334        if (lroot) then
     335          write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE    #### '
     336          write(*,*) ' #### Releaseheight is forced to 0 - 20km      #### '
     337          write(*,*) ' #### Release is performed above ground lev    #### '
     338        end if
    337339         WETBKDEP=.true.
    338340         allocate(xscav_frac1(maxpart,maxspec))
    339341     case (4)  ! 4 .. dry deposition in outputfield
    340342         ind_rel = 4
    341          write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE    #### '
    342          write(*,*) ' #### Releaseheight is forced to 0 - 2*href    #### '
    343          write(*,*) ' #### Release is performed above ground lev    #### '
     343         if (lroot) then
     344           write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE    #### '
     345           write(*,*) ' #### Releaseheight is forced to 0 - 2*href    #### '
     346           write(*,*) ' #### Release is performed above ground lev    #### '
     347         end if
    344348         DRYBKDEP=.true.
    345349         allocate(xscav_frac1(maxpart,maxspec))
  • src/readwind_ecmwf_mpi.f90

    r61e07ba r20963b1  
    103103
    104104  integer :: isec1(56),isec2(22+nxmax+nymax)
    105   real(kind=4) :: zsec4(jpunp)
    106   real(kind=4) :: xaux,yaux
    107   real(kind=8) :: xauxin,yauxin
    108   real,parameter :: eps=1.e-4
    109   real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
    110   real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
     105  real(sp) :: zsec4(jpunp)
     106  real(sp) :: xaux,yaux
     107  real(dp) :: xauxin,yauxin
     108  real(sp),parameter :: eps=1.e-4
     109  real(sp) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
     110  real(sp) :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
    111111
    112112  logical :: hflswitch,strswitch !,readcloud
     
    123123!ZHG test the grib fields that have lcwc without using them
    124124!  readcloud=.false.
    125 !hg end
     125
    126126  levdiff2=nlev_ec-nwz+1
    127127  iumax=0
  • src/releaseparticles_mpi.f90

    r16b61a5 r20963b1  
    216216              xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
    217217                   *timecorrect(k)/average_timecorrect
    218 !             write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k
     218              if (DRYBKDEP.or.WETBKDEP) then ! if there is no scavenging in wetdepo it will be set to 0
     219!              if ( henry(k).gt.0 .or. &
     220!                   crain_aero(k).gt.0. .or. csnow_aero(k).gt.0. .or. &
     221!                   ccn_aero(k).gt.0. .or. in_aero(k).gt.0. )  then
     222                xscav_frac1(ipart,k)=-1.
     223               endif
    219224! Assign certain properties to particle
    220225!**************************************
     
    366371!Af ind_rel is defined in readcommand.f
    367372
    368             if (ind_rel .eq. 1) then
     373            if ((ind_rel .eq. 1).or.(ind_rel .eq. 3).or.(ind_rel .eq. 4)) then
    369374
    370375! Interpolate the air density
  • src/timemanager_mpi.f90

    ra9cf4b1 r20963b1  
    115115  integer :: ip
    116116  integer :: loutnext,loutstart,loutend
    117   integer :: ix,jy,ldeltat,itage,nage
     117  integer :: ix,jy,ldeltat,itage,nage,idummy
    118118  integer :: i_nan=0,ii_nan,total_nan_intl=0  !added by mc to check instability in CBL scheme
    119119  integer :: numpart_tot_mpi ! for summing particles on all processes
    120   real :: outnum,weight,prob(maxspec)
    121   real :: decfact
     120  real :: outnum,weight,prob(maxspec), prob_rec(maxspec), decfact,wetscav
    122121
    123122  real(sp) :: gridtotalunc
     
    125124       & drydeposit(maxspec)=0_dep_prec
    126125  real :: xold,yold,zold,xmassfract
     126  real :: grfraction(3)
    127127  real, parameter :: e_inv = 1.0/exp(1.0)
    128128
     
    159159!CGZ-lifetime: set lifetime to 0
    160160
     161  if (.not.lusekerneloutput) write(*,*) 'Not using the kernel'
     162  if (turboff) write(*,*) 'Turbulence switched off'
    161163
    162164
     
    708710        zold=ztra1(j)
    709711
     712   
     713  ! RECEPTOR: dry/wet depovel
     714  !****************************
     715  ! Before the particle is moved
     716  ! the calculation of the scavenged mass shall only be done once after release
     717  ! xscav_frac1 was initialised with a negative value
     718
     719      if  (DRYBKDEP) then
     720       do ks=1,nspec
     721         if  ((xscav_frac1(j,ks).lt.0)) then
     722            call get_vdep_prob(itime,xtra1(j),ytra1(j),ztra1(j),prob_rec)
     723            if (DRYDEPSPEC(ks)) then        ! dry deposition
     724               xscav_frac1(j,ks)=prob_rec(ks)
     725             else
     726                xmass1(j,ks)=0.
     727                xscav_frac1(j,ks)=0.
     728             endif
     729         endif
     730        enddo
     731       endif
     732
     733       if (WETBKDEP) then
     734       do ks=1,nspec
     735         if  ((xscav_frac1(j,ks).lt.0)) then
     736            call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy,wetscav)
     737            if (wetscav.gt.0) then
     738                xscav_frac1(j,ks)=wetscav* &
     739                       (zpoint2(npoint(j))-zpoint1(npoint(j)))*grfraction(1)
     740            else
     741                xmass1(j,ks)=0.
     742                xscav_frac1(j,ks)=0.
     743            endif
     744         endif
     745        enddo
     746       endif
     747
    710748! Integrate Lagevin equation for lsynctime seconds
    711749!*************************************************
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG