From a117bfbbc2dd6d9deb961ebe92224db9199dc4fb Mon Sep 17 00:00:00 2001 From: ronesy <rlt@nilu.no> Date: Mon, 27 May 2024 15:40:00 +0200 Subject: [PATCH] Removed need to specify receptor for each time stamp when all receptors have fixed location, means RECEPTORS file can be given with no timestamp in which case receptor output is given for each lrecoutstep time. Changed reading of satellite receptors to read in only those for current day. Fixed bug in interpolation of reagent fields (CL_fields) between months. Former-commit-id: 1e07841f9d8c44180b8ca8b999d4dd99c1c918dd --- src/FLEXPART.f90 | 8 +- src/binary_output_mod.f90 | 4 +- src/chemistry_mod.f90 | 118 +++++++--- src/com_mod.f90 | 4 +- src/emissions_mod.f90 | 23 +- src/gitversion.txt | 1 - src/initdomain_mod.f90 | 29 +-- src/readoptions_mod.f90 | 24 +- src/receptor_mod.f90 | 43 ++-- src/receptor_netcdf_mod.f90 | 436 +++++++++++++++++++++++------------- src/timemanager_mod.f90 | 19 +- src/totals_mod.f90 | 6 +- 12 files changed, 460 insertions(+), 255 deletions(-) delete mode 100644 src/gitversion.txt diff --git a/src/FLEXPART.f90 b/src/FLEXPART.f90 index 7a20ec48..ec735eb0 100644 --- a/src/FLEXPART.f90 +++ b/src/FLEXPART.f90 @@ -44,7 +44,7 @@ program flexpart real :: s_timemanager character(len=256) :: & inline_options ! pathfile, flexversion, arg2 - character(len=256) :: gitversion_tmp="8d9c680 Mon Apr 22 15:37:39 2024 +0200" + character(len=256) :: gitversion_tmp="undefined" ! Keeping track of the total running time of FLEXPART, printed out at the end. !***************************************************************************** @@ -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 53df70e9..7f4118c6 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 3bd7f844..45e08927 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 ed717eff..3ea51658 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 de85615a..60656c8f 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/gitversion.txt b/src/gitversion.txt deleted file mode 100644 index bc95db4e..00000000 --- a/src/gitversion.txt +++ /dev/null @@ -1 +0,0 @@ -8d9c680 Mon Apr 22 15:37:39 2024 +0200 diff --git a/src/initdomain_mod.f90 b/src/initdomain_mod.f90 index 8e994760..059f29ff 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 f71bb33a..424247a0 100644 --- a/src/readoptions_mod.f90 +++ b/src/readoptions_mod.f90 @@ -1948,6 +1948,8 @@ subroutine readreceptors lon = -999. lat = -999. + time = -999. + lrecregular = .false. ! try namelist input read(unitreceptor,receptors,iostat=ios) @@ -1969,7 +1971,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 @@ -1992,15 +1995,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 @@ -2012,6 +2020,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 0883b60b..44fa3a4e 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 1c7de4ad..7d311f51 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 931477c1..e7666b68 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 65491241..bf6af51f 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 -- GitLab