diff --git a/src/chemistry_mod.f90 b/src/chemistry_mod.f90
index 533c5daa74eba42138fc4310ecabe01dce8331dd..19abd5f283707a0c3066716ce6f2500710a5072e 100644
--- a/src/chemistry_mod.f90
+++ b/src/chemistry_mod.f90
@@ -485,6 +485,7 @@ module chemistry_mod
use particle_mod, only: count, part, mass
use point_mod, only: xlon0, ylat0, dx, dy
use windfields_mod, only: xresoln, yresoln, xln, yln, xrn, yrn, height, tt, nz
+ use omp_lib
implicit none
@@ -492,7 +493,7 @@ module chemistry_mod
integer :: ii,i,j,ks,ix,jy,kz,jrx,jry,nr
integer :: ngrid,interp_time,n,indz
integer :: jjjjmmdd,hhmmss,mm,hh,m1,m2
- integer :: clx,cly,clz,clzm
+ integer :: clx,cly,clz,clzm,ithread
real, dimension(nzCL) :: altCLtop
real, dimension(2) :: cl_tmp
real :: xlon,ylat
@@ -504,6 +505,9 @@ module chemistry_mod
real, parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
real(kind=dp) :: jul
real :: lonjrx,latjry
+#ifdef _OPENMP
+ real(kind=dp), allocatable, dimension(:,:,:) :: chem_loss_tmp
+#endif
! use middle of synchronisation time step
interp_time=nint(itime+0.5*lsynctime)
@@ -514,14 +518,24 @@ module chemistry_mod
! initialization
chem_loss(:,:)=0d0
+#ifdef _OPENMP
+ allocate( chem_loss_tmp(nreagent,nspec,numthreads) )
+ chem_loss_tmp(:,:,:) = 0d0
+#endif
! Loop over particles
!*****************************************
!$OMP PARALLEL &
!$OMP PRIVATE(ii,jpart,ngrid,j,xtn,ytn,ix,jy, &
!$OMP xlon,ylat,clx,cly,clz,clzm,kz,altCLtop,dz1,dz2,dzz,nr,i, &
-!$OMP cl_cur,indz,temp,ks,clrate,restmass,clreacted) &
-!$OMP REDUCTION(+:chem_loss)
+!$OMP cl_cur,indz,temp,ks,clrate,restmass,clreacted,ithread)
+! !$OMP REDUCTION(+:chem_loss)
+
+#ifdef _OPENMP
+ ithread = OMP_GET_THREAD_NUM()+1 ! Starts with 1
+#else
+ ithread = 1
+#endif
!$OMP DO
do ii=1,count%alive
@@ -562,8 +576,8 @@ module chemistry_mod
endif
! get position in the chem field
! assumes chem field dimensions given as grid midpoints
- clx=int((xlon-(lonCL(1)-0.5*dxCL))/dxCL)+1
- cly=int((ylat-(latCL(1)-0.5*dyCL))/dyCL)+1
+ clx=min(nxCL, int((xlon-(lonCL(1)-0.5*dxCL))/dxCL)+1)
+ cly=min(nyCL, int((ylat-(latCL(1)-0.5*dyCL))/dyCL)+1)
! get the level of the chem field for the particle
! z is the z-coord of the trajectory above model orography in metres
@@ -626,7 +640,11 @@ module chemistry_mod
clreacted=mass(jpart,ks)
mass(jpart,ks)=0.
endif
+#ifdef _OPENMP
+ chem_loss_tmp(nr,ks,ithread)=chem_loss_tmp(nr,ks,ithread)+real(clreacted,kind=dp)
+#else
chem_loss(nr,ks)=chem_loss(nr,ks)+real(clreacted,kind=dp)
+#endif
endif
end do ! nspec
endif
@@ -637,6 +655,13 @@ module chemistry_mod
!$OMP END DO
!$OMP END PARALLEL
+
+#ifdef _OPENMP
+ do ithread=1,numthreads
+ chem_loss(:,:) = chem_loss(:,:)+chem_loss_tmp(:,:,ithread)
+ end do
+ deallocate( chem_loss_tmp )
+#endif
end subroutine chemreaction
diff --git a/src/emissions_mod.f90 b/src/emissions_mod.f90
index 43fb5c79af4d3ec0a96faf3da5024d4d2a81cd3a..e0414708c42f3532dc6bd3f719848fc23ead5932 100644
--- a/src/emissions_mod.f90
+++ b/src/emissions_mod.f90
@@ -52,12 +52,13 @@ module emissions_mod
!*****************************************************************************
use windfields_mod, only: hmix
+ use omp_lib
implicit none
integer :: itime
real :: xlon,ylat
- integer :: ix, jy, ixp, jyp, ii, ks, em_ix, em_jy
+ integer :: ix, jy, ixp, jyp, ii, ks, em_ix, em_jy, ithread
real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
real :: em_dt1, em_dt2, em_dtt
real, dimension(2) :: hm
@@ -68,6 +69,10 @@ module emissions_mod
character(len=20) :: frmt
real :: em_frac
real, allocatable, dimension(:,:,:) :: em_neg
+#ifdef _OPENMP
+ real, allocatable, dimension(:,:) :: tot_em_up_tmp
+ real, allocatable, dimension(:,:,:,:) :: mass_field_tmp, em_neg_tmp
+#endif
real :: f_m
logical :: lexist
integer, parameter :: unittest=120
@@ -86,29 +91,49 @@ module emissions_mod
dt2=float(memtime(2)-itime)
dtt=1./(dt1+dt2)
- tot_em_up(:) = 0.
-
! estimate mass in PBL from particle positions
!**********************************************
- mass_field(:,:,:) = 0.
allocate( em_neg(nspec-1,nxem,nyem) )
+#ifdef _OPENMP
+ allocate( mass_field_tmp(nspec,nxem,nyem,numthreads))
+ allocate( em_neg_tmp(nspec-1,nxem,nyem,numthreads))
+ allocate( tot_em_up_tmp(nspec,numthreads) )
+#endif
+
+!$OMP PARALLEL
+!$OMP WORKSHARE
+ mass_field_tmp(:,:,:,:) = 0.
+ em_neg_tmp(:,:,:,:) = 0.
+ tot_em_up_tmp(:,:) = 0.
+ tot_em_up(:) = 0.
+ mass_field(:,:,:) = 0.
em_neg(:,:,:) = 0.
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
!$OMP PARALLEL &
!$OMP PRIVATE(ii,xlon,ylat,em_ix,em_jy,ix,jy,ixp,jyp,ddx,ddy, &
-!$OMP rddx,rddy,p1,p2,p3,p4,mm,hm,hmixi,ks,em_cur,tmp) &
-!$OMP SHARED(em_cond)
+!$OMP rddx,rddy,p1,p2,p3,p4,mm,hm,hmixi,ks,ithread) &
+!$OMP SHARED(em_cond,mass_field_tmp)
+
+#ifdef _OPENMP
+ ithread = OMP_GET_THREAD_NUM()+1 ! Starts with 1
+#else
+ ithread = 1
+#endif
+
+!$OMP DO
-!$OMP DO REDUCTION(+:mass_field)
+! REDUCTION(+:mass_field)
do ii=1,count%alive ! loop over all particles
xlon=xlon0+part(ii)%xlon*dx
ylat=ylat0+part(ii)%ylat*dy
! assume emission dimensions given as grid midpoints
- em_ix=int((xlon-(lonem(1)-0.5*dxem))/dxem)+1
- em_jy=int((ylat-(latem(1)-0.5*dyem))/dyem)+1
+ em_ix=min(nxem, int((xlon-(lonem(1)-0.5*dxem))/dxem)+1)
+ em_jy=min(nyem, int((ylat-(latem(1)-0.5*dyem))/dyem)+1)
!! testing
if (ii.lt.20) print*, 'lonem, lon, latem, lat = ',lonem(em_ix),xlon,latem(em_jy),ylat
@@ -141,19 +166,42 @@ module emissions_mod
em_cond(ii) = part(ii)%z.le.hmixi
if (em_cond(ii)) then
- mass_field(1:nspec,em_ix, em_jy)=mass_field(1:nspec, em_ix, em_jy) + &
- mass(ii,1:nspec)
+#ifdef _OPENMP
+ mass_field_tmp(1:nspec,em_ix,em_jy,ithread)= &
+ mass_field_tmp(1:nspec,em_ix,em_jy,ithread) + &
+ mass(ii,1:nspec)
+#else
+ mass_field(1:nspec,em_ix,em_jy)=mass_field(1:nspec,em_ix,em_jy) + &
+ mass(ii,1:nspec)
+#endif
endif
end do ! end of particle loop
!$OMP END DO
+!$OMP END PARALLEL
+ ! Manual reduction of mass_field
+ do ithread=1,numthreads
+ mass_field(:,:,:) = mass_field(:,:,:)+mass_field_tmp(:,:,:,ithread)
+ end do
f_m = exp(-1.*real(lsynctime)/tau_ipm)
! estimate emissions for each particle
!**************************************
-!$OMP DO REDUCTION(+:tot_em_up,em_neg)
+!$OMP PARALLEL &
+!$OMP PRIVATE(ii,xlon,ylat,em_ix,em_jy, &
+!$OMP ks,em_cur,tmp,ithread) &
+!$OMP SHARED(em_cond,tot_em_up_tmp,em_neg_tmp)
+
+#ifdef _OPENMP
+ ithread = OMP_GET_THREAD_NUM()+1 ! Starts with 1
+#else
+ ithread = 1
+#endif
+
+!$OMP DO
+! REDUCTION(+:tot_em_up,em_neg)
do ii=1,count%alive ! loop over particles
if (.not.em_cond(ii)) cycle ! skip particles not in PBL
@@ -161,8 +209,8 @@ module emissions_mod
xlon=xlon0+part(ii)%xlon*dx
ylat=ylat0+part(ii)%ylat*dy
! assume emission dimensions given as grid midpoints
- em_ix=int((xlon-(lonem(1)-0.5*dxem))/dxem)+1
- em_jy=int((ylat-(latem(1)-0.5*dyem))/dyem)+1
+ em_ix=min(nxem, int((xlon-(lonem(1)-0.5*dxem))/dxem)+1)
+ em_jy=min(nyem, int((ylat-(latem(1)-0.5*dyem))/dyem)+1)
! loop over species
! skip species 1 as it is always air tracer with no emission
@@ -185,18 +233,23 @@ module emissions_mod
if (-1.*tmp.gt.mass(ii,ks)) then
tmp = tmp + mass(ii,ks)
! subtract mass from atmosphere by setting it to zero below
- tot_em_up(ks) = tot_em_up(ks) - real(mass(ii,ks),kind=dp)
+ tot_em_up_tmp(ks,ithread) = tot_em_up_tmp(ks,ithread) - real(mass(ii,ks),kind=dp)
mass(ii,ks) = 0.
- ! add remaining uptake to em_neg
+ ! add remaining uptake to em_neg
+#ifdef _OPENMP
+ em_neg_tmp(ks-1,em_ix,em_jy,ithread) = em_neg_tmp(ks-1,em_ix,em_jy,ithread) + &
+ tmp/mass(ii,1)*mass_field(1,em_ix,em_jy)
+#else
em_neg(ks-1, em_ix, em_jy) = em_neg(ks-1, em_ix, em_jy) + &
tmp/mass(ii,1)*mass_field(1,em_ix,em_jy)
+#endif
else
mass(ii,ks)=mass(ii,ks)+tmp
- tot_em_up(ks) = tot_em_up(ks) + real(tmp,kind=dp)
+ tot_em_up_tmp(ks,ithread) = tot_em_up_tmp(ks,ithread) + real(tmp,kind=dp)
endif
else
mass(ii,ks)=mass(ii,ks)+tmp
- tot_em_up(ks) = tot_em_up(ks) + real(tmp,kind=dp)
+ tot_em_up_tmp(ks,ithread) = tot_em_up_tmp(ks,ithread) + real(tmp,kind=dp)
endif ! negative emissions
end do ! nspec
@@ -226,6 +279,16 @@ module emissions_mod
!$OMP END DO
!$OMP END PARALLEL
+#ifdef _OPENMP
+ do ithread=1,numthreads
+ em_neg(:,:,:) = em_neg(:,:,:)+em_neg_tmp(:,:,:,ithread)
+ tot_em_up(:) = tot_em_up(:) + tot_em_up_tmp(:,ithread)
+ end do
+ deallocate( mass_field_tmp,em_neg_tmp,tot_em_up_tmp )
+#else
+ tot_em_up(:) = tot_em_up_tmp(:,1)
+#endif
+
! update for negative emissions
em_res(:,:,:) = em_res(:,:,:)+em_neg(:,:,:)
deallocate(em_neg)
diff --git a/src/netcdf_output_mod.f90 b/src/netcdf_output_mod.f90
index 8dd9dffbbfbf2478b9c276b1f9a67301428800fa..e3d36b663f46681931db049f5d8f49adc705bd28 100644
--- a/src/netcdf_output_mod.f90
+++ b/src/netcdf_output_mod.f90
@@ -2055,7 +2055,7 @@ subroutine writeheader_partoutput_vars(np,ncid,totpart,timeDimID,partDimID,latDi
call write_to_file(ncid,trim(partopt(np)%short_name),nf90_float,(/ timeDimID,partDimID /), &
varid,(/ 1,totpart /),'m/s',.false.,'settling_velocity_average','settling velocity averaged')
case ('MA') ! Mass
- if (mdomainfill.ge.1) then
+ if ((mdomainfill.ge.1).and.(nspec.eq.1)) then
call nf90_err(nf90_def_var(ncid=ncid, name=trim(partopt(np)%short_name), xtype=nf90_float, &
dimids=1, varid=varid))
call nf90_err(nf90_put_att(ncid, varid, 'units', 'kg'))
@@ -2073,7 +2073,7 @@ subroutine writeheader_partoutput_vars(np,ncid,totpart,timeDimID,partDimID,latDi
end do
endif
case ('ma') ! Mass averaged
- if (mdomainfill.ge.1) then
+ if ((mdomainfill.ge.1).and.(nspec.eq.1)) then
call nf90_err(nf90_def_var(ncid=ncid, name=trim(partopt(np)%short_name), xtype=nf90_float, dimids=1, varid=varid))
call nf90_err(nf90_put_att(ncid, varid, 'units', 'kg'))
call nf90_err(nf90_put_att(ncid, varid, '_FillValue', fillval))
@@ -2364,7 +2364,7 @@ subroutine partoutput_netcdf(itime,field,np,imass,ncid)
endif
else if (partopt(np)%name.eq.'MA') then
- if ((mdomainfill.ge.1).and.(imass.eq.1)) then
+ if ((mdomainfill.ge.1).and.(imass.eq.1).and.(nspec.eq.1)) then
if (mass_written.eqv..false.) then
call nf90_err(nf90_inq_varid(ncid=ncid,name=trim(partopt(np)%short_name),varid=tempIDend))
call nf90_err(nf90_put_var(ncid=ncid,varid=tempIDend,values=field(1)))
@@ -2376,7 +2376,7 @@ subroutine partoutput_netcdf(itime,field,np,imass,ncid)
call nf90_err(nf90_put_var(ncid,tempIDend,field, (/ tpointer_part,1 /),(/ 1,count%allocated /)))
endif
else if (partopt(np)%name.eq.'ma') then
- if ((mdomainfill.ge.1).and.(imass.eq.1)) then
+ if ((mdomainfill.ge.1).and.(imass.eq.1).and.(nspec.eq.1)) then
if (mass_written.eqv..false.) then
call nf90_err(nf90_inq_varid(ncid=ncid,name=trim(partopt(np)%short_name),varid=tempIDend))
call nf90_err(nf90_put_var(ncid,tempIDend,field, (/ tpointer_part,1 /),(/ 1,count%allocated /)))
@@ -2477,26 +2477,26 @@ subroutine readpartpositions_netcdf(ibtime,ibdate)
! And give them the correct positions
! Longitude
- call nf90_err(nf90_inq_varid(ncid=ncidend,name='longitude',varid=tempIDend))
+ call nf90_err(nf90_inq_varid(ncid=ncidend,name='lon',varid=tempIDend))
call nf90_err(nf90_get_var(ncid=ncidend,varid=tempIDend,values=part(:)%xlon, &
start=(/ tlen, 1 /),count=(/ 1, plen /)))
part(:)%xlon=(part(:)%xlon-xlon0)/dx
! Latitude
- call nf90_err(nf90_inq_varid(ncid=ncidend,name='latitude',varid=tempIDend))
+ call nf90_err(nf90_inq_varid(ncid=ncidend,name='lat',varid=tempIDend))
call nf90_err(nf90_get_var(ncid=ncidend,varid=tempIDend,values=part(:)%ylat, &
start=(/ tlen, 1 /),count=(/ 1, plen /)))
part(:)%ylat=(part(:)%ylat-ylat0)/dx
! Height
- call nf90_err(nf90_inq_varid(ncid=ncidend,name='height',varid=tempIDend))
+ call nf90_err(nf90_inq_varid(ncid=ncidend,name='z',varid=tempIDend))
call nf90_err(nf90_get_var(ncid=ncidend,varid=tempIDend,values=part(:)%z, &
start=(/ tlen, 1 /),count=(/ 1, plen /)))
! Mass
! allocate(mass_temp(count%allocated), stat=stat)
! if (stat.ne.0) error stop "Could not allocate mass_temp"
- if (mdomainfill.eq.0) then
+ if ((mdomainfill.eq.0).or.(nspec.gt.1)) then
do j=1,nspec
write(anspec, '(i3.3)') j
- call nf90_err(nf90_inq_varid(ncid=ncidend,name='mass'//anspec,varid=tempIDend))
+ call nf90_err(nf90_inq_varid(ncid=ncidend,name='m'//anspec,varid=tempIDend))
call nf90_err(nf90_get_var(ncid=ncidend,varid=tempIDend,values=mass(:,j), &
start=(/ tlen, 1 /),count=(/ 1, plen /)))
! do i=1,count%allocated
diff --git a/src/timemanager_mod.f90 b/src/timemanager_mod.f90
index 72b77400ae2468e5add3b6f61c77f987a52a632e..0871a411abd55c83c7e86a6587dfa4dbdaee78c6 100644
--- a/src/timemanager_mod.f90
+++ b/src/timemanager_mod.f90
@@ -704,8 +704,17 @@ subroutine timemanager
#ifdef USE_NCF
if ((mdomainfill.eq.1).and.(llcmoutput)) then
+
+ if (.not.allocated(tot_mass)) then
+ allocate( tot_mass(nspec) )
+ tot_mass(:)=0.
+ endif
+
do ks=1,nspec
- tot_mass(ks)=sum(real(mass(1:count%alive,ks),kind=dp))
+ do i=1,count%alive
+ j=count%ialive(i)
+ tot_mass(ks)=tot_mass(ks)+real(mass(j,ks),kind=dp)
+ end do
end do
call totals_write(itime)
endif