diff --git a/src/FLEXPART.f90 b/src/FLEXPART.f90
index d1bdc211420b495c3fe24f05612ae058673187af..ec735eb06f5f63bdc451d7e6fc66d9db4f37d2d7 100644
--- a/src/FLEXPART.f90
+++ b/src/FLEXPART.f90
@@ -177,7 +177,7 @@ subroutine read_options_and_initialise_flexpart
#ifdef USE_NCF
use chemistry_mod, only: readreagents
use totals_mod
- use receptor_netcdf_mod, only: readreceptors_satellite, receptorout_init, satelliteout_init
+ use receptor_netcdf_mod, only: read_satellite_info, receptorout_init
#endif
use receptor_mod, only: alloc_receptor
@@ -301,7 +301,7 @@ subroutine read_options_and_initialise_flexpart
numsatreceptor=0
nlayermax=1
#ifdef USE_NCF
- call readreceptors_satellite
+ call read_satellite_info
#endif
call readreceptors
@@ -366,11 +366,9 @@ subroutine read_options_and_initialise_flexpart
if (lnetcdfout.eq.1) then
#ifdef USE_NCF
call receptorout_init
- call satelliteout_init
#endif
else
call receptorout_init_binary
- call satelliteout_init_binary
endif
if ((iout.eq.4).or.(iout.eq.5)) call openouttraj ! CHECK ETA
diff --git a/src/binary_output_mod.f90 b/src/binary_output_mod.f90
index b19fa1dd13c07d1bc92421626b8c3688fd32c3ba..49d23a486b145b1a55569ad2da207433d32e3074 100644
--- a/src/binary_output_mod.f90
+++ b/src/binary_output_mod.f90
@@ -708,7 +708,7 @@ subroutine satelliteout_init_binary
else
open(unitoutsatellite,file=path(2)(1:length(2))//'satellite_pptv', &
form='unformatted',err=998)
- write(unitoutsatellite) numsatreceptor
+! write(unitoutsatellite) numsatreceptor
if (llcmoutput) then
write(unitoutsatellite) nspec-1 ! first species is mass of air
else
@@ -799,7 +799,7 @@ subroutine write_satellite_binary(crec,cunc,nnrec,xkrec,lonrec,latrec,altrec,tim
real, dimension(maxrecsample,nlayermax) :: nnrec, xkrec, altrec
real, dimension(maxrecsample) :: lonrec, latrec
integer, dimension(maxrecsample) :: timerec
- character(len=16), dimension(maxrecsample) :: namerec
+ character(len=24), dimension(maxrecsample) :: namerec
if (llcmoutput) then
ks_start=2
diff --git a/src/chemistry_mod.f90 b/src/chemistry_mod.f90
index 3bd7f844567dcbe5c1e9f5a5ff829359533018e8..45e08927129bca3ca1191b72a2507a03b65ff925 100644
--- a/src/chemistry_mod.f90
+++ b/src/chemistry_mod.f90
@@ -155,6 +155,7 @@ module chemistry_mod
integer :: itime
integer :: jjjjmmdd, hhmmss, mm, eomday
integer :: nr, memid
+ real(kind=dp) :: jd, jdmid
character(len=2) :: amonth
character(len=256):: CL_name, jr_name
logical :: lexist
@@ -177,30 +178,35 @@ module chemistry_mod
! Current time is after 2nd chem field
!*************************************
+ write(*,*) 'Updating CL fields... '
+
+ memCLtime(1)=memCLtime(2) ! time in sec
+ CL_time(1)=CL_time(2) ! julian date
+ memid=2
+
+ ! Get date/time of next chem field
+ !*********************************
+ ! assumes fields are monthly
+
+ call caldate(CL_time(1), jjjjmmdd, hhmmss)
+ eomday=calceomday(jjjjmmdd/100)
+ memCLtime(2)=memCLtime(1)+ldirect*eomday*24*3600 ! time in sec
+ CL_time(2)=CL_time(1)+real(ldirect*eomday,kind=dp) ! julian date
+ write(*,*) 'getchemfield: memCLtime = ',memCLtime(1), memCLtime(2)
+ write(*,*) 'getchemfield: CL_time = ',CL_time(1), CL_time(2)
+ call caldate(CL_time(2), jjjjmmdd,hhmmss)
+ mm=int((jjjjmmdd-(jjjjmmdd/10000)*10000)/100)
+ write(amonth,'(I2.2)') mm
+
do nr=1,nreagent
write(*,*) 'Updating CL field for: ',trim(reagents(nr))
CL_field(:,:,:,nr,1)=CL_field(:,:,:,nr,2)
- memCLtime(1)=memCLtime(2) ! time in sec
- CL_time(1)=CL_time(2) ! julian date
- memid=2
! Read new chem field and store in 2nd position
!**********************************************
- ! julian date of next chem field assuming monthly fields
- call caldate(CL_time(1), jjjjmmdd, hhmmss)
- eomday=calceomday(jjjjmmdd/100)
- memCLtime(2)=memCLtime(1)+ldirect*eomday*24*3600 ! time in sec
- CL_time(2)=CL_time(1)+real(ldirect*eomday,kind=dp) ! julian date
- !! test
- write(*,*) 'memCLtime = ',memCLtime(1), memCLtime(2)
- write(*,*) 'CL_time = ',CL_time(1), CL_time(2)
- !!
- call caldate(CL_time(2), jjjjmmdd,hhmmss)
- mm=int((jjjjmmdd-(jjjjmmdd/10000)*10000)/100)
- write(amonth,'(I2.2)') mm
write(CL_name,'(A)') trim(reag_path(nr))//trim(reagents(nr))//'_'//amonth//'.nc'
inquire(file=CL_name,exist=lexist)
if (lexist) then
@@ -239,22 +245,63 @@ module chemistry_mod
! No chem field in memory that can be used
!******************************************
- do nr=1,nreagent
-
- write(*,*) 'Reading two new CL fields for: ',trim(reagents(nr))
-
- ! read both fields into memory
- do memid=1,2
- if (memid.eq.1) then
- CL_time(memid)=bdate+real(ldirect*itime,kind=dp)/86400._dp
+ write(*,*) 'Reading two new CL fields...'
+
+ ! Get date/time of both chem fields
+ !**********************************
+ ! assumes fields are monthly
+
+ do memid=1,2
+ if (memid.eq.1) then
+ jd=bdate+real(ldirect*itime,kind=dp)/86400._dp
+ call caldate(jd, jjjjmmdd, hhmmss)
+ ! middle of month day
+ jdmid=juldate(int(jjjjmmdd/100)*100+15,0)
+ !! testing
+ print*, 'getchemfield: jjjjmmdd, jjjjmmdd_mid = ',jjjjmmdd, int(jjjjmmdd/100)*100+15
+ ! julian date of fields
+ if (ldirect.gt.0) then
+ if (jd.ge.jdmid) then
+ ! use current month
+ CL_time(memid)=jdmid
+ else
+ ! use last month
+ eomday=calceomday(jjjjmmdd/100)
+ CL_time(memid)=jdmid-real(eomday,kind=dp)
+ endif
else
- CL_time(memid)=CL_time(memid-1)+real(ldirect*eomday,kind=dp)
+ if (jd.ge.jdmid) then
+ ! use next month
+ CL_time(memid)=jdmid+real(eomday,kind=dp)
+ else
+ ! use current month
+ CL_time(memid)=jdmid
+ endif
endif
- memCLtime(memid)=int((CL_time(memid)-bdate)*86400._dp)*ldirect
- call caldate(CL_time(memid), jjjjmmdd, hhmmss)
- mm=int((jjjjmmdd-(jjjjmmdd/10000)*10000)/100)
+ !! testing
+ print*, 'getchemfield: memid, jd, jdmid, CL_time = ',memid,jd,jdmid,CL_time(1)
+ else
+ call caldate(jd, jjjjmmdd, hhmmss)
eomday=calceomday(jjjjmmdd/100)
- write(amonth,'(I2.2)') mm
+ CL_time(memid)=CL_time(memid-1)+real(ldirect*eomday,kind=dp)
+ endif
+ ! time of field in seconds from start
+ memCLtime(memid)=int((CL_time(memid)-bdate)*86400._dp)
+
+ write(*,*) 'getchemfield: memid, memCLtime = ',memCLtime(memid)
+ write(*,*) 'getchemfield: memid, CL_time = ',CL_time(memid)
+
+ call caldate(CL_time(memid), jjjjmmdd, hhmmss)
+ mm=int((jjjjmmdd-(jjjjmmdd/10000)*10000)/100)
+ write(amonth,'(I2.2)') mm
+
+ do nr=1,nreagent
+
+ write(*,*) 'Reading two new CL fields for: ',trim(reagents(nr))
+
+ ! Read new chem field
+ !********************
+
write(CL_name,'(A)') trim(reag_path(nr))//trim(reagents(nr))//'_'//amonth//'.nc'
inquire(file=CL_name,exist=lexist)
if (lexist) then
@@ -266,6 +313,10 @@ module chemistry_mod
write(*,*) '#### '//trim(CL_name)//' ####'
error stop
endif
+
+ ! Read average jrates
+ !********************
+
if (reag_hourly(nr).gt.0) then
! Read average jrates into memory
write(jr_name,'(A)') trim(reag_path(nr))//'jrate_average.nc'
@@ -280,10 +331,11 @@ module chemistry_mod
error stop
endif
endif ! reag_hourly
- end do ! memid
- end do ! nreagent
+ end do ! nreagent
+ end do ! memid
+
endif ! update hourly fields
@@ -418,7 +470,7 @@ module chemistry_mod
! current hour of simulation
curhour=itime/3600
- print*, 'getchemhourly: curhour, curCLhour = ',curhour, curCLhour
+ write(*,*) 'getchemhourly: curhour, curCLhour = ',curhour, curCLhour
jul=bdate+real(itime,kind=dp)/86400._dp
call caldate(jul,jjjjmmdd,hhmmss)
@@ -434,7 +486,7 @@ module chemistry_mod
dt2=float(memCLtime(2)-interp_time)
dtt=1./(dt1+dt2)
!! testing
-! print*, 'getchemhourly: dt1, dt2, dtt = ',dt1,dt2,dtt
+ print*, 'getchemhourly: dt1, dt2, dtt = ',dt1,dt2,dtt
do nr=1,nreagent
write(*,*) 'Interpolating fields for reagent: ',reagents(nr)
clfield_cur(:,:,:,nr)=(dt2*CL_field(:,:,:,nr,1) + dt1*CL_field(:,:,:,nr,2))*dtt
@@ -707,7 +759,7 @@ module chemistry_mod
dzz=1./(dz1+dz2)
!! testing
- if(ii.lt.100) print*, 'chemreaction: nr, clx, cly, clz, clzm = ',nr, clx, cly, clz, clzm
+! if(ii.lt.100) print*, 'chemreaction: nr, clx, cly, clz, clzm = ',nr, clx, cly, clz, clzm
! chem reagent for particle time and location
cl_cur=(dz2*clfield_cur(clx,cly,clzm,nr) + &
diff --git a/src/com_mod.f90 b/src/com_mod.f90
index 45eb215616aca1f1c80c5dfde3c7433928cb2900..37a337a996c7a7b63de80453e3c13cff1071fd41 100644
--- a/src/com_mod.f90
+++ b/src/com_mod.f90
@@ -420,6 +420,7 @@ module com_mod
character(len=16), allocatable, dimension(:) :: receptorname
integer :: cpointer(maxrecsample)
integer :: numreceptor, numcurrec
+ logical :: lrecregular
! satellite receptors
real, allocatable, dimension(:) :: xsatellite,ysatellite
@@ -428,7 +429,7 @@ module com_mod
real, allocatable, dimension(:,:) :: zsatellite
real, allocatable, dimension(:,:,:) :: csatellite, csatuncert
real, allocatable, dimension(:,:) :: xksatellite, nnsatellite
- character(len=16), allocatable, dimension(:) :: satellitename
+ character(len=24), allocatable, dimension(:) :: satellitename
integer :: numsatreceptor, nlayermax, numsatellite, numcursat
integer, allocatable, dimension(:) :: nnsatlayer
integer :: csatpointer(maxrecsample)
@@ -438,6 +439,7 @@ module com_mod
! receptorarea area of 1*1 grid cell at receptor point
! numreceptor number of receptors (non-satellite)
! numcurrec number of receptors in current time interval (updated each time interval)
+ ! lrecregular logical if receptor output should be at regular intervals (and not according to RECEPTORS namelist)
! numsatreceptor number of satellite receptors (aka. retrievals)
! numcursat number of satellite receptors in current time interval (updated each time interval)
! numsatellite number of satellite instruments
diff --git a/src/emissions_mod.f90 b/src/emissions_mod.f90
index de85615a2b0266a8e52e43327b1e942ad2624d93..60656c8f5f9ae57fcfc018fc90de75a166204a5d 100644
--- a/src/emissions_mod.f90
+++ b/src/emissions_mod.f90
@@ -340,7 +340,7 @@ module emissions_mod
implicit none
integer :: itime
- real(kind=dp) :: julstart
+ real(kind=dp) :: julstart, jd, jdmid
integer :: jjjjmmdd, hhmmss, dd, mm, yyyy
integer :: nn, ks, eomday, memid
character(len=4) :: ayear
@@ -448,14 +448,29 @@ module emissions_mod
do memid=1,2
if (memid.eq.1) then
- em_time(memid)=bdate+real(ldirect*itime,kind=dp)/86400._dp
+ jd=bdate+real(ldirect*itime,kind=dp)/86400._dp
+ call caldate(jd, jjjjmmdd, hhmmss)
+ ! middle of month day
+ jdmid=juldate(int(jjjjmmdd/100)*100+15,0)
+ if (jd.ge.jdmid) then
+ ! use current month
+ em_time(memid)=jdmid
+ else
+ ! use last month
+ eomday=calceomday(jjjjmmdd/100)
+ em_time(memid)=jdmid-real(eomday,kind=dp)
+ endif
else
+ call caldate(jd, jjjjmmdd, hhmmss)
+ eomday=calceomday(jjjjmmdd/100)
em_time(memid)=em_time(memid-1)+real(ldirect*eomday,kind=dp)
endif
- em_memtime(memid)=int((em_time(memid)-bdate)*86400._dp)*ldirect
+ em_memtime(memid)=int((em_time(memid)-bdate)*86400._dp)
+
+ write(*,*) 'getemissions: memid, em_time = ',memid, em_time(memid)
+ write(*,*) 'getemissions: memid, em_memtime = ',memid, em_memtime(memid)
call caldate(em_time(memid), jjjjmmdd, hhmmss)
- eomday=calceomday(jjjjmmdd/100)
yyyy=jjjjmmdd/10000
mm=(jjjjmmdd-yyyy*10000)/100
dd=jjjjmmdd-yyyy*10000-mm*100
diff --git a/src/initdomain_mod.f90 b/src/initdomain_mod.f90
index 8e994760ff19cf183ac4a7efa5947fe6a02df678..059f29ff052936a025ab5eb9eb9e8d3165f86fde 100644
--- a/src/initdomain_mod.f90
+++ b/src/initdomain_mod.f90
@@ -485,7 +485,7 @@ module initdomain_mod
integer :: j,ix,jy,kz,ncolumn,numparttot,iterminate,stat
real :: ylat,ylatp,ylatm,hzone
- real :: cosfactm,cosfactp,deltacol,dz1,dz2,dz,pnew,fractus
+ real :: cosfactm,cosfactp,deltacol,dz1,dz2,dz,pnew,pnew_temp,fractus
real,parameter :: pih=pi/180.
real :: colmasstotal,zposition
real :: hgt_tmp
@@ -752,21 +752,22 @@ module initdomain_mod
! poles), distribute the particles randomly. When only few particles are
! left distribute them randomly
- if ((ncolumn.gt.20)) then !.and.(ncolumn-j.gt.20)) then
+ if ((ncolumn.gt.20).and.(ncolumn-j.gt.20)) then
+ pnew_temp=pnew-ran1(idummy,0)*deltacol
pnew=pnew-deltacol
-! else if (ncolumn.gt.20) then
-! pnew=pnew-ran1(idummy,0)*(pnew-pp(nz))
+ else if (ncolumn.gt.20) then
+ pnew_temp=pnew-ran1(idummy,0)*(pnew-pp(nz))
else
- pnew=pp(1)-ran1(idummy,0)*(pp(1)-pp(nz))
+ pnew_temp=pp(1)-ran1(idummy,0)*(pp(1)-pp(nz))
endif
- pnew=min(pp(1),pnew)
- pnew=max(pp(nz)+eps,pnew)
+ pnew_temp=min(pp(1),pnew_temp)
+ pnew_temp=max(pp(nz)+eps,pnew_temp)
! find vertical layer
do kz=1,nz-1
- if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew)) then
- dz1=log(pp(kz))-log(pnew)
- dz2=log(pnew)-log(pp(kz+1))
+ if ((pp(kz).ge.pnew_temp).and.(pp(kz+1).lt.pnew_temp)) then
+ dz1=log(pp(kz))-log(pnew_temp)
+ dz2=log(pnew_temp)-log(pp(kz+1))
dz=1./(dz1+dz2)
! Assign particle position
@@ -885,9 +886,9 @@ module initdomain_mod
ixm=min(ixm,nxini(indxn))
jym=min(jym,nyini(indxn))
!! testing
- if (jj.eq.1.and.jy.lt.5.and.ix.lt.5) then
- print*, 'init_domainfill: lonini, xl, latini, yl = ',lonini(ixm,indxn),xl,latini(jym,indxn),yl
- endif !!
+! if (jj.eq.1.and.jy.lt.5.and.ix.lt.5) then
+! print*, 'init_domainfill: lonini, xl, latini, yl = ',lonini(ixm,indxn),xl,latini(jym,indxn),yl
+! endif !!
! Get vertical position in gridini
if (any(altini(:,indxn).gt.0)) then
! vertical coordinate in metres above ground
@@ -907,7 +908,7 @@ module initdomain_mod
dz1*gridini(ixm,jym,indzp,indxn))*dz
else if (any(prsini(:,:,:,indxn).gt.0)) then
! vertical coordinate in pressure (Pa)
- presspart=pnew
+ presspart=pnew_temp
indzm=nzini(indxn)-1
indzp=nzini(indxn)
do ii=2,nzini(indxn)
diff --git a/src/readoptions_mod.f90 b/src/readoptions_mod.f90
index 5e92c357aee095c389f2e874f746cdb919f6acc0..d82bcee0f51f38f261d130331a936ede4cd379e3 100644
--- a/src/readoptions_mod.f90
+++ b/src/readoptions_mod.f90
@@ -1973,6 +1973,8 @@ subroutine readreceptors
lon = -999.
lat = -999.
+ time = -999.
+ lrecregular = .false.
! try namelist input
read(unitreceptor,receptors,iostat=ios)
@@ -1994,7 +1996,8 @@ subroutine readreceptors
lon=-999.9
read(unitreceptor,receptors,iostat=ios)
if ((lon.lt.-900).or.(ios.ne.0)) exit
- if ((time.lt.bdate).or.(time.ge.edate)) cycle ! skip receptors not in simulation window
+ ! skip receptors for which a timestamp is given but are not in simulation window
+ if ((time.ne.-999.).and.((time.lt.bdate).or.(time.ge.edate))) cycle
j=j+1
end do
numreceptor=j
@@ -2017,15 +2020,20 @@ subroutine readreceptors
lon=-999.9
read(unitreceptor,receptors,iostat=ios)
if ((lon.lt.-900).or.(ios.ne.0)) exit ! read error
- if ((time.lt.bdate).or.(time.ge.edate)) cycle ! skip receptors not in simulation window
+ ! skip receptors for which a timestamp is given but are not in simulation window
+ if ((time.ne.-999.).and.((time.lt.bdate).or.(time.ge.edate))) cycle
j=j+1
receptorname(j)=receptor
xreceptor(j)=(lon-xlon0)/dx ! transform to grid coordinates
yreceptor(j)=(lat-ylat0)/dy
zreceptor(j)=alt
- treceptor(j)=int((time-bdate)*24.*3600.) ! time in sec
- ! round to nearest 10 seconds
- treceptor(j)=nint(real(treceptor(j))/10.)*10
+ if (time.ne.-999.) then
+ treceptor(j)=int((time-bdate)*24.*3600.) ! time in sec
+ ! round to nearest 10 seconds
+ treceptor(j)=nint(real(treceptor(j))/10.)*10
+ else
+ treceptor(j)=-999
+ endif
xm=r_earth*cos(lat*pi/180.)*dx/180.*pi
ym=r_earth*dy/180.*pi
receptorarea(j)=xm*ym
@@ -2037,6 +2045,12 @@ subroutine readreceptors
endif
+ ! if not timestamp given in namelist assume regular output
+ ! according to COMMAND file settings
+ if (.not.any(treceptor.ne.-999)) then
+ lrecregular=.true.
+ endif
+
!! testing
! write(*,*) 'readreceptors: '
! do j=1,numreceptor
diff --git a/src/receptor_mod.f90 b/src/receptor_mod.f90
index 0883b60b3e4f754c018e1a63d606033665ee4d77..44fa3a4ec47dd92935c38929af4805580ace82fa 100644
--- a/src/receptor_mod.f90
+++ b/src/receptor_mod.f90
@@ -58,17 +58,6 @@ module receptor_mod
xkreceptor(:)=0.
endif
- if (numsatreceptor.gt.0) then
- allocate( csatellite(nspec,nlayermax,maxrecsample) )
- allocate( csatuncert(nspec,nlayermax,maxrecsample) )
- allocate( nnsatellite(nlayermax,maxrecsample) )
- allocate( xksatellite(nlayermax,maxrecsample) )
- csatellite(:,:,:)=0.
- csatuncert(:,:,:)=0.
- nnsatellite(:,:)=0.
- xksatellite(:,:)=0.
- endif
-
end subroutine alloc_receptor
@@ -106,11 +95,13 @@ module receptor_mod
integer, dimension(:), allocatable :: timerec
real, dimension(:,:,:), allocatable :: crec, cunc
character(len=16), dimension(:), allocatable :: namerec
+ character(len=24), dimension(:), allocatable :: namesatrec
real, dimension(:,:,:), allocatable :: nnrec_omp, xkrec_omp, altrec_omp
real, dimension(:,:), allocatable :: lonrec_omp, latrec_omp
integer, dimension(:,:), allocatable :: timerec_omp
real, dimension(:,:,:,:), allocatable :: crec_omp, cunc_omp
character(len=16), dimension(:,:), allocatable :: namerec_omp
+ character(len=24), dimension(:,:), allocatable :: namesatrec_omp
integer, dimension(:), allocatable :: nr_omp
real, parameter :: weightair=28.97
integer, parameter :: unitoutrecdates=109
@@ -200,6 +191,7 @@ module receptor_mod
latrec(maxrecsample),&
altrec(maxrecsample,nlayermax),&
namerec(maxrecsample),&
+ namesatrec(maxrecsample),&
timerec(maxrecsample),&
crec(nspec,maxrecsample,nlayermax),&
cunc(nspec,maxrecsample,nlayermax))
@@ -209,6 +201,7 @@ module receptor_mod
latrec_omp(maxrecsample,nthreads),&
altrec_omp(maxrecsample,nlayermax,nthreads),&
namerec_omp(maxrecsample,nthreads),&
+ namesatrec_omp(maxrecsample,nthreads),&
timerec_omp(maxrecsample,nthreads),&
crec_omp(nspec,maxrecsample,nlayermax,nthreads),&
cunc_omp(nspec,maxrecsample,nlayermax,nthreads))
@@ -351,7 +344,11 @@ module receptor_mod
lonrec_omp(nr,thread)=xreceptor(n)*dx+xlon0
latrec_omp(nr,thread)=yreceptor(n)*dy+ylat0
altrec_omp(nr,1,thread)=zreceptor(n)
- timerec_omp(nr,thread)=treceptor(n)
+ if ( lrecregular ) then
+ timerec_omp(nr,thread)=itime
+ else
+ timerec_omp(nr,thread)=treceptor(n)
+ endif
namerec_omp(nr,thread)=receptorname(n)
end do ! numcurrec
@@ -536,7 +533,7 @@ module receptor_mod
latrec_omp(nr,thread)=ysatellite(n)*dy+ylat0
altrec_omp(nr,:,thread)=zsatellite(:,n)
timerec_omp(nr,thread)=tsatellite(n)
- namerec_omp(nr,thread)=satellitename(n)
+ namesatrec_omp(nr,thread)=satellitename(n)
end do ! numcursat
!$OMP END DO
@@ -556,7 +553,7 @@ module receptor_mod
latrec(nr:nr+nr_omp(ithread)-1)=latrec_omp(1:nr_omp(ithread),ithread)
altrec(nr:nr+nr_omp(ithread)-1,:)=altrec_omp(1:nr_omp(ithread),:,ithread)
timerec(nr:nr+nr_omp(ithread)-1)=timerec_omp(1:nr_omp(ithread),ithread)
- namerec(nr:nr+nr_omp(ithread)-1)=namerec_omp(1:nr_omp(ithread),ithread)
+ namesatrec(nr:nr+nr_omp(ithread)-1)=namesatrec_omp(1:nr_omp(ithread),ithread)
nr=nr+nr_omp(ithread)
!! testing
! print*, 'receptor_mod: satellites, nr = ',nr
@@ -569,18 +566,18 @@ module receptor_mod
if (nr.gt.0) then
if (lnetcdfout.eq.1) then
#if USE_NCF
- call write_satellite_netcdf(crec,cunc,nnrec,xkrec,lonrec,latrec,altrec,timerec,namerec,nr)
+ call write_satellite_netcdf(crec,cunc,nnrec,xkrec,lonrec,latrec,altrec,timerec,namesatrec,nr)
#endif
else
- call write_satellite_binary(crec,cunc,nnrec,xkrec,lonrec,latrec,altrec,timerec,namerec,nr)
+ call write_satellite_binary(crec,cunc,nnrec,xkrec,lonrec,latrec,altrec,timerec,namesatrec,nr)
endif
endif
! advance output index
spointer=spointer+nr
!! testing
-! print*, 'receptor_mod: satellites, nr_omp = ',nr_omp
-! print*, 'receptor_mod: spointer = ',spointer
+ print*, 'receptor_mod: satellites, nr = ',nr
+ print*, 'receptor_mod: spointer = ',spointer
! close files
if (lnetcdfout.eq.1) then
@@ -597,8 +594,8 @@ module receptor_mod
deallocate(densityoutrecept)
deallocate(crec,cunc,nnrec,xkrec)
deallocate(crec_omp,cunc_omp,nnrec_omp,xkrec_omp)
- deallocate(lonrec,latrec,altrec,timerec,namerec)
- deallocate(lonrec_omp,latrec_omp,altrec_omp,timerec_omp,namerec_omp)
+ deallocate(lonrec,latrec,altrec,timerec,namerec,namesatrec)
+ deallocate(lonrec_omp,latrec_omp,altrec_omp,timerec_omp,namesatrec_omp)
deallocate(nr_omp)
! End of receptor output
@@ -620,7 +617,7 @@ module receptor_mod
crecuncert(:,:)=0.
nnreceptor(:)=0.
xkreceptor(:)=0.
- if (numsatreceptor.gt.0) then
+ if (numsatellite.gt.0) then
csatellite(:,:,:)=0.
csatuncert(:,:,:)=0.
nnsatellite(:,:)=0.
@@ -718,8 +715,8 @@ module receptor_mod
do n=1,numreceptor
- if ((treceptor(n).lt.lrecoutstart).or. &
- (treceptor(n).ge.lrecoutend)) cycle ! skip if not in current sampling time interval
+ if ((.not.lrecregular).and.((treceptor(n).lt.lrecoutstart).or. &
+ (treceptor(n).ge.lrecoutend))) cycle ! skip if not in current sampling time interval
! update pointer
k=k+1
diff --git a/src/receptor_netcdf_mod.f90 b/src/receptor_netcdf_mod.f90
index 1c7de4ad36a5a1a7ff214a659df0283419abc434..7d311f51fa71cc5c2d9e249dad4ba4dfc752be7e 100644
--- a/src/receptor_netcdf_mod.f90
+++ b/src/receptor_netcdf_mod.f90
@@ -32,15 +32,42 @@ module receptor_netcdf_mod
integer :: nnvar_id, xkvar_id
! satellite receptors
integer :: spointer
+ integer :: memsatday
character(len=8), dimension(:), allocatable :: sat_name
+ character(len=256), dimension(:), allocatable :: sat_file, sat_path
integer :: ncsat_id
- integer :: satrecdim_id, sataltdim_id, satlayerdim_id
+ integer :: satrecdim_id, ncharsatdim_id, sataltdim_id, satlayerdim_id
integer :: sattimevar_id, satrecvar_id, satlonvar_id, satlatvar_id, sataltvar_id, satnamevar_id
integer, dimension(:), allocatable :: satvar_id, satuncvar_id
integer :: satnnvar_id, satxkvar_id
contains
+
+ subroutine alloc_satellite
+
+ !*****************************************************************************
+ ! *
+ ! This routine allocates variables for satellite concentrations *
+ ! *
+ !*****************************************************************************
+
+ implicit none
+
+ if (numsatellite.gt.0) then
+ allocate( csatellite(nspec,nlayermax,maxrecsample) )
+ allocate( csatuncert(nspec,nlayermax,maxrecsample) )
+ allocate( nnsatellite(nlayermax,maxrecsample) )
+ allocate( xksatellite(nlayermax,maxrecsample) )
+ csatellite(:,:,:)=0.
+ csatuncert(:,:,:)=0.
+ nnsatellite(:,:)=0.
+ xksatellite(:,:)=0.
+ endif
+
+ end subroutine alloc_satellite
+
+
subroutine receptorout_init
!*****************************************************************************
@@ -105,7 +132,7 @@ module receptor_netcdf_mod
! define dimensions
!******************
- call nf90_err( nf90_def_dim(nc_id, "rec", numreceptor, recdim_id) )
+ call nf90_err( nf90_def_dim(nc_id, "rec", nf90_unlimited, recdim_id) )
call nf90_err( nf90_def_dim(nc_id, 'nchar', 16, nchardim_id) )
! define variables
@@ -221,7 +248,7 @@ module receptor_netcdf_mod
character(10) :: time
character(5) :: zone
- if (numsatreceptor.eq.0) then
+ if (numsatellite.eq.0) then
return
endif
@@ -258,10 +285,10 @@ module receptor_netcdf_mod
! define dimensions
!******************
- call nf90_err( nf90_def_dim(ncsat_id, "rec", numsatreceptor, satrecdim_id) )
+ call nf90_err( nf90_def_dim(ncsat_id, "rec", nf90_unlimited, satrecdim_id) )
call nf90_err( nf90_def_dim(ncsat_id, "layer", (nlayermax-1), satlayerdim_id) )
call nf90_err( nf90_def_dim(ncsat_id, "level", nlayermax, sataltdim_id) )
- call nf90_err( nf90_def_dim(ncsat_id, 'nchar', 16, nchardim_id) )
+ call nf90_err( nf90_def_dim(ncsat_id, 'nchar', 24, ncharsatdim_id) )
! define variables
!*****************
@@ -277,7 +304,7 @@ module receptor_netcdf_mod
call nf90_err( nf90_put_att(ncsat_id, sattimevar_id, 'calendar', 'proleptic_gregorian') )
! receptor names
- call nf90_err( nf90_def_var(ncsat_id, "receptorname", nf90_char, (/ nchardim_id, satrecdim_id /), &
+ call nf90_err( nf90_def_var(ncsat_id, "receptorname", nf90_char, (/ ncharsatdim_id, satrecdim_id /), &
satnamevar_id) )
call nf90_err( nf90_put_att(ncsat_id, satnamevar_id, "longname", "receptor name") )
@@ -297,19 +324,19 @@ module receptor_netcdf_mod
call nf90_err( nf90_put_att(ncsat_id, satlatvar_id, "units", "degrees_north") )
! receptor altitude
- call nf90_err( nf90_def_var(ncsat_id, "alt", nf90_float, (/ satrecdim_id, sataltdim_id /), sataltvar_id) )
+ call nf90_err( nf90_def_var(ncsat_id, "alt", nf90_float, (/ sataltdim_id, satrecdim_id /), sataltvar_id) )
call nf90_err( nf90_put_att(ncsat_id, sataltvar_id, "longname", "receptor altitude of levels") )
call nf90_err( nf90_put_att(ncsat_id, sataltvar_id, "units", "meters") )
! species specific variables
do ks=ks_start,nspec
! mixing ratio output for each layer of retrieval
- call nf90_err( nf90_def_var(ncsat_id, trim(species(ks)), nf90_float, (/ satrecdim_id, satlayerdim_id /), &
+ call nf90_err( nf90_def_var(ncsat_id, trim(species(ks)), nf90_float, (/ satlayerdim_id, satrecdim_id /), &
satvar_id(ks)) )
call nf90_err( nf90_put_att(ncsat_id, satvar_id(ks), "units", "pptv") )
call nf90_err( nf90_put_att(ncsat_id, satvar_id(ks), "longname", "mean VMR") )
! uncertainty output for each layer of retrieval
- call nf90_err( nf90_def_var(ncsat_id, trim(species(ks))//"_uncert", nf90_float, (/ satrecdim_id, satlayerdim_id /), &
+ call nf90_err( nf90_def_var(ncsat_id, trim(species(ks))//"_uncert", nf90_float, (/ satlayerdim_id, satrecdim_id /), &
satuncvar_id(ks)) )
call nf90_err( nf90_put_att(ncsat_id, satuncvar_id(ks), "units", "pptv") )
call nf90_err( nf90_put_att(ncsat_id, satuncvar_id(ks), "longname", "uncertainty VMR") )
@@ -317,11 +344,11 @@ module receptor_netcdf_mod
! not species specific variables
! number of particles in receptor output
- call nf90_err( nf90_def_var(ncsat_id,"npart", nf90_float, (/ satrecdim_id, satlayerdim_id /), satnnvar_id) )
+ call nf90_err( nf90_def_var(ncsat_id,"npart", nf90_float, (/ satlayerdim_id, satrecdim_id /), satnnvar_id) )
call nf90_err( nf90_put_att(ncsat_id, satnnvar_id, "units", "counts") )
call nf90_err( nf90_put_att(ncsat_id, satnnvar_id, "longname","number of particles at receptor") )
! average kernel weight at receptor
- call nf90_err( nf90_def_var(ncsat_id,"kernel", nf90_float, (/ satrecdim_id, satlayerdim_id /), satxkvar_id) )
+ call nf90_err( nf90_def_var(ncsat_id,"kernel", nf90_float, (/ satlayerdim_id, satrecdim_id /), satxkvar_id) )
call nf90_err( nf90_put_att(ncsat_id, satxkvar_id, "units", "") )
call nf90_err( nf90_put_att(ncsat_id, satxkvar_id, "longname", "average kernel weight at receptor") )
@@ -471,7 +498,7 @@ module receptor_netcdf_mod
endif
!! testing
- print*, 'write_receptor_netcdf: rpointer, nrec = ',rpointer, nrec
+! print*, 'write_receptor_netcdf: rpointer, nrec = ',rpointer, nrec
! species specific
do ks=ks_start,nspec
@@ -516,7 +543,7 @@ module receptor_netcdf_mod
real, dimension(maxrecsample,nlayermax) :: nnrec, xkrec, altrec
real, dimension(maxrecsample) :: lonrec, latrec
integer, dimension(maxrecsample) :: timerec
- character(len=16), dimension(maxrecsample) :: namerec
+ character(len=24), dimension(maxrecsample) :: namerec
if (llcmoutput) then
ks_start=2
@@ -525,17 +552,19 @@ module receptor_netcdf_mod
endif
!! testing
- print*, 'write_satellite_netcdf: spointer, nrec = ',spointer, nrec
+! print*, 'write_satellite_netcdf: spointer, nrec = ',spointer, nrec
! species specific output
do ks = ks_start,nspec
- call nf90_err( nf90_put_var(ncsat_id,satvar_id(ks),crec(ks,1:nrec,1:nlayermax-1),(/spointer,1/),(/nrec,nlayermax-1/) ) )
- call nf90_err( nf90_put_var(ncsat_id,satuncvar_id(ks),cunc(ks,1:nrec,1:nlayermax-1),(/spointer,1/),(/nrec,nlayermax-1/) ) )
+ call nf90_err( nf90_put_var(ncsat_id,satvar_id(ks),transpose(crec(ks,1:nrec,1:nlayermax-1)),&
+ (/1,spointer/),(/nlayermax-1,nrec/) ) )
+ call nf90_err( nf90_put_var(ncsat_id,satuncvar_id(ks),transpose(cunc(ks,1:nrec,1:nlayermax-1)),&
+ (/1,spointer/),(/nlayermax-1,nrec/) ) )
end do
! non-species specific
! receptor name
- call nf90_err( nf90_put_var(ncsat_id, satnamevar_id, namerec(1:nrec), (/1,spointer/), (/16,nrec/) ) )
+ call nf90_err( nf90_put_var(ncsat_id, satnamevar_id, namerec(1:nrec), (/1,spointer/), (/24,nrec/) ) )
! receptor time
call nf90_err( nf90_put_var(ncsat_id, sattimevar_id, timerec(1:nrec), (/spointer/), (/nrec/) ) )
! receptor longitude
@@ -543,11 +572,13 @@ module receptor_netcdf_mod
! receptor latitude
call nf90_err( nf90_put_var(ncsat_id, satlatvar_id, latrec(1:nrec), (/spointer/), (/nrec/) ) )
! receptor altitude
- call nf90_err( nf90_put_var(ncsat_id, sataltvar_id, altrec(1:nrec,:), (/spointer,1/), (/nrec,(nlayermax)/) ) )
+ call nf90_err( nf90_put_var(ncsat_id, sataltvar_id, transpose(altrec(1:nrec,:)), (/1,spointer/), (/nlayermax,nrec/) ) )
! average number of particles all timesteps for each receptor
- call nf90_err( nf90_put_var(ncsat_id, satnnvar_id, nnrec(1:nrec,1:nlayermax-1), (/spointer,1/), (/nrec,nlayermax-1/) ) )
+ call nf90_err( nf90_put_var(ncsat_id, satnnvar_id, transpose(nnrec(1:nrec,1:nlayermax-1)),&
+ (/1,spointer/), (/nlayermax-1,nrec/) ) )
! average kernel all timesteps
- call nf90_err( nf90_put_var(ncsat_id, satxkvar_id, xkrec(1:nrec,1:nlayermax-1), (/spointer,1/), (/nrec,nlayermax-1/) ) )
+ call nf90_err( nf90_put_var(ncsat_id, satxkvar_id, transpose(xkrec(1:nrec,1:nlayermax-1)), &
+ (/1,spointer/), (/nlayermax-1,nrec/) ) )
end subroutine write_satellite_netcdf
@@ -571,12 +602,12 @@ module receptor_netcdf_mod
!*****************************************************************************
- subroutine readreceptors_satellite
+ subroutine read_satellite_info
!*****************************************************************************
! *
- ! This routine reads the satellite retrieval information for which *
- ! mixing ratios should be modelled *
+ ! This routine reads the satellite information for which satellite *
+ ! retrievals should be read *
! *
! Author: R. Thompson *
! October 2023 *
@@ -585,31 +616,16 @@ module receptor_netcdf_mod
implicit none
- character(len=256) :: file_name, ppath, pfile, strtmp1, strtmp2
- character(len=256), dimension(:), allocatable :: sat_file, sat_path
- character(len=8) :: anretr, psatname
- character(len=4) :: ayear
- character(len=2) :: amonth, aday
- integer :: jjjjmmdd, hhmmss, yyyy, mm, dd
- integer :: nc_id, dimid, varid
- real :: xm, ym
- real(kind=dp) :: jul, jd
- real, allocatable, dimension(:,:) :: xpts, ypts, zpt1, zpt2
- integer, allocatable, dimension(:) :: npt, sdate, stime
- integer :: readerror, writeerror, stat
- integer :: j, nn, nr, nretr, nlayer, tmp_numsat
+ character(len=256) :: ppath, pfile
+ character(len=8) :: psatname
+ integer :: readerror, writeerror
+ integer :: j
integer,parameter :: unitreceptorout=2
logical :: lexist
- real, allocatable, dimension(:) :: tmp_xsat, tmp_ysat, tmp_satarea
- real, allocatable, dimension(:,:) :: tmp_zsat
- integer, allocatable, dimension(:) :: tmp_tsat
- character(len=16), allocatable, dimension(:) :: tmp_satname
! declare namelist
namelist /satellites/ psatname, ppath, pfile
- numsatreceptor=0
-
! For backward runs, do not allow receptor output
!************************************************
@@ -623,7 +639,7 @@ module receptor_netcdf_mod
open(unitreceptor,file=path(1)(1:length(1))//'SATELLITES',form='formatted',status='old',iostat=readerror)
if (readerror.ne.0) then
- write(*,*) 'FLEXPART WARNING readreceptors_satellite: no satellite file found'
+ write(*,*) 'FLEXPART WARNING read_satellite_info: no satellite file found'
return
endif
@@ -632,7 +648,7 @@ module receptor_netcdf_mod
close(unitreceptor)
if (readerror.ne.0) then
- write(*,*) ' #### FLEXPART ERROR in readreceptors_satellite: #### '
+ write(*,*) ' #### FLEXPART ERROR in read_satellite_info: #### '
write(*,*) ' #### error in namelist input #### '
error stop
endif ! only namelist input
@@ -643,7 +659,7 @@ module receptor_netcdf_mod
open(unitreceptorout,file=path(2)(1:length(2))//'SATELLITES.namelist',&
&access='append',status='replace',iostat=writeerror)
if (writeerror.ne.0) then
- write(*,*) ' #### FLEXPART ERROR readreceptors_satellite: cannot open file #### '
+ write(*,*) ' #### FLEXPART ERROR read_satellite_info: cannot open file #### '
write(*,*) ' #### '//trim(path(2)(1:length(2)))//'SATELLITES.namelist #### '
error stop
endif
@@ -652,6 +668,7 @@ module receptor_netcdf_mod
! Get number of satellites
!*************************
+ numsatellite=0
open (unitreceptor,file=trim(path(1))//'SATELLITES',status='old',iostat=readerror)
j=0
do while (readerror.eq.0)
@@ -678,7 +695,7 @@ module receptor_netcdf_mod
read(unitreceptor,satellites,iostat=readerror)
if (readerror.ne.0) exit
j=j+1
- write(*,*) 'readreceptors_satellite: psatname, ppath = ',trim(psatname),', ',trim(ppath)
+ write(*,*) 'read_satellite_info: psatname, ppath = ',trim(psatname),', ',trim(ppath)
sat_name(j)=trim(psatname)
sat_path(j)=trim(ppath)
sat_file(j)=trim(pfile)
@@ -690,39 +707,133 @@ module receptor_netcdf_mod
close(unitreceptor)
close(unitreceptorout)
+ end subroutine read_satellite_info
+
+
+ subroutine readreceptors_satellite(itime)
+
+ !*****************************************************************************
+ ! *
+ ! This routine reads the satellite retrieval information for which *
+ ! mixing ratios should be modelled *
+ ! *
+ ! Author: R. Thompson *
+ ! October 2023 *
+ ! *
+ !*****************************************************************************
+
+ use binary_output_mod, only: satelliteout_init_binary
+
+ implicit none
+
+ integer :: itime
+ character(len=256) :: file_name
+ character(len=6) :: anretr
+ character(len=8) :: adate
+ integer :: jjjjmmdd, hhmmss, yyyy, mm, dd
+ integer :: nc_id, dimid, varid
+ real :: xm, ym
+ real(kind=dp) :: jul, jd, curday
+ real, allocatable, dimension(:,:) :: xpts, ypts, zpt1, zpt2
+ integer, allocatable, dimension(:) :: npt, sdate, stime
+ integer :: stat
+ integer :: j, nn, nr, nretr, nlayer, tmp_numsat
+ logical :: lexist
+ real, allocatable, dimension(:) :: tmp_xsat, tmp_ysat, tmp_satarea
+ real, allocatable, dimension(:,:) :: tmp_zsat
+ integer, allocatable, dimension(:) :: tmp_tsat
+ character(len=24), allocatable, dimension(:) :: tmp_satname
+
+ ! For backward runs, do not allow receptor output
+ !************************************************
+
+ if (ldirect.lt.0) then
+ return
+ endif
+
+ ! If no satellites do nothing
+ !****************************
+
+ if (numsatellite.eq.0) then
+ return
+ endif
+
+ ! Check if retrievals already in memory for current day
+ !******************************************************
+
+ ! If itime is a full day do not update retrievals
+ print*, 'readreceptors_satellite: mod = ',mod(real(itime),86400.)
+ if ((itime.ne.0).and.(mod(real(itime),86400.).eq.0)) then
+ return
+ endif
+
+ curday=bdate+real(itime,kind=dp)/86400._dp
+ curday=floor(curday)
+ print*, 'readreceptors_satellite: curday = ',curday
+ print*, 'readreceptors_satellite: memsatday = ',memsatday
+
+ if (memsatday.eq.curday ) then
+ ! The retrievals for the current day are in memory -> don't do anything
+ return
+ endif
+ memsatday=curday
+
! Get number of satellite receptors
!**********************************
- ! Loop over satellites and days
- jd=bdate
+ ! Loop over satellites
+ jd=curday
tmp_numsat=0
- do while (jd.le.edate)
- do j=1,numsatellite
- ! get filename for current day
- call caldate(jd, jjjjmmdd, hhmmss)
- file_name=sat_file(j)
- call getfilename(jjjjmmdd, file_name)
- write(*,*) 'readreceptors_satellite: file_name = ',file_name
- inquire(file=trim(sat_path(j))//'/'//trim(file_name),exist=lexist)
- if (.not.lexist) then
- write(*,*) 'readreceptors_satellite: no retrievals file for '//ayear//amonth//aday
- cycle
- endif
- ! open file
- call nf90_err( nf90_open(trim(sat_path(j))//'/'//trim(file_name), nf90_nowrite, nc_id) )
- ! read dimensions
- call nf90_err( nf90_inq_dimid(nc_id, 'retrieval', dimid) )
- call nf90_err( nf90_inquire_dimension(nc_id, dimid, len=nretr) )
- call nf90_err( nf90_inq_dimid(nc_id, 'nlayer', dimid) )
- call nf90_err( nf90_inquire_dimension(nc_id, dimid, len=nlayer) )
- tmp_numsat=tmp_numsat+nretr
+ do j=1,numsatellite
+ ! get filename for current day
+ call caldate(jd, jjjjmmdd, hhmmss)
+ write(adate,'(I8.8)') jjjjmmdd
+ file_name=sat_file(j)
+ call getfilename(jjjjmmdd, file_name)
+ write(*,*) 'readreceptors_satellite: file_name = ',file_name
+ inquire(file=trim(sat_path(j))//'/'//trim(file_name),exist=lexist)
+ if (.not.lexist) then
+ write(*,*) 'readreceptors_satellite: no retrievals file for '//adate
+ cycle
+ endif
+ ! open file
+ call nf90_err( nf90_open(trim(sat_path(j))//'/'//trim(file_name), nf90_nowrite, nc_id) )
+ ! read dimensions
+ call nf90_err( nf90_inq_dimid(nc_id, 'retrieval', dimid) )
+ call nf90_err( nf90_inquire_dimension(nc_id, dimid, len=nretr) )
+ call nf90_err( nf90_inq_dimid(nc_id, 'nlayer', dimid) )
+ call nf90_err( nf90_inquire_dimension(nc_id, dimid, len=nlayer) )
+ tmp_numsat=tmp_numsat+nretr
+ ! Define nlayermax only for first day (must not change by day)
+ if (jd.eq.bdate) then
nlayermax=max(nlayermax,nlayer)
nnsatlayer(j)=nlayer
- call nf90_err(nf90_close(nc_id))
- end do
- jd=jd+1d0
+ endif
+ call nf90_err(nf90_close(nc_id))
end do
- nlayermax=nlayermax+1 ! for levels
+ if (jd.eq.bdate) then
+ nlayermax=nlayermax+1 ! for levels
+ endif
+ print*, 'readreceptors_satellite: nlayermax = ',nlayermax
+
+ ! Initialize satellite output
+ !****************************
+
+ ! Only do once
+ print*, 'readreceptors_satellite: jd, bdate = ',jd, bdate
+ if (jd.eq.bdate) then
+ print*, 'readreceptors_satellite: allocating satellite variables'
+ call alloc_satellite
+ if (lnetcdfout.eq.1) then
+#ifdef USE_NCF
+ print*, 'readreceptors_satellite: initialising output file'
+ call satelliteout_init
+#endif
+ else
+ print*, 'readreceptors_satellite: initialising binary output file'
+ call satelliteout_init_binary
+ endif
+ endif
! Allocate temporary arrays
!**************************
@@ -731,94 +842,101 @@ module receptor_netcdf_mod
tmp_tsat(tmp_numsat),tmp_satarea(tmp_numsat),&
tmp_zsat(nlayermax,tmp_numsat),tmp_satname(tmp_numsat))
+ ! Deallocate existing satellite arrays
+ !*************************************
+
+ if (allocated(xsatellite)) then
+ deallocate(xsatellite,ysatellite,&
+ tsatellite,satellitearea,&
+ zsatellite,satellitename)
+ endif
+
! Read satellite retrievals info
!*******************************
- ! Loop over days in simulation
numsatreceptor=0
- jd=bdate
- do while (jd.le.edate)
-
- do j=1,numsatellite
-
- ! get filename for current day
- ! assumes same format as output from prep_satellite
- call caldate(jd, jjjjmmdd, hhmmss)
- file_name=sat_file(j)
- call getfilename(jjjjmmdd, file_name)
- write(*,*) 'readreceptors_satellite: file_name = ',file_name
- inquire(file=trim(sat_path(j))//'/'//trim(file_name),exist=lexist)
- if (.not.lexist) then
- write(*,*) 'readreceptors_satellite: no retrievals file for '//ayear//amonth//aday
- cycle
- endif
+ jd=curday
+ call caldate(jd, jjjjmmdd, hhmmss)
+ write(adate,'(I8.8)') jjjjmmdd
+ print*, 'readreceptors_satellite: adate = ',adate
+
+ do j=1,numsatellite
+
+ ! get filename for current day
+ ! assumes same format as output from prep_satellite
+ call caldate(jd, jjjjmmdd, hhmmss)
+ file_name=sat_file(j)
+ call getfilename(jjjjmmdd, file_name)
+ write(*,*) 'readreceptors_satellite: file_name = ',file_name
+ inquire(file=trim(sat_path(j))//'/'//trim(file_name),exist=lexist)
+ if (.not.lexist) then
+ write(*,*) 'readreceptors_satellite: no retrievals file for '//adate
+ cycle
+ endif
- ! open file
- call nf90_err( nf90_open(trim(sat_path(j))//'/'//trim(file_name), nf90_nowrite, nc_id) )
-
- ! read dimensions
- call nf90_err( nf90_inq_dimid(nc_id, 'retrieval', dimid) )
- call nf90_err( nf90_inquire_dimension(nc_id, dimid, len=nretr) )
- call nf90_err( nf90_inq_dimid(nc_id, 'nlayer', dimid) )
- call nf90_err( nf90_inquire_dimension(nc_id, dimid, len=nlayer) )
-
- ! allocate temporary variables
- allocate(sdate(nretr),stat=stat)
- if (stat.ne.0) write(*,*)'ERROR: could not allocate sdate'
- allocate(stime(nretr),stat=stat)
- if (stat.ne.0) write(*,*)'ERROR: could not allocate stime'
- allocate(xpts(nretr,4),stat=stat)
- if (stat.ne.0) write(*,*)'ERROR: could not allocate xpts'
- allocate(ypts(nretr,4),stat=stat)
- if (stat.ne.0) write(*,*)'ERROR: could not allocate ypts'
- allocate(zpt1(nretr,nlayer),stat=stat)
- if (stat.ne.0) write(*,*)'ERROR: could not allocate zpt1'
- allocate(zpt2(nretr,nlayer),stat=stat)
- if (stat.ne.0) write(*,*)'ERROR: could not allocate zpt2'
-
- ! read coordinate variables
- call nf90_err( nf90_inq_varid(nc_id,'idate',varid) )
- call nf90_err( nf90_get_var(nc_id,varid,sdate) )
- call nf90_err( nf90_inq_varid(nc_id,'itime',varid) )
- call nf90_err( nf90_get_var(nc_id,varid,stime) )
- call nf90_err( nf90_inq_varid(nc_id,'xpoints',varid) )
- call nf90_err( nf90_get_var(nc_id,varid,xpts) )
- call nf90_err( nf90_inq_varid(nc_id,'ypoints',varid) )
- call nf90_err( nf90_get_var(nc_id,varid,ypts) )
- call nf90_err( nf90_inq_varid(nc_id,'zpoint1',varid) )
- call nf90_err( nf90_get_var(nc_id,varid,zpt1) )
- call nf90_err( nf90_inq_varid(nc_id,'zpoint2',varid) )
- call nf90_err( nf90_get_var(nc_id,varid,zpt2) )
- call nf90_err( nf90_close(nc_id) )
-
- ! write to coordinates receptor variables
- do nr=1,nretr
- ! time of retrieval from bdate in seconds
- jul=juldate(sdate(nr),stime(nr))
- if ((jul.lt.bdate).or.(jul.gt.edate)) cycle
- numsatreceptor=numsatreceptor+1
- write(anretr,'(I8.8)') nr
- tmp_satname(numsatreceptor)=trim(sat_name(j))//'_'//anretr
- tmp_tsat(numsatreceptor)=int((jul-bdate)*24.*3600.) ! time in sec
- ! transform to grid coordinates
- tmp_xsat(numsatreceptor)=(sum(xpts(nr,:))/4.-xlon0)/dx
- tmp_ysat(numsatreceptor)=(sum(ypts(nr,:))/4.-ylat0)/dy
- ! vertical coordinates layer boundaries in Pa
- tmp_zsat(1:nlayer,numsatreceptor)=zpt1(nr,:)
- tmp_zsat(nlayer+1,numsatreceptor)=zpt2(nr,nlayer)
- ! area for mixing ratio calc (used if ind_samp = -1)
- xm=r_earth*cos((sum(ypts(nr,:))/4.)*pi/180.)*dx/180.*pi
- ym=r_earth*dy/180.*pi
- tmp_satarea(numsatreceptor)=xm*ym
- end do ! nretr
-
- deallocate(sdate, stime, xpts, ypts, zpt1, zpt2)
-
- end do ! numsatellite
-
- jd=jd+1d0
-
- end do ! jd
+ ! open file
+ call nf90_err( nf90_open(trim(sat_path(j))//'/'//trim(file_name), nf90_nowrite, nc_id) )
+
+ ! read dimensions
+ call nf90_err( nf90_inq_dimid(nc_id, 'retrieval', dimid) )
+ call nf90_err( nf90_inquire_dimension(nc_id, dimid, len=nretr) )
+ call nf90_err( nf90_inq_dimid(nc_id, 'nlayer', dimid) )
+ call nf90_err( nf90_inquire_dimension(nc_id, dimid, len=nlayer) )
+ print*, 'readreceptors_satellite: nretr = ',nretr
+
+ ! allocate temporary variables
+ allocate(sdate(nretr),stat=stat)
+ if (stat.ne.0) write(*,*)'ERROR: could not allocate sdate'
+ allocate(stime(nretr),stat=stat)
+ if (stat.ne.0) write(*,*)'ERROR: could not allocate stime'
+ allocate(xpts(nretr,4),stat=stat)
+ if (stat.ne.0) write(*,*)'ERROR: could not allocate xpts'
+ allocate(ypts(nretr,4),stat=stat)
+ if (stat.ne.0) write(*,*)'ERROR: could not allocate ypts'
+ allocate(zpt1(nretr,nlayer),stat=stat)
+ if (stat.ne.0) write(*,*)'ERROR: could not allocate zpt1'
+ allocate(zpt2(nretr,nlayer),stat=stat)
+ if (stat.ne.0) write(*,*)'ERROR: could not allocate zpt2'
+
+ ! read coordinate variables
+ call nf90_err( nf90_inq_varid(nc_id,'idate',varid) )
+ call nf90_err( nf90_get_var(nc_id,varid,sdate) )
+ call nf90_err( nf90_inq_varid(nc_id,'itime',varid) )
+ call nf90_err( nf90_get_var(nc_id,varid,stime) )
+ call nf90_err( nf90_inq_varid(nc_id,'xpoints',varid) )
+ call nf90_err( nf90_get_var(nc_id,varid,xpts) )
+ call nf90_err( nf90_inq_varid(nc_id,'ypoints',varid) )
+ call nf90_err( nf90_get_var(nc_id,varid,ypts) )
+ call nf90_err( nf90_inq_varid(nc_id,'zpoint1',varid) )
+ call nf90_err( nf90_get_var(nc_id,varid,zpt1) )
+ call nf90_err( nf90_inq_varid(nc_id,'zpoint2',varid) )
+ call nf90_err( nf90_get_var(nc_id,varid,zpt2) )
+ call nf90_err( nf90_close(nc_id) )
+
+ ! write to coordinates receptor variables
+ do nr=1,nretr
+ jul=juldate(sdate(nr),stime(nr))
+ ! skip retrievals not for current day
+ if (floor(jul).ne.curday) cycle
+ numsatreceptor=numsatreceptor+1
+ write(anretr,'(I6.6)') nr
+ tmp_satname(numsatreceptor)=trim(sat_name(j))//'_'//adate//'_'//anretr
+ tmp_tsat(numsatreceptor)=int((jul-bdate)*24.*3600.) ! time in sec
+ ! transform to grid coordinates
+ tmp_xsat(numsatreceptor)=(sum(xpts(nr,:))/4.-xlon0)/dx
+ tmp_ysat(numsatreceptor)=(sum(ypts(nr,:))/4.-ylat0)/dy
+ ! vertical coordinates layer boundaries in Pa
+ tmp_zsat(1:nlayer,numsatreceptor)=zpt1(nr,:)
+ tmp_zsat(nlayer+1,numsatreceptor)=zpt2(nr,nlayer)
+ ! area for mixing ratio calc (used if ind_samp = -1)
+ xm=r_earth*cos((sum(ypts(nr,:))/4.)*pi/180.)*dx/180.*pi
+ ym=r_earth*dy/180.*pi
+ tmp_satarea(numsatreceptor)=xm*ym
+ end do ! nretr
+
+ deallocate(sdate, stime, xpts, ypts, zpt1, zpt2)
+
+ end do ! numsatellite
write(*,*) 'readreceptors_satellite: numsatreceptor = ',numsatreceptor
@@ -837,7 +955,13 @@ module receptor_netcdf_mod
satellitename=tmp_satname(1:numsatreceptor)
deallocate(tmp_xsat,tmp_ysat,tmp_tsat,tmp_satarea,tmp_zsat,tmp_satname)
- deallocate(sat_path,sat_file)
+
+ ! Transform vertical coordinates of satellite receptors
+ !******************************************************
+
+ if (numsatreceptor.gt.0) then
+ call verttransform_satellite
+ endif
end subroutine readreceptors_satellite
diff --git a/src/timemanager_mod.f90 b/src/timemanager_mod.f90
index 931477c1a4ea865d8a873fda1ef3350303e7f288..e7666b688d1c379777df04e925a574a52e7af7ea 100644
--- a/src/timemanager_mod.f90
+++ b/src/timemanager_mod.f90
@@ -108,7 +108,7 @@ subroutine timemanager
#ifdef USE_NCF
use chemistry_mod
use initdomain_mod
- use receptor_netcdf_mod, only: verttransform_satellite
+ use receptor_netcdf_mod, only: readreceptors_satellite, verttransform_satellite
use emissions_mod
use totals_mod
#endif
@@ -149,7 +149,7 @@ subroutine timemanager
drytmp ! dry deposition related
logical :: itsopen
real, dimension(maxrecsample) :: recoutnum ! number of samples for receptor calculation
- real, dimension(nlayermax,maxrecsample) :: recoutnumsat ! number of samples for satellite receptor calculation
+ real, allocatable, dimension(:,:) :: recoutnumsat ! number of samples for satellite receptor calculation
! First output for time 0
!************************
@@ -162,7 +162,7 @@ subroutine timemanager
lrecoutnext=lrecoutstep/2
outnum=0.
recoutnum(:)=0.
- recoutnumsat(:,:)=0.
+! recoutnumsat(:,:)=0.
endif
loutstart=loutnext-loutaver/2
loutend=loutnext+loutaver/2
@@ -283,13 +283,16 @@ subroutine timemanager
endif
#endif
- ! Transform vertical coordinates of satellite receptors
- !******************************************************
+ ! Read satellite receptors
+ !*************************
+
#ifdef USE_NCF
- if ((numsatreceptor.gt.0).and.(itime.eq.0)) then
- call verttransform_satellite
- endif
+ call readreceptors_satellite(itime)
#endif
+ if (.not.allocated(recoutnumsat)) then
+ allocate(recoutnumsat(nlayermax,maxrecsample))
+ recoutnumsat(:,:)=0.
+ endif
! Release particles
!******************
diff --git a/src/totals_mod.f90 b/src/totals_mod.f90
index 65491241618ec983564247c16d62c09f6017b0ab..bf6af51f64ebcbc7969c6119a3b3b808ef5ac164 100644
--- a/src/totals_mod.f90
+++ b/src/totals_mod.f90
@@ -96,7 +96,7 @@ module totals_mod
call nf90_err( nf90_def_dim(nc_id, "species", nspec, specdim_id) )
call nf90_err( nf90_def_dim(nc_id, "reagents", nreagent, reagdim_id) )
call nf90_err( nf90_def_dim(nc_id, "time", nf90_unlimited, timedim_id) )
- call nf90_err( nf90_def_dim(nc_id, "nchar", 18, nchardim_id) )
+ call nf90_err( nf90_def_dim(nc_id, "nchar", 10, nchardim_id) )
! define variables
!*****************
@@ -145,7 +145,8 @@ module totals_mod
! write species info
do ks=1,nspec
- call nf90_err( nf90_put_var(nc_id, spec_id, species(ks), (/1,ks/), (/18,1/)) )
+ print*, 'totals: species = ',trim(species(ks))
+ call nf90_err( nf90_put_var(nc_id, spec_id, species(ks), (/1,ks/), (/10,1/)) )
end do
! close file
@@ -171,7 +172,6 @@ module totals_mod
integer :: var_id
! open file
- print*, 'fn_totals = ',fn_totals
call nf90_err( nf90_open(trim(fn_totals), nf90_write, nc_id) )
! get length of time dimension -> increase index by one to write new data