...
 
Commits (18)
&SPECIES_PARAMS
PSPECIES="C-CO2 ",
PWEIGHTMOLAR= 12.0000000 ,
/
GRIB1| Level| flexpart | flexpart | flexpart |GRIB2|GRIB2|GRIB2|GRIB2|
Param| Type | Name | Units | Description |Discp|Catgy|Param|Level|
-----+------+----------+----------+------------------------+-----------------------+
130 | 109 | TT | K | Temperature | 0 | 0 | 0 | 105 |
131 | 109 | UU | m s-1 | U | 0 | 2 | 2 | 105 |
132 | 109 | VV | m s-1 | V | 0 | 2 | 3 | 105 |
133 | 109 | QV | kg kg-1 | Specific humidity | 0 | 1 | 0 | 105 |
134 | 1 | PS | Pa | Surface pressure | 0 | 3 | 0 | 1 |
135 | 109 | ETADOT | Pa s-1 | Vertical velocity | 0 | 2 | 32 | 105 |
141 | 1 | SD | m | Snow depth | 0 | 1 | 11 | 1 |
151 | 1 | MSL | Pa | Sea-level Pressure | 0 | 3 | 0 | 101 |
164 | 1 | TCC | (0-1) | Total cloud cover | 192 | 128 | 164 | 1 |
165 | 1 | U10 | m s-1 | 10m U | 0 | 2 | 2 | 103 |
166 | 1 | V10 | m s-1 | 10m V | 0 | 2 | 3 | 103 |
167 | 1 | T2 | K | 2m Temperature | 0 | 0 | 0 | 103 |
168 | 1 | TD2 | K | 2m Dewpoint | 0 | 0 | 6 | 103 |
142 | 1 | LSPREC | m | Large scale precip | 192 | 128 | 142 | 1 |
143 | 1 | CONVPREC | m | Convective precip | 0 | 1 | 10 | 1 |
146 | 1 | SHF | J m-2 | Sensible heat flux | 0 | 0 | 11 | 1 |
176 | 1 | SR | J m-2 | Solar radiation | 0 | 4 | 9 | 1 |
180 | 1 | EWSS | N m-2 s | EW turb surf stress | 0 | 2 | 38 | 1 |
181 | 1 | NSSS | N m-2 s | NS turb surf stress | 0 | 2 | 37 | 1 |
129 | 1 | ORO | m2 s2??? | Surface geopotential | 0 | 3 | 4 | 1 |
160 | 1 | EXCESSORO| none | Std dev of orography | 192 | 128 | 160 | 1 |
172 | 1 | LSM | 0/1 Flag | Land/Sea flag | 2 | 0 | 0 | 1 |
246 | 109 | CLWC | kg/kg | Cloud Liquid Water | 0 | 1 | 83 | 105 |
247 | 109 | CIWC | kg/kg | Cloud Ice Water | 0 | 1 | 84 | 105 |
201031| 109 | QC | kg/kg | Cloud Total Water | 192 | 201 | 31 | 105 |
-----+------+----------+----------+------------------------+-----------------------+
GRIB1| Level| flexpart | flexpart | flexpart |GRIB2|GRIB2|GRIB2|GRIB2|
Param| Type | Name | Units | Description |Discp|Catgy|Param|Level|
-----+------+----------+----------+------------------------+-----------------------+
11 | 100 | TT | K | Temperature | 0 | 0 | 0 | 100 |
33 | 100 | UU | m s-1 | U | 0 | 2 | 2 | 100 |
34 | 100 | VV | m s-1 | V | 0 | 2 | 3 | 100 |
52 | 100 | RH | % | Relative humidity | 0 | 1 | 1 | 100 |
1 | 1 | PS | Pa | Surface pressure | 0 | 3 | 0 | 1 |
39 | 100 | WW | Pa s-1 | Vertical velocity | 0 | 2 | 8 | 100 |
65 | 1 | SD | kg m-2 | Snow depth water equiv | 0 | 1 | 13 | 1 |
2 | 102 | SLP | Pa | Mean sea level press | 0 | 3 | 1 | 101 |
33 | 105 | U10 | m s-1 | 10m U | 0 | 2 | 2 | 103 |
34 | 105 | V10 | m s-1 | 10m V | 0 | 2 | 3 | 103 |
11 | 105 | T2 | K | 2m Temperature | 0 | 0 | 0 | 103 |
* | * | TD2 | K | 2m Dewpoint | 0 | 0 | 6 | 103 |
71 | 244 | TCC | % | Total cloud cover | 0 | 6 | 1 | 244 |
62 | 1 | LSPREC | kg m-2 | Total precip | 0 | 1 | 7 | 1 |
63 | 1 | CONVPREC | kg m-2 | Convective precip | 0 | 1 | 196 | 1 |
7 | 1 | ORO | m | Orography/geopotential | 0 | 3 | 5 | 1 |
81 | 1 | LSM | 0/1 Flag | Land/Sea flag | 2 | 0 | 0 | 1 |
221 | 1 | HMIX | m | Mixing height (BLH) | 0 | 3 | 196 | 1 |
52 | 105 | RH2 | % | Relative humidity | 0 | 1 | 1 | 103 |
11 | 107 | TSIG1 | K | Temperature (sigma) | 0 | 0 | 0 | 104 |
33 | 107 | USIG1 | m s-1 | U wind (sigma) | 0 | 2 | 2 | 104 |
34 | 107 | VSIG1 | m s-1 | V wind (sigma) | 0 | 2 | 3 | 104 |
-----+------+----------+----------+------------------------+-----------------------+
......@@ -138,7 +138,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)
......
......@@ -150,7 +150,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)
......
......@@ -221,6 +221,11 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
endif
ixp=ix+1
jyp=jy+1
! fix (Espen 02.05.2017)
if(jyp>=nymax) then
write(*,*) 'WARNING: advance.f90 jyp>nymax'
jyp=jyp-1
endif
! Determine the lower left corner and its distance to the current position
......
......@@ -32,7 +32,7 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! 21 May 1995 *
! *
! ------------------------------------------------------------------ *
! Petra Seibert, Feb 2000: *
! Petra Seibert, Feb 2000: *
! convection scheme: *
! new variables in call to richardson *
! *
......
......@@ -120,6 +120,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
......@@ -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
......@@ -652,6 +658,7 @@ module com_mod
real :: xreceptor(maxreceptor),yreceptor(maxreceptor)
real :: receptorarea(maxreceptor)
real :: creceptor(maxreceptor,maxspec)
real, allocatable, dimension(:,:) :: creceptor0
character(len=16) :: receptorname(maxreceptor)
integer :: numreceptor
......
......@@ -101,6 +101,11 @@ subroutine conccalc(itime,weight)
jy=int(ytra1(i))
ixp=ix+1
jyp=jy+1
! fix (Espen 02.05.2017)
if(jyp>=nymax) then
write(*,*) 'WARNING: conccalc.f90 jyp>nymax'
jyp=jyp-1
endif
ddx=xtra1(i)-real(ix)
ddy=ytra1(i)-real(jy)
rddx=1.-ddx
......@@ -132,7 +137,7 @@ subroutine conccalc(itime,weight)
! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
!*****************************************************************************
do ind=indz,indzp
rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,memind(2)) &
rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
+p2*rho(ixp,jy ,ind,2) &
+p3*rho(ix ,jyp,ind,2) &
+p4*rho(ixp,jyp,ind,2)
......
......@@ -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
......@@ -362,7 +375,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
if ((ldirect.eq.1).and.(WETDEP)) then
do jy=0,numygrid-1
do ix=0,numxgrid-1
!oncentraion greater zero
!concentraion greater zero
if (wetgrid(ix,jy).gt.smallnum) then
if (sp_zer.eqv..true.) then ! first non zero value
sp_count_i=sp_count_i+1
......@@ -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
!*************************
......
This diff is collapsed.
This diff is collapsed.
......@@ -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
......
......@@ -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
......@@ -459,6 +472,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
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)
endif ! concentration output
......@@ -503,6 +517,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
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)
......@@ -595,6 +610,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 +679,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
!*************************
......
......@@ -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,18 +161,27 @@ 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
do i=1,numreceptor
xl=xreceptor(i)
yl=yreceptor(i)
iix=max(min(nint(xl),nxmin1),0)
jjy=max(min(nint(yl),nymin1),0)
densityoutrecept(i)=rho(iix,jjy,1,2)
end do
do i=1,numreceptor
xl=xreceptor(i)
yl=yreceptor(i)
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
......
......@@ -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
......@@ -203,12 +207,35 @@ subroutine getfields(itime,nstop,metdata_format)
end do
60 indmin=indj
if (WETBKDEP) then
call writeprecip(itime,memind(1))
endif
if (WETBKDEP) then
call writeprecip(itime,memind(1))
endif
endif
! 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))
if (lwindinterv.gt.idiffmax) nstop=3
......
......@@ -41,16 +41,25 @@ subroutine gethourlyOH(itime)
implicit none
integer :: itime
integer :: ix,jy,kz,m1,m2
integer :: ix,jy,kz,m1,m2,d1,d2,m,dd,ndays
integer :: ijx,jjy
integer :: jjjjmmdd,hhmmss
real :: sza,jrate,photo_O1D,zenithangle
real :: sza,jrate,photo_O1D,zenithangle,OH_tmp,jrate_tmp
real(kind=dp) :: jul1,jul2
integer, dimension(12) :: monthdays
monthdays=(/31,28,31,30,31,30,31,31,30,31,30,31/)
! print*, 'itime: ',itime
! print*, 'memOHtime(1):',memOHtime(1)
! print*, 'memOHtime(2):',memOHtime(2)
!! testing OH
! open(1000,file=path(2)(1:length(2))//'OH_50N_0E.txt',action='write',status='old',access='append')
! open(999,file=path(2)(1:length(2))//'jrate_50N_0E.txt',action='write',status='old',access='append')
! open(998,file=path(2)(1:length(2))//'OH_50S_0E.txt',action='write',status='old',access='append')
! open(997,file=path(2)(1:length(2))//'jrate_50S_0E.txt',action='write',status='old',access='append')
! Check hourly OH field is available for the current time step
!**************************************************************
......@@ -79,9 +88,13 @@ subroutine gethourlyOH(itime)
jul2=bdate+memOHtime(2)/86400._dp ! date for next hour
call caldate(jul2,jjjjmmdd,hhmmss)
m2=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
dd=(jjjjmmdd-(jjjjmmdd/100)*100)
ndays=monthdays(m2)
! print*, 'jul2:',jul2
! print*, 'm2:',m2
! print*, 'dd:',dd
! print*, 'ndays:',ndays
do kz=1,nzOH
do jy=1,nyOH
......@@ -92,19 +105,33 @@ subroutine gethourlyOH(itime)
sza=zenithangle(latOH(jy),lonOH(ix),jul2)
! calculate J(O1D) (jrate)
jrate=photo_O1D(sza)
! interpolate monthly fields to current day
if(dd.ge.(ndays/2)) then
m=min((m2+1),12)
OH_tmp=OH_field(ix,jy,kz,m2)+(OH_field(ix,jy,kz,m)-OH_field(ix,jy,kz,m2))*real(dd-ndays/2)/real(ndays)
jrate_tmp=jrate_average(ijx,jjy,m2)+(jrate_average(ijx,jjy,m)-jrate_average(ijx,jjy,m2))*real(dd-ndays/2)/real(ndays)
else
m=max((m2-1),1)
OH_tmp=OH_field(ix,jy,kz,m)+(OH_field(ix,jy,kz,m2)-OH_field(ix,jy,kz,m))*real(dd+ndays/2)/real(ndays)
jrate_tmp=jrate_average(ijx,jjy,m)+(jrate_average(ijx,jjy,m2)-jrate_average(ijx,jjy,m))*real(dd+ndays/2)/real(ndays)
endif
! apply hourly correction to OH
if(jrate_average(ijx,jjy,m2).gt.0.) then
OH_hourly(ix,jy,kz,2)=OH_field(ix,jy,kz,m2)*jrate/jrate_average(ijx,jjy,m2)
if(jrate_tmp.gt.0.) then
OH_hourly(ix,jy,kz,2)=OH_tmp*jrate/jrate_tmp
else
OH_hourly(ix,jy,kz,2)=0.
endif
!! for testing !!
! if(jy.eq.36.and.ix.eq.36.and.kz.eq.1) then
! write(999,fmt='(F6.3)') jrate/jrate_average(ijx,jjy,m2)
! endif
! if(jy.eq.11.and.ix.eq.36.and.kz.eq.1) then
! write(998,fmt='(F6.3)') jrate/jrate_average(ijx,jjy,m2)
! endif
!! testing OH
! if(jy.eq.36.and.ix.eq.37.and.kz.eq.1) then
! write(999,fmt='(E12.6)') OH_hourly(ix,jy,kz,2)
! write(1000,fmt='(E12.6)') OH_tmp
! write(999,fmt='(E12.6)') jrate/jrate_tmp
! endif
! if(jy.eq.11.and.ix.eq.37.and.kz.eq.1) then
! write(998,fmt='(E12.6)') OH_hourly(ix,jy,kz,2)
! write(998,fmt='(E12.6)') OH_tmp
! write(997,fmt='(E12.6)') jrate/jrate_tmp
! endif
end do
end do
end do
......@@ -118,15 +145,18 @@ subroutine gethourlyOH(itime)
call caldate(jul1,jjjjmmdd,hhmmss)
m1=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
memOHtime(1)=0.
d1=(jjjjmmdd-(jjjjmmdd/100)*100)
jul2=bdate+ldirect*real(1./24.,kind=dp) ! date for next hour
call caldate(jul2,jjjjmmdd,hhmmss)
m2=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
memOHtime(2)=ldirect*3600.
d2=(jjjjmmdd-(jjjjmmdd/100)*100)
! print*, 'jul1:',jul1
! print*, 'jul2:',jul2
! print*, 'm1,m2:',m1,m2
! print*, 'd1,d2:',d1,d2
do kz=1,nzOH
do jy=1,nyOH
......@@ -137,9 +167,18 @@ subroutine gethourlyOH(itime)
sza=zenithangle(latOH(jy),lonOH(ix),jul1)
! calculate J(O1D) (jrate), beginning
jrate=photo_O1D(sza)
! interpolate monthly fields to current day
ndays=monthdays(m1)
if(d1.ge.(ndays/2)) then
m=min((m1+1),12)
OH_tmp=OH_field(ix,jy,kz,m1)+(OH_field(ix,jy,kz,m)-OH_field(ix,jy,kz,m1))*real(d1-ndays/2)/real(ndays)
else
m=max((m1-1),1)
OH_tmp=OH_field(ix,jy,kz,m)+(OH_field(ix,jy,kz,m1)-OH_field(ix,jy,kz,m))*real(d1+ndays/2)/real(ndays)
endif
! apply hourly correction to OH
if(jrate_average(ijx,jjy,m1).gt.0.) then
OH_hourly(ix,jy,kz,1)=OH_field(ix,jy,kz,m1)*jrate/jrate_average(ijx,jjy,m1)
OH_hourly(ix,jy,kz,1)=OH_tmp*jrate/jrate_average(ijx,jjy,m1)
else
OH_hourly(ix,jy,kz,1)=0.
endif
......@@ -147,9 +186,18 @@ subroutine gethourlyOH(itime)
sza=zenithangle(latOH(jy),lonOH(ix),jul2)
! calculate J(O1D) (jrate), after 1-hour
jrate=photo_O1D(sza)
! interpolate monthly fields to current day
ndays=monthdays(m2)
if(d2.ge.(ndays/2)) then
m=min((m2+1),12)
OH_tmp=OH_field(ix,jy,kz,m2)+(OH_field(ix,jy,kz,m)-OH_field(ix,jy,kz,m2))*real(d2-ndays/2)/real(ndays)
else
m=max((m2-1),1)
OH_tmp=OH_field(ix,jy,kz,m)+(OH_field(ix,jy,kz,m2)-OH_field(ix,jy,kz,m))*real(d2+ndays/2)/real(ndays)
endif
! apply hourly correction to OH
if(jrate_average(ijx,jjy,m2).gt.0.) then
OH_hourly(ix,jy,kz,2)=OH_field(ix,jy,kz,m2)*jrate/jrate_average(ijx,jjy,m2)
OH_hourly(ix,jy,kz,2)=OH_tmp*jrate/jrate_average(ijx,jjy,m2)
else
OH_hourly(ix,jy,kz,2)=0.
endif
......@@ -159,5 +207,10 @@ subroutine gethourlyOH(itime)
endif
!! testing OH
! close(999)
! close(998)
end subroutine gethourlyOH
This diff is collapsed.
This diff is collapsed.
......@@ -42,6 +42,11 @@ subroutine gridcheck_gfs
! Marian Harustak, 12.5.2017 *
! - Renamed routine from gridcheck to gridcheck_gfs *
! *
! Implementation of the Vtables approach *
! D. Morton, D. Arnold 17.11.2017 *
! - Inclusion of specific code and usage of class_vtable *
! *
! *
!**********************************************************************
! *
! DESCRIPTION: *
......@@ -78,6 +83,7 @@ subroutine gridcheck_gfs
use com_mod
use conv_mod
use cmapf_mod, only: stlmbr,stcm2p
use class_vtable
implicit none
......@@ -101,7 +107,10 @@ subroutine gridcheck_gfs
! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
integer :: isec1(8),isec2(3)
!!!! DJM - isec1 no longer user - integer :: isec1(8),isec2(3)
integer :: isec2(3)
real(kind=4) :: zsec4(jpunp)
character(len=1) :: opt
......@@ -109,6 +118,25 @@ subroutine gridcheck_gfs
character(len=24) :: gribErrorMsg = 'Error reading grib file'
character(len=20) :: gribFunction = 'gridcheckwind_gfs'
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! Vtable related variables
!
! Path to Vtable - current implementation assumes it's in cwd, named
! "Vtable"
! ESO: Changed to use default Vtable file in options directory
! CHARACTER(LEN=255), PARAMETER :: VTABLE_PATH = "Vtable"
character(LEN=255) :: VTABLE_PATH
character(LEN=15) :: fpname ! stores FLEXPART name for curr grib mesg.
type(Vtable) :: my_vtable ! unallocated
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! DJM
integer current_grib_level ! this was isec1(8) in previous versions
if (numbnests.ge.1) then
write(*,*) ' ###########################################'
write(*,*) ' FLEXPART ERROR SUBROUTINE GRIDCHECK:'
......@@ -125,6 +153,38 @@ subroutine gridcheck_gfs
else
ifn=numbwf
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! Vtable code
VTABLE_PATH = path(1)(1:length(1))//'Vtables/Vtable.gfs'
PRINT *, 'Loading Vtable: ', VTABLE_PATH
call vtable_load_by_name(VTABLE_PATH, my_vtable)
!! Debugging tool
! PRINT *, 'Dump of Vtable...'
! call vtable_dump_records(my_vtable)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!! VTABLE code
! This is diagnostic/debugging code, and will normally be commented out.
! It's purpose is to look at the provided grib file and produce an
! inventory of the FP-related messages, relative to the Vtable that's
! already been open.
! CALL vtable_gribfile_inventory(path(3)(1:length(3)) // trim(wfname(ifn)), &
!& my_vtable)
!!!!!!!!!!!!!!!!!!! VTABLE code
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! OPENING OF DATA FILE (GRIB CODE)
!
......@@ -148,66 +208,61 @@ subroutine gridcheck_gfs
goto 999 ! ERROR DETECTED
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!! VTABLE code
! Get the fpname
fpname = vtable_get_fpname(igrib, my_vtable)
!print *, 'fpname: ', trim(fpname)
!!!!!!!!!!!!!!!!!!! VTABLE code
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!first see if we read GRIB1 or GRIB2
call grib_get_int(igrib,'editionNumber',gribVer,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
if (gribVer.eq.1) then ! GRIB Edition 1
!read the grib1 identifiers
call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'level',isec1(8),iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'level',current_grib_level,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
!get the size and data of the values array
call grib_get_real4_array(igrib,'values',zsec4,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
!!! Added by DJM 2017-08-04 - if this is GRIB1 we assume that
!!! level units are hPa and need to be multiplied by 100 for Pa
!!! We only do this for the UU field, since that's the only 3D field
!!! used in this routine
current_grib_level = current_grib_level*100.0
else ! GRIB Edition 2
!read the grib2 identifiers
call grib_get_int(igrib,'discipline',discipl,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'parameterCategory',parCat,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'parameterNumber',parNum,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
valSurf,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
!convert to grib1 identifiers
isec1(6)=-1
isec1(7)=-1
isec1(8)=-1
if ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
isec1(6)=33 ! indicatorOfParameter
isec1(7)=100 ! indicatorOfTypeOfLevel
isec1(8)=valSurf/100 ! level, convert to hPa
elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
isec1(6)=7 ! indicatorOfParameter
isec1(7)=1 ! indicatorOfTypeOfLevel
isec1(8)=0
elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
.and.(discipl.eq.2)) then ! LSM
isec1(6)=81 ! indicatorOfParameter
isec1(7)=1 ! indicatorOfTypeOfLevel
isec1(8)=0
endif
call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
current_grib_level,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
if (isec1(6).ne.-1) then
! get the size and data of the values array
call grib_get_real4_array(igrib,'values',zsec4,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
endif
!!! Added by DJM 2017-08-04 - if this is GRIB2 we assume that
!!! level units are Pa and don't need to be modified
endif ! gribVer
IF (TRIM(fpname) .NE. 'NOFP') THEN
!get the size and data of the values array
call grib_get_real4_array(igrib,'values',zsec4,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
END IF
if(ifield.eq.1) then
!get the required fields from section 2
......@@ -231,8 +286,13 @@ subroutine gridcheck_gfs
yaux2in,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
! Fix for flexpart.eu ticket #48
if (xaux2in.lt.0) xaux2in = 359.0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!! DJM ARTIFICIAL CHANGE FOR NCEP GRIB1 - change value from -1 to 359
!!!!!!!! See flexpart.eu CTBTO project Ticket #112
!!!!!!!! Also see flexpart.eu mainstream Ticket #48
if (xaux2in .lt. 0) xaux2in = 359.0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
xaux1=xaux1in
xaux2=xaux2in
......@@ -316,9 +376,12 @@ subroutine gridcheck_gfs
! NCEP ISOBARIC LEVELS
!*********************
if((isec1(6).eq.33).and.(isec1(7).eq.100)) then ! check for U wind
!!!! DJM - orig - if((isec1(6).eq.33).and.(isec1(7).eq.100)) then ! check for U wind
IF (TRIM(fpname) .EQ. 'UU') THEN ! check for U wind
iumax=iumax+1
pres(iumax)=real(isec1(8))*100.0
!!!! DJM - orig - pres(iumax)=real(isec1(8))*100.0
pres(iumax) = REAL(current_grib_level)
endif
......@@ -334,7 +397,8 @@ subroutine gridcheck_gfs
! NCEP TERRAIN
!*************
if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
!!!!!! DJM - orig - if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
IF (TRIM(fpname) .EQ. 'ORO') THEN
do jy=0,ny-1
do ix=0,nxfield-1
help=zsec4(nxfield*(ny-jy-1)+ix+1)
......@@ -352,7 +416,8 @@ subroutine gridcheck_gfs
! NCEP LAND SEA MASK
!*******************
if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
!!!!!! DJM - orig - if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
IF (TRIM(fpname) .EQ. 'LSM') THEN
do jy=0,ny-1
do ix=0,nxfield-1
help=zsec4(nxfield*(ny-jy-1)+ix+1)
......
This diff is collapsed.
!**********************************************************************
! 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 *
! *
!*****************************************************************************
! *
! Variables: *
! ncells number of cells with non-zero concentrations *
! sparse .true. if in sparse matrix format, else .false. *
! *
!*****************************************************************************
use unc_mod
use point_mod
use outg_mod
use par_mod
use com_mod
implicit none
integer :: itime,i,ix,jy,kz,ks,kp,sp_count_i,sp_count_r
integer :: jjjjmmdd, ihmmss
real(kind=dp) :: jul
real :: sp_fact,fact_recept
real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
logical :: sp_zer,lexist
logical,save :: listart=.true.
logical,save,allocatable,dimension(:) :: listartrel
character :: adate*8,atime*6
character :: areldate*8,areltime*6
character(len=3) :: anspec
if(listart) then
allocate(listartrel(maxpointspec_act))
listartrel(:)=.true.
endif
print*, 'listartrel = ',listartrel
!*********************************************************************
! Determine the standard deviation of the mean concentration or mixing
! ratio (uncertainty of the output) and the dry and wet deposition
!*********************************************************************
do ks=1,nspec
write(anspec,'(i3.3)') ks
do kp=1,maxpointspec_act
! calculate date of release
jul=bdate+real(ireleasestart(kp),kind=dp)/86400._dp ! this is the current day
call caldate(jul,jjjjmmdd,ihmmss)
write(areldate,'(i8.8)') jjjjmmdd
write(areltime,'(i6.6)') ihmmss
print*, areldate//areltime
! calculate date of field
jul=bdate+real(itime,kind=dp)/86400._dp
call caldate(jul,jjjjmmdd,ihmmss)
write(adate,'(i8.8)') jjjjmmdd
write(atime,'(i6.6)') ihmmss
print*, adate//atime
inquire(file=path(2)(1:length(2))//'grid_initial_'//areldate// &
areltime//'_'//anspec,exist=lexist)
if(lexist.and..not.listartrel(kp)) then
! open and append to existing file
open(97,file=path(2)(1:length(2))//'grid_initial_'//areldate// &
areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
else
! open new file
open(97,file=path(2)(1:length(2))//'grid_initial_'//areldate// &
areltime//'_'//anspec,form='unformatted',status='replace',action='write')
endif
write(97) jjjjmmdd
write(97) ihmmss
listartrel(kp)=.false.
if (ind_rel.eq.1) then
fact_recept=rho_rel(kp)
else
fact_recept=1.
endif
!*******************************************************************
! Generate output: may be in concentration (ng/m3) or in mixing
! ratio (ppt) or both
! Output the position and the values alternated multiplied by
! 1 or -1, first line is number of values, number of positions
! For backward simulations, the unit is seconds, stored in grid_time
!*******************************************************************
! Write out dummy "wet and dry deposition" fields, to keep same format
! as for concentration output
! sp_count_i=0
! sp_count_r=0
! write(97) sp_count_i
! write(97) (sparse_dump_i(i),i=1,sp_count_i)
! write(97) sp_count_r
! write(97) (sparse_dump_r(i),i=1,sp_count_r)
! write(97) sp_count_i
! write(97) (sparse_dump_i(i),i=1,sp_count_i)
! write(97) sp_count_r
! write(97) (sparse_dump_r(i),i=1,sp_count_r)
! Write out sensitivity to initial conditions
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 (init_cond(ix,jy,kz,ks,kp).gt.smallnum) then
if (sp_zer.eqv..true.) then ! first non zero value
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* &
init_cond(ix,jy,kz,ks,kp)/xmass(kp,ks)*fact_recept
else ! concentration is zero
sp_zer=.true.
endif
end do
end do
end do
write(97) sp_count_i
write(97) (sparse_dump_i(i),i=1,sp_count_i)
write(97) sp_count_r
write(97) (sparse_dump_r(i),i=1,sp_count_r)
close(97)
end do
end do
! reset listart
if (listart) then
listart=.false.
endif
end subroutine initial_cond_output_inversion
......@@ -36,6 +36,8 @@ SHELL = /bin/bash
#
################################################################################
#FPPFLAGS+=-DUSE_MPIINPLACE
## PROGRAMS
# Unified executable names
# The same executable is used for both ECMWF and GFS metdata
......@@ -50,24 +52,24 @@ FLEXPART-MPI-DBG = DBG_FLEXPART_MPI
FLEXPART-SERIAL = FLEXPART
ifeq ($(gcc), 4.9)
# Compiled libraries under user ~flexpart, gfortran v4.9
ROOT_DIR = /homevip/flexpart/
F90 = ${ROOT_DIR}/gcc-4.9.1/bin/gfortran
MPIF90 = ${ROOT_DIR}/bin/mpifort
F90 = gfortran
MPIF90 = mpif90
INCPATH1 = ${ROOT_DIR}/gcc-4.9.1/include
INCPATH2 = ${ROOT_DIR}/include
LIBPATH1 = ${ROOT_DIR}/lib
else
# Compiled libraries under user ~flexpart, gfortran v5.4
ROOT_DIR = /homevip/flexpart/
F90 = /usr/bin/gfortran
MPIF90 = /usr/bin/mpifort
else
ROOT_DIR = /homevip/flexpart/
F90 = /usr/bin/gfortran
MPIF90 = /usr/bin/mpif90