*
L. Bakels 2022: This module contains most output related subroutines * *
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt ! SPDX-License-Identifier: GPL-3.0-or-later !***************************************************************************** ! * ! L. Bakels 2022: This module contains most output related subroutines * ! * !***************************************************************************** module output_mod use com_mod use par_mod use date_mod #ifdef USE_NCF use netcdf_output_mod #endif use binary_output_mod use txt_output_mod implicit none contains subroutine init_output(itime,filesize) implicit none integer, intent(in) :: itime real, intent(inout) :: filesize #ifdef USE_NCF real(kind=dp) :: & jul integer :: & jjjjmmdd,ihmmss,i #endif ! Writing header information to either binary or NetCDF format if (itime.eq.itime_init) then if (iout.ne.0) then ! No gridded output for iout=0 if (lnetcdfout.eq.1) then #ifdef USE_NCF call writeheader_netcdf(lnest=.false.) if (nested_output.eq.1) call writeheader_netcdf(lnest=.true.) #endif else if ((ipin.ne.1).or.(ipin.ne.4)) then ! Not necessary for restart call writeheader_bin !if (nested_output.eq.1) call writeheader_nest if ((nested_output.eq.1).and.(sfc_only.ne.1)) call writeheader_bin_nest if ((nested_output.eq.1).and.(sfc_only.eq.1)) call writeheader_bin_sfc_nest if ((nested_output.ne.1).and.(sfc_only.eq.1)) call writeheader_bin_sfc endif endif ! iout.ne.0 ! FLEXPART 9.2 ticket ?? write header in ASCII format if ((ipin.ne.1).or.(ipin.ne.4)) call writeheader_txt ! NetCDF only: Create file for storing initial particle positions. #ifdef USE_NCF if (ipout.ge.1) then if (itime_init.ne.0) then jul=bdate+real(itime,kind=dp)/86400._dp call caldate(jul,jjjjmmdd,ihmmss) endif ! if ((mdomainfill.eq.0).and.(ipin.le.1)) then ! if (itime_init.ne.0) then ! if (ldirect.eq.1) then ! call create_particles_initialoutput(ihmmss,jjjjmmdd,ibtime,ibdate) ! else ! call create_particles_initialoutput(ihmmss,jjjjmmdd,ietime,iedate) ! endif ! else if (ldirect.eq.1) then ! call create_particles_initialoutput(ibtime,ibdate,ibtime,ibdate) ! else ! call create_particles_initialoutput(ietime,iedate,ietime,iedate) ! endif ! endif ! Create header files for files that store the particle dump output if (itime_init.ne.0) then if (ldirect.eq.1) then call writeheader_partoutput(ihmmss,jjjjmmdd,ibtime,ibdate) else call writeheader_partoutput(ihmmss,jjjjmmdd,ietime,iedate) endif else if (ldirect.eq.1) then call writeheader_partoutput(ibtime,ibdate,ibtime,ibdate) else call writeheader_partoutput(ietime,iedate,ietime,iedate) endif endif #endif ! In case the particle output file is becoming larger than the maximum set ! in par_mod, create a new one while keeping track of the filesize. ! Also if a new restart file is created. else if ((mod(itime,ipoutfac*loutstep).eq.0).and.(ipout.ge.1)) then #ifdef USE_NCF if ((filesize.ge.maxfilesize).or. & ((loutrestart.ne.-1).and.(mod(itime,loutrestart).eq.0))) then jul=bdate+real(itime,kind=dp)/86400._dp call caldate(jul,jjjjmmdd,ihmmss) if (ldirect.eq.1) then call writeheader_partoutput(ihmmss,jjjjmmdd,ibtime,ibdate) else call writeheader_partoutput(ihmmss,jjjjmmdd,ietime,iedate) endif filesize = 0. endif do i=1,numpoint filesize = filesize + npart(i)*13.*4./1000000. end do #endif endif end subroutine init_output subroutine finalise_output(itime) ! Complete the calculation of initial conditions for particles not yet terminated use particle_mod implicit none integer, intent(in) :: itime integer :: i,j,ithread if (linit_cond.ge.1) then do i=1,count%alive j=count%ialive(i) call initcond_calc(itime,j,1) end do #ifdef _OPENMP do ithread=1,numthreads init_cond(:,:,:,:,:)=init_cond(:,:,:,:,:)+init_cond_omp(:,:,:,:,:,ithread) end do #endif endif if (ipout.eq.2) call output_particles(itime)!,active_per_rel) ! dump particle positions if (linit_cond.ge.1) then if(linversionout.eq.1) then call initcond_output_inv(itime) ! dump initial cond. field else call initcond_output(itime) ! dump initial cond. fielf endif endif end subroutine finalise_output subroutine output_particles(itime,initial_output) ! i !***************************************************************************** ! * ! Dump all particle positions * ! Author: A. Stohl * ! * ! 12 March 1999 * ! * ! Changes L. Bakels, 2021 * ! Output is chosen by the fields set in PARTOPTIONS * ! Binary output is no longer supported. If required, function can be * ! added below at "Put binary function here" * !***************************************************************************** ! * ! Variables: * ! * !***************************************************************************** use interpol_mod #ifdef ETA use coord_ecmwf_mod #endif use particle_mod #ifdef USE_NCF use netcdf use netcdf_output_mod, only: partoutput_netcdf,open_partoutput_file, & close_partoutput_file,partinitpointer1 use omp_lib, only: OMP_GET_THREAD_NUM #endif implicit none integer,intent(in) :: itime logical,optional,intent(in) :: initial_output logical :: init_out,lskip integer :: i,j,m,jjjjmmdd,ihmmss,np,ns,i_av real(kind=dp) :: jul real :: tmp(2) character :: adate*8,atime*6 real :: dummy(2) real :: masstemp(count%allocated,nspec),masstemp_av(count%allocated,nspec) real :: wetdepotemp(count%allocated,nspec),drydepotemp(count%allocated,nspec) real :: output(num_partopt, count%allocated) real :: cartxyz(3) logical :: cartxyz_comp #ifdef USE_NCF integer :: ncid,ncid_tmp #else error stop 'NETCDF missing! Please compile with netcdf if you want the particle dump.' #endif #ifdef USE_NCF if (present(initial_output)) then init_out=initial_output else init_out=.false. endif !$OMP PARALLEL PRIVATE(i,j,m,tmp,ns,i_av,cartxyz_comp,cartxyz,np,lskip) ! Some variables needed for temporal interpolation !************************************************* call find_time_vars(itime) !$OMP DO do i=1,count%allocated ! LB: Loop over all particles, including terminated ones because of ! averages that could still be available. There should be a better way. !Initialise fields lskip=.false. if (.not. part(i)%spawned) lskip=.true. ! Not spawned yet if ((.not. init_out) .and. (part(i)%tstart.eq.itime)) lskip=.true. ! No information avail yet for new parts if (((.not. part(i)%alive).and.(abs(part(i)%tend-itime).ge.ipoutfac*loutstep)) .or. & (init_out .and. (i.lt.partinitpointer1-1 .or. (part(i)%alive .eqv. .false.) ))) lskip=.true. ! no particles that have been dead for longer than a write interval if (lskip) then output(:,i) = -1 masstemp(i,:) = -1 masstemp_av(i,:) = -1 if (wetdep) wetdepotemp(i,:) = -1 if (drydep) drydepotemp(i,:) = -1 cycle endif !***************************************************************************** ! Interpolate several variables (PV, specific humidity, etc.) to particle position !***************************************************************************** ! Where in the grid? Stereographic (ngrid<0) or nested (ngrid>0) !*************************************************************** call find_ngrid(real(part(i)%xlon),real(part(i)%ylat)) call find_grid_indices(real(part(i)%xlon),real(part(i)%ylat)) call find_grid_distances(real(part(i)%xlon),real(part(i)%ylat)) ! First set dz1out from interpol_mod to -1 so it only is calculated once per particle !************************************************************************************ dz1out=-1 cartxyz_comp=.false. do np=1,num_partopt if (.not. partopt(np)%print) cycle ! Only compute when field should be printed i_av = partopt(np)%i_average if (init_out.and.(i_av.ne.0)) cycle ! no averages for initial particle output if ((i_av.ne.0).and.(part(i)%ntime.eq.0)) cycle ! no averages for freshly spawned particles select case (partopt(np)%name) case ('LO') output(np,i)=xlon0+real(part(i)%xlon)*dx cycle case ('LA') output(np,i)=ylat0+real(part(i)%ylat)*dy cycle case ('TO') ! Topography if (topo_written) cycle if (ngrid.le.0) then call hor_interpol(oro,output(np,i)) else call hor_interpol_nest(oron,output(np,i)) endif cycle case ('TR') ! Tropopause if (ngrid.le.0) then do m=1,2 call hor_interpol(tropopause,tmp(m),1,memind(m),1) end do else do m=1,2 call hor_interpol_nest(tropopausen,tmp(m),1,memind(m),1) end do endif call temporal_interpolation(tmp(1),tmp(2),output(np,i)) cycle case ('HM') ! PBL height if (ngrid.le.0) then do m=1,2 call hor_interpol(hmix,tmp(m),1,memind(m),1) end do else do m=1,2 call hor_interpol_nest(hmixn,tmp(m),1,memind(m),1) end do endif call temporal_interpolation(tmp(1),tmp(2),output(np,i)) cycle case ('ZZ') ! Height #ifdef ETA call update_zeta_to_z(itime, i) ! Convert eta z coordinate to meters if necessary #endif output(np,i)=real(part(i)%z) cycle ! case ('UU') ! Longitudinal velocity ! output(np,i)=part(i)%vel%u !This would be preferred, but not implemented yet ! cycle case ('VS') ! Settling velocity output(np,i)=part(i)%settling cycle case ('MA') ! Mass if (mass_written) cycle masstemp(i,:)=mass(i,:) cycle case ('ma') ! Mass averaged do ns=1,nspec masstemp_av(i,ns)=val_av(i, i_av+(ns-1))/part(i)%ntime end do cycle case ('WD') ! Wet deposition if (wetdep) then wetdepotemp(i,:)=wetdeposit(i,:) endif cycle case ('DD') ! dry deposition if (drydep) then drydepotemp(i,:)=drydeposit(i,:) endif cycle case ('lo') if (.not. cartxyz_comp) then cartxyz(1) = part(i)%cartx_av/part(i)%ntime cartxyz(2) = part(i)%carty_av/part(i)%ntime cartxyz(3) = part(i)%cartz_av/part(i)%ntime cartxyz_comp=.true. endif output(np,i) = atan2(cartxyz(1),-1.*cartxyz(2))/pi180 if (output(np,i).gt.360.) output(np,i)=output(np,i)-360. if (output(np,i).lt.0.) output(np,i)=output(np,i)+360. cycle case ('la') if (.not. cartxyz_comp) then cartxyz(1) = part(i)%cartx_av/part(i)%ntime cartxyz(2) = part(i)%carty_av/part(i)%ntime cartxyz(3) = part(i)%cartz_av/part(i)%ntime cartxyz_comp=.true. endif output(np,i) = atan2(cartxyz(3),sqrt(cartxyz(1)*cartxyz(1)+ & cartxyz(2)*cartxyz(2)))/pi180 case default if (.not. partopt(np)%average) then call interpol_partoutput_val(partopt(np)%name,output(np,i),i) else output(np,i) = val_av(i,i_av)/part(i)%ntime endif end select end do ! Reset dz1out !************* dz1out=-1 cartxyz_comp=.false. if ((.not. init_out).and.(n_average.gt.0)) then val_av(i,:) = 0. part(i)%ntime = 0. part(i)%cartx_av = 0. part(i)%carty_av = 0. part(i)%cartz_av = 0. endif end do !$OMP END DO !$OMP END PARALLEL if ((.not. init_out).and.(numpart.gt.0)) then do np=1,num_partopt if (.not. partopt(np)%print) cycle if (partopt(np)%name.eq.'MA') then write(*,*) partopt(np)%long_name, masstemp(1,:) else if (partopt(np)%name.eq.'ma') then write(*,*) partopt(np)%long_name, masstemp_av(1,:) else if (partopt(np)%name.eq.'WD'.and.wetdep) then write(*,*) partopt(np)%long_name, wetdepotemp(1,:) else if (partopt(np)%name.eq.'DD'.and.drydep) then write(*,*) partopt(np)%long_name, drydepotemp(1,:) else write(*,*) partopt(np)%long_name, output(np,1) endif end do write(*,*) 'Alive: ', count%alive, 'Total spawned: ', count%spawned, 'Terminated: ', count%terminated endif ! Determine current calendar date, needed for the file name !********************************************************** jul=bdate+real(itime,kind=dp)/86400._dp call caldate(jul,jjjjmmdd,ihmmss) write(adate,'(i8.8)') jjjjmmdd write(atime,'(i6.6)') ihmmss j=1 ! if (lnetcdfout.eq.1) then ! open output file if (init_out) then call open_partinit_file(ncid) else if (.not. lpartoutputperfield) then call open_partoutput_file(ncid) ! First allocate the time and particle dimensions within the netcdf file else do np=1,num_partopt if (.not. partopt(np)%print) cycle call open_partoutput_file(partopt(np)%ncid,np) end do endif call update_partoutput_pointers(itime, ncid) !ppointer_part = count%allocated endif ! Fill the fields in parallel if (numpart.gt.0) then ! OpenMP output does not work on all systems depending on how they are set-up ! !$OMP PARALLEL PRIVATE(np,ns,ncid_tmp) ! !$OMP DO SCHEDULE(dynamic) do np=1,num_partopt if (.not. partopt(np)%print) cycle if (init_out.and.(partopt(np)%i_average.ne.0)) cycle ! no averages for initial particle output if (lpartoutputperfield.and. (.not. init_out)) then ncid_tmp=partopt(np)%ncid else ncid_tmp = ncid endif if (partopt(np)%name.eq.'MA') then do ns=1,nspec if (init_out) then call partinit_netcdf(masstemp(:,ns),'MA',ns,ncid_tmp) else call partoutput_netcdf(itime,masstemp(:,ns),np,ns,ncid_tmp) endif end do else if (partopt(np)%name.eq.'ma') then do ns=1,nspec call partoutput_netcdf(itime,masstemp_av(:,ns),np,ns,ncid_tmp) end do else if ((.not. init_out).and.(partopt(np)%name.eq.'WD').and.wetdep) then do ns=1,nspec call partoutput_netcdf(itime,wetdepotemp(:,ns),np,ns,ncid_tmp) end do else if ((.not. init_out).and.(partopt(np)%name.eq.'DD').and.drydep) then do ns=1,nspec call partoutput_netcdf(itime,drydepotemp(:,ns),np,ns,ncid_tmp) end do else if (init_out) then call partinit_netcdf(output(np,:),partopt(np)%name,j,ncid_tmp) else call partoutput_netcdf(itime,output(np,:),np,j,ncid_tmp) endif endif if (lpartoutputperfield) call close_partoutput_file(ncid_tmp) end do ! !$OMP END DO ! !$OMP END PARALLEL endif if (.not. lpartoutputperfield) call close_partoutput_file(ncid) if (mdomainfill.ge.1 .and. (.not. init_out)) then mass_written=.true. ! needs to be reduced within openmp loop topo_written=.true. ! same endif !else ! Put binary function here !endif #else ! Put binary function here #endif end subroutine output_particles subroutine output_conc(itime,loutstart,loutend,loutnext,outnum) use unc_mod use outgrid_mod use par_mod use com_mod use binary_output_mod implicit none integer,intent(in) :: & itime ! time index integer,intent(inout) :: & loutstart,loutend, & ! concentration calculation starting and ending time loutnext real,intent(inout) :: & outnum ! concentration calculation sample number real(sp) :: & gridtotalunc ! concentration calculation related real(dep_prec) :: & wetgridtotalunc, & ! concentration calculation related drygridtotalunc ! concentration calculation related real :: & weight ! concentration calculation sample weight ! Is the time within the computation interval, if not, return !************************************************************ if ((ldirect*itime.lt.ldirect*loutstart).or.(ldirect*itime.gt.ldirect*loutend)) then return endif ! If we are exactly at the start or end of the concentration averaging interval, ! give only half the weight to this sample !***************************************************************************** if (mod(itime-loutstart,loutsample).eq.0) then if ((itime.eq.loutstart).or.(itime.eq.loutend)) then weight=0.5 else weight=1.0 endif outnum=outnum+weight if (iout.ne.0) call conccalc(itime,weight) endif ! If no grid is to be written to file, return (LB) !************************************************* if (iout.eq.0) then if (itime.ne.loutend) return loutnext=loutnext+loutstep loutstart=loutnext-loutaver/2 loutend=loutnext+loutaver/2 if (itime.eq.loutstart) then weight=0.5 outnum=outnum+weight endif return endif ! If it is not time yet to write outputs, return !*********************************************** if ((itime.ne.loutend).or.(outnum.le.0).or.(itime.eq.itime_init)) then return endif ! Output and reinitialization of grid ! If necessary, first sample of new grid is also taken !***************************************************** if ((iout.le.3.).or.(iout.eq.5)) then if (linversionout.eq.0) then ! regular output format if (lnetcdfout.eq.1) then #ifdef USE_NCF call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc) #endif else call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc) endif else ! inversion output - one file each release if (lnetcdfout.eq.1) then write(*,*) 'FLEXPART ERROR: netcdf output not avaiable yet for inversion format' error stop else call concoutput_inversion(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc) endif endif if (nested_output .eq. 1) then if (linversionout.eq.0) then ! regular output format if (lnetcdfout.eq.1) then #ifdef USE_NCF call concoutput_nest_netcdf(itime,outnum) #endif else call concoutput_nest(itime,outnum) endif else ! inversion output if (lnetcdfout.eq.1) then write(*,*) 'FLEXPART ERROR: netcdf output not avaiable yet for inversion format' error stop else call concoutput_inversion_nest(itime,outnum) endif endif endif outnum=0. endif write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc 45 format(i13,' Seconds simulated: ',i13, ' Particles: Uncertainty: ',3f7.3) loutnext=loutnext+loutstep loutstart=loutnext-loutaver/2 loutend=loutnext+loutaver/2 if (itime.eq.loutstart) then weight=0.5 outnum=outnum+weight call conccalc(itime,weight) endif end subroutine output_conc subroutine conccalc(itime,weight) ! i i !***************************************************************************** ! * ! Calculation of the concentrations on a regular grid using volume * ! sampling * ! * ! Author: A. Stohl * ! * ! 24 May 1996 * ! * ! April 2000: Update to calculate age spectra * ! Bug fix to avoid negative conc. at the domain boundaries, * ! as suggested by Petra Seibert * ! * ! 2 July 2002: re-order if-statements in order to optimize CPU time * ! * ! 2021, LB: OpenMP parallelisation * ! * !***************************************************************************** ! * ! Variables: * ! nspeciesdim = nspec for forward runs, 1 for backward runs * ! * !***************************************************************************** use unc_mod use outgrid_mod use par_mod use com_mod use omp_lib, only: OMP_GET_THREAD_NUM use interpol_mod, only: interpol_density #ifdef ETA use coord_ecmwf_mod #endif use particle_mod implicit none integer,intent(in) :: itime real,intent(in) :: weight integer :: itage,i,j,kz,ks,n,nage,inage,thread,ithread integer :: nrelpointer integer :: ix,jy,ixp,jyp real :: ddx,ddy real :: hx,hy,hz,hxyz,xd,yd,zd,xkern,r2,c(maxspec) real :: rhoi real :: xl,yl,wx,wy,w real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150. ! integer xscav_count ! For forward simulations, make a loop over the number of species; ! for backward simulations, make an additional loop over the ! releasepoints !*************************************************************************** ! xscav_count=0 #ifdef _OPENMP call omp_set_num_threads(numthreads_grid) #endif !$OMP PARALLEL PRIVATE(i,itage,nage,inage,rhoi,nrelpointer,kz,xl,yl,ks,wx,wy,w,thread,ddx,ddy, & !$OMP ix,jy,ixp,jyp) #if (defined _OPENMP) thread = OMP_GET_THREAD_NUM()+1 ! Starts with 1 #else thread = 1 #endif !$OMP DO do j=1,count%alive i=count%ialive(j) ! Determine age class of the particle itage=abs(itime-part(i)%tstart) nage=1 if (lagespectra.eq.1) then do inage=1,nageclass nage=inage if (itage.lt.lage(nage)) exit end do endif ! if (xscav_frac1(i,1).lt.0) xscav_count=xscav_count+1 !************************************************************************ ! For special runs, interpolate the air density to the particle position !************************************************************************ !AF IND_SOURCE switches between different units for concentrations at the source !Af NOTE that in backward simulations the release of particles takes place !Af at the receptor and the sampling at the source. !Af 1="mass" !Af 2="mass mixing ratio" !Af IND_RECEPTOR switches between different units for concentrations at the receptor !Af 1="mass" !Af 2="mass mixing ratio" !Af switches for the conccalcfile: !AF IND_SAMP = 0 : xmass * 1 !Af IND_SAMP = -1 : xmass / rho !Af ind_samp is defined in readcommand.f !************************************************************************ if ( ind_samp .eq. -1 ) then #ifdef ETA call update_zeta_to_z(itime,i) #endif call interpol_density(itime,i,rhoi) elseif (ind_samp.eq.0) then rhoi = 1. endif !**************************************************************************** ! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy !**************************************************************************** ! For backward simulations, look from which release point the particle comes from ! For domain-filling trajectory option, npoint contains a consecutive particle ! number, not the release point information. Therefore, nrelpointer is set to 1 ! for the domain-filling option. !***************************************************************************** if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then nrelpointer=1 else nrelpointer=part(i)%npoint endif do kz=1,numzgrid ! determine height of cell if (outheight(kz).gt.part(i)%z) exit end do if (kz.le.numzgrid) then ! inside output domain !******************************** ! Do everything for mother domain !******************************** xl=(real(part(i)%xlon)*dx+xoutshift)/dxout yl=(real(part(i)%ylat)*dy+youtshift)/dyout ix=int(xl) if (xl.lt.0.) ix=ix-1 jy=int(yl) if (yl.lt.0.) jy=jy-1 ! For particles aged less than 3 hours, attribute particle mass to grid cell ! it resides in rather than use the kernel, in order to avoid its smoothing effect. ! For older particles, use the uniform kernel. ! If a particle is close to the domain boundary, do not use the kernel either. !***************************************************************************** if ((.not.lusekerneloutput).or.(itage.lt.10800).or. & (xl.lt.0.5).or.(yl.lt.0.5).or. & (xl.gt.real(numxgrid-1)-0.5).or. & (yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. & (jy.le.numygrid-1)) then if (DRYBKDEP.or.WETBKDEP) then do ks=1,nspec #ifdef _OPENMP gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0) #else gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0) #endif end do else if (lparticlecountoutput) then do ks=1,nspec #ifdef _OPENMP gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+1 #else gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+1 #endif end do else if (llcmoutput) then ! special case LCM output use mass ratio species to airtracer ! species 1 is always airtracer do ks=2,nspec #ifdef _OPENMP gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/mass(i,1)*weight #else gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/mass(i,1)*weight #endif end do #ifdef _OPENMP gridcnt_omp(ix,jy,kz,thread)=gridcnt_omp(ix,jy,kz,thread)+weight #else gridcnt(ix,jy,kz)=gridcnt(ix,jy,kz)+weight #endif else do ks=1,nspec #ifdef _OPENMP gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight #else gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight #endif end do end if ! llcmoutput end if endif endif else ! attribution via uniform kernel ddx=xl-real(ix) ! distance to left cell border ddy=yl-real(jy) ! distance to lower cell border if (ddx.gt.0.5) then ixp=ix+1 wx=1.5-ddx else ixp=ix-1 wx=0.5+ddx endif if (ddy.gt.0.5) then jyp=jy+1 wy=1.5-ddy else jyp=jy-1 wy=0.5+ddy endif ! Determine mass fractions for four grid points !********************************************** if ((ix.ge.0).and.(ix.le.numxgrid-1)) then if ((jy.ge.0).and.(jy.le.numygrid-1)) then w=wx*wy if (DRYBKDEP.or.WETBKDEP) then do ks=1,nspec #ifdef _OPENMP gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0) #else gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0) #endif end do else if (llcmoutput) then ! special case CTM output use mass ratio species to airtracer ! species 1 is always airtracer do ks=2,nspec #ifdef _OPENMP gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/mass(i,1)*weight*w #else gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/mass(i,1)*weight*w #endif end do #ifdef _OPENMP gridcnt_omp(ix,jy,kz,thread)=gridcnt_omp(ix,jy,kz,thread)+w*weight #else gridcnt(ix,jy,kz)=gridcnt(ix,jy,kz)+w*weight #endif else do ks=1,nspec #ifdef _OPENMP gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w #else gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w #endif end do endif ! llcmoutput endif endif if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then w=wx*(1.-wy) if (DRYBKDEP.or.WETBKDEP) then do ks=1,nspec #ifdef _OPENMP gridunc_omp(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #else gridunc(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #endif end do else if (llcmoutput) then ! special case CTM output use mass ratio species to airtracer ! species 1 is always airtracer do ks=2,nspec #ifdef _OPENMP gridunc_omp(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/mass(i,1)*weight*w #else gridunc(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/mass(i,1)*weight*w #endif end do #ifdef _OPENMP gridcnt_omp(ix,jyp,kz,thread)=gridcnt_omp(ix,jyp,kz,thread)+w*weight #else gridcnt(ix,jyp,kz)=gridcnt(ix,jyp,kz)+w*weight #endif else do ks=1,nspec #ifdef _OPENMP gridunc_omp(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w #else gridunc(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w #endif end do endif ! llcmoutput endif endif endif !ix ge 0 if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then w=(1.-wx)*(1.-wy) if (DRYBKDEP.or.WETBKDEP) then do ks=1,nspec #ifdef _OPENMP gridunc_omp(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0) #else gridunc(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0) #endif end do else if (llcmoutput) then ! special case CTM output use mass ratio species to airtracer ! species 1 is always airtracer do ks=2,nspec #ifdef _OPENMP gridunc_omp(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/mass(i,1)*weight*w #else gridunc(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/mass(i,1)*weight*w #endif end do #ifdef _OPENMP gridcnt_omp(ixp,jyp,kz,thread)=gridcnt_omp(ixp,jyp,kz,thread)+w*weight #else gridcnt(ixp,jyp,kz)=gridcnt(ixp,jyp,kz)+w*weight #endif else do ks=1,nspec #ifdef _OPENMP gridunc_omp(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w #else gridunc(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w #endif end do endif ! llcmoutput endif endif if ((jy.ge.0).and.(jy.le.numygrid-1)) then w=(1.-wx)*wy if (DRYBKDEP.or.WETBKDEP) then do ks=1,nspec #ifdef _OPENMP gridunc_omp(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #else gridunc(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #endif end do else if (llcmoutput) then ! special case CTM output use mass ratio species to airtracer ! species 1 is always airtracer do ks=2,nspec #ifdef _OPENMP gridunc_omp(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/mass(i,1)*weight*w #else gridunc(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/mass(i,1)*weight*w #endif end do #ifdef _OPENMP gridcnt_omp(ixp,jy,kz,thread)=gridcnt_omp(ixp,jy,kz,thread)+w*weight #else gridcnt(ixp,jy,kz)=gridcnt(ixp,jy,kz)+w*weight #endif else do ks=1,nspec #ifdef _OPENMP gridunc_omp(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & gridunc_omp(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w #else gridunc(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & gridunc(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w #endif end do endif ! llcmoutput endif endif endif !ixp ge 0 endif !************************************ ! Do everything for the nested domain !************************************ if (nested_output.eq.1) then xl=(real(part(i)%xlon)*dx+xoutshiftn)/dxoutn yl=(real(part(i)%ylat)*dy+youtshiftn)/dyoutn ix=int(xl) if (xl.lt.0.) ix=ix-1 jy=int(yl) if (yl.lt.0.) jy=jy-1 ! For particles aged less than 3 hours, attribute particle mass to grid cell ! it resides in rather than use the kernel, in order to avoid its smoothing effect. ! For older particles, use the uniform kernel. ! If a particle is close to the domain boundary, do not use the kernel either. !***************************************************************************** if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. & (xl.gt.real(numxgridn-1)-0.5).or. & (yl.gt.real(numygridn-1)-0.5).or.((.not.lusekerneloutput))) then ! no kernel, direct attribution to grid cell if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. & (jy.le.numygridn-1)) then if (DRYBKDEP.or.WETBKDEP) then do ks=1,nspec #ifdef _OPENMP griduncn_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & griduncn_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0) #else griduncn(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & griduncn(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0) #endif end do else if (lparticlecountoutput) then do ks=1,nspec #ifdef _OPENMP griduncn_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & griduncn_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+1 #else griduncn(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & griduncn(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+1 #endif end do else do ks=1,nspec #ifdef _OPENMP griduncn_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & griduncn_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight #else griduncn(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & griduncn(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight #endif end do endif endif endif else ! attribution via uniform kernel ddx=xl-real(ix) ! distance to left cell border ddy=yl-real(jy) ! distance to lower cell border if (ddx.gt.0.5) then ixp=ix+1 wx=1.5-ddx else ixp=ix-1 wx=0.5+ddx endif if (ddy.gt.0.5) then jyp=jy+1 wy=1.5-ddy else jyp=jy-1 wy=0.5+ddy endif ! Determine mass fractions for four grid points !********************************************** if ((ix.ge.0).and.(ix.le.numxgridn-1)) then if ((jy.ge.0).and.(jy.le.numygridn-1)) then w=wx*wy if (DRYBKDEP.or.WETBKDEP) then do ks=1,nspec #ifdef _OPENMP griduncn_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & griduncn_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #else griduncn(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & griduncn(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #endif end do else do ks=1,nspec #ifdef _OPENMP griduncn_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & griduncn_omp(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w #else griduncn(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & griduncn(ix,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w #endif end do endif endif if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then w=wx*(1.-wy) if (DRYBKDEP.or.WETBKDEP) then do ks=1,nspec #ifdef _OPENMP griduncn_omp(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & griduncn_omp(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #else griduncn(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)= & griduncn(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #endif end do else do ks=1,nspec #ifdef _OPENMP griduncn_omp(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & griduncn_omp(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w #else griduncn(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)= & griduncn(ix,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w #endif end do endif endif endif if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then w=(1.-wx)*(1.-wy) if (DRYBKDEP.or.WETBKDEP) then do ks=1,nspec #ifdef _OPENMP griduncn_omp(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & griduncn_omp(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #else griduncn(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)= & griduncn(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #endif end do else do ks=1,nspec #ifdef _OPENMP griduncn_omp(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & griduncn_omp(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w #else griduncn(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)= & griduncn(ixp,jyp,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w #endif end do endif endif if ((jy.ge.0).and.(jy.le.numygridn-1)) then w=(1.-wx)*wy if (DRYBKDEP.or.WETBKDEP) then do ks=1,nspec #ifdef _OPENMP griduncn_omp(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & griduncn_omp(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #else griduncn(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & griduncn(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0) #endif end do else do ks=1,nspec #ifdef _OPENMP griduncn_omp(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)= & griduncn_omp(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage,thread)+ & mass(i,ks)/rhoi*weight*w #else griduncn(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage)= & griduncn(ixp,jy,kz,ks,nrelpointer,part(i)%nclass,nage)+ & mass(i,ks)/rhoi*weight*w #endif end do endif endif endif endif endif endif end do !$OMP END DO !$OMP END PARALLEL #ifdef _OPENMP call omp_set_num_threads(numthreads) #endif ! Reduction of gridunc and griduncn #ifdef _OPENMP do ithread=1,numthreads_grid gridunc(:,:,:,:,:,:,:)=gridunc(:,:,:,:,:,:,:)+gridunc_omp(:,:,:,:,:,:,:,ithread) gridcnt(:,:,:)=gridcnt(:,:,:)+gridcnt_omp(:,:,:,ithread) gridunc_omp(:,:,:,:,:,:,:,ithread)=0. gridcnt_omp(:,:,:,ithread)=0. end do if (nested_output.eq.1) then do ithread=1,numthreads_grid griduncn(:,:,:,:,:,:,:)=griduncn(:,:,:,:,:,:,:)+griduncn_omp(:,:,:,:,:,:,:,ithread) griduncn_omp(:,:,:,:,:,:,:,ithread)=0. end do endif #endif end subroutine conccalc subroutine partpos_avg(itime,j) !********************************************************************** ! This subroutine averages particle quantities, to be used for particle ! dump (in partoutput.f90). Averaging is done over output interval. ! Author: A. Stohl ! Changes L Bakels: ! - Computing fields defined in PARTOPTIONS !********************************************************************** use par_mod use com_mod use interpol_mod #ifdef ETA use coord_ecmwf_mod #endif implicit none integer,intent(in) :: itime,j integer :: np,i_av,ns,m real :: xlon,ylat,x,y,z real :: hm(2) real :: output real :: tr(2)!,energy logical :: cart_comp if (ipout.eq.0) return ! No need to compute averages since there is no particle output if (n_average.eq.0) return if (.not. part(j)%alive) return if (part(j)%nstop) return ! If particle is to be killed, averages cannot be computed ! Some variables needed for temporal interpolation !************************************************* call find_time_vars(itime) xlon=xlon0+real(part(j)%xlon)*dx ylat=ylat0+real(part(j)%ylat)*dy !***************************************************************************** ! Interpolate several variables (PV, specific humidity, etc.) to particle position !***************************************************************************** ! Where in the grid? Stereographic (ngrid<0) or nested (ngrid>0) !*************************************************************** call find_ngrid(real(part(j)%xlon),real(part(j)%ylat)) call find_grid_indices(real(part(j)%xlon),real(part(j)%ylat)) call find_grid_distances(real(part(j)%xlon),real(part(j)%ylat)) ! First set dz1out from interpol_mod to -1 so it only is calculated once per particle !************************************************************************************ part(j)%ntime=part(j)%ntime + 1 dz1out=-1 cart_comp=.false. do np=1,num_partopt if ((.not. partopt(np)%print) .or. (.not. partopt(np)%average)) cycle i_av = partopt(np)%i_average select case (partopt(np)%name) case ('to') if (ngrid.le.0) then call hor_interpol(oro,output) else call hor_interpol_nest(oron,output) endif val_av(j,i_av)=val_av(j,i_av)+output case ('tr') if (ngrid.le.0) then do m=1,2 call hor_interpol(tropopause,tr(m),1,memind(m),1) end do else do m=1,2 call hor_interpol_nest(tropopausen,tr(m),1,memind(m),1) end do endif call temporal_interpolation(tr(1),tr(2),output) val_av(j,i_av)=val_av(j,i_av)+output case ('hm') if (ngrid.le.0) then do m=1,2 call hor_interpol(hmix,hm(m),1,memind(m),1) end do else do m=1,2 call hor_interpol_nest(hmixn,hm(m),1,memind(m),1) end do endif call temporal_interpolation(hm(1),hm(2),output) val_av(j,i_av)=val_av(j,i_av)+output case ('lo') if (.not. cart_comp) then ! Calculate Cartesian 3D coordinates suitable for averaging !********************************************************** xlon=xlon*pi180 ylat=ylat*pi180 x = cos(ylat)*sin(xlon) y = -1.*cos(ylat)*cos(xlon) z = sin(ylat) part(j)%cartx_av=part(j)%cartx_av+x part(j)%carty_av=part(j)%carty_av+y part(j)%cartz_av=part(j)%cartz_av+z cart_comp=.true. endif case ('la') if (.not. cart_comp) then ! Calculate Cartesian 3D coordinates suitable for averaging !********************************************************** xlon=xlon*pi180 ylat=ylat*pi180 x = cos(ylat)*sin(xlon) y = -1.*cos(ylat)*cos(xlon) z = sin(ylat) part(j)%cartx_av=part(j)%cartx_av+x part(j)%carty_av=part(j)%carty_av+y part(j)%cartz_av=part(j)%cartz_av+z cart_comp=.true. endif case ('zz') ! Convert eta z coordinate to meters if necessary. Can be moved to output only !************************************************ #ifdef ETA call update_zeta_to_z(itime,j) #endif val_av(j,i_av)=val_av(j,i_av)+real(part(j)%z) case ('ma') do ns=1,nspec val_av(j,i_av+(ns-1))=val_av(j,i_av+(ns-1))+mass(j,ns) end do case ('vs') val_av(j,i_av)=val_av(j,i_av)+part(j)%settling case default call interpol_partoutput_val(partopt(np)%name,output,j) val_av(j,i_av)=val_av(j,i_av)+output end select end do ! Reset dz1out !************* dz1out=-1 cart_comp=.false. return end subroutine partpos_avg end module output_mod