Skip to content

F
flexpart
Commit 2eefa58a authored by Espen Sollum's avatar Espen Sollum
Browse files
Options

Added Ronas changes for inversion output

parent f963113e
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 2135 additions and 42 deletions
src/FLEXPART.f90
View file @ 2eefa58a
......@@ -136,7 +136,7 @@ program flexpart
print*,'nxshift=',nxshift
write(*,*) 'call readpaths'
endif
call readpaths(pathfile)
call readpaths
if (verbosity.gt.1) then !show clock info
!print*,'length(4)',length(4)
......@@ -451,7 +451,9 @@ program flexpart
print*,'call timemanager'
endif
if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
call timemanager(metdata_format)
if (verbosity.gt.0) then
! NIK 16.02.2005
......@@ -466,7 +468,6 @@ program flexpart
write(*,*) '**********************************************'
endif
end do
write (*,*) 'timemanager> call wetdepo'
endif
write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
......
src/FLEXPART_MPI.f90
View file @ 2eefa58a
......@@ -145,7 +145,7 @@ program flexpart
if (verbosity.gt.0) then
write(*,*) 'call readpaths'
endif
call readpaths(pathfile)
call readpaths
if (verbosity.gt.1) then !show clock info
!print*,'length(4)',length(4)
......
src/com_mod.f90
View file @ 2eefa58a
......@@ -123,6 +123,9 @@ module com_mod
integer :: lnetcdfout
! lnetcdfout 1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input)
integer :: linversionout
! linversionout 1 for one grid_time output file for each release containing all timesteps
integer :: nageclass,lage(maxageclass)
! nageclass number of ageclasses for the age spectra calculation
......@@ -175,7 +178,7 @@ module com_mod
real :: vset(maxspec,ni),schmi(maxspec,ni),fract(maxspec,ni)
real :: ri(5,numclass),rac(5,numclass),rcl(maxspec,5,numclass)
real :: rgs(maxspec,5,numclass),rlu(maxspec,5,numclass)
real :: rm(maxspec),dryvel(maxspec)
real :: rm(maxspec),dryvel(maxspec),kao(maxspec)
real :: ohcconst(maxspec),ohdconst(maxspec),ohnconst(maxspec)
real :: area_hour(maxspec,24),point_hour(maxspec,24)
......@@ -360,7 +363,9 @@ module com_mod
real :: clwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !liquid [kg/kg]
real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !ice [kg/kg]
real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !combined [m3/m3]
! RLT add pressure and dry air density
real :: prs(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: rho_dry(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: pv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: rho(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: drhodz(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
......@@ -382,6 +387,7 @@ module com_mod
! uu,vv,ww [m/2] wind components in x,y and z direction
! uupol,vvpol [m/s] wind components in polar stereographic projection
! tt [K] temperature data
! prs air pressure
! qv specific humidity data
! pv (pvu) potential vorticity
! rho [kg/m3] air density
......
src/concoutput.f90
View file @ 2eefa58a
......@@ -71,6 +71,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,
......@@ -105,6 +108,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
......@@ -202,6 +206,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
......@@ -213,8 +220,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
......@@ -352,7 +365,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! Concentration output
!*********************
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
! Wet deposition
sp_count_i=0
......@@ -613,6 +626,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
......@@ -639,7 +695,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
!*************************
......
src/concoutput_inversion.f90 0 → 100644
View file @ 2eefa58a
This diff is collapsed.
src/concoutput_inversion_nest.f90 0 → 100644
View file @ 2eefa58a
This diff is collapsed.
src/concoutput_nest.f90
View file @ 2eefa58a
......@@ -69,6 +69,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,
......@@ -95,6 +98,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
......@@ -163,6 +167,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
......@@ -174,8 +181,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
......@@ -550,6 +563,47 @@ subroutine concoutput_nest(itime,outnum)
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
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
......
src/concoutput_surf.f90
View file @ 2eefa58a
......@@ -71,6 +71,9 @@ subroutine concoutput_surf(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,
......@@ -100,6 +103,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
logical :: sp_zer
character :: adate*8,atime*6
character(len=3) :: anspec
logical :: lexist
if (verbosity.eq.1) then
......@@ -179,6 +183,9 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
jjy=max(min(nint(yl),nymin1),0)
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
rho(iix,jjy,kzz-1,2)*dz2)/dz
! RLT
densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
end do
end do
end do
......@@ -189,8 +196,14 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
iix=max(min(nint(xl),nxmin1),0)
jjy=max(min(nint(yl),nymin1),0)
densityoutrecept(i)=rho(iix,jjy,1,2)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,2)
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
......@@ -595,6 +608,49 @@ subroutine concoutput_surf(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,1
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
......@@ -621,7 +677,23 @@ subroutine concoutput_surf(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
!*************************
......
src/concoutput_surf_nest.f90
View file @ 2eefa58a
......@@ -69,6 +69,9 @@ subroutine concoutput_surf_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,
......@@ -95,7 +98,7 @@ subroutine concoutput_surf_nest(itime,outnum)
logical :: sp_zer
character :: adate*8,atime*6
character(len=3) :: anspec
logical :: lexist
! Determine current calendar date, needed for the file name
!**********************************************************
......@@ -158,6 +161,9 @@ subroutine concoutput_surf_nest(itime,outnum)
jjy=max(min(nint(yl),nymin1),0)
densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
rho(iix,jjy,kzz-1,2)*dz2)/dz
! RLT
densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
end do
end do
end do
......@@ -168,8 +174,14 @@ subroutine concoutput_surf_nest(itime,outnum)
iix=max(min(nint(xl),nxmin1),0)
jjy=max(min(nint(yl),nymin1),0)
densityoutrecept(i)=rho(iix,jjy,1,2)
! RLT
densitydryrecept(i)=rho_dry(iix,jjy,1,2)
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
......@@ -316,8 +328,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgrid) sp_count_r
! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
! Dry deposition
sp_count_i=0
......@@ -353,8 +365,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgrid) sp_count_r
! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
......@@ -398,8 +410,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgrid) sp_count_r
! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
else
! write full vertical resolution
......@@ -439,8 +451,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgrid) sp_count_r
write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgrid) sp_count_r
! write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
endif ! surf_only
......@@ -486,8 +498,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgridppt) sp_count_r
! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! Dry deposition
......@@ -525,8 +537,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgridppt) sp_count_r
! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! Mixing ratios
......@@ -569,8 +581,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgridppt) sp_count_r
! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
else
! write full vertical resolution
......@@ -610,8 +622,8 @@ subroutine concoutput_surf_nest(itime,outnum)
write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
write(unitoutgridppt) sp_count_r
write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
! write(unitoutgridppt) sp_count_r
! write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
endif ! surf_only
endif ! output for ppt
......@@ -624,6 +636,48 @@ subroutine concoutput_surf_nest(itime,outnum)
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,1
do jy=0,numygridn-1
do ix=0,numxgridn-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*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
......
src/getfields.f90
View file @ 2eefa58a
......@@ -86,6 +86,10 @@ subroutine getfields(itime,nstop,metdata_format)
real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
! RLT added partial pressure water vapor
real :: pwater(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
integer :: kz, ix
character(len=100) :: rowfmt
integer :: indmin = 1
......@@ -152,9 +156,9 @@ subroutine getfields(itime,nstop,metdata_format)
end do
40 indmin=indj
if (WETBKDEP) then
call writeprecip(itime,memind(1))
endif
if (WETBKDEP) then
call writeprecip(itime,memind(1))
endif
else
......@@ -203,11 +207,34 @@ subroutine getfields(itime,nstop,metdata_format)
end do
60 indmin=indj
if (WETBKDEP) then
call writeprecip(itime,memind(1))
endif
endif
if (WETBKDEP) then
call writeprecip(itime,memind(1))
endif
end if
! RLT calculate dry air density
pwater=qv*prs/((r_air/r_water)*(1.-qv)+qv)
rho_dry=(prs-pwater)/(r_air*tt)
! test density
! write(rowfmt,'(A,I6,A)') '(',nymax,'(E11.4,1X))'
! if(itime.eq.0) then
! open(500,file=path(2)(1:length(2))//'rho_dry.txt',status='replace',action='write')
! do kz=1,nzmax
! do ix=1,nxmax
! write(500,fmt=rowfmt) rho_dry(ix,:,kz,1)
! end do
! end do
! close(500)
! open(500,file=path(2)(1:length(2))//'rho.txt',status='replace',action='write')
! do kz=1,nzmax
! do ix=1,nxmax
! write(500,fmt=rowfmt) rho(ix,:,kz,1)
! end do
! end do
! close(500)
! endif
lwindinterv=abs(memtime(2)-memtime(1))
......
src/initial_cond_output_inversion.f90 0 → 100644
View file @ 2eefa58a
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! 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 <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine initial_cond_output_inversion(itime)
! i
!*****************************************************************************
! *
! Output of the initial condition sensitivity field. *
! *
! Author: A. Stohl *
! *
! 24 May 1995 *
! *
! 13 April 1999, Major update: if output size is smaller, dump output *
! in sparse matrix format; additional output of *
! uncertainty *
! *
! 05 April 2000, Major update: output of age classes; output for backward*
! runs is time spent in grid cell times total mass of *
! species. *
! *
! 17 February 2002, Appropriate dimensions for backward and forward runs *
! are now specified in file par_mod *
! *
! June 2006, write grid in sparse matrix with a single write command *
! in order to save disk space *
! *
! 2008 new sparse matrix format *
! *
!*****************************************************************************
!