From f36a8890df53b9183821dc03e0128b780e856e25 Mon Sep 17 00:00:00 2001
From: ronesy <rlt@nilu.no>
Date: Mon, 8 Jul 2024 16:11:37 +0200
Subject: [PATCH] Bug fix for prior profile reading from tropomi data and
completion of the routine to read OCO2 data
---
prep_satellite/get_bremen.f90 | 37 +++++++++++++++-------
prep_satellite/get_oco2.f90 | 22 ++++++-------
prep_satellite/get_tropomi.f90 | 56 +++++++++++++++++++++++++++-------
3 files changed, 82 insertions(+), 33 deletions(-)
diff --git a/prep_satellite/get_bremen.f90 b/prep_satellite/get_bremen.f90
index e428896..a9c8e23 100644
--- a/prep_satellite/get_bremen.f90
+++ b/prep_satellite/get_bremen.f90
@@ -235,27 +235,27 @@ subroutine get_bremen(settings)
! select retrievals and extract needed info
if( n.eq.1 ) then
- allocate( xpoints(maxretr,ncoord) )
+ allocate( xpoints(maxretr,ncoord), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate xpoints'
- allocate( ypoints(maxretr,ncoord) )
+ allocate( ypoints(maxretr,ncoord), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate ypoints'
- allocate( zpoint1(maxretr,nlayer) )
+ allocate( zpoint1(maxretr,nlayer), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate zpoint1'
- allocate( zpoint2(maxretr,nlayer) )
+ allocate( zpoint2(maxretr,nlayer), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate zpoint2'
- allocate( idate(maxretr) )
+ allocate( idate(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate idate'
- allocate( itime(maxretr) )
+ allocate( itime(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate itime'
- allocate( cak_out(maxretr,nlayer) )
+ allocate( cak_out(maxretr,nlayer), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate cak_out'
- allocate( vapri_out(maxretr,nlayer) )
+ allocate( vapri_out(maxretr,nlayer), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vapri_out'
- allocate( vdryair_out(maxretr,nlayer) )
+ allocate( vdryair_out(maxretr,nlayer), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vdryair_out'
- allocate( vmr_out(maxretr) )
+ allocate( vmr_out(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vmr_out'
- allocate( verr_out(maxretr) )
+ allocate( verr_out(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate verr_out'
nretr = 0
endif
@@ -352,6 +352,21 @@ subroutine get_bremen(settings)
vapri_out, cpri, cakpri, cak_out, cdryair, vdryair_out, nretr, ncoord, nlayer)
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
! -----------------
diff --git a/prep_satellite/get_oco2.f90 b/prep_satellite/get_oco2.f90
index e206ef5..9f3c00c 100644
--- a/prep_satellite/get_oco2.f90
+++ b/prep_satellite/get_oco2.f90
@@ -262,27 +262,27 @@ subroutine get_oco2(settings)
! select retrievals and extract needed info
if( n.eq.1 ) then
- allocate( xpoints(maxretr,ncoord) )
+ allocate( xpoints(maxretr,ncoord), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate xpoints'
- allocate( ypoints(maxretr,ncoord) )
+ allocate( ypoints(maxretr,ncoord), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate ypoints'
- allocate( zpoint1(maxretr,nlevel) )
+ allocate( zpoint1(maxretr,nlevel), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate zpoint1'
- allocate( zpoint2(maxretr,nlevel) )
+ allocate( zpoint2(maxretr,nlevel), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate zpoint2'
- allocate( idate(maxretr) )
+ allocate( idate(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate idate'
- allocate( itime(maxretr) )
+ allocate( itime(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate itime'
- allocate( cak_out(maxretr,nlevel) )
+ allocate( cak_out(maxretr,nlevel), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate cak_out'
- allocate( vapri_out(maxretr,nlevel) )
+ allocate( vapri_out(maxretr,nlevel), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vapri_out'
- allocate( vdryair_out(maxretr,nlevel) )
+ allocate( vdryair_out(maxretr,nlevel), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vdryair_out'
- allocate( vmr_out(maxretr) )
+ allocate( vmr_out(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vmr_out'
- allocate( verr_out(maxretr) )
+ allocate( verr_out(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate verr_out'
nretr = 0
diff --git a/prep_satellite/get_tropomi.f90 b/prep_satellite/get_tropomi.f90
index bf8477d..0462e7a 100644
--- a/prep_satellite/get_tropomi.f90
+++ b/prep_satellite/get_tropomi.f90
@@ -69,6 +69,7 @@ subroutine get_tropomi(settings)
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 :: cpri, cakpri, cdryair
integer, dimension(:), allocatable :: idate, itime, time
@@ -350,27 +351,27 @@ subroutine get_tropomi(settings)
! select retrievals and extract needed info
if( n.eq.1 ) then
- allocate( xpoints(maxretr,ncoord) )
+ allocate( xpoints(maxretr,ncoord), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate xpoints'
- allocate( ypoints(maxretr,ncoord) )
+ allocate( ypoints(maxretr,ncoord), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate ypoints'
- allocate( zpoint1(maxretr,nlayer) )
+ allocate( zpoint1(maxretr,nlayer), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate zpoint1'
- allocate( zpoint2(maxretr,nlayer) )
+ allocate( zpoint2(maxretr,nlayer), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate zpoint2'
- allocate( idate(maxretr) )
+ allocate( idate(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate idate'
- allocate( itime(maxretr) )
+ allocate( itime(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate itime'
- allocate( cak_out(maxretr,nlayer) )
+ allocate( cak_out(maxretr,nlayer), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate cak_out'
- allocate( vapri_out(maxretr,nlayer) )
+ allocate( vapri_out(maxretr,nlayer), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vapri_out'
- allocate( vdryair_out(maxretr,nlayer) )
+ allocate( vdryair_out(maxretr,nlayer), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vdryair_out'
- allocate( vmr_out(maxretr) )
+ allocate( vmr_out(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate vmr_out'
- allocate( verr_out(maxretr) )
+ allocate( verr_out(maxretr), stat=ierr )
if ( ierr.ne.0 ) write(logid,*) 'ERROR: cannot allocate verr_out'
cak_out(:,:)=0.
vapri_out(:,:)=0.
@@ -452,6 +453,26 @@ subroutine get_tropomi(settings)
end do ! nread
+ ! 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)
+
print*, 'itime(1), itime(nretr) = ',itime(1),itime(nretr)
! calculate column prior (units mol/m2)
@@ -476,6 +497,19 @@ subroutine get_tropomi(settings)
vapri_out, cpri, cakpri, cak_out, cdryair, vdryair_out, nretr, ncoord, nlayer)
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
! -----------------
--
GitLab