From 3c1c46fd1a4e2e5c000f23a9c203312d17546a03 Mon Sep 17 00:00:00 2001
From: ronesy <rlt@nilu.no>
Date: Fri, 20 Dec 2024 12:05:17 +0100
Subject: [PATCH] Fix to average.f90 to include new retrieval OCO-2 and
S5P-BLND and inclusion of routines and SETTINGS for S5P-BLND
---
prep_satellite/README_satellites.txt | 50 +++
prep_satellite/SETTINGS_bremen | 2 +
prep_satellite/SETTINGS_tropomi | 2 +
prep_satellite/SETTINGS_tropomi_blnd | 159 ++++++++
prep_satellite/average.f90 | 63 ++--
prep_satellite/get_bremen.f90 | 62 ++--
prep_satellite/get_oco2.f90 | 37 +-
prep_satellite/get_tropomi.f90 | 109 ++++--
prep_satellite/get_tropomi_blnd.f90 | 518 +++++++++++++++++++++++++++
prep_satellite/main.f90 | 3 +
prep_satellite/makefile | 1 +
prep_satellite/mod_save.f90 | 57 +--
prep_satellite/mod_settings.f90 | 14 +
13 files changed, 942 insertions(+), 135 deletions(-)
create mode 100644 prep_satellite/README_satellites.txt
create mode 100644 prep_satellite/SETTINGS_tropomi_blnd
create mode 100644 prep_satellite/get_tropomi_blnd.f90
diff --git a/prep_satellite/README_satellites.txt b/prep_satellite/README_satellites.txt
new file mode 100644
index 0000000..5441efc
--- /dev/null
+++ b/prep_satellite/README_satellites.txt
@@ -0,0 +1,50 @@
+================================================================
+
+ FLEXPART/FLEXINVERT Satellites Pre-processor
+
+================================================================
+
+Description:
+ prep_satellite is a pre-processor for satellite retrievals
+ to be run before running FLEXPART/FLEXINVERT. It contains
+ SETTINGS files and get_[satellite_id].f90 routines to read satellite
+ retrieval data from the original format and extract the
+ required variables needed for FLEXPART/FLEXINVERT, perform
+ any averaging to super-observations, and write the data into
+ two output files for each day:
+ i) releases_[satellite_id]_YYYYMMDD.nc
+ ii) retrieval_[satellite_id]_YYYYMMDD.nc
+ The first file is needed for FLEXPART to calculate total
+ column SRRs, the second file contains the column average
+ mixing ratios and is needed for FLEXINVERT
+
+Currently available retrievals:
+ 1) XCH4 TROPOMI Official product provided by SRON
+ routine: get_tropomi.f90
+ satellite id in SETTINGS: tropomi
+ 2) XCH4 TROPOMI WFMD retrieval product provided by University
+ of Bremen
+ routine: get_bremen.f90
+ satellite id in SETTINGS: bremen
+ 3) XCH4 TROPOMI-GOSAT blended product provided by Harvard
+ University
+ routine: get_tropomi_blnd.f90
+ satellite id in SETTINGS: s5p_blnd
+ 4) XCO2 OCO-2 Official "Lite" product provided by NASA
+ routine: get_oco2.f90
+ satellite id in SETTINGS: oco2
+
+To add new retrievals:
+ i) create a new get_[satellite_id].f90 file
+ ii) edit main.f90 and average.f90 to include this
+ new retrieval type
+ iii) add get_[satellite_id].f90 to the makefile
+
+Usage:
+ i) Edit the SETTINGS file
+ ii) Edit sbatch_satellite.sh
+ and run this to prepare the output files
+ iii) Edit run_flexpart.sh
+ and run this to run FLEXPART
+
+
diff --git a/prep_satellite/SETTINGS_bremen b/prep_satellite/SETTINGS_bremen
index 4788ef4..fe3121b 100644
--- a/prep_satellite/SETTINGS_bremen
+++ b/prep_satellite/SETTINGS_bremen
@@ -69,6 +69,8 @@ pdel_name: pressure_levels
dryair_name: pressure_weight
# land fraction
landfrac_name: land_fraction
+# albedo SWIR
+albedo_name: apparent_albedo
# COMMAND
diff --git a/prep_satellite/SETTINGS_tropomi b/prep_satellite/SETTINGS_tropomi
index d1142ff..6cb97e6 100644
--- a/prep_satellite/SETTINGS_tropomi
+++ b/prep_satellite/SETTINGS_tropomi
@@ -74,6 +74,8 @@ psurf_name: PRODUCT/SUPPORT_DATA/INPUT_DATA/surface_pressure
pdel_name: PRODUCT/SUPPORT_DATA/INPUT_DATA/pressure_interval
# surface classification
surfclass_name: PRODUCT/SUPPORT_DATA/INPUT_DATA/surface_classification
+# albedo SWIR
+albedo_name: PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/surface_albedo_SWIR
# COMMAND
diff --git a/prep_satellite/SETTINGS_tropomi_blnd b/prep_satellite/SETTINGS_tropomi_blnd
new file mode 100644
index 0000000..57171f6
--- /dev/null
+++ b/prep_satellite/SETTINGS_tropomi_blnd
@@ -0,0 +1,159 @@
+# ==================================================
+#
+# SATELLITE DATA SETTINGS
+#
+# comment lines start with '#'
+# each parameter line starts with 'parameter name:'
+#
+# ==================================================
+
+# PATHS AND FILES
+
+# use scratch (0: no, 1: yes)
+lscratch: 1
+# satellite instrument/retrieval (one of: tropomi, bremen, s5p_blnd, oco2)
+satellite: s5p_blnd
+# name to include in output files
+satname: s5pblnd_ch4
+# flexpart source path
+path_flexpart: /mypath/FLEXPART_v104/
+# path where to write options folders (folders are created daily)
+path_options: /mypath/flexpart_input/
+# path where to write flexpart output
+path_output: /mypath/flexpart_output/
+# path to OH fields (if use OH chemistry)
+path_ohfield: /mypath/OH_fields/
+# path to satellite data
+path_obs: /mypath/satellite_data/
+# suffix satellite files
+suffix: ESACCI
+# path where to write observation output
+path_obsout: /mypath/satellite_obs_out/
+# AVAILABLE files
+file_avail: /mypath/AVAILABLE
+file_availnest: /mypath/AVAILABLE
+# windfields to use (determines lon of domain for flexpart run) (oper, era5)
+windfield: era5
+
+# SATELLITE VARIABLES
+
+# select only retrievals only land (no sunglint)
+landonly: .true.
+# number measurements dimension
+nmeas_name: nobs
+# number of lat/lon coordinates for each pixel
+ncoord_name: corner
+# number of retrieval levels
+nlevel_name: layer
+# time variable
+time_name: time_utc
+# quality flag
+qa_name: qa_value
+# quality cutoff value (0 = bad, 1 = good, values < qa_cutoff discarded)
+qa_cutoff: 0.99
+# mixing ratio (ppbv)
+vmr_name: methane_mixing_ratio_blended
+# mixing ratio precision (ppbv)
+verr_name: methane_mixing_ratio_precision
+# latitude bounds
+lat_name: latitude_bounds
+# longitude bounds
+lon_name: longitude_bounds
+# column averaging kernel
+ak_name: column_averaging_kernel
+# mixing ratio apriori (mol/m2)
+apri_name: methane_profile_apriori
+# pressure intervals (Pa)
+pdel_name: pressure_interval
+# surface pressure (Pa)
+psurf_name: surface_pressure
+# dry air subcolumns (mol/m2)
+dryair_name: dry_air_subcolumns
+# surface classification
+surfclass_name: surface_classification
+# albedo SWIR
+albedo_name: surface_albedo_SWIR
+
+# COMMAND
+
+# start date of simulation (yyyymmdd)
+datei: 20200826
+# end date of simulation (yyyymmdd)
+datef: 20200831
+# output rate (secs)
+outrate: 86400
+# time average of output (secs)
+outaverage: 86400
+# sampling rate of output (secs)
+outsample: 1200
+# (mass unit = 1, mixing ratio unit = 2)
+ind_source: 1
+# (mass unit = 1, mixing ratio unit = 2)
+ind_receptor: 2
+# use nested output (0 = no, 1 = yes)
+lnested: 1
+# output sensitivity to initial conditions (for backward runs, 0 = none, 1 = mass unit, 2 = mass mixing ratio)
+linit_cond: 2
+
+# OUTGRID
+
+# NOTE: domain specified will determine where retrievals are selected
+
+# longitude of lower left corner of output grid
+# for global use -180 (adjusts automatically for windfield)
+outlonleft: -180.00
+# latitude of lower left corner of output grid
+outlatlower: -90.00
+# number of longitudinal grid cells in output grid
+numxgrid: 180
+# number of latitudinal grid cells in output grid
+numygrid: 90
+# longitudinal resolution in output grid
+dxout: 2.0
+# latitudinal resolution in outout grid
+dyout: 2.0
+# longitude of lower left corner of nested grid (used if nested_output = 1)
+outlonnest: -15.0
+# latitude of lower left corner of nested grid (used if nested_output = 1)
+outlatnest: 33.0
+# number of longitudinal grid cells in nested grid (used if nested_output = 1)
+numxnest: 100
+# number of latitudinal grid cells in nested grid (used if nested_output = 1)
+numynest: 80
+# longitudinal resolution in nested grid (used if nested_output = 1)
+dxoutnest: 0.5
+# latitudinal resolution in nested grid (used if nested_output = 1)
+dyoutnest: 0.5
+# comma-separated list of vertical levels in output grid (meters)
+zlevel: 250, 750, 1500, 2500, 3500, 4500, 5500, 6500, 7500, 8500, 9500, 11000, 13000, 15000, 17000, 19000, 21250, 23750, 26250, 28750, 32500, 37500, 45000, 60000
+
+# RELEASES
+
+# species (see: /options/SPECIES/spec_overview)
+species: CH4
+# ageclass (sec)
+ageclass: 864000
+# number of particles (no longer used in flexpart!)
+npart: 2000
+# mass of particles (no longer used in flexpart!)
+mass: 1000
+
+# AVERAGING
+
+# average retrievals etc. (true/false)
+avg_pixels: .true.
+# If averaging need to set following (else ignored):
+# use predefined grid from file (true/false)
+usegridfile: .false.
+# if usegridfile is true specify file:
+filegrid:
+# min grid resolution (degrees)
+dmin: 0.25
+# number of resolution steps
+nsteps: 2
+# cutoff standard deviation
+cutoff: 10.
+# min number of retrievals in grid box in order to keep it
+nmin: 1
+
+
diff --git a/prep_satellite/average.f90 b/prep_satellite/average.f90
index 335590a..166b173 100644
--- a/prep_satellite/average.f90
+++ b/prep_satellite/average.f90
@@ -39,8 +39,9 @@
!!
!---------------------------------------------------------------------------------------
-subroutine average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, vmr_out, verr_out,&
- vapri_out, cpri, cakpri, cak_out, cdryair, vdryair_out, nretr, ncoord, nlayer)
+subroutine average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, vmr_out, verr_out, &
+ vapri_out, cpri, cakpri, cak_out, cdryair, vdryair_out, albedo_out, nretr, &
+ ncoord, nlayer, numbox)
use mod_settings
use mod_var
@@ -53,21 +54,22 @@ subroutine average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, v
real, dimension(maxretr,ncoord), intent(in out) :: xpoints, ypoints
real, dimension(maxretr,nlayer), intent(in out) :: zpoint1, zpoint2
real, dimension(maxretr,nlayer), intent(in out) :: vapri_out, vdryair_out, cak_out
- real, dimension(maxretr), intent(in out) :: vmr_out, verr_out
+ real, dimension(maxretr), intent(in out) :: vmr_out, verr_out, albedo_out
real, dimension(maxretr), intent(in out) :: cdryair, cpri, cakpri
- integer, intent(in out) :: nretr
+ integer, intent(in) :: nretr
integer, intent(in) :: ncoord, nlayer
+ integer, intent(out) :: numbox
character(max_name_len) :: line
real :: llx, lly, urx, ury, xav, yav
real, dimension(settings%nsteps) :: stepres
real :: res, box_mean, box_sd, box_err, rwork, area
- integer :: nbox, nbox_new, num, numbox, numx, numy
+ integer :: nbox, nbox_new, num, numx, numy
integer :: ix, jy, ixy, n, nn, nr, ns, nb, nl, ierr
integer :: jjjjmmdd, hhmiss
integer, dimension(:), allocatable :: retr_nbox, iretr, iretr_keep, mark_box, idate_mean, itime_mean
real, dimension(:), allocatable :: llx_box, lly_box, res_box, warea
- real, dimension(:), allocatable :: vmr_mean, vmr_sd, cpri_mean, cdryair_mean, cakpri_mean
+ real, dimension(:), allocatable :: vmr_mean, vmr_sd, cpri_mean, cdryair_mean, cakpri_mean, albedo_mean
real, dimension(:,:), allocatable :: cak_mean, vdryair_mean, xpts_box, ypts_box, zpt1_box, zpt2_box
integer, dimension(4) :: num_subbox
real(kind=8) :: jdmean
@@ -106,9 +108,7 @@ subroutine average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, v
allocate( iretr(nretr) )
allocate( iretr_keep(nretr) )
res = stepres(settings%nsteps)
-! numx = (urx - llx)/res
- numx = (urx - llx)/res + 1 ! to cover whole domain
-! numy = (ury - lly)/res
+ numx = (urx - llx)/res + 1 ! +1 to cover whole domain
numy = (ury - lly)/res + 1
nbox = numx*numy
if ( nbox.gt.maxbox ) then
@@ -197,7 +197,6 @@ subroutine average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, v
num_subbox(n) = num_subbox(n) + 1
endif
end do ! n
-! endif ! nb
end do ! nn
print*, 'num_subbox = ',num_subbox(:)
if ( all(num_subbox.gt.0) ) then
@@ -319,6 +318,7 @@ subroutine average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, v
allocate( cpri_mean(nbox) )
allocate( cdryair_mean(nbox) )
allocate( cakpri_mean(nbox) )
+ allocate( albedo_mean(nbox) )
allocate( cak_mean(nbox,nlayer) )
allocate( vdryair_mean(nbox,nlayer) )
allocate( xpts_box(nbox,ncoord) )
@@ -402,9 +402,13 @@ subroutine average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, v
! mean column prior convolved with averaging kernel
rwork = 0.
do n = 1, num
- if ( settings%satellite.eq.'tropomi' ) then
+ if ( settings%satellite.eq.'tropomi'.or. &
+ settings%satellite.eq.'s5p_blnd' ) then
+ ! retrieval products with AK and prior profile in mol/m2
rwork = rwork + warea(n)*dot_product(cak_out(iretr(n),:),vapri_out(iretr(n),:))
- else if (settings%satellite.eq.'bremen' ) then
+ else if (settings%satellite.eq.'bremen'.or. &
+ settings%satellite.eq.'oco2' ) then
+ ! retrieval products with AK and prior profile in mole fraction
do nl = 1, nlayer
rwork = rwork + warea(n)*vdryair_out(iretr(n),nl)*vapri_out(iretr(n),nl)*cak_out(iretr(n),nl)
end do
@@ -419,6 +423,8 @@ subroutine average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, v
cdryair_mean(numbox) = dot_product(warea(1:num),cdryair(iretr(1:num)))
! mean dry air concentration
vdryair_mean(numbox,:) = matmul(warea(1:num),vdryair_out(iretr(1:num),:))
+ ! mean albedo
+ albedo_mean(numbox) = dot_product(warea(1:num),albedo_out(iretr(1:num)))
! box coordinates
xpts_box(numbox,(/1,4/)) = llx_box(nb)
xpts_box(numbox,(/2,3/)) = llx_box(nb) + res_box(nb)
@@ -436,27 +442,26 @@ subroutine average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, v
idate_mean(numbox) = jjjjmmdd
itime_mean(numbox) = hhmiss
end do
- print*, 'zpt1_box(1:10,:) = ',zpt1_box(1:10,:)
- nretr = numbox
- write(logid,*) 'Number of averaged retrievals = ',nretr
- print*, 'nretr = ',nretr
- vmr_out(1:nretr) = vmr_mean(1:nretr)
- verr_out(1:nretr) = vmr_sd(1:nretr)
- cdryair(1:nretr) = cdryair_mean(1:nretr)
- vdryair_out(1:nretr,:) = vdryair_mean(1:nretr,:)
- cak_out(1:nretr,:) = cak_mean(1:nretr,:)
- cpri(1:nretr) = cpri_mean(1:nretr)
- cakpri(1:nretr) = cakpri_mean(1:nretr)
- xpoints(1:nretr,:) = xpts_box(1:nretr,:)
- ypoints(1:nretr,:) = ypts_box(1:nretr,:)
- zpoint1(1:nretr,:) = zpt1_box(1:nretr,:)
- zpoint2(1:nretr,:) = zpt2_box(1:nretr,:)
- idate(1:nretr) = idate_mean(1:nretr)
- itime(1:nretr) = itime_mean(1:nretr)
+ write(logid,*) 'Number of averaged retrievals = ',numbox
+ vmr_out(1:numbox) = vmr_mean(1:numbox)
+ verr_out(1:numbox) = vmr_sd(1:numbox)
+ cdryair(1:numbox) = cdryair_mean(1:numbox)
+ albedo_out(1:numbox) = albedo_mean(1:numbox)
+ vdryair_out(1:numbox,:) = vdryair_mean(1:numbox,:)
+ cak_out(1:numbox,:) = cak_mean(1:numbox,:)
+ cpri(1:numbox) = cpri_mean(1:numbox)
+ cakpri(1:numbox) = cakpri_mean(1:numbox)
+ xpoints(1:numbox,:) = xpts_box(1:numbox,:)
+ ypoints(1:numbox,:) = ypts_box(1:numbox,:)
+ zpoint1(1:numbox,:) = zpt1_box(1:numbox,:)
+ zpoint2(1:numbox,:) = zpt2_box(1:numbox,:)
+ idate(1:numbox) = idate_mean(1:numbox)
+ itime(1:numbox) = itime_mean(1:numbox)
deallocate( vmr_mean )
deallocate( vmr_sd )
deallocate( cdryair_mean )
deallocate( vdryair_mean )
+ deallocate( albedo_mean )
deallocate( cak_mean )
deallocate( cpri_mean )
deallocate( cakpri_mean )
diff --git a/prep_satellite/get_bremen.f90 b/prep_satellite/get_bremen.f90
index 5592546..4328353 100644
--- a/prep_satellite/get_bremen.f90
+++ b/prep_satellite/get_bremen.f90
@@ -53,7 +53,7 @@ subroutine get_bremen(settings)
character(len=3) :: aspec
character(len=6), dimension(4) :: units
integer :: i, ierr, n, nf, nl, nm, np, nt, nbits, iswater
- integer :: nfiles, nread, nretr, npx_in_scan
+ integer :: nfiles, nread, nretr, npx_in_scan, numbox
real(kind=8) :: jdatei, jdatef, jd, days_since_ref, jdateref, jdmeas
real :: sec_since_midnight
integer :: jjjjmmdd, hhmiss, hh, mi, ss
@@ -61,12 +61,12 @@ subroutine get_bremen(settings)
integer, parameter :: lsynctime = 600
integer :: ntime, nmeas, npixel, nlevel, ncoord, nlayer
real, dimension(:,:), allocatable :: lat_bnds, lon_bnds
- real, dimension(:), allocatable :: vmr, verr, time, landfrac, qa
+ real, dimension(:), allocatable :: vmr, verr, time, landfrac, qa, albedo
real, dimension(:,:), allocatable :: pres, vapri, vdryair, cak
real, dimension(:,:), allocatable :: xpoints, ypoints
real, dimension(:,:), allocatable :: zpoint1, zpoint2
real, dimension(:,:), allocatable :: cak_out, vapri_out, vdryair_out
- real, dimension(:), allocatable :: vmr_out, verr_out
+ real, dimension(:), allocatable :: vmr_out, verr_out, albedo_out
real, dimension(:), allocatable :: cpri, cakpri, cdryair
integer, dimension(:), allocatable :: idate, itime
integer :: specnum, nchar
@@ -159,8 +159,8 @@ subroutine get_bremen(settings)
cycle read_loop
endif
- ! loop over files
- ! ---------------
+ ! loop over files for current day
+ ! -------------------------------
do n = 1, nread
@@ -228,8 +228,13 @@ subroutine get_bremen(settings)
call check( nf90_inq_varid(ncid,trim(settings%dryair_name),varid) )
call check( nf90_get_var(ncid,varid,vdryair) )
allocate( landfrac(nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate landfrac'
call check( nf90_inq_varid(ncid,trim(settings%landfrac_name),varid) )
call check( nf90_get_var(ncid,varid,landfrac) )
+ allocate( albedo(nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate albedo'
+ call check( nf90_inq_varid(ncid,trim(settings%albedo_name),varid) )
+ call check( nf90_get_var(ncid,varid,albedo) )
call check( nf90_close(ncid))
@@ -257,6 +262,8 @@ subroutine get_bremen(settings)
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vmr_out'
allocate( verr_out(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate verr_out'
+ allocate( albedo_out(maxretr), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate albedo_out'
nretr = 0
endif
do nm = 1, nmeas
@@ -306,6 +313,7 @@ subroutine get_bremen(settings)
vmr_out(nretr) = vmr(nm)
! unlike official retrieval no factor of 2 necessary for error
verr_out(nretr) = verr(nm)
+ albedo_out(nretr) = albedo(nm)
end do
print*, 'range(pres) = ',minval(pres),maxval(pres)
print*, 'range(zpoint1) = ',minval(zpoint1),maxval(zpoint1)
@@ -323,22 +331,20 @@ subroutine get_bremen(settings)
if ( allocated(vmr) ) deallocate(vmr)
if ( allocated(verr) ) deallocate(verr)
if ( allocated(landfrac) ) deallocate(landfrac)
+ if ( allocated(albedo) ) deallocate(albedo)
end do ! nread
if ( nretr.eq.0 ) go to 10
! calculate column prior (units ppbv)
- allocate( cpri(nretr), stat=ierr )
+ allocate( cpri(maxretr), stat=ierr )
do nm = 1, nretr
cpri(nm) = dot_product(vapri_out(nm,:), vdryair_out(nm,:))
end do
- !! test
- print*, 'vdryair_out(1:10,:) = ',vdryair_out(1:10,:)
-
! calculate column prior convolved with averaging kernel (units ppbv)
- allocate( cakpri(nretr), stat=ierr )
+ allocate( cakpri(maxretr), stat=ierr )
do nm = 1, nretr
do nl = 1, nlayer
cakpri(nm) = cakpri(nm) + vdryair_out(nm,nl)*vapri_out(nm,nl)*cak_out(nm,nl)
@@ -346,38 +352,29 @@ subroutine get_bremen(settings)
end do
! column dry-air not needed but only as dummy variable
- allocate( cdryair(nretr) )
+ allocate( cdryair(maxretr) )
cdryair(:) = 1.
! average pixels
! --------------
if ( settings%avg_pixels ) then
- call average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, vmr_out, verr_out,&
- vapri_out, cpri, cakpri, cak_out, cdryair, vdryair_out, nretr, ncoord, nlayer)
+ ! nretr is number of individual observations
+ ! after averaging numbox is number of super observations
+ call average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, vmr_out, verr_out, &
+ vapri_out, cpri, cakpri, cak_out, cdryair, vdryair_out, albedo_out, nretr, &
+ ncoord, nlayer, numbox)
+ else
+ ! no averaging
+ numbox = nretr
endif
- !! testing
-! open(500,file='bremen_cak.txt',action='write',status='replace')
-! do nm = 1, nretr
-! write(500,fmt='(20(F7.5,1X))') cak_out(nm,:)
-! end do
-! close(500)
-! open(500,file='bremen_cpri.txt',action='write',status='replace')
-! do nm = 1, nretr
-! write(500,fmt='(20(F7.2,1X))') vapri_out(nm,:)
-! end do
-! close(500)
-! stop
- !!
-
-
! write output data
! -----------------
! adjust itime
! needed so that with flexpart rounding to lsynctime all retrievals are still for current day
- do nm = 1, nretr
+ do nm = 1, numbox
if ( itime(nm).gt.(235959-lsynctime*100/60) ) then
itime(nm) = 235950-lsynctime*100/60
endif
@@ -402,12 +399,12 @@ subroutine get_bremen(settings)
units(4) = 'ppbv'
! write releases data for flexpart
- call write_releases(settings,nretr,nlayer,ncoord,xpoints,ypoints,zpoint1,zpoint2,&
+ call write_releases(settings,numbox,nlayer,ncoord,xpoints,ypoints,zpoint1,zpoint2,&
idate,itime,cak_out,vdryair_out,specnum,units)
! write retrieval data for flexinvert
- call write_retrieval(settings,nretr,ncoord,xpoints,ypoints,idate,itime,&
- cpri,cakpri,vmr_out,verr_out,cdryair,units)
+ call write_retrieval(settings,numbox,ncoord,xpoints,ypoints,idate,itime,&
+ cpri,cakpri,vmr_out,verr_out,cdryair,albedo_out,units)
! copy standard input files to options folder
write(aspec,fmt='(I3.3)') specnum
@@ -437,6 +434,7 @@ subroutine get_bremen(settings)
if ( allocated(cpri) ) deallocate(cpri)
if ( allocated(cakpri) ) deallocate(cakpri)
if ( allocated(cdryair) ) deallocate(cdryair)
+ if ( allocated(albedo_out) ) deallocate(albedo_out)
! iterate over days
jd = jd + 1d0
diff --git a/prep_satellite/get_oco2.f90 b/prep_satellite/get_oco2.f90
index 9f3c00c..8ea646a 100644
--- a/prep_satellite/get_oco2.f90
+++ b/prep_satellite/get_oco2.f90
@@ -56,7 +56,7 @@ subroutine get_oco2(settings)
character(len=3) :: aspec
character(len=6), dimension(4) :: units
integer :: i, ierr, n, nf, nl, ln, nm, np, nt
- integer :: nfiles, nread, nretr
+ integer :: nfiles, nread, nretr, numbox
real(kind=8) :: jdatei, jdatef, jd, jdateref, jdmeas
integer :: jjjjmmdd, jjmmdd, hhmiss
integer :: hh, mi, ss
@@ -65,12 +65,12 @@ subroutine get_oco2(settings)
integer :: ntime, nmeas, nlevel, ncoord
real(kind=8), dimension(:), allocatable :: time
real, dimension(:,:), allocatable :: lat_bnds, lon_bnds
- real, dimension(:), allocatable :: vmr, verr, landfrac, qa
+ real, dimension(:), allocatable :: vmr, verr, landfrac, qa
real, dimension(:,:), allocatable :: pres, vapri, vdryair, cak
real, dimension(:,:), allocatable :: xpoints, ypoints
real, dimension(:,:), allocatable :: zpoint1, zpoint2
real, dimension(:,:), allocatable :: cak_out, vapri_out, vdryair_out
- real, dimension(:), allocatable :: vmr_out, verr_out
+ real, dimension(:), allocatable :: vmr_out, verr_out, albedo_out ! albedo is a dummy variable
real, dimension(:), allocatable :: cpri, cakpri, cdryair
integer, dimension(:), allocatable :: idate, itime, oper_mode
integer :: specnum, nchar
@@ -284,6 +284,9 @@ subroutine get_oco2(settings)
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vmr_out'
allocate( verr_out(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate verr_out'
+ allocate( albedo_out(maxretr), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate albedo_out'
+ albedo_out(:) = 0. ! dummy variable
nretr = 0
! checking
@@ -358,16 +361,13 @@ subroutine get_oco2(settings)
if ( nretr.eq.0 ) go to 10
! calculate column prior (units ppbv)
- allocate( cpri(nretr), stat=ierr )
+ allocate( cpri(maxretr), stat=ierr )
do nm = 1, nretr
cpri(nm) = dot_product(vapri_out(nm,:), vdryair_out(nm,:))
end do
- !! test
- print*, 'vdryair_out(1:10,:) = ',vdryair_out(1:10,:)
-
! calculate column prior convolved with averaging kernel (units ppbv)
- allocate( cakpri(nretr), stat=ierr )
+ allocate( cakpri(maxretr), stat=ierr )
do nm = 1, nretr
do nl = 1, nlevel
cakpri(nm) = cakpri(nm) + vdryair_out(nm,nl)*vapri_out(nm,nl)*cak_out(nm,nl)
@@ -375,15 +375,21 @@ subroutine get_oco2(settings)
end do
! column dry-air not needed but only as dummy variable
- allocate( cdryair(nretr) )
+ allocate( cdryair(maxretr) )
cdryair(:) = 1.
! average pixels
! --------------
if ( settings%avg_pixels ) then
- call average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, vmr_out, verr_out,&
- vapri_out, cpri, cakpri, cak_out, cdryair, vdryair_out, nretr, ncoord, nlevel)
+ ! nretr is number of individual observations
+ ! after averaging numbox is number of super observations
+ call average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, vmr_out, verr_out, &
+ vapri_out, cpri, cakpri, cak_out, cdryair, vdryair_out, albedo_out, nretr, &
+ ncoord, nlevel, numbox)
+ else
+ ! no averaging
+ numbox = nretr
endif
! write output data
@@ -391,7 +397,7 @@ subroutine get_oco2(settings)
! adjust itime
! needed so that with flexpart rounding to lsynctime all retrievals are still for current day
- do nm = 1, nretr
+ do nm = 1, numbox
if ( itime(nm).gt.(235959-lsynctime*100/60) ) then
itime(nm) = 235950-lsynctime*100/60
endif
@@ -416,12 +422,12 @@ subroutine get_oco2(settings)
units(4) = 'ppmv'
! write releases data for flexpart
- call write_releases(settings,nretr,nlevel,ncoord,xpoints,ypoints,zpoint1,zpoint2,&
+ call write_releases(settings,numbox,nlevel,ncoord,xpoints,ypoints,zpoint1,zpoint2,&
idate,itime,cak_out,vdryair_out,specnum,units)
! write retrieval data for flexinvert
- call write_retrieval(settings,nretr,ncoord,xpoints,ypoints,idate,itime,&
- cpri,cakpri,vmr_out,verr_out,cdryair,units)
+ call write_retrieval(settings,numbox,ncoord,xpoints,ypoints,idate,itime,&
+ cpri,cakpri,vmr_out,verr_out,cdryair,albedo_out,units)
! copy standard input files to options folder
write(aspec,fmt='(I3.3)') specnum
@@ -451,6 +457,7 @@ subroutine get_oco2(settings)
if ( allocated(cpri) ) deallocate(cpri)
if ( allocated(cakpri) ) deallocate(cakpri)
if ( allocated(cdryair) ) deallocate(cdryair)
+ if ( allocated(albedo_out) ) deallocate(albedo_out)
! iterate over days
jd = jd + 1d0
diff --git a/prep_satellite/get_tropomi.f90 b/prep_satellite/get_tropomi.f90
index 850d7c6..d20aefa 100644
--- a/prep_satellite/get_tropomi.f90
+++ b/prep_satellite/get_tropomi.f90
@@ -56,7 +56,7 @@ subroutine get_tropomi(settings)
character(len=2) :: amon, aday
character(len=6), dimension(4) :: units
integer :: i, ierr, n, nf, nl, nm, np, nt, nbits, iswater
- integer :: nfiles, nread, nretr
+ integer :: nfiles, nread, nretr, numbox
real(kind=8) :: jdatei, jdatef, jd, ndays, jref, jdmeas
integer :: jjjjmmdd, hhmiss, mm, dd, hh, mi, ss
integer, parameter :: refdate = 20100101
@@ -64,13 +64,13 @@ subroutine get_tropomi(settings)
integer, parameter :: lsynctime = 600
integer :: ntime, nmeas, npixel, nlevel, ncoord, nlayer
real, dimension(:,:,:,:), allocatable :: lat_bnds, lon_bnds
- real, dimension(:,:,:), allocatable :: psurf, pdel, vmr, verr, qa
+ real, dimension(:,:,:), allocatable :: psurf, pdel, vmr, verr, qa, albedo
real, dimension(:,:,:,:), allocatable :: vapri, vdryair, cak
real, dimension(:,:), allocatable :: xpoints, ypoints
real, dimension(:,:), allocatable :: zpoint1, zpoint2
real, dimension(:,:), allocatable :: cak_out, vapri_out, vdryair_out
real, dimension(:,:), allocatable :: cak_tmp, vapri_tmp, vdryair_tmp
- real, dimension(:), allocatable :: vmr_out, verr_out
+ real, dimension(:), allocatable :: vmr_out, verr_out, albedo_out
real, dimension(:), allocatable :: cpri, cakpri, cdryair
integer, dimension(:), allocatable :: idate, itime, time
integer, dimension(:,:), allocatable :: timedel
@@ -80,8 +80,10 @@ subroutine get_tropomi(settings)
integer :: dimid, varid
real :: nest_llx, nest_lly, nest_urx, nest_ury
real :: llx, lly, urx, ury
- real :: sclfact, offset
+ real :: sclfact, offset
real, parameter :: bignum=huge(0.)
+ real :: const, alb, factor_corr
+ real, dimension(11) :: poly
write(logid,*) 'largest number that can be handled (used to define NA) = ',bignum
@@ -162,6 +164,12 @@ subroutine get_tropomi(settings)
write(aday,fmt='(I2.2)') dd
nread = 0
do n = 1, nfiles
+ ! corrupted files have "xxx" at start of name
+ i=len_trim(settings%path_obs)
+ if ( filelist(n)((i+1):(i+3)).eq.'xxx' ) then
+ print*, 'skipping file: ',filelist(n)
+ cycle
+ endif
do i = 1, len_trim(filelist(n))-18
if ( filelist(n)(i:(i+7)).eq.adate ) then
nread = nread + 1
@@ -176,8 +184,8 @@ subroutine get_tropomi(settings)
cycle read_loop
endif
- ! loop over files
- ! ---------------
+ ! loop over files for current day
+ ! --------------------------------
do n = 1, nread
@@ -302,6 +310,16 @@ subroutine get_tropomi(settings)
call check( nf90_inq_varid(det_grpid,trim(string),varid) )
call check( nf90_get_var(det_grpid,varid,cak) )
write(logid,*) 'read column averaging kernel'
+ allocate( albedo(npixel,nmeas,ntime), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate albedo'
+ string = settings%albedo_name
+ call split(string,'/',grpname)
+ call split(string,'/',suppdata_grpname)
+ call split(string,'/',details_grpname)
+ call check( nf90_inq_ncid(supp_grpid,trim(details_grpname),det_grpid) )
+ call check( nf90_inq_varid(det_grpid,trim(string),varid) )
+ call check( nf90_get_var(det_grpid,varid,albedo) )
+ write(logid,*) 'read surface albedo SWIR'
! product variables
allocate( qa(npixel,nmeas,ntime), stat=ierr )
@@ -373,10 +391,13 @@ subroutine get_tropomi(settings)
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vmr_out'
allocate( verr_out(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate verr_out'
+ allocate( albedo_out(maxretr), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate albedo_out'
cak_out(:,:)=0.
vapri_out(:,:)=0.
vmr_out(:)=0.
verr_out(:)=0.
+ albedo_out(:)=0.
nretr = 0
endif
do nt = 1, ntime
@@ -431,6 +452,7 @@ subroutine get_tropomi(settings)
zpoint2(nretr,nl) = zpoint2(nretr,nl-1) - pdel(np,nm,nt)
end do
cak_out(nretr,:) = cak(:,np,nm,nt)
+ albedo_out(nretr) = albedo(np,nm,nt)
vapri_out(nretr,:) = vapri(:,np,nm,nt)
vdryair_out(nretr,:) = vdryair(:,np,nm,nt)
vmr_out(nretr) = vmr(np,nm,nt)
@@ -450,6 +472,7 @@ subroutine get_tropomi(settings)
if ( allocated(vapri) ) deallocate(vapri)
if ( allocated(vdryair) ) deallocate(vdryair)
if ( allocated(cak) ) deallocate(cak)
+ if ( allocated(albedo) ) deallocate(albedo)
if ( allocated(timedel) ) deallocate(timedel)
if ( allocated(time) ) deallocate(time)
if ( allocated(qa) ) deallocate(qa)
@@ -479,49 +502,70 @@ subroutine get_tropomi(settings)
vdryair_out = vdryair_tmp
deallocate( cak_tmp, vapri_tmp, vdryair_tmp)
- print*, 'itime(1), itime(nretr) = ',itime(1),itime(nretr)
-
! calculate column prior (units mol/m2)
- allocate( cpri(nretr), stat=ierr )
- cpri = sum(vapri_out(1:nretr,:), dim=2)
+ allocate( cpri(maxretr), stat=ierr )
+ cpri(1:nretr) = sum(vapri_out(1:nretr,:), dim=2)
! calculate column prior convolved with averaging kernel (units mol/m2)
- allocate( cakpri(nretr), stat=ierr )
+ allocate( cakpri(maxretr), stat=ierr )
do nm = 1, nretr
cakpri(nm) = dot_product(vapri_out(nm,:), cak_out(nm,:))
end do
+ ! apply albedo correction used in Official retrievals to cak, cpri and cakpri
+ ! this is polynomial function order of 11 (Soumyajit Mandal, SRON)
+! const = 1.0282413517102085
+! poly(1) = -0.4717177565799524
+! poly(2) = 2.1491179428159324
+! poly(3) = 16.952582059603163
+! poly(4) = -250.66075995913351
+! poly(5) = 1307.5191440027886
+! poly(6) =-3791.6452849323673
+! poly(7) =6783.218182991253
+! poly(8) =-7689.387517819326
+! poly(9) = 5404.118866385074
+! poly(10) =-2154.928509865762
+! poly(11) =373.4579485937264
+! if ( settings%albedo_corr ) then
+! open(100,file=trim(settings%path_obsout)//'aldbedo_corr_'//adate//'.txt', &
+! status='replace',action='write')
+! do nm = 1, nretr
+! alb = albedo_out(nm)
+! factor_corr = const + poly(1)*alb + poly(2)*alb**2 + poly(3)*alb**3 + poly(4)*alb**4 + &
+! poly(5)*alb**5 + poly(6)*alb**6 + poly(7)*alb**7 + poly(8)*alb**8 + &
+! poly(9)*alb**9 + poly(10)*alb**10 + poly(11)*alb**11
+! cakpri(nm) = factor_corr*cakpri(nm)
+! cpri(nm) = factor_corr*cpri(nm)
+! cak_out(nm,:) = factor_corr*cak_out(nm,:)
+! write(100,'(F6.4,1X,F6.4)') factor_corr, alb
+! end do
+! close(100)
+! endif
+
! calculate column dry air (units mol/m2)
- allocate( cdryair(nretr), stat=ierr )
- cdryair = sum(vdryair_out(1:nretr,:), dim=2)
+ allocate( cdryair(maxretr), stat=ierr )
+ cdryair(1:nretr) = sum(vdryair_out(1:nretr,:), dim=2)
! average pixels
! --------------
if ( settings%avg_pixels ) then
- call average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, vmr_out, verr_out,&
- vapri_out, cpri, cakpri, cak_out, cdryair, vdryair_out, nretr, ncoord, nlayer)
+ ! nretr is number of individual observations
+ ! after averaging numbox is number of super observations
+ call average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, vmr_out, verr_out, &
+ vapri_out, cpri, cakpri, cak_out, cdryair, vdryair_out, albedo_out, nretr, &
+ ncoord, nlayer, numbox)
+ else
+ ! no averaging
+ numbox = nretr
endif
- !! testing
-! open(500,file='tropomi_cak.txt',action='write',status='replace')
-! do nm = 1, nretr
-! write(500,fmt='(12(F7.5,1X))') cak_out(nm,:)
-! end do
-! close(500)
-! open(500,file='tropomi_cpri.txt',action='write',status='replace')
-! do nm = 1, nretr
-! write(500,fmt='(12(F7.2,1X))') vapri_out(nm,:)*1.e9/vdryair_out(nm,:)
-! end do
-! close(500)
-! !!
-
! write output data
! -----------------
! adjust itime
! needed so that with flexpart rounding to lsynctime all retrievals are still for current day
- do nm = 1, nretr
+ do nm = 1, numbox
if ( itime(nm).gt.(235959-lsynctime*100/60) ) then
itime(nm) = 235950-lsynctime*100/60
endif
@@ -546,12 +590,12 @@ subroutine get_tropomi(settings)
units(4) = 'ppbv'
! write releases data for flexpart
- call write_releases(settings,nretr,nlayer,ncoord,xpoints,ypoints,zpoint1,zpoint2,&
+ call write_releases(settings,numbox,nlayer,ncoord,xpoints,ypoints,zpoint1,zpoint2,&
idate,itime,cak_out,vdryair_out,specnum,units)
! write retrieval data for flexinvert
- call write_retrieval(settings,nretr,ncoord,xpoints,ypoints,idate,itime,&
- cpri,cakpri,vmr_out,verr_out,cdryair,units)
+ call write_retrieval(settings,numbox,ncoord,xpoints,ypoints,idate,itime,&
+ cpri,cakpri,vmr_out,verr_out,cdryair,albedo_out,units)
! copy standard input files to options folder
write(aspec,fmt='(I3.3)') specnum
@@ -574,6 +618,7 @@ subroutine get_tropomi(settings)
if ( allocated(vapri_out) ) deallocate(vapri_out)
if ( allocated(vdryair_out) ) deallocate(vdryair_out)
if ( allocated(vmr_out) ) deallocate(vmr_out)
+ if ( allocated(albedo_out) ) deallocate(albedo_out)
if ( allocated(verr_out) ) deallocate(verr_out)
if ( allocated(cpri) ) deallocate(cpri)
if ( allocated(cakpri) ) deallocate(cakpri)
diff --git a/prep_satellite/get_tropomi_blnd.f90 b/prep_satellite/get_tropomi_blnd.f90
new file mode 100644
index 0000000..d935e02
--- /dev/null
+++ b/prep_satellite/get_tropomi_blnd.f90
@@ -0,0 +1,518 @@
+!---------------------------------------------------------------------------------------
+! PREP_SATELLITE: get_tropomi_blnd
+!---------------------------------------------------------------------------------------
+! FLEXINVERT is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! FLEXINVERT is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 2022, Rona Thompson
+!---------------------------------------------------------------------------------------
+!
+!> get_tropomi_blnd
+!!
+!! Purpose: Reads observations and auxilliary data from the specified satellite and
+!! and writes the output into files that can be directly used by FLEXPART
+!! and, for inversions, by FLEXINVERT
+!!
+!! Interface:
+!!
+!! Inputs
+!! settings - settings data structure
+!!
+!! Outputs
+!!
+!! Externals
+!!
+!---------------------------------------------------------------------------------------
+
+subroutine get_tropomi_blnd(settings)
+
+ use mod_settings
+ use mod_var
+ use mod_dates
+ use mod_save
+ use netcdf
+
+ implicit none
+
+ type (settings_t), intent(in) :: settings
+
+ character(max_name_len), dimension(:), allocatable :: filelist, filereadlist
+ character(max_name_len) :: line, filename
+ character(len=27), dimension(:), allocatable :: timestr
+ character(len=8) :: adate
+ character(len=6) :: atime
+ character(len=3) :: aspec
+ character(len=6), dimension(4) :: units
+ integer :: i, ierr, n, nf, nl, nm, np, nt, nbits, iswater
+ integer :: nfiles, nread, nretr, npx_in_scan, numbox
+ real(kind=8) :: jdatei, jdatef, jd, days_since_ref, jdateref, jdmeas
+ real :: sec_since_midnight
+ integer :: jjjjmmdd, hhmiss, jjjj, mm, dd, hh, mi, ss
+ integer, parameter :: refdate = 19700101
+ integer, parameter :: lsynctime = 600
+ integer :: ntime, nmeas, ncoord, nlayer
+ real, dimension(:,:), allocatable :: lat_bnds, lon_bnds
+ real, dimension(:), allocatable :: vmr, verr, time, qa, albedo
+ real, dimension(:,:), allocatable :: vapri, vdryair, cak
+ real, dimension(:), allocatable :: psurf, pdel
+ real, dimension(:,:), allocatable :: xpoints, ypoints
+ real, dimension(:,:), allocatable :: zpoint1, zpoint2
+ real, dimension(:,:), allocatable :: cak_out, vapri_out, vdryair_out
+ real, dimension(:,:), allocatable :: cak_tmp, vapri_tmp, vdryair_tmp
+ real, dimension(:), allocatable :: vmr_out, verr_out, albedo_out
+ real, dimension(:), allocatable :: cpri, cakpri, cdryair
+ integer, dimension(:), allocatable :: idate, itime, surfclass
+ character(len=8) :: surfclass_bin, qa_bin
+ integer :: specnum, nchar
+ integer :: ncid, dimid, varid
+ real :: nest_llx, nest_lly, nest_urx, nest_ury
+ real :: llx, lly, urx, ury
+ integer :: nargs, cnt
+ character(max_name_len), dimension(2) :: args
+ character(len=28) :: str
+ real :: sclfact, offset
+ character(len=2) :: anread
+
+ ! get species number
+ ! ------------------
+
+ open(100,file=trim(settings%path_flexpart)//'options/SPECIES/spec_overview',status='old',action='read',iostat=ierr)
+ nchar = len_trim(settings%species)
+ do while ( ierr.eq.0 )
+ read (100, fmt='(A)', iostat=ierr) line
+ if ( line(len_trim(line)-nchar+1:len_trim(line)) == trim(settings%species) ) ierr = 1
+ end do
+ read(line(9:11),*) specnum
+ write(logid,*) 'species: ',settings%species
+ write(logid,*) 'specnum: ',specnum
+ write(logid,*) 'qa cutoff: ',settings%qa_cutoff
+ close(100)
+
+ ! define domains
+ ! --------------------
+
+ llx = settings%outlonleft
+ lly = settings%outlatlower
+ urx = llx + settings%numxgrid*settings%dxout
+ ury = lly + settings%numygrid*settings%dyout
+
+ if ( settings%lnested.eq.1 ) then
+ nest_llx = settings%outlonnest
+ nest_lly = settings%outlatnest
+ nest_urx = nest_llx + settings%numxnest*settings%dxoutnest
+ nest_ury = nest_lly + settings%numynest*settings%dyoutnest
+ endif
+ print*, 'llx,lly,urx,ury = ',llx,lly,urx,ury
+ if (settings%lnested.eq.1) print*, 'nest_llx,nest_lly,nest_urx,nest_ury = ',nest_llx,nest_lly,nest_urx,nest_ury
+
+ ! list satellite files
+ ! ----------------------
+
+ call system('find '//trim(settings%path_obs)//' -type f | grep '//trim(settings%suffix)//&
+ ' | wc -l > '//trim(settings%path_obsout)//'obsfiles.txt')
+ call system('find '//trim(settings%path_obs)//' -type f | grep '//trim(settings%suffix)//&
+ ' >> '//trim(settings%path_obsout)//'obsfiles.txt')
+
+ open(100,file=trim(settings%path_obsout)//'obsfiles.txt',action='read',status='old',iostat=ierr)
+ if ( ierr.ne. 0 ) then
+ write(logid,*) 'ERROR: cannot open obsfiles.txt'
+ stop
+ endif
+ read(100,*,iostat=ierr) nfiles
+ allocate ( filelist(nfiles), stat = ierr )
+ allocate ( filereadlist(nfiles), stat = ierr )
+ do nf = 1, nfiles
+ read(100,fmt='(A)',iostat=ierr) filelist(nf)
+ if (ierr.ne.0) exit
+ end do
+ close(100)
+
+ ! loop over days
+ ! --------------
+
+ jdatei = juldate(settings%datei, 0)
+ jdatef = juldate(settings%datef, 0)
+ jd = jdatei
+
+ read_loop: do
+
+ if ( jd.gt.jdatef ) exit read_loop
+
+ ! search files for current date
+ call caldate(jd, jjjjmmdd, hhmiss)
+ write(adate,fmt='(I8.8)') jjjjmmdd
+ write(logid,*) 'current day: ',adate
+ nread = 0
+ do n = 1, nfiles
+ do i = 1, len_trim(filelist(n))-18
+ if ( filelist(n)(i:(i+7)).eq.adate ) then
+ nread = nread + 1
+ filereadlist(nread) = filelist(n)
+ exit
+ endif
+ end do
+ end do
+ write(logid,*) 'number of files for current day: ',nread
+ if ( nread.eq.0 ) then
+ jd = jd + 1d0
+ cycle read_loop
+ endif
+
+ ! loop over files for current day
+ ! -------------------------------
+
+ do n = 1, nread
+
+ ! open file
+ write(logid,*) 'reading file: ',filereadlist(n)
+ call check( nf90_open(trim(filereadlist(n)),nf90_NOWRITE,ncid) )
+
+ ! dimension variables
+ call check( nf90_inq_dimid(ncid,trim(settings%nmeas_name),dimid) )
+ call check( nf90_inquire_dimension(ncid,dimid,len=nmeas) )
+ write(logid,*) 'number of retrievals: ',nmeas
+ if ( nmeas.eq.0 ) then
+ write(logid,*) 'file contains no retrievals moving to next file'
+ cycle
+ endif
+ call check( nf90_inq_dimid(ncid,trim(settings%ncoord_name),dimid) )
+ call check( nf90_inquire_dimension(ncid,dimid,len=ncoord) )
+ write(logid,*) 'number of geocoordinates: ',ncoord
+ call check( nf90_inq_dimid(ncid,trim(settings%nlevel_name),dimid) )
+ call check( nf90_inquire_dimension(ncid,dimid,len=nlayer) )
+ write(logid,*) 'number of vertical layers: ',nlayer
+
+ ! geolocation variables
+ allocate( lat_bnds(ncoord,nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate lat_bnds'
+ call check( nf90_inq_varid(ncid,trim(settings%lat_name),varid) )
+ call check( nf90_get_var(ncid,varid,lat_bnds) )
+ allocate( lon_bnds(ncoord,nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate lon_bnds'
+ call check( nf90_inq_varid(ncid,trim(settings%lon_name),varid) )
+ call check( nf90_get_var(ncid,varid,lon_bnds) )
+
+ ! product variables
+ allocate( timestr(nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate timestr'
+ ! cannot read variable type string in fortran
+! call check( nf90_inq_varid(ncid,trim(settings%time_name),varid) )
+! call check( nf90_get_var(ncid,varid,timestr) )
+ ! work around -> extract time variable to ascii
+ call system("ncdump -v "//trim(settings%time_name)//" "//trim(filereadlist(n))//&
+ " | sed -e '1,/data:/d' -e '$d' | tail -n +2 > "//trim(settings%path_obsout)//"time.txt")
+ open(100,file=trim(settings%path_obsout)//"time.txt",action='read',status='old',iostat=ierr)
+ if (ierr.ne.0) then
+ write(logid,*) 'ERROR: cannot open time.txt'
+ stop
+ endif
+ cnt=0
+ do while (ierr.eq.0)
+ read (100, fmt='(A)', iostat=ierr) line
+ if (ierr.ne.0) exit
+ ! remove pesky characters in first line
+ if (cnt==0) then
+ line = line(12:len(line))
+ endif
+ call parse_string(line, ',', args, nargs)
+ do i = 1, nargs
+ cnt=cnt+1
+ str=trim(args(i))
+ timestr(cnt)=str(2:28)
+ end do
+ end do
+ close(100)
+ call system("rm -f "//trim(settings%path_obsout)//"time.txt")
+ allocate( psurf(nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate psurf'
+ call check( nf90_inq_varid(ncid,trim(settings%psurf_name),varid) )
+ call check( nf90_get_var(ncid,varid,psurf) )
+ allocate( pdel(nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate pdel'
+ call check( nf90_inq_varid(ncid,trim(settings%pdel_name),varid) )
+ call check( nf90_get_var(ncid,varid,pdel) )
+ allocate( qa(nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate qa'
+ call check( nf90_inq_varid(ncid,trim(settings%qa_name),varid) )
+ call check( nf90_get_var(ncid,varid,qa) )
+ call check( nf90_get_att(ncid,varid,'scale_factor',sclfact) )
+ call check( nf90_get_att(ncid,varid,'add_offset',offset) )
+ qa = qa*sclfact + offset
+ print*, 'qa = ',qa(1:20)
+ allocate( vmr(nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vmr'
+ call check( nf90_inq_varid(ncid,trim(settings%vmr_name),varid) )
+ call check( nf90_get_var(ncid,varid,vmr) )
+ allocate( verr(nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate verr'
+ call check( nf90_inq_varid(ncid,trim(settings%verr_name),varid) )
+ call check( nf90_get_var(ncid,varid,verr) )
+
+ ! detailed results variables
+ allocate( cak(nlayer,nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate cak'
+ call check( nf90_inq_varid(ncid,trim(settings%ak_name),varid) )
+ call check( nf90_get_var(ncid,varid,cak) )
+ allocate( vapri(nlayer,nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vapri'
+ call check( nf90_inq_varid(ncid,trim(settings%apri_name),varid) )
+ call check( nf90_get_var(ncid,varid,vapri) )
+ allocate( vdryair(nlayer,nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vdryair'
+ call check( nf90_inq_varid(ncid,trim(settings%dryair_name),varid) )
+ call check( nf90_get_var(ncid,varid,vdryair) )
+ allocate( surfclass(nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate surfclass'
+ call check( nf90_inq_varid(ncid,trim(settings%surfclass_name),varid) )
+ call check( nf90_get_var(ncid,varid,surfclass) )
+ allocate( albedo(nmeas), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate albedo'
+ call check( nf90_inq_varid(ncid,trim(settings%albedo_name),varid) )
+ call check( nf90_get_var(ncid,varid,albedo) )
+
+ call check( nf90_close(ncid))
+
+ ! select retrievals and extract needed info
+ if( n.eq.1 ) then
+ allocate( xpoints(maxretr,ncoord), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate xpoints'
+ allocate( ypoints(maxretr,ncoord), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate ypoints'
+ allocate( zpoint1(maxretr,nlayer), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate zpoint1'
+ allocate( zpoint2(maxretr,nlayer), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate zpoint2'
+ allocate( idate(maxretr), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate idate'
+ allocate( itime(maxretr), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate itime'
+ allocate( cak_out(maxretr,nlayer), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate cak_out'
+ allocate( vapri_out(maxretr,nlayer), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vapri_out'
+ allocate( vdryair_out(maxretr,nlayer), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vdryair_out'
+ allocate( vmr_out(maxretr), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vmr_out'
+ allocate( verr_out(maxretr), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate verr_out'
+ allocate( albedo_out(maxretr), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate albedo_out'
+ nretr = 0
+ endif
+ print*, 'nmeas = ',nmeas
+ do nm = 1, nmeas
+ ! exclude retrievals above quality cutoff
+ if ( qa(nm).lt.settings%qa_cutoff ) cycle
+ write(surfclass_bin, fmt='(B0)') surfclass(nm)
+ nbits = len_trim(surfclass_bin)
+ read(surfclass_bin(nbits:nbits),*) iswater
+ if(nm.lt.10) print*, 'surfclass_bin, nbits, iswater = ',surfclass_bin, nbits, iswater
+ ! if landonly exclude pixels over water
+ if ( settings%landonly.and.iswater.eq.1 ) cycle
+ if ( settings%lnested.eq.1 ) then
+ ! if nested domain exclude pixels outside
+ if ( any(lon_bnds(:,nm).lt.nest_llx).or. &
+ any(lon_bnds(:,nm).gt.nest_urx).or. &
+ any(lat_bnds(:,nm).lt.nest_lly).or. &
+ any(lat_bnds(:,nm).gt.nest_ury) ) cycle
+ else
+ ! if not nested exclude pixels outside global domain
+ if ( any(lon_bnds(:,nm).lt.llx).or. &
+ any(lon_bnds(:,nm).gt.urx).or. &
+ any(lat_bnds(:,nm).lt.lly).or. &
+ any(lat_bnds(:,nm).gt.ury) ) cycle
+ endif
+ read(timestr(nm)(1:4),'(I4)') jjjj
+ read(timestr(nm)(6:7),'(I2)') mm
+ read(timestr(nm)(9:10),'(I2)') dd
+ read(timestr(nm)(12:13),'(I2)') hh
+ read(timestr(nm)(15:16),'(I2)') mi
+ read(timestr(nm)(18:19),'(I2)') ss
+ hhmiss=hh*10000+mi*100+ss
+ jjjjmmdd=jjjj*10000+mm*100+dd
+ jdmeas=juldate(jjjjmmdd,hhmiss)
+ ! exclude retrievals if not for current day
+ if ( floor(jdmeas).ne.jd ) cycle
+ ! check date and time
+ if (isnan(real(hhmiss)).or.(jjjjmmdd.lt.0)) then
+ write(logid,*) 'WARNING: problem with itime, idate for file, nmeas, jdmeas = ',trim(filereadlist(n)), nm, jdmeas
+ cycle
+ endif
+ nretr = nretr + 1
+ if (nretr.lt.10) print*, 'jjjjmmdd, hhmiss = ',jjjjmmdd, hhmiss
+ idate(nretr) = jjjjmmdd
+ itime(nretr) = hhmiss
+ ! corners ordered anticlockwise from lower left
+ xpoints(nretr,:) = lon_bnds(:,nm)
+ ypoints(nretr,:) = lat_bnds(:,nm)
+ ! pressure levels are equipressure apart
+ zpoint1(nretr,1) = psurf(nm)
+ zpoint2(nretr,1) = psurf(nm) - pdel(nm)
+ do nl = 2, nlayer
+ zpoint1(nretr,nl) = zpoint1(nretr,nl-1) - pdel(nm)
+ zpoint2(nretr,nl) = zpoint2(nretr,nl-1) - pdel(nm)
+ end do
+ cak_out(nretr,:) = cak(:,nm)
+ vapri_out(nretr,:) = vapri(:,nm)
+ vdryair_out(nretr,:) = vdryair(:,nm)
+ vmr_out(nretr) = vmr(nm)
+ verr_out(nretr) = verr(nm)
+ albedo_out(nretr) = albedo(nm)
+ end do
+
+ write(logid,*) 'number of retrievals: ',nretr
+
+ if ( allocated(lat_bnds) ) deallocate(lat_bnds)
+ if ( allocated(lon_bnds) ) deallocate(lon_bnds)
+ if ( allocated(vapri) ) deallocate(vapri)
+ if ( allocated(vdryair) ) deallocate(vdryair)
+ if ( allocated(cak) ) deallocate(cak)
+ if ( allocated(timestr) ) deallocate(timestr)
+ if ( allocated(psurf) ) deallocate(psurf)
+ if ( allocated(pdel) ) deallocate(pdel)
+ if ( allocated(qa) ) deallocate(qa)
+ if ( allocated(vmr) ) deallocate(vmr)
+ if ( allocated(verr) ) deallocate(verr)
+ if ( allocated(surfclass) ) deallocate(surfclass)
+ if ( allocated(albedo) ) deallocate(albedo)
+
+ end do ! nread
+
+ if ( nretr.eq.0 ) go to 10
+
+ ! reverse direction of vertical coordinate
+ ! raw data ordered from top of atmosphere to surface
+ ! which is opposite to the pressure levels calculated (zpoint)
+ allocate( cak_tmp(maxretr,nlayer), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate cak_tmp'
+ allocate( vapri_tmp(maxretr,nlayer), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vapri_tmp'
+ allocate( vdryair_tmp(maxretr,nlayer), stat=ierr )
+ if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vdryair_tmp'
+ do nl = 1, nlayer
+ cak_tmp(:,nl) = cak_out(:,(nlayer-nl+1))
+ vapri_tmp(:,nl) = vapri_out(:,(nlayer-nl+1))
+ vdryair_tmp(:,nl) = vdryair_out(:,(nlayer-nl+1))
+ end do
+ cak_out = cak_tmp
+ vapri_out = vapri_tmp
+ vdryair_out = vdryair_tmp
+ deallocate( cak_tmp, vapri_tmp, vdryair_tmp)
+
+ ! calculate column prior (units mol/m2)
+ allocate( cpri(maxretr), stat=ierr )
+ cpri(1:nretr) = sum(vapri_out(1:nretr,:), dim=2)
+
+ ! calculate column prior convolved with averaging kernel (units mol/m2)
+ allocate( cakpri(maxretr), stat=ierr )
+ do nm = 1, nretr
+ cakpri(nm) = dot_product(vapri_out(nm,:), cak_out(nm,:))
+ end do
+
+ ! calculate column dry air (units mol/m2)
+ allocate( cdryair(maxretr), stat=ierr )
+ cdryair(1:nretr) = sum(vdryair_out(1:nretr,:), dim=2)
+
+ ! average pixels
+ ! --------------
+
+ if ( settings%avg_pixels ) then
+ ! nretr is number of individual observations
+ ! after averaging numbox is number of super observations
+ call average(settings, idate, itime, xpoints, ypoints, zpoint1, zpoint2, vmr_out, verr_out, &
+ vapri_out, cpri, cakpri, cak_out, cdryair, vdryair_out, albedo_out, nretr, &
+ ncoord, nlayer, numbox)
+ else
+ ! no averaging
+ numbox = nretr
+ endif
+
+ ! write output data
+ ! -----------------
+
+ ! adjust itime
+ ! needed so that with flexpart rounding to lsynctime all retrievals are still for current day
+ do nm = 1, numbox
+ if ( itime(nm).gt.(235959-lsynctime*100/60) ) then
+ itime(nm) = 235950-lsynctime*100/60
+ endif
+ end do
+
+ ! write command file
+ call prep_command(settings,jd)
+
+ ! write outgrid file
+ call prep_outgrid(settings,jd)
+
+ ! write ageclasses file
+ call prep_ageclass(settings,jd)
+
+ ! write pathnames file
+ call prep_pathnames(settings,jd)
+
+ ! units of output data: dryair, cdryair, cpri, vmr
+ units(1) = 'mol/m2'
+ units(2) = 'mol/m2'
+ units(3) = 'mol/m2'
+ units(4) = 'ppbv'
+
+ ! write releases data for flexpart
+ call write_releases(settings,numbox,nlayer,ncoord,xpoints,ypoints,zpoint1,zpoint2,&
+ idate,itime,cak_out,vdryair_out,specnum,units)
+
+ ! write retrieval data for flexinvert
+ call write_retrieval(settings,numbox,ncoord,xpoints,ypoints,idate,itime,&
+ cpri,cakpri,vmr_out,verr_out,cdryair,albedo_out,units)
+
+ ! copy standard input files to options folder
+ write(aspec,fmt='(I3.3)') specnum
+ call system('mkdir -p '//trim(settings%path_options)//adate//'/options/SPECIES/')
+ filename = trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//aspec
+ call system('cp '//trim(filename)//' '//trim(settings%path_options)//adate//'/options/SPECIES/')
+ call system('cp '//trim(settings%path_flexpart)//'options/*.dat'//' '//trim(settings%path_options)//adate//'/options/')
+ call system('cp '//trim(settings%path_flexpart)//'options/*.t'//' '//trim(settings%path_options)//adate//'/options/')
+
+ ! no retrievals on this day
+10 continue
+
+ ! prepare for next day
+ ! --------------------
+
+ if ( allocated(xpoints) ) deallocate(xpoints)
+ if ( allocated(ypoints) ) deallocate(ypoints)
+ if ( allocated(zpoint1) ) deallocate(zpoint1)
+ if ( allocated(zpoint2) ) deallocate(zpoint2)
+ if ( allocated(idate) ) deallocate(idate)
+ if ( allocated(itime) ) deallocate(itime)
+ if ( allocated(cak_out) ) deallocate(cak_out)
+ if ( allocated(vapri_out) ) deallocate(vapri_out)
+ if ( allocated(vdryair_out) ) deallocate(vdryair_out)
+ if ( allocated(vmr_out) ) deallocate(vmr_out)
+ if ( allocated(verr_out) ) deallocate(verr_out)
+ if ( allocated(cpri) ) deallocate(cpri)
+ if ( allocated(cakpri) ) deallocate(cakpri)
+ if ( allocated(cdryair) ) deallocate(cdryair)
+ if ( allocated(albedo_out) ) deallocate(albedo_out)
+
+ ! iterate over days
+ jd = jd + 1d0
+
+ end do read_loop
+
+ if ( allocated(filelist) ) deallocate(filelist)
+ if ( allocated(filereadlist) ) deallocate(filereadlist)
+
+
+end subroutine get_tropomi_blnd
+
diff --git a/prep_satellite/main.f90 b/prep_satellite/main.f90
index 60685ff..72c1c63 100644
--- a/prep_satellite/main.f90
+++ b/prep_satellite/main.f90
@@ -73,6 +73,9 @@ program main
else if ( settings%satellite.eq.'bremen' ) then
! bremen retrieval for tropomi
call get_bremen(settings)
+ else if ( settings%satellite.eq.'s5p_blnd' ) then
+ ! tropomi-gosat blended product
+ call get_tropomi_blnd(settings)
else if ( settings%satellite.eq.'oco2' ) then
call get_oco2(settings)
else
diff --git a/prep_satellite/makefile b/prep_satellite/makefile
index e54f331..7e2d44e 100644
--- a/prep_satellite/makefile
+++ b/prep_satellite/makefile
@@ -28,6 +28,7 @@ OBJECTS = gridarea.o \
get_tropomi.o \
get_bremen.o \
get_oco2.o \
+ get_tropomi_blnd.o \
main.o
%.o: %.mod
diff --git a/prep_satellite/mod_save.f90 b/prep_satellite/mod_save.f90
index f3cfaf5..6e9e2ac 100644
--- a/prep_satellite/mod_save.f90
+++ b/prep_satellite/mod_save.f90
@@ -67,11 +67,11 @@ module mod_save
!
! --------------------------------------------------
- subroutine write_releases(settings,nretr,nlayer,ncoord,xpoints,ypoints,zpoint1,zpoint2,&
+ subroutine write_releases(settings,numbox,nlayer,ncoord,xpoints,ypoints,zpoint1,zpoint2,&
idate,itime,cak_out,vdryair_out,specnum,units)
type (settings_t), intent(in) :: settings
- integer, intent(in) :: nretr, nlayer, ncoord
+ integer, intent(in) :: numbox, nlayer, ncoord
real, dimension(maxretr,ncoord), intent(in) :: xpoints, ypoints
real, dimension(maxretr,nlayer), intent(in) :: zpoint1, zpoint2
integer, dimension(maxretr), intent(in) :: idate, itime
@@ -95,7 +95,7 @@ module mod_save
call check( nf90_create(trim(pathout)//trim(filename), nf90_clobber, ncid) )
! define dimensions and variables
- call check( nf90_def_dim(ncid, 'retrieval', nretr, nrid) )
+ call check( nf90_def_dim(ncid, 'retrieval', numbox, nrid) )
call check( nf90_def_dim(ncid, 'nlayer', nlayer, nlid) )
call check( nf90_def_dim(ncid, 'ncoord', ncoord, ncoid) )
call check( nf90_def_var(ncid, 'xpoints', nf90_real, (/nrid, ncoid/), xptid) )
@@ -122,14 +122,14 @@ module mod_save
call check( nf90_enddef(ncid) )
! insert variables
- call check( nf90_put_var(ncid, xptid, xpoints(1:nretr,:)) )
- call check( nf90_put_var(ncid, yptid, ypoints(1:nretr,:)) )
- call check( nf90_put_var(ncid, zpt1id, zpoint1(1:nretr,:)) )
- call check( nf90_put_var(ncid, zpt2id, zpoint2(1:nretr,:)) )
- call check( nf90_put_var(ncid, dateid, idate(1:nretr)) )
- call check( nf90_put_var(ncid, timeid, itime(1:nretr)) )
- call check( nf90_put_var(ncid, cakid, cak_out(1:nretr,:)) )
- call check( nf90_put_var(ncid, dairid, vdryair_out(1:nretr,:)) )
+ call check( nf90_put_var(ncid, xptid, xpoints(1:numbox,:)) )
+ call check( nf90_put_var(ncid, yptid, ypoints(1:numbox,:)) )
+ call check( nf90_put_var(ncid, zpt1id, zpoint1(1:numbox,:)) )
+ call check( nf90_put_var(ncid, zpt2id, zpoint2(1:numbox,:)) )
+ call check( nf90_put_var(ncid, dateid, idate(1:numbox)) )
+ call check( nf90_put_var(ncid, timeid, itime(1:numbox)) )
+ call check( nf90_put_var(ncid, cakid, cak_out(1:numbox,:)) )
+ call check( nf90_put_var(ncid, dairid, vdryair_out(1:numbox,:)) )
! close file
call check( nf90_close(ncid) )
@@ -154,22 +154,22 @@ module mod_save
!
! --------------------------------------------------
- subroutine write_retrieval(settings,nretr,ncoord,xpoints,ypoints,idate,itime,&
- cpri,cakpri,vmr_out,verr_out,cdryair,units)
+ subroutine write_retrieval(settings,numbox,ncoord,xpoints,ypoints,idate,itime,&
+ cpri,cakpri,vmr_out,verr_out,cdryair,albedo_out,units)
type (settings_t), intent(in) :: settings
- integer, intent(in) :: nretr, ncoord
+ integer, intent(in) :: numbox, ncoord
real, dimension(maxretr,ncoord), intent(in) :: xpoints, ypoints
integer, dimension(maxretr), intent(in) :: idate, itime
- real, dimension(maxretr), intent(in) :: vmr_out, verr_out
- real, dimension(nretr), intent(in) :: cpri, cakpri, cdryair
+ real, dimension(maxretr), intent(in) :: vmr_out, verr_out, albedo_out
+ real, dimension(maxretr), intent(in) :: cpri, cakpri, cdryair
character(len=6), dimension(4), intent(in) :: units
character(max_name_len) :: filename
character(len=8) :: adate
integer :: ncid, nrid, nlid, ncoid
integer :: xptid, yptid, dateid, timeid
- integer :: cprid, cakprid, vmrid, verrid, cdairid
+ integer :: cprid, cakprid, vmrid, verrid, albid, cdairid
write(adate,fmt='(I8.8)') idate(1)
filename = 'retrieval_'//trim(settings%satname)//'_'//adate//'.nc'
@@ -179,7 +179,7 @@ module mod_save
call check( nf90_create(trim(settings%path_obsout)//trim(filename), nf90_clobber, ncid) )
! define dimensions and variables
- call check( nf90_def_dim(ncid, 'retrieval', nretr, nrid) )
+ call check( nf90_def_dim(ncid, 'retrieval', numbox, nrid) )
call check( nf90_def_dim(ncid, 'ncoord', ncoord, ncoid) )
call check( nf90_def_var(ncid, 'xpoints', nf90_real, (/nrid, ncoid/), xptid) )
call check( nf90_put_att(ncid, xptid, 'units', 'degrees_east') )
@@ -197,6 +197,8 @@ module mod_save
call check( nf90_put_att(ncid, vmrid, 'units', trim(units(4))) )
call check( nf90_def_var(ncid, 'vmr_err', nf90_real, nrid, verrid) )
call check( nf90_put_att(ncid, verrid, 'units', trim(units(4))) )
+ call check( nf90_def_var(ncid, 'albedo_swir', nf90_real, nrid, albid) )
+ call check( nf90_put_att(ncid, albid, 'units', 'none') )
call check( nf90_def_var(ncid, 'cdryair', nf90_real, nrid, cdairid) )
call check( nf90_put_att(ncid, cdairid, 'units', trim(units(2))) )
@@ -204,15 +206,16 @@ module mod_save
call check( nf90_enddef(ncid) )
! insert variables
- call check( nf90_put_var(ncid, xptid, xpoints(1:nretr,:)) )
- call check( nf90_put_var(ncid, yptid, ypoints(1:nretr,:)) )
- call check( nf90_put_var(ncid, dateid, idate(1:nretr)) )
- call check( nf90_put_var(ncid, timeid, itime(1:nretr)) )
- call check( nf90_put_var(ncid, cprid, cpri) )
- call check( nf90_put_var(ncid, cakprid, cakpri) )
- call check( nf90_put_var(ncid, vmrid, vmr_out(1:nretr)) )
- call check( nf90_put_var(ncid, verrid, verr_out(1:nretr)) )
- call check( nf90_put_var(ncid, cdairid, cdryair) )
+ call check( nf90_put_var(ncid, xptid, xpoints(1:numbox,:)) )
+ call check( nf90_put_var(ncid, yptid, ypoints(1:numbox,:)) )
+ call check( nf90_put_var(ncid, dateid, idate(1:numbox)) )
+ call check( nf90_put_var(ncid, timeid, itime(1:numbox)) )
+ call check( nf90_put_var(ncid, cprid, cpri(1:numbox)) )
+ call check( nf90_put_var(ncid, cakprid, cakpri(1:numbox)) )
+ call check( nf90_put_var(ncid, vmrid, vmr_out(1:numbox)) )
+ call check( nf90_put_var(ncid, verrid, verr_out(1:numbox)) )
+ call check( nf90_put_var(ncid, albid, albedo_out(1:numbox)) )
+ call check( nf90_put_var(ncid, cdairid, cdryair(1:numbox)) )
! close file
call check( nf90_close(ncid) )
diff --git a/prep_satellite/mod_settings.f90 b/prep_satellite/mod_settings.f90
index e03234e..96ed844 100644
--- a/prep_satellite/mod_settings.f90
+++ b/prep_satellite/mod_settings.f90
@@ -36,6 +36,7 @@ module mod_settings
type :: settings_t
+ integer :: lscratch
character(len=max_name_len) :: satellite
character(len=max_name_len) :: satname
character(len=max_path_len) :: path_flexpart
@@ -49,6 +50,7 @@ module mod_settings
character(len=max_path_len) :: file_availnest
character(len=max_name_len) :: windfield
+ logical :: albedo_corr
logical :: landonly
integer :: opmode_sel
character(len=max_name_len) :: nmeas_name
@@ -71,6 +73,7 @@ module mod_settings
character(len=max_name_len) :: pdel_name
character(len=max_name_len) :: surfclass_name
character(len=max_name_len) :: landfrac_name
+ character(len=max_name_len) :: albedo_name
integer :: datei
integer :: datef
@@ -310,7 +313,9 @@ module mod_settings
settings%avg_pixels = .false.
settings%usegridfile = .false.
settings%landonly = .false.
+ settings%albedo_corr = .false.
settings%opmode_sel = -1
+ settings%lscratch = 1 ! use scratch drive
! open file
open (100, file = trim (filename), status = 'old', iostat=ierr)
@@ -336,6 +341,9 @@ module mod_settings
if ( len (trim (line)) .eq. 0 ) cycle read_loop
! Path and file settings
+ identifier = "lscratch:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) settings%lscratch = int(cn)
identifier = "satellite:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%satellite = cc
@@ -374,6 +382,9 @@ module mod_settings
if ( match ) settings%windfield = cc
! satellite variables
+ identifier = "albedo_corr:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) settings%albedo_corr = cl
identifier = "landonly:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%landonly = cl
@@ -440,6 +451,9 @@ module mod_settings
identifier = "landfrac_name:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%landfrac_name = cc
+ identifier = "albedo_name:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) settings%albedo_name = cc
! COMMAND settings
identifier = "datei:"
--
GitLab