Commit f3054ea9 authored by Espen Sollum's avatar Espen Sollum
Browse files

Changed from grib_api to eccodes. MPI: implemented linversionout=1; fixed...

Changed from grib_api to eccodes. MPI: implemented linversionout=1; fixed calculation of grid_initial fields.
parent a756649a
./options_test/
./output_test/
/
/xnilu_wrk/flex_wrk/WIND_FIELDS/AVAILABLE_ECMWF_OPER_fields_global
./options/
./output/
./preprocess/flex_extract/work/
./AVAILABLE
......@@ -82,8 +82,6 @@ program flexpart
flexversion='Version '//trim(flexversion_major)//'.4 (2019-07-23)'
verbosity=0
write(*,*) "flexversion: ", flexversion
! Read the pathnames where input/output files are stored
!*******************************************************
......
......@@ -76,6 +76,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
integer :: sp_count_i,sp_count_r
real :: sp_fact
real :: outnum,densityoutrecept(maxreceptor),xl,yl
! RLT
real :: densitydryrecept(maxreceptor)
real :: factor_dryrecept(maxreceptor)
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
......@@ -110,6 +113,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! mind eso: added to ensure identical results between 2&3-fields versions
character(LEN=8),save :: file_stat='REPLACE'
logical :: ldates_file
logical :: lexist
integer :: ierr
character(LEN=100) :: dates_char
......@@ -211,6 +215,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! rho(iix,jjy,kzz-1,2)*dz2)/dz
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
rho(iix,jjy,kzz-1,mind)*dz2)/dz
! RLT
densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz
end do
end do
end do
......@@ -222,8 +229,14 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
jjy=max(min(nint(yl),nymin1),0)
!densityoutrecept(i)=rho(iix,jjy,1,2)
densityoutrecept(i)=rho(iix,jjy,1,mind)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
end do
! RLT
! conversion factor for output relative to dry air
factor_drygrid=densityoutgrid/densitydrygrid
factor_dryrecept=densityoutrecept/densitydryrecept
! Output is different for forward and backward simulations
do kz=1,numzgrid
......@@ -276,11 +289,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
endif
write(unitoutgrid) itime
endif
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
atime//'_'//anspec,form='unformatted')
write(unitoutgridppt) itime
endif
endif ! if deposition output
......@@ -624,6 +635,49 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do
! RLT Aug 2017
! Write out conversion factor for dry air
inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
if (lexist) then
! open and append
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
status='old',action='write',access='append')
else
! create new
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
status='new',action='write')
endif
sp_count_i=0
sp_count_r=0
sp_fact=-1.
sp_zer=.true.
do kz=1,numzgrid
do jy=0,numygrid-1
do ix=0,numxgrid-1
if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
if (sp_zer.eqv..true.) then ! first value not equal to one
sp_count_i=sp_count_i+1
sparse_dump_i(sp_count_i)= &
ix+jy*numxgrid+kz*numxgrid*numygrid
sp_zer=.false.
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
sparse_dump_r(sp_count_r)= &
sp_fact*factor_drygrid(ix,jy,kz)
else ! factor is one
sp_zer=.true.
endif
end do
end do
end do
write(unitoutfactor) sp_count_i
write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutfactor) sp_count_r
write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
close(unitoutfactor)
if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
wetgridtotal
......@@ -650,7 +704,23 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do
endif
! RLT Aug 2017
! Write out conversion factor for dry air
if (numreceptor.gt.0) then
inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
if (lexist) then
! open and append
open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
status='old',action='write',access='append')
else
! create new
open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
status='new',action='write')
endif
write(unitoutfactor) itime
write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
close(unitoutfactor)
endif
! Reinitialization of grid
!*************************
......
......@@ -605,7 +605,6 @@ subroutine concoutput_nest(itime,outnum)
write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
close(unitoutfactor)
! Reinitialization of grid
!*************************
......
......@@ -57,7 +57,6 @@ subroutine concoutput_nest(itime,outnum)
! *
!*****************************************************************************
use unc_mod
use point_mod
use outg_mod
......@@ -73,6 +72,9 @@ subroutine concoutput_nest(itime,outnum)
integer :: sp_count_i,sp_count_r
real :: sp_fact
real :: outnum,densityoutrecept(maxreceptor),xl,yl
! RLT
real :: densitydryrecept(maxreceptor)
real :: factor_dryrecept(maxreceptor)
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
......@@ -99,6 +101,7 @@ subroutine concoutput_nest(itime,outnum)
logical :: sp_zer
character :: adate*8,atime*6
character(len=3) :: anspec
logical :: lexist
integer :: mind
! mind eso:added to ensure identical results between 2&3-fields versions
......@@ -175,6 +178,9 @@ subroutine concoutput_nest(itime,outnum)
! rho(iix,jjy,kzz-1,2)*dz2)/dz
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
rho(iix,jjy,kzz-1,mind)*dz2)/dz
! RLT
densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz
end do
end do
end do
......@@ -186,8 +192,14 @@ subroutine concoutput_nest(itime,outnum)
jjy=max(min(nint(yl),nymin1),0)
!densityoutrecept(i)=rho(iix,jjy,1,2)
densityoutrecept(i)=rho(iix,jjy,1,mind)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
end do
! RLT
! conversion factor for output relative to dry air
factor_drygrid=densityoutgrid/densitydrygrid
factor_dryrecept=densityoutrecept/densitydryrecept
! Output is different for forward and backward simulations
do kz=1,numzgrid
......@@ -210,6 +222,16 @@ subroutine concoutput_nest(itime,outnum)
do ks=1,nspec
write(anspec,'(i3.3)') ks
if (DRYBKDEP.or.WETBKDEP) then !scavdep output
if (DRYBKDEP) &
open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_nest_'//adate// &
atime//'_'//anspec,form='unformatted')
if (WETBKDEP) &
open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_nest_'//adate// &
atime//'_'//anspec,form='unformatted')
write(unitoutgrid) itime
else
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
if (ldirect.eq.1) then
open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_nest_' &
......@@ -222,6 +244,7 @@ subroutine concoutput_nest(itime,outnum)
endif
write(unitoutgrid) itime
endif
endif
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_nest_' &
......@@ -552,29 +575,71 @@ subroutine concoutput_nest(itime,outnum)
end do
! Reinitialization of grid
!*************************
do ks=1,nspec
do kp=1,maxpointspec_act
do i=1,numreceptor
creceptor(i,ks)=0.
end do
! RLT Aug 2017
! Write out conversion factor for dry air
inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
if (lexist) then
! open and append
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
status='old',action='write',access='append')
else
! create new
open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
status='new',action='write')
endif
sp_count_i=0
sp_count_r=0
sp_fact=-1.
sp_zer=.true.
do kz=1,numzgrid
do jy=0,numygridn-1
do ix=0,numxgridn-1
do l=1,nclassunc
do nage=1,nageclass
do kz=1,numzgrid
griduncn(ix,jy,kz,ks,kp,l,nage)=0.
end do
end do
end do
end do
if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
if (sp_zer.eqv..true.) then ! first value not equal to one
sp_count_i=sp_count_i+1
sparse_dump_i(sp_count_i)= &
ix+jy*numxgridn+kz*numxgridn*numygridn
sp_zer=.false.
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
sparse_dump_r(sp_count_r)= &
sp_fact*factor_drygrid(ix,jy,kz)
else ! factor is one
sp_zer=.true.
endif
end do
end do
end do
write(unitoutfactor) sp_count_i
write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutfactor) sp_count_r
write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
close(unitoutfactor)
! Reinitialization of grid
!*************************
! do ks=1,nspec
! do kp=1,maxpointspec_act
! do i=1,numreceptor
! creceptor(i,ks)=0.
! end do
! do jy=0,numygridn-1
! do ix=0,numxgridn-1
! do l=1,nclassunc
! do nage=1,nageclass
! do kz=1,numzgrid
! griduncn(ix,jy,kz,ks,kp,l,nage)=0.
! end do
! end do
! end do
! end do
! end do
! end do
! end do
creceptor(:,:)=0.
griduncn(:,:,:,:,:,:,:)=0.
if (mp_measure_time) call mpif_mtime('iotime',1)
! if (mp_measure_time) then
......
......@@ -62,14 +62,16 @@ ifeq ($(gcc), 4.9)
LIBPATH1 = ${ROOT_DIR}/lib
else
# Compiled libraries under user ~flexpart, gfortran v5.4
ROOT_DIR = /homevip/flexpart/
# ROOT_DIR = /homevip/flexpart/
F90 = /usr/bin/gfortran
MPIF90 = /usr/bin/mpifort
INCPATH1 = ${ROOT_DIR}/gcc-5.4.0/include
# INCPATH1 = ${ROOT_DIR}/gcc-5.4.0/include
# INCPATH2 = /usr/include
# LIBPATH1 = ${ROOT_DIR}/gcc-5.4.0/lib
INCPATH1 = /usr/include
INCPATH2 = /usr/include
LIBPATH1 = ${ROOT_DIR}/gcc-5.4.0/lib
endif
......@@ -91,8 +93,7 @@ O_LEV = 0 # [0,1,2,3,g,s,fast]
O_LEV_DBG = g # [0,g]
## LIBRARIES
#LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff
LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper $(NCOPT)
LIBS = -leccodes -leccodes_f90 -lm $(NCOPT)
FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(NCOPT) $(FUSER) #-Warray-bounds -fcheck=all # -march=native
......@@ -127,10 +128,10 @@ OBJECTS_SERIAL = \
boundcond_domainfill.o \
redist.o \
concoutput_surf.o concoutput_surf_nest.o \
concoutput_inversion_nest.o \
concoutput_inversion.o \
getfields.o \
readwind_ecmwf.o
readwind_ecmwf.o \
initial_cond_output.o \
initial_cond_output_inversion.o
## For MPI version
OBJECTS_MPI = releaseparticles_mpi.o partoutput_mpi.o \
......@@ -144,11 +145,14 @@ OBJECTS_MPI = releaseparticles_mpi.o partoutput_mpi.o \
redist_mpi.o \
concoutput_surf_mpi.o concoutput_surf_nest_mpi.o \
getfields_mpi.o \
readwind_ecmwf_mpi.o
readwind_ecmwf_mpi.o \
initial_cond_output_mpi.o \
initial_cond_output_inversion_mpi.o
OBJECTS_NCF = netcdf_output_mod.o
OBJECTS = \
initial_cond_calc.o \
advance.o initialize.o \
writeheader.o writeheader_txt.o \
partpos_average.o writeprecip.o \
......@@ -201,10 +205,12 @@ writeheader_nest.o writeheader_nest_surf.o \
wetdepokernel_nest.o \
drydepokernel_nest.o zenithangle.o \
ohreaction.o getvdep_nests.o \
initial_cond_calc.o initial_cond_output.o initial_cond_output_inversion.o \
dynamic_viscosity.o get_settling.o \
initialize_cbl_vel.o re_initialize_particle.o \
cbl.o
cbl.o \
concoutput_inversion_nest.o \
concoutput_inversion.o \
ifeq ($(ncf), yes)
OBJECTS := $(OBJECTS) $(OBJECTS_NCF)
......@@ -324,9 +330,13 @@ hanna_short.o: com_mod.o hanna_mod.o par_mod.o
init_domainfill.o: com_mod.o par_mod.o point_mod.o random_mod.o
init_domainfill_mpi.o: com_mod.o mpi_mod.o par_mod.o point_mod.o random_mod.o
initial_cond_calc.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
initial_cond_calc_mpi.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
initial_cond_output.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
initial_cond_output_mpi.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mpi_mod.o
initial_cond_output_inversion.o: com_mod.o outg_mod.o par_mod.o point_mod.o \
unc_mod.o
initial_cond_output_inversion_mpi.o: com_mod.o outg_mod.o par_mod.o point_mod.o \
unc_mod.o
initialize.o: com_mod.o hanna_mod.o interpol_mod.o par_mod.o random_mod.o
initialize_cbl_vel.o: com_mod.o par_mod.o random_mod.o
interpol_all.o: com_mod.o hanna_mod.o interpol_mod.o par_mod.o
......
......@@ -2597,6 +2597,58 @@ contains
601 end subroutine mpif_tm_reduce_grid_nest
subroutine mpif_tm_reduce_initcond
!***********************************************************************
! Collect init_cond to PID 0, adding from all processes.
!
!
!***********************************************************************
use com_mod
use unc_mod
use par_mod
implicit none
integer :: grid_size
!**********************************************************************
grid_size=numxgrid*numygrid*numzgrid*maxspec* &
& maxpointspec_act
! Time for MPI communications
if (mp_measure_time) call mpif_mtime('commtime',0)
#ifdef USE_MPIINPLACE
! Using in-place reduction
if (lroot) then
call MPI_Reduce(MPI_IN_PLACE, init_cond, grid_size, mp_sp, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
if (mp_ierr /= 0) goto 600
else
call MPI_Reduce(init_cond, 0, grid_size, mp_sp, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
if (mp_ierr /= 0) goto 600
end if
#else
call MPI_Reduce(init_cond, init_cond0, grid_size, mp_sp, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
if (mp_ierr /= 0) goto 600
! if (lroot) init_cond = init_cond0
#endif
if (mp_measure_time) call mpif_mtime('commtime',1)
goto 601
600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
stop
601 end subroutine mpif_tm_reduce_initcond
subroutine mpif_mtime(ident,imode)
!***********************************************************************
! Measure CPU/Wall time in various parts of the code
......
......@@ -314,6 +314,14 @@ subroutine outgrid_init
allocate(init_cond(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
maxpointspec_act),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate init_cond'
if (mpi_mode.gt.0) then
if (lroot) then
allocate(init_cond0(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec, &
maxpointspec_act),stat=stat)
else
allocate(init_cond0(1,1,1,1,1))
end if
end if
endif
!************************
......
......@@ -107,11 +107,12 @@ subroutine timemanager(metdata_format)
implicit none
integer(selected_int_kind(16)) :: idummy,idummy2
integer :: metdata_format
integer :: j,ks,kp,l,n,itime=0,nstop,nstop1
! integer :: ksp
integer :: loutnext,loutstart,loutend
integer :: ix,jy,ldeltat,itage,nage,idummy,idummy2
integer :: ix,jy,ldeltat,itage,nage
integer :: i_nan=0,ii_nan,total_nan_intl=0 !added by mc to check instability in CBL scheme
real :: outnum,weight,prob_rec(maxspec),prob(maxspec),decfact,wetscav
! real :: uap(maxpart),ucp(maxpart),uzp(maxpart)
......@@ -396,7 +397,7 @@ subroutine timemanager(metdata_format)
else
call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
endif
else
else ! surf only
if (verbosity.eq.1) then
print*,'call concoutput_surf '
call system_clock(count_clock)
......
<
......@@ -108,13 +108,14 @@ subroutine timemanager(metdata_format)
implicit none
integer(selected_int_kind(16)) :: idummy,idummy2
integer :: metdata_format
logical :: reqv_state=.false. ! .true. if waiting for a MPI_Irecv to complete
integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0
! integer :: ksp
integer :: ip,irec
integer :: loutnext,loutstart,loutend
integer :: ix,jy,ldeltat,itage,nage,idummy
integer :: ix,jy,ldeltat,itage,nage
integer :: i_nan=0,ii_nan,total_nan_intl=0 !added by mc to check instability in CBL scheme
integer :: numpart_tot_mpi ! for summing particles on all processes
real :: outnum,weight,prob(maxspec), prob_rec(maxspec), decfact,wetscav
......@@ -176,7 +177,7 @@ subroutine timemanager(metdata_format)
! changed by Petra Seibert 9/02
!********************************************************************
if (mp_dbg_mode) write(*,*) 'myid, itime: ',mp_pid,itime
! if (mp_dbg_mode) write(*,*) 'myid, itime: ',mp_pid,itime
if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
if (verbosity.gt.0) then
......@@ -297,7 +298,7 @@ subroutine timemanager(metdata_format)
!*********************************************************************
if (lmpreader.and.lmp_use_reader) then
if (itime.lt.ideltas*ldirect) then
if (abs(itime).lt.ideltas*ldirect) then
cycle
else
goto 999
......@@ -497,16 +498,20 @@ subroutine timemanager(metdata_format)
gridunc(:,:,:,:,:,:,:)=0.
creceptor(:,:)=0.
end if
else
else ! surf only
if (lroot) then
if (lnetcdfout.eq.1) then
#ifdef USE_NCF
call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,&
&drygridtotalunc)
#endif
else
if (linversionout.eq.1) then
call concoutput_inversion(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
else
call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
end if
end if
else
! zero arrays on non-root processes
gridunc(:,:,:,:,:,:,:)=0.
......@@ -533,7 +538,19 @@ subroutine timemanager(metdata_format)
end if
else
if(linversionout.eq.1) then
if (lroot) then
call concoutput_inversion_nest(itime,outnum)
else
griduncn(:,:,:,:,:,:,:)=0.
end if
else
if (lroot) then
call concoutput_surf_nest(itime,outnum)
else
griduncn(:,:,:,:,:,:,:)=0.
end if
end if
end if
else
#ifdef USE_NCF
......@@ -741,7 +758,7 @@ subroutine timemanager(metdata_format)
if (WETBKDEP) then
do ks=1,nspec
if ((xscav_frac1(j,ks).lt.0)) then
call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy,wetscav)
call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy2,wetscav)
if (wetscav.gt.0) then
xscav_frac1(j,ks)=wetscav* &
(zpoint2(npoint(j))-zpoint1(npoint(j)))*grfraction(1)
......@@ -894,14 +911,12 @@ subroutine timemanager(metdata_format)
! Complete the calculation of initial conditions for particles not yet terminated
!*****************************************************************************
! eso :TODO: this not implemented yet (transfer particles to PID 0 or rewrite)
! the tools to do this are already in mpi_mod.f90
if (lroot) then
do j=1,numpart
if (linit_cond.ge.1) call initial_cond_calc(itime,j)
end do
end if
! Transfer sum of init_cond field to root process, for output
call mpif_tm_reduce_initcond
if (ipout.eq.2) then
! MPI process 0 creates the file, the other processes append to it
......@@ -914,8 +929,14 @@ subroutine timemanager(metdata_format)
end do
end if
! eso :TODO: MPI