From c1c9fd9e6eb55d1d614f30b599b49d7b08fdd4ae Mon Sep 17 00:00:00 2001
From: ronesy <rlt@nilu.no>
Date: Mon, 24 Apr 2023 17:35:58 +0200
Subject: [PATCH] Updated calculation of initial mixing ratios and harmonized
reading of initial mixing ratio fields
---
source/get_initconc.f90 | 409 +++++++++++++++++++++++++++++++------
source/init_cini.f90 | 186 ++++++++---------
source/init_cini_month.f90 | 346 +++++++++++++++++++++++++++++++
source/init_co2.f90 | 8 +-
source/init_ghg.f90 | 26 +--
source/initialize.f90 | 17 +-
source/interpz.f90 | 108 +++++++---
source/main.f90 | 3 +
source/makefile | 8 +-
source/mod_perturb.f90 | 1 -
source/mod_save.f90 | 66 +++---
source/mod_settings.f90 | 50 ++++-
source/mod_var.f90 | 12 +-
source/read_orog.f90 | 109 ++++++++++
source/simulate.f90 | 10 +-
15 files changed, 1071 insertions(+), 288 deletions(-)
create mode 100644 source/init_cini_month.f90
create mode 100644 source/read_orog.f90
diff --git a/source/get_initconc.f90 b/source/get_initconc.f90
index 9bc7ab1..d06ac3e 100644
--- a/source/get_initconc.f90
+++ b/source/get_initconc.f90
@@ -25,8 +25,10 @@
!! Interface:
!!
!! Inputs
+!! files - file data structure
+!! config - configuration settings data structure
!! filename - name of 4D concentration file
-!! varname - variable name
+!! ntime - time dimension of concentration field
!! concini - concentrations on flexpart grid (zero array on input)
!!
!! Externals
@@ -34,34 +36,46 @@
!! interpz
!! check
!!
+!! To do:
+!!
+!! Add temperature to EGG4 and CAMS-N2O files for more accurate geopotential
+!! heights
+!!
!---------------------------------------------------------------------------------------
-subroutine get_initconc(filename, varname, concini)
+subroutine get_initconc(files, config, filename, ntime, concini)
use mod_var
+ use mod_settings
+ use mod_tools, only : sort
use mod_save, only : check
use netcdf
implicit none
- character(len=max_path_len), intent(in) :: filename
- character(len=max_name_len), intent(in) :: varname
- real, dimension(nxgrid,nygrid,nzgrid), intent(in out) :: concini
-
- character(len=*), parameter :: lonname='longitude'
- character(len=*), parameter :: latname='latitude'
- character(len=*), parameter :: prsname='pressure'
- character(len=*), parameter :: timename='time'
- integer :: ncid, xid, yid, zid, tid, varid
- integer :: kz
- integer :: nx, ny, nz, nt
- integer :: nx_out, ny_out
- logical :: midpoints
- real, dimension(:), allocatable :: lon, lat, prs, geoh
- real, dimension(nxgrid) :: lon_out
- real, dimension(nygrid) :: lat_out
- real, dimension(:,:,:), allocatable :: work, workout, concout
- real, dimension(:,:,:,:), allocatable :: conc
+ type (files_t), intent(in) :: files
+ type (config_t), intent(in) :: config
+ character(len=max_path_len), intent(in) :: filename
+ integer, intent(in) :: ntime
+ real, dimension(nxgrid,nygrid,nzgrid,ntime), intent(in out) :: concini
+
+ character(len=max_name_len) :: lonname, latname, varname, prsname, timename
+ logical :: midpoints
+ integer :: ncid, xid, yid, zid, tid, varid
+ integer :: n, kz, status, iproduct, zkind
+ integer :: nx, ny, nz, nt
+ integer :: nx_out, ny_out
+ real, dimension(:), allocatable :: lon, lat
+ real(kind=8), dimension(:), allocatable :: lono, lato, alto
+ real, dimension(:), allocatable :: ap, bp, mpres
+ integer, dimension(:), allocatable :: indx, indy, indz
+ real, dimension(nxgrid) :: lon_out
+ real, dimension(nygrid) :: lat_out
+ real, dimension(:,:), allocatable :: rpres, scalh
+ real, dimension(:,:,:), allocatable :: work1, work2, work3, work4
+ real, dimension(:,:,:), allocatable :: geoh, geoho, pres, surfp, pw
+ real, dimension(:,:,:,:), allocatable :: conc, conco, humid, temp
+ real :: sclfact, offset
! initialize
! ----------
@@ -72,6 +86,56 @@ subroutine get_initconc(filename, varname, concini)
nx_out = nxgrid
ny_out = nygrid
+ ! which initial concentration field product
+ ! 1 = CAMS CH4, 2 = EGG4, 3 = CAMS N2O, 4 = FP-CTM, 5 = NOAA interpolated, 6 = FP-CTM (monthly)
+ iproduct = int(files%init_type)
+ varname = files%varname_init
+
+ select case ( iproduct)
+ case (1)
+ ! CAMS CH4
+ lonname = 'longitude'
+ latname = 'latitude'
+ timename = 'time'
+ prsname = 'level'
+ case (2)
+ ! EGG4
+ lonname = 'longitude'
+ latname = 'latitude'
+ timename = 'time'
+ prsname = 'level'
+ case (3)
+ ! CAMS N2O
+ lonname = 'longitude'
+ latname = 'latitude'
+ timename = 'time'
+ prsname = 'level'
+ case (4)
+ ! FP-CTM
+ lonname = 'longitude'
+ latname = 'latitude'
+ timename = 'Date'
+ prsname = 'lev'
+ case (5)
+ ! NOAA
+ lonname = 'longitude'
+ latname = 'latitude'
+ timename = 'time'
+ prsname = 'pressure'
+ case (6)
+ ! FP-CTM (monthly)
+ lonname = 'longitude'
+ latname = 'latitude'
+ timename = 'none'
+ prsname = 'lev'
+ end select
+
+ ! from here on FP-CTM daily and monthly treated the same
+ if ( iproduct.eq.6 ) iproduct = 4
+
+ print*, 'get_initconc: lonname, latname, prsname =',lonname,latname,prsname
+
+
! read netcdf file
! ----------------
@@ -79,70 +143,285 @@ subroutine get_initconc(filename, varname, concini)
! read dimensions
call check( nf90_open(trim(filename),nf90_NOWRITE,ncid) )
- call check( nf90_inq_dimid(ncid,lonname,xid) )
- call check( nf90_inq_dimid(ncid,latname,yid) )
- call check( nf90_inq_dimid(ncid,prsname,zid) )
- call check( nf90_inq_dimid(ncid,timename,tid) )
+ call check( nf90_inq_dimid(ncid,trim(lonname),xid) )
call check( nf90_inquire_dimension(ncid,xid,len=nx) )
+ call check( nf90_inq_dimid(ncid,trim(latname),yid) )
call check( nf90_inquire_dimension(ncid,yid,len=ny) )
+ call check( nf90_inq_dimid(ncid,trim(prsname),zid) )
call check( nf90_inquire_dimension(ncid,zid,len=nz) )
- call check( nf90_inquire_dimension(ncid,tid,len=nt) )
-
- ! mean so time dimension should be 1
- if( nt.ne.1 ) then
- write(logid,*) 'ERROR get_initconc: time dimension should have length of 1 in '//trim(filename)
- stop
+ if ( trim(timename).ne.'none' ) then
+ call check( nf90_inq_dimid(ncid,trim(timename),tid) )
+ call check( nf90_inquire_dimension(ncid,tid,len=nt) )
+ else
+ nt = 1
endif
+ if ( ntime.ne.nt ) write(logid,*) 'ERROR get_initconc: inconsistent time dimension'
+ ! read dimension variables
allocate( lon(nx) )
allocate( lat(ny) )
- allocate( prs(nz) )
+ call check( nf90_inq_varid(ncid,trim(lonname),varid) )
+ call check( nf90_get_var(ncid,varid,lon) )
+ call check( nf90_inq_varid(ncid,trim(latname),varid) )
+ call check( nf90_get_var(ncid,varid,lat) )
+
+ ! read supporting variables
+ select case (iproduct)
+ case (1)
+ ! CAMS CH4
+ ! hybrid pressure coordinates
+ allocate( ap(nz) )
+ allocate( bp(nz) )
+ call check( nf90_inq_varid(ncid,'hyam',varid) )
+ call check( nf90_get_var(ncid,varid,ap) )
+ call check( nf90_inq_varid(ncid,'hybm',varid) )
+ call check( nf90_get_var(ncid,varid,bp) )
+ ! read temperature
+ allocate( temp(nx,ny,nz,nt) )
+ call check( nf90_inq_varid(ncid,'T',varid) )
+ call check( nf90_get_var(ncid,varid,temp) )
+ call check( nf90_get_att(ncid, varid, 'scale_factor', sclfact) )
+ call check( nf90_get_att(ncid, varid, 'add_offset', offset) )
+ temp = temp*sclfact + offset
+ ! read surface pressure
+ allocate( surfp(nx,ny,nt) )
+ call check( nf90_inq_varid(ncid,'ps',varid) )
+ call check( nf90_get_var(ncid,varid,surfp) )
+ call check( nf90_get_att(ncid, varid, 'scale_factor', sclfact) )
+ call check( nf90_get_att(ncid, varid, 'add_offset', offset) )
+ surfp = surfp*sclfact + offset
+ ! mixing ratio is in humid air -> read specific humidity
+ allocate( humid(nx,ny,nz,nt) )
+ call check( nf90_inq_varid(ncid,'Q',varid) )
+ call check( nf90_get_var(ncid,varid,humid) )
+ call check( nf90_get_att(ncid, varid, 'scale_factor', sclfact) )
+ call check( nf90_get_att(ncid, varid, 'add_offset', offset) )
+ humid = humid*sclfact + offset
+ case (2)
+ ! EGG4
+ ! vertical coordinate is pressure
+ allocate( mpres(nz) )
+ call check( nf90_inq_varid(ncid,trim(prsname),varid) )
+ call check( nf90_get_var(ncid,varid,mpres) )
+ ! convert to unit of Pa
+ mpres = mpres*100.
+ case (3)
+ ! CAMS N2O
+ ! vertical coordinate is pressure
+ allocate( mpres(nz) )
+ call check( nf90_inq_varid(ncid,trim(prsname),varid) )
+ call check( nf90_get_var(ncid,varid,mpres) )
+ case (4)
+ ! FP-CTM
+ allocate( mpres(nz) )
+ call check( nf90_inq_varid(ncid,trim(prsname),varid) )
+ call check( nf90_get_var(ncid,varid,mpres) )
+ case (5)
+ ! NOAA
+ allocate( mpres(nz) )
+ call check( nf90_inq_varid(ncid,trim(prsname),varid) )
+ call check( nf90_get_var(ncid,varid,mpres) )
+ end select
+
+ ! read mixing ratio
allocate( conc(nx,ny,nz,nt) )
- allocate( work(nx,ny,1) )
- allocate( workout(nxgrid,nygrid,1) )
- allocate( concout(nxgrid,nygrid,nz) )
- allocate( geoh(nz) )
-
- ! read variables
- call check( nf90_inq_varid(ncid,lonname,xid) )
- call check( nf90_get_var(ncid,xid,lon) )
- call check( nf90_inq_varid(ncid,latname,yid) )
- call check( nf90_get_var(ncid,yid,lat) )
- call check( nf90_inq_varid(ncid,prsname,zid) )
- call check( nf90_get_var(ncid,zid,prs) )
call check( nf90_inq_varid(ncid,trim(varname),varid) )
call check( nf90_get_var(ncid,varid,conc(:,:,:,:)) )
-
+ ! check if file contains scalar and offset
+ status = nf90_inquire_attribute(ncid,varid,'scale_factor')
+ if ( status.eq.nf90_noerr ) then
+ call check( nf90_get_att(ncid, varid, 'scale_factor', sclfact) )
+ call check( nf90_get_att(ncid, varid, 'add_offset', offset) )
+ conc = conc*sclfact + offset
+ endif
call check( nf90_close(ncid) )
- ! calculate geopotential height for each layer
- do kz = 1, nz
- geoh(kz) = (gasc*ts/gc/mmol)*log(psurf/prs(kz))
- enddo
+ ! adjust mixing ratios
+ ! --------------------
- ! interpolate to flexpart grid
- ! ----------------------------
-
- ! interpolate horizontally
- midpoints = .false.
- do kz = 1, nz
- work(:,:,1) = conc(:,:,kz,1)
- ! gbl_lon and gbl_lat are lon and lat of global flexpart domain
- call convertgrid(midpoints, nx, ny, lon, lat, work, nx_out, ny_out, lon_out, lat_out, workout)
- concout(:,:,kz) = workout(:,:,1)
- end do
+ select case ( iproduct )
+ case (1)
+ ! CAMS CH4
+ ! mole fraction of humid air -> convert to dry air
+ ! calculate mid layer pressure (mean over nt steps)
+ allocate( pres(nx,ny,nz) )
+ allocate( pw(nx,ny,nz) )
+ do n = 1, nz
+ pres(:,:,n) = ap(n) + bp(n)*sum(surfp(:,:,:),dim=3)/real(nt)
+ end do
+ ! partial pressure water from specific humidity
+ pw = sum(humid,dim=4)*pres/real(nt)/(0.622 + 0.378*(sum(humid,dim=4)/real(nt)))
+ do n = 1, nt
+ conc(:,:,:,n) = conc(:,:,:,n)*pres/(pres - pw)
+ end do
+ print*, 'pres = ',pres(nx/2,ny/2,:)
+ print*, 'range corr dry air = ',minval(pres/(pres-pw)),maxval(pres/(pres-pw))
+ deallocate( humid )
+ deallocate( pw )
+ case (2)
+ ! EGG4
+ ! adjust units from mass mixing ratio to mole fraction
+ conc = conc*1.e9*mmair/config%molarmass
+ end select
+
+ ! calculate geopotential height
+ ! -----------------------------
+
+ allocate( geoh(nx,ny,nz) )
+
+ select case (iproduct)
+ case (1)
+ ! CAMS CH4
+ allocate( rpres(nx,ny) )
+ allocate( scalh(nx,ny) )
+ do n = 1, nz
+ rpres = sum(surfp(:,:,:),dim=3)/real(nt)/pres(:,:,n)
+ scalh = gasc*sum(temp(:,:,n,:),dim=3)/real(nt)/mmol/gc
+ geoh(:,:,n) = scalh*log(rpres)
+ print*, 'range(scalh) = ',minval(scalh), maxval(scalh)
+ end do
+ deallocate( rpres )
+ deallocate( scalh )
+ deallocate( surfp )
+ zkind = 0
+ case (2)
+ ! EGG4
+ do n = 1, nz
+ ! at present temp not included so use estimate
+ geoh(:,:,n) = (gasc*ts/gc/mmol)*log(psurf/mpres(n))
+ end do
+ zkind = 0
+ case (3)
+ ! CAMS N2O
+ do n = 1, nz
+ ! at present temp not included so use estimate
+ geoh(:,:,n) = (gasc*ts/gc/mmol)*log(psurf/mpres(n))
+ end do
+ zkind = 0
+ case (4)
+ ! FP-CTM
+ ! vertical coordinate is meters above ground
+ do n = 1, nz
+ geoh(:,:,n) = mpres(n)
+ end do
+ zkind = 1
+ case (5)
+ ! NOAA
+ do n = 1, nz
+ ! at present temp not included so use estimate
+ geoh(:,:,n) = (gasc*ts/gc/mmol)*log(psurf/mpres(n))
+ end do
+ zkind = 0
+ end select
+
+ if ( allocated(mpres) ) deallocate( mpres )
+ if ( allocated(pres) ) deallocate( pres )
+ if ( allocated(temp) ) deallocate( temp )
+
+ ! adjust lat, lon and alt of grid
+ ! -------------------------------
- ! interpolate vertically
- call interpz(nz, geoh, concout, concini)
+ select case ( iproduct )
+ case (1)
+ ! CAMS CH4
+ midpoints = .true.
+ case (2)
+ ! EGG4
+ ! adjust longitude to be -180 to 180
+ where( lon.gt.180. ) lon(:) = lon(:) - 360.
+ ! longitude given as eastern boundary
+ lon(:) = lon(:) - 0.5*(lon(2) - lon(1))
+ ! file contains extra latitude
+ ! latitude given from 90 to -90
+ do n = 1, ny
+ lat(n) = 90. - real(n-1)*(180./real(ny-1)) - 0.5*180./real(ny-1)
+ end do
+ ny = ny - 1
+ midpoints = .true.
+ case (3)
+ ! CAMS N2O
+ ! file contains extra longitude
+ nx = nx - 1
+ ! latitude given from 90 to -90
+ do n = 1, ny
+ lat(n) = 90. - real(n)*(180./real(ny))
+ end do
+ midpoints = .false.
+ case (4)
+ ! FP-CTM
+ ! file contains extra longitude
+ nx = nx - 1
+ print*, 'get_initconc: nx = ',nx
+ lon(:) = lon(:) + 0.5*(lon(2) - lon(1))
+ midpoints = .true.
+ case (5)
+ ! NOAA
+ midpoints = .false.
+ end select
+ ! sort from south to north and west to east
+ print*, 'get_initconc: lon = ',lon
+ print*, 'get_iniconc: lat = ',lat
+ allocate( indy(ny) )
+ allocate( indx(nx) )
+ allocate( lato(ny) )
+ allocate( lono(nx) )
+ allocate( conco(nx,ny,nz,nt) )
+ allocate( geoho(nx,ny,nz) )
+ lato = dble(lat(1:ny))
+ lono = dble(lon(1:nx))
+ call sort(ny, lato, indy)
+ call sort(nx, lono, indx)
+ ! sort from lowest to highest level
+ allocate( indz(nz) )
+ allocate( alto(nz) )
+ alto = dble(geoh(nx/2,ny/2,:))
+ call sort(nz, alto, indz)
+ print*, 'get_initconc: geoh = ',geoh(nx/2,ny/2,:)
+ print*, 'get_initconc: indx = ',indx
+ print*, 'get_initconc: indy = ',indy
+ print*, 'get_initconc: lono = ',lono
+ print*, 'get_initconc: lato = ',lato
+ print*, 'get_initconc: indz = ',indz
+ conco = conc(indx,indy,indz,:)
+ geoho = geoh(indx,indy,indz)
+ deallocate( indx )
+ deallocate( indy )
deallocate( lon )
deallocate( lat )
- deallocate( prs )
+
+ ! interpolate to flexpart grid
+ ! ----------------------------
+
+ allocate( work1(nx,ny,nz) )
+ allocate( work2(nx,ny,nzgrid) )
+ allocate( work3(nx,ny,1) )
+ allocate( work4(nxgrid,nygrid,1) )
+
+ do n = 1, nt
+ work1(:,:,:) = conco(:,:,:,n)
+ ! interpolate vertically
+ call interpz(nx, ny, nz, zkind, midpoints, sngl(lono), sngl(lato), geoho, work1, work2)
+ ! interpolate horizontally
+ do kz = 1, nzgrid
+ work3(:,:,1) = work2(:,:,kz)
+ ! lon_out and lat_out are for global flexpart domain
+ call convertgrid(midpoints, nx, ny, sngl(lono), sngl(lato), work3, nx_out, ny_out, lon_out, lat_out, work4)
+ concini(:,:,kz,n) = work4(:,:,1)
+ end do
+ end do
+
+ deallocate( lono )
+ deallocate( lato )
+ deallocate( alto )
deallocate( conc )
- deallocate( work )
- deallocate( workout )
- deallocate( concout )
deallocate( geoh )
+ deallocate( conco )
+ deallocate( geoho )
+ deallocate( work1 )
+ deallocate( work2 )
+ deallocate( work3 )
+ deallocate( work4 )
end subroutine get_initconc
diff --git a/source/init_cini.f90 b/source/init_cini.f90
index 93088e6..19a493d 100644
--- a/source/init_cini.f90
+++ b/source/init_cini.f90
@@ -52,37 +52,55 @@ subroutine init_cini(files, config, obs)
type (config_t), intent (in) :: config
type (obs_t), intent (in out) :: obs
- character(max_path_len) :: path_flexrec, filename
- character(max_name_len) :: varname
+ character(max_path_len) :: path_flexrec, pathname, filename
character(len=8) :: adate, areldate, areltime
+ character(len=6) :: ayrmon
+ character(len=4) :: ayear
+ character(len=2) :: amonth
logical :: lexist
- integer :: lastjjjjmm, prejjjjmm, projjjjmm
- integer :: jjjjmm, jjjjmmdd, hhmiss, mm, eomday
- integer :: ix, jy, kz, i, n
+
+ integer :: hfreq
+ integer :: jjjjmm, jjjjmmdd, hhmiss, eomday
+ integer :: lasttime, startmonth
integer :: month, prevmonth
- real :: conc
- real(kind=8) :: jdmid, jdpre, jdlast
+ integer :: ix, jy, kz, i, n, nt, nz, nhours, ntime, ierr
+ real(kind=8) :: jdi
real(kind=8), dimension(nobs) :: jdates
real, dimension(nobs,ncini) :: cini
integer, dimension(nobs) :: ind
real, dimension(nxgrid,nygrid,nzgrid) :: gridini, gridini_tmp
- real, dimension(nxgrid,nygrid,nzgrid) :: concini, preconcini, proconcini
+ real, dimension(:,:,:,:), allocatable :: concini
real :: bndx, bndy, delx, dely
integer :: numx, numy, xshift
integer :: numpoint, release_nr, nr_init
integer :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2) :: releases
+ integer :: nretr
+ character(len=recname_len) :: anretr
+
+
+ ! initialize
+ ! ----------
+ ! which initial concentration field product
+ ! 1 = CAMS CH4, 2 = EGG4, 3 = CAMS N2O, 4 = FP-CTM (daily)
+ select case ( files%init_type )
+ case (1)
+ hfreq = 24
+ case (2)
+ hfreq = 6
+ case (3)
+ hfreq = 3
+ case (4)
+ hfreq = 24
+ end select
! loop over observations
! ----------------------
- lastjjjjmm = 0
- concini(:,:,:) = 0.
- preconcini(:,:,:) = 0.
- proconcini(:,:,:) = 0.
cini(:,:) = 0.
+ lasttime = 0
prevmonth = 0
! sort obs chronologically for efficiency
@@ -95,9 +113,12 @@ subroutine init_cini(files, config, obs)
gridini(:,:,:) = 0.
nr_init = 0
+ print*, 'obs%recs(i) = ',obs%recs(ind(i))
+
! check month of observation
call caldate(jdates(i), jjjjmmdd, hhmiss)
write(adate,'(I6.6)') jjjjmmdd/100
+ write(areldate,'(I8.8)') jjjjmmdd
month = jjjjmmdd/100
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
@@ -131,7 +152,6 @@ subroutine init_cini(files, config, obs)
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
-! print*, 'init_cini: filename = ',filename
! read first flexpart grid_init files
inquire ( file=trim(filename), exist=lexist )
if ( lexist ) then
@@ -148,7 +168,6 @@ subroutine init_cini(files, config, obs)
write(areltime,'(I6.6)') hhmiss
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
-! print*, 'init_cini: filename next = ',filename
! sum the grid with previous initial grids
call read_init(filename, gridini_tmp)
gridini = gridini + gridini_tmp
@@ -187,114 +206,73 @@ subroutine init_cini(files, config, obs)
endif
! calculate termination time of trajectory
- call caldate(jdates(i) - trajdays, jjjjmmdd, hhmiss)
+ jdi = jdates(i) - trajdays
+ call caldate(jdi, jjjjmmdd, hhmiss)
jjjjmm = jjjjmmdd/100
-
- ! preceding and proceding months
- mm = jjjjmm - (jjjjmm/100)*100
- if( mm.gt.1 ) then
- prejjjjmm = jjjjmm - 1
- else
- prejjjjmm = (jjjjmm/100 - 1)*100 + 12
- endif
- if( mm.lt.12 ) then
- projjjjmm = jjjjmm + 1
- else
- projjjjmm = (jjjjmm/100 + 1)*100 + 1
- endif
- jdmid = juldate(jjjjmm*100+14, 0)
- jdpre = juldate(prejjjjmm*100+14, 0)
- jdlast = juldate(lastjjjjmm*100+14, 0)
- eomday = calceomday(prejjjjmm)
-
- ! get new concentrations if not same month as previous observation
- if( jjjjmm.ne.lastjjjjmm ) then
- if( (lastjjjjmm.ne.0).and.(jdmid.eq.(jdlast+dble(eomday))) ) then
- ! this is proceeding month
- preconcini = concini
- concini = proconcini
- ! get proceeding concentration field
- write(adate,'(I6.6)') projjjjmm
- filename = str_replace(files%file_initconc, 'YYYYMM', adate)
- filename = trim(files%path_initconc)//trim(filename)
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- varname = files%varname_init
- call get_initconc(filename, varname, proconcini)
- else
- write(logid,*) 'WARNING: cannot find ',filename
- write(logid,*) 'WARNING: next background file missing, copy present into next to avoid 0 value'
- proconcini = concini
- endif
- else
- ! get current concentration field
- write(adate,'(I6.6)') jjjjmm
- filename = str_replace(files%file_initconc, 'YYYYMM', adate)
- filename = trim(files%path_initconc)//trim(filename)
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- varname = files%varname_init
- call get_initconc(filename, varname, concini)
- else
- write(logid,*) 'WARNING: cannot find ',filename
- go to 10
- endif
- ! get previous concentration field
- write(adate,'(I6.6)') prejjjjmm
- filename = str_replace(files%file_initconc, 'YYYYMM', adate)
- filename = trim(files%path_initconc)//trim(filename)
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- varname = files%varname_init
- call get_initconc(filename, varname, preconcini)
- else
- write(logid,*) 'WARNING: cannot find ',filename
- go to 10
- endif
- ! get next concentration field
- write(adate,'(I6.6)') projjjjmm
- filename = str_replace(files%file_initconc, 'YYYYMM', adate)
- filename = trim(files%path_initconc)//trim(filename)
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- varname = files%varname_init
- call get_initconc(filename, varname, proconcini)
- else
- write(logid,*) 'WARNING: cannot find ',filename
- write(logid,*) 'WARNING: next background file missing, copy present into next to avoid 0 value'
- proconcini = concini
- endif
- endif
+ ! read concentration field only if not already in memory
+ if ( jjjjmm.ne.lasttime ) then
+ ! allocate concini knowing data freq (from hfreq)
+ if (allocated( concini ) ) deallocate( concini )
+ eomday = calceomday(jjjjmm)
+ ntime = eomday*24/hfreq
+ allocate( concini(nxgrid,nygrid,nzgrid,ntime), stat = ierr )
+ if ( ierr.ne.0 ) then
+ write(logid,*) 'ERROR init_cini: not enough memory'
+ stop
+ endif
+ concini(:,:,:,:) = 0.
+ ! read new concentration field
+ write(ayear,'(I4.4)') jjjjmm/100
+ write(ayrmon,'(I6.6)') jjjjmm
+ write(amonth,'(I2.2)') jjjjmm-(jjjjmm/100)*100
+ if ( index(files%path_initconc, 'YYYY').ne.0 ) then
+ pathname = str_replace(files%path_initconc, 'YYYY', ayear)
+ else
+ pathname = files%path_initconc
+ endif
+ filename = str_replace(files%file_initconc, 'YYYY', ayear)
+ filename = str_replace(filename, 'MM', amonth)
+ filename = trim(pathname)//trim(filename)
+ print*, 'init_cini: filename = ',filename
+ inquire(file=trim(filename),exist=lexist)
+ if(lexist) then
+ call get_initconc(files, config, filename, ntime, concini)
+ else
+ write(logid,*) 'ERROR init_cini: cannot find file ',filename
+ stop
+ endif
endif
! calculate contribution of initial concentration to receptor
+ ! index time dimension of concini
+ startmonth = (jjjjmmdd/100)*100 + 1
+ nhours = int((jdi - juldate(startmonth, 0))*24d0)
+ nt = nhours/hfreq + 1
+ print*, 'init_cini: startmonth = ',startmonth
+ print*, 'init_cini: nhours = ',nhours
+ print*, 'init_cini: nt = ',nt
do kz = 1, nzgrid
+ nz = 0
+ do while(.true.)
+ nz = nz + 1
+ if ( outheight(kz).lt.cini_alt(nz) ) exit
+ end do
do jy = 1, nygrid
n = 0
do while(.true.)
n = n + 1
if ( gbl_lat(jy).lt.cini_lat(n) ) exit
end do
-! print*, 'init_cini: gbl_lat(jy), n = ',gbl_lat(jy), n
do ix = 1, nxgrid
- ! interpolate between concentrations
- if( (jdates(i) - trajdays).lt.jdmid ) then
- conc = preconcini(ix,jy,kz) + (concini(ix,jy,kz) - preconcini(ix,jy,kz)) * &
- real((jdates(i) - trajdays) - jdpre, kind=4)/30.
- else
- conc = concini(ix,jy,kz) + (proconcini(ix,jy,kz) - concini(ix,jy,kz)) * &
- real((jdates(i) - trajdays) - jdmid, kind=4)/30.
- endif
- cini(i,n) = cini(i,n) + gridini(ix,jy,kz)*conc
+ cini(i,(nz-1)*nycini+n) = cini(i,(nz-1)*nycini+n) + gridini(ix,jy,kz)*concini(ix,jy,kz,nt)
end do
end do
end do
10 continue
- lastjjjjmm = jjjjmm
- prevmonth = month
+ lasttime = jjjjmm
end do
@@ -303,6 +281,8 @@ subroutine init_cini(files, config, obs)
obs%cini(ind(i),:) = cini(i,:)
end do
+ if (allocated( concini )) deallocate( concini )
+
end subroutine init_cini
diff --git a/source/init_cini_month.f90 b/source/init_cini_month.f90
new file mode 100644
index 0000000..d709da6
--- /dev/null
+++ b/source/init_cini_month.f90
@@ -0,0 +1,346 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT: init_cini_month
+!---------------------------------------------------------------------------------------
+! 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 2017, Rona Thompson
+!---------------------------------------------------------------------------------------
+!
+!> init_cini_month
+!!
+!! Purpose: Calculates the initial concentration for each observation using the
+!! 4D initial concentration fields and the flexpart grid_initial files.
+!!
+!! Interface:
+!!
+!! Inputs
+!! files - file data structure
+!! config - configuration settings data structure
+!! obs - observation data structure
+!!
+!! Externals
+!! caldate
+!! read_init
+!! get_initconc
+!!
+!---------------------------------------------------------------------------------------
+
+subroutine init_cini_month(files, config, obs)
+
+ use mod_flexpart
+ use mod_settings
+ use mod_strings
+ use mod_var
+ use mod_dates
+ use mod_obs
+ use mod_tools
+
+ implicit none
+
+ type (files_t), intent (in) :: files
+ type (config_t), intent (in) :: config
+ type (obs_t), intent (in out) :: obs
+
+ character(max_path_len) :: path_flexrec, filename, pathname
+ character(max_name_len) :: varname
+ character(len=8) :: adate, areldate, areltime
+ character(len=6) :: ayear
+ logical :: lexist
+ integer :: lastjjjjmm, prejjjjmm, projjjjmm
+ integer :: jjjjmm, jjjjmmdd, hhmiss, mm, eomday
+ integer :: ix, jy, kz, i, n, nz, ntime
+ integer :: month, prevmonth
+ real :: conc
+ real(kind=8) :: jdmid, jdpre, jdlast
+ real(kind=8), dimension(nobs) :: jdates
+ real, dimension(nobs,ncini) :: cini
+ integer, dimension(nobs) :: ind
+ real, dimension(nxgrid,nygrid,nzgrid) :: gridini, gridini_tmp
+ real, dimension(nxgrid,nygrid,nzgrid,1):: concini, preconcini, proconcini
+
+ real :: bndx, bndy, delx, dely
+ integer :: numx, numy, xshift
+ integer :: numpoint, release_nr, nr_init
+ integer :: ibdate, ibtime
+ real(kind=8), dimension(maxpoint,2) :: releases
+ integer :: nretr, ierr
+ character(len=6) :: anretr
+
+ ! loop over observations
+ ! ----------------------
+
+ lastjjjjmm = 0
+ concini(:,:,:,:) = 0.
+ preconcini(:,:,:,:) = 0.
+ proconcini(:,:,:,:) = 0.
+ cini(:,:) = 0.
+ prevmonth = 0
+ ntime = 1
+
+ ! sort obs chronologically for efficiency
+ jdates = obs%jdate(:)
+ call sort(nobs,jdates,ind)
+
+ do i = 1, nobs
+
+ ! reset values
+ gridini(:,:,:) = 0.
+ nr_init = 0
+
+ ! check month of observation
+ call caldate(jdates(i), jjjjmmdd, hhmiss)
+ write(adate,'(I6.6)') jjjjmmdd/100
+ write(areldate,'(I8.8)') jjjjmmdd
+ month = jjjjmmdd/100
+ path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
+
+ if ( config%average_fp ) then
+
+ ! match and, if needed, average grid_initial files to observation timestamp
+
+ ! get header and release information related to current observation location
+ ! (if not already read and same for previous observation)
+ filename = trim(path_flexrec)//'header'
+ if ( i.eq.1 .or. month.ne.prevmonth ) then
+ write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
+ call read_header(filename, numpoint, ibdate, ibtime, releases, &
+ numx, numy, bndx, bndy, delx, dely, xshift)
+ else
+ if(trim(obs%recs(ind(i))).ne.trim(obs%recs(ind(i-1))))then
+ write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
+ call read_header(filename, numpoint, ibdate, ibtime, releases, &
+ numx, numy, bndx, bndy, delx, dely, xshift)
+ endif
+ endif
+ ! find the first release number corresponding to the current observation
+ release_nr=1
+ do while ( int(releases(release_nr,1)*1.e4,kind=8).lt.int(jdates(i)*1.e4,kind=8) .and. release_nr.lt.numpoint )
+ release_nr = release_nr+1
+ end do
+ ! define first filename
+ call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
+ ! round seconds
+ hhmiss = (hhmiss/100)*100
+ write(areldate,'(I8.8)') jjjjmmdd
+ write(areltime,'(I6.6)') hhmiss
+ filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
+ ! read first flexpart grid_init files
+ inquire ( file=trim(filename), exist=lexist )
+ if ( lexist ) then
+ call read_init(filename, gridini_tmp)
+ gridini = gridini + gridini_tmp
+ nr_init = 1
+ ! current release_nr was already read, start with next
+ release_nr = release_nr + 1
+ do while ( int(releases(release_nr,2)*1.e4,kind=8).le.int((jdates(i)+obs%avetime(ind(i)))*1.e4,kind=8) .and. release_nr.lt.numpoint )
+ call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
+ write(adate,'(I6.6)') jjjjmmdd/100
+ write(areldate,'(I8.8)') jjjjmmdd
+ hhmiss=int(hhmiss/100)*100
+ write(areltime,'(I6.6)') hhmiss
+ path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
+ filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
+ ! sum the grid with previous initial grids
+ call read_init(filename, gridini_tmp)
+ gridini = gridini + gridini_tmp
+ release_nr = release_nr+1
+ nr_init = nr_init + 1
+ end do ! finished summing grid_initial
+ ! convert from ppt to fraction of one
+ gridini = gridini*1.e-12
+ gridini = gridini/real(nr_init)
+ else
+ ! filename not found
+ write(logid,*) 'WARNING: cannot find '//trim(filename)
+ go to 10
+ endif
+
+ else
+
+ ! timestamps of grid_init files match observations
+
+ ! define file name
+ write(areldate,'(I8.8)') jjjjmmdd
+ write(areltime,'(I6.6)') hhmiss
+ filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
+ print*, 'init_cini_month: filename = ',filename
+ ! read flexpart grid_init file
+ inquire(file=trim(filename),exist=lexist)
+ if(lexist) then
+ call read_init(filename, gridini)
+ ! convert from ppt to fraction of one
+ gridini = gridini*1.e-12
+ else
+ write(logid,*) 'WARNING: cannot find ',filename
+ go to 10
+ endif
+
+ endif
+
+ ! calculate termination time of trajectory
+ call caldate(jdates(i) - trajdays, jjjjmmdd, hhmiss)
+ jjjjmm = jjjjmmdd/100
+
+ ! preceding and proceding months
+ mm = jjjjmm - (jjjjmm/100)*100
+ if( mm.gt.1 ) then
+ prejjjjmm = jjjjmm - 1
+ else
+ prejjjjmm = (jjjjmm/100 - 1)*100 + 12
+ endif
+ if( mm.lt.12 ) then
+ projjjjmm = jjjjmm + 1
+ else
+ projjjjmm = (jjjjmm/100 + 1)*100 + 1
+ endif
+
+ jdmid = juldate(jjjjmm*100+14, 0)
+ jdpre = juldate(prejjjjmm*100+14, 0)
+ jdlast = juldate(lastjjjjmm*100+14, 0)
+ eomday = calceomday(prejjjjmm)
+
+ ! get new concentrations if not same month as previous observation
+ if( jjjjmm.ne.lastjjjjmm ) then
+ if( (lastjjjjmm.ne.0).and.(jdmid.eq.(jdlast+dble(eomday))) ) then
+ ! this is proceeding month
+ preconcini = concini
+ concini = proconcini
+ ! get proceeding concentration field
+ write(adate,'(I6.6)') projjjjmm
+ write(ayear,'(I4.4)') projjjjmm/100
+ select case ( files%init_type )
+ case (4)
+ pathname = str_replace(files%path_initconc, 'YYYY', ayear)
+ filename = str_replace(files%file_initconc, 'YYYYMM', adate)
+ filename = trim(pathname)//trim(filename)
+ case (5)
+ filename = str_replace(files%file_initconc, 'YYYYMM', adate)
+ filename = trim(files%path_initconc)//trim(filename)
+ end select
+ inquire(file=trim(filename),exist=lexist)
+ if(lexist) then
+ varname = files%varname_init
+ call get_initconc(files, config, filename, ntime, proconcini)
+ else
+ write(logid,*) 'WARNING: cannot find ',filename
+ write(logid,*) 'WARNING: next background file missing, copy present into next to avoid 0 value'
+ proconcini = concini
+ endif
+ else
+ ! get current concentration field
+ write(adate,'(I6.6)') jjjjmm
+ write(ayear,'(I4.4)') jjjjmm/100
+ select case ( files%init_type )
+ case (4)
+ pathname = str_replace(files%path_initconc, 'YYYY', ayear)
+ filename = str_replace(files%file_initconc, 'YYYYMM', adate)
+ filename = trim(pathname)//trim(filename)
+ case (5)
+ filename = str_replace(files%file_initconc, 'YYYYMM', adate)
+ filename = trim(files%path_initconc)//trim(filename)
+ end select
+ inquire(file=trim(filename),exist=lexist)
+ if(lexist) then
+ varname = files%varname_init
+ call get_initconc(files, config, filename, ntime, concini)
+ else
+ write(logid,*) 'WARNING: cannot find ',filename
+ go to 10
+ endif
+ ! get previous concentration field
+ write(adate,'(I6.6)') prejjjjmm
+ write(ayear,'(I4.4)') prejjjjmm/100
+ select case ( files%init_type )
+ case (4)
+ pathname = str_replace(files%path_initconc, 'YYYY', ayear)
+ filename = str_replace(files%file_initconc, 'YYYYMM', adate)
+ filename = trim(pathname)//trim(filename)
+ case (5)
+ filename = str_replace(files%file_initconc, 'YYYYMM', adate)
+ filename = trim(files%path_initconc)//trim(filename)
+ end select
+ inquire(file=trim(filename),exist=lexist)
+ if(lexist) then
+ varname = files%varname_init
+ call get_initconc(files, config, filename, ntime, preconcini)
+ else
+ write(logid,*) 'WARNING: cannot find ',filename
+ go to 10
+ endif
+ ! get next concentration field
+ write(adate,'(I6.6)') projjjjmm
+ write(ayear,'(I4.4)') projjjjmm/100
+ select case ( files%init_type )
+ case (4)
+ pathname = str_replace(files%path_initconc, 'YYYY', ayear)
+ filename = str_replace(files%file_initconc, 'YYYYMM', adate)
+ filename = trim(pathname)//trim(filename)
+ case (5)
+ filename = str_replace(files%file_initconc, 'YYYYMM', adate)
+ filename = trim(files%path_initconc)//trim(filename)
+ end select
+ inquire(file=trim(filename),exist=lexist)
+ if(lexist) then
+ varname = files%varname_init
+ call get_initconc(files, config, filename, ntime, proconcini)
+ else
+ write(logid,*) 'WARNING: cannot find ',filename
+ write(logid,*) 'WARNING: next background file missing, copy present into next to avoid 0 value'
+ proconcini = concini
+ endif
+ endif
+ endif
+
+ ! calculate contribution of initial concentration to receptor
+ do kz = 1, nzgrid
+ nz = 0
+ do while(.true.)
+ nz = nz + 1
+ if ( outheight(kz).lt.cini_alt(nz) ) exit
+ end do
+ do jy = 1, nygrid
+ n = 0
+ do while(.true.)
+ n = n + 1
+ if ( gbl_lat(jy).lt.cini_lat(n) ) exit
+ end do
+ do ix = 1, nxgrid
+ ! interpolate between concentrations
+ if( (jdates(i) - trajdays).lt.jdmid ) then
+ conc = preconcini(ix,jy,kz,1) + (concini(ix,jy,kz,1) - preconcini(ix,jy,kz,1)) * &
+ real((jdates(i) - trajdays) - jdpre, kind=4)/30.
+ else
+ conc = concini(ix,jy,kz,1) + (proconcini(ix,jy,kz,1) - concini(ix,jy,kz,1)) * &
+ real((jdates(i) - trajdays) - jdmid, kind=4)/30.
+ endif
+ cini(i,(nz-1)*nycini+n) = cini(i,(nz-1)*nycini+n) + gridini(ix,jy,kz)*conc
+ end do
+ end do
+ end do
+
+10 continue
+
+ lastjjjjmm = jjjjmm
+ prevmonth = month
+
+ end do
+
+ ! reorder cini for observations
+ do i = 1, nobs
+ obs%cini(ind(i),:) = cini(i,:)
+ end do
+
+end subroutine init_cini_month
+
+
diff --git a/source/init_co2.f90 b/source/init_co2.f90
index 34aa406..987645b 100644
--- a/source/init_co2.f90
+++ b/source/init_co2.f90
@@ -686,9 +686,11 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
! initial concentrations from termination of trajectories
! only do for a fresh run
if ( config%restart.eq.0 ) then
- call init_cini(files, config, obs)
-! call init_cams(files, config, obs)
-! call init_fpctm_hour(files, config, obs)
+ if ( files%init_type.lt.5 ) then
+ call init_cini(files, config, obs)
+ else
+ call init_cini_month(files, config, obs)
+ endif
endif
! add random perturbation for MC
diff --git a/source/init_ghg.f90 b/source/init_ghg.f90
index 7071bde..c652bbb 100644
--- a/source/init_ghg.f90
+++ b/source/init_ghg.f90
@@ -101,6 +101,7 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
call alloc_fluxes(config, fluxes)
fluxes%flx(:,:,:) = 0.
+ fluxes%flxocn(:,:,:) = 0.
! loop over prior flux files
@@ -303,11 +304,6 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
if ( config%lognormal ) then
! lognormal distribution
-! ! replace values <= 0 with smallnum since need to transform to log
-! do n = 1, ntstate
-! where( flx_box(:,n).lt.smallnum ) flx_box(:,n) = smallnum
-! states%px0((n-1)*nbox+1:n*nbox) = flx_box(:,n)
-! end do
! optimize scalar z = ln(x/xb)
! prior state vector z is zero
states%px0(1:npvar) = 0.
@@ -320,18 +316,12 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! assign prior initial mixing ratio scalars
if ( config%opt_cini ) then
if ( config%lognormal ) then
-! states%px0(npvar+1:nvar) = 1.
states%px0(npvar+1:nvar) = 0.
else
states%px0(npvar+1:nvar) = 1.
endif
endif
- ! lognormal transform to log
-! if ( config%lognormal ) then
-! states%px0(:) = log(states%px0)
-! endif
-
! time stamp of flux variables
states%xtime(:) = fluxes%time
@@ -343,10 +333,6 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
if ( config%lognormal ) then
! lognormal
-! ! errors are for log transformed state variables
-! xerr = sum(flx_box(:,:),dim=2)/real(ntstate)
-! xerr = log(xerr)
-! xerr = config%flxerr*xerr
! uncertainty of a scalar of fluxes -> use constant uncertainty
xerr(:) = config%flxerr
else
@@ -410,9 +396,13 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! initial concentrations from termination of trajectories
! only do for a fresh run
if ( config%restart.eq.0 ) then
-! call init_cams(files, config, obs)
- call init_cini(files, config, obs)
-! call init_fpctm(files, config, obs)
+ if ( files%init_type.lt.5 ) then
+ ! daily or higher frequency
+ call init_cini(files, config, obs)
+ else
+ ! monthly frequency
+ call init_cini_month(files, config, obs)
+ endif
endif
! add random perturbation for MC
diff --git a/source/initialize.f90 b/source/initialize.f90
index 5f682a5..39e640d 100644
--- a/source/initialize.f90
+++ b/source/initialize.f90
@@ -212,31 +212,20 @@ subroutine initialize(files, config)
! statres is time interval over which state variables are averaged
! statres_hr is the sub-daily time interval for state variables
! ndt is number of time steps per 24h in prior NEE
- ! flxtres is temporal resolution of low frequency fluxes (biomass burning and ocean)
! ntflx is number of time steps of low frequency fluxes
! ntstate is number of averaged time invervals for state variables
statres = config%statres
statres_hr = config%statres_hr
ndt = int(24./statres_hr)
-! flxtres = aint(365./real(config%nt_flx))
-! ntflx = nint(real(nday)/flxtres)
-! ntstate = nint(real(nday)/statres)
ntstate = int(floor(real(nday)/statres))
else
! statres is time resolution of state variables
- ! flxtres is temporal resolution of fluxes (same as state variables for non-CO2 tracers)
! ntflx is number of state time invervals
! ntstate is also number of state time invervals
statres = config%statres
-! flxtres = config%statres
-! ntflx = nint(real(nday)/statres)
ntstate = nint(real(nday)/statres)
ndt = 1
endif
-! if ( ntflx.lt.1 ) then
-! write(logid,*) 'ERROR initialize: number flux time steps < 1 -> increase time interval'
-! stop
-! endif
if ( ntstate.lt.1 ) then
write(logid,*) 'ERROR initialize: state time steps < 1 -> increase time interval'
stop
@@ -323,12 +312,14 @@ subroutine initialize(files, config)
! -----------------------
maxiter = config%maxiter
+ ncini = nycini*nzcini
+
+ ! write info to log file
+ ! ----------------------
write(logid,*) 'nday: ',nday
write(logid,*) 'statres: ',statres
-! write(logid,*) 'flxtres: ',flxtres
write(logid,*) 'ntstate: ',ntstate
-! write(logid,*) 'ntflx: ',ntflx
if ( trim(config%spec).eq.'co2' ) then
write(logid,*) 'nd_nee: ',nd_nee
write(logid,*) 'ndt: ',ndt
diff --git a/source/interpz.f90 b/source/interpz.f90
index cbb7add..52bdfb2 100644
--- a/source/interpz.f90
+++ b/source/interpz.f90
@@ -24,50 +24,106 @@
!! Interface:
!!
!! Inputs
+!! nx - number of grid cells longitude
+!! ny - number of grid cells latitude
!! nz - number of vertical layers in input
-!! geoh - geopotential heights of layers in input
+!! zkind - indicator of height reference (0 = sea level, 1 = ground level)
+!! geoh - geopotential height mid-layer in input (metres above sea level)
!! gridin - input data
!! gridout - output data (zero array on input)
!!
+!! Note: outheight is for top of layers in metres above surface so need to adjust
+!! to metres above sea level using orography
+!!
!---------------------------------------------------------------------------------------
-subroutine interpz(nz, geoh, gridin, gridout)
+subroutine interpz(nx, ny, nz, zkind, midpoints, lon, lat, geoh, gridin, gridout)
use mod_var
implicit none
- integer, intent(in) :: nz
- real, dimension(nz), intent(in) :: geoh
- real, dimension(nxgrid,nygrid,nz), intent(in) :: gridin
- real, dimension(nxgrid,nygrid,nzgrid), intent(in out) :: gridout
+ integer, intent(in) :: nx, ny, nz, zkind
+ logical, intent(in) :: midpoints
+ real, dimension(nx), intent(in) :: lon
+ real, dimension(ny), intent(in) :: lat
+ real, dimension(nx,ny,nz), intent(in) :: geoh
+ real, dimension(nx,ny,nz), intent(in) :: gridin
+ real, dimension(nx,ny,nzgrid), intent(in out) :: gridout
+
+ integer :: ix, jy, kz
+ integer :: nn, ii, jj
+ real :: hadj, dz1, dz2, ddx, ddy
+ real, dimension(nzgrid) :: halfh
+ !*** TESTING
+! real, dimension(nx,ny,nzgrid) :: htest
- integer :: ix, jy, kz
- integer :: z, z1, z2
- real :: ht
+ do nn = 1, nzgrid
+ if ( nn.eq.1 ) then
+ halfh(nn) = 0.5*outheight(nn)
+ else
+ halfh(nn) = 0.5*(outheight(nn-1) + outheight(nn))
+ endif
+ end do
+! print*, 'outheight = ',outheight
+! print*, 'halfh = ',halfh
- do kz = 1, nzgrid
- do jy = 1, nygrid
- do ix = 1, nxgrid
- ht = outheight(kz)
- z = minloc(abs(geoh(:)-ht),dim=1,mask=abs(geoh(:)-ht).eq.minval(abs(geoh(:)-ht)))
- if( geoh(z).lt.ht ) then
- z1 = z
- z2 = min(z+1,nz)
+ ! flexpart coordinates are western and southern boundary
+ if ( midpoints ) then
+ ddy = dy/2.
+ ddx = dx/2.
+ else
+ ddy = 0.
+ ddx = 0.
+ endif
+! print*, 'ddx, ddy = ',ddx,ddy
+
+ do jy = 1, ny
+ jj = minloc(gbl_lat,dim=1,mask=abs(gbl_lat+ddy-lat(jy)).eq.minval(abs(gbl_lat+ddy-lat(jy))))
+ do ix = 1, nx
+ ii = minloc(gbl_lon,dim=1,mask=abs(gbl_lon+ddx-lon(ix)).eq.minval(abs(gbl_lon+ddx-lon(ix))))
+ if ( zkind.eq.0 ) then
+ ! input is relative to sea level
+ hadj = max(orog(ii,jj),0.)
+ else
+ ! input is relative to ground level
+ hadj = 0.
+ endif
+! if (jy.eq.ny/2.and.ix.eq.nx/2) print*, 'gbl_lat, lat, gbl_lon, lon = ',gbl_lat(jj), lat(jy), gbl_lon(ii), lon(ix)
+! if (jy.eq.ny/2.and.ix.eq.nx/2) print*, 'hadj = ',hadj
+ do nn = 1, nzgrid
+ !*** TESTING
+ ! htest(ix,jy,nn) = halfh(nn)+hadj
+ do kz = 1, nz
+ if ( geoh(ix,jy,kz).gt.(halfh(nn)+hadj) ) exit
+ end do
+! if(jy.eq.ny/2.and.ix.eq.nx/2) print*, 'kz, halfh(nn)+hadj, geoh(kz) = ',kz, halfh(nn)+hadj, geoh(ix,jy,kz)
+ if ( kz.gt.nz ) then
+ gridout(ix,jy,nn) = gridin(ix,jy,nz)
+ else if ( kz.gt.1 ) then
+ dz1 = (halfh(nn) + hadj) - geoh(ix,jy,kz-1)
+ dz2 = geoh(ix,jy,kz) - (halfh(nn) + hadj)
+ gridout(ix,jy,nn) = (gridin(ix,jy,kz-1)*dz2+gridin(ix,jy,kz)*dz1)/(dz1+dz2)
+ ! if (jy.eq.ny/2.and.ix.eq.nx/2) print*, 'w1, w2 = ',dz1/(dz1+dz2), dz2/(dz1+dz2)
else
- z1 = max(z-1,1)
- z2 = z
- endif
- if( z1.eq.z2 ) then
- gridout(ix,jy,kz) = gridin(ix,jy,z)
- else
- gridout(ix,jy,kz) = gridin(ix,jy,z1)+(gridin(ix,jy,z2)-gridin(ix,jy,z1)) * &
- (ht-geoh(z1))/(geoh(z2)-geoh(z1))
- endif
+ gridout(ix,jy,nn) = gridin(ix,jy,kz)
+ endif
end do
end do
end do
+ !*** TESTING
+! open(200,file='conc_profile_in.txt',action='write',status='replace')
+! do nn = 1, nz
+! write(200,fmt='(F8.2,1X,F8.2)') geoh(ix/2,jy/2,nn), gridin(ix/2,jy/2,nn)
+! end do
+! close(200)
+! open(200,file='conc_profile_out.txt',action='write',status='replace')
+! do nn = 1, nzgrid
+! write(200,fmt='(F8.2,1X,F8.2)') htest(ix/2,jy/2,nn), gridout(ix/2,jy/2,nn)
+! end do
+! close(200)
+
end subroutine interpz
diff --git a/source/main.f90 b/source/main.f90
index 640bd6c..232cca8 100644
--- a/source/main.f90
+++ b/source/main.f90
@@ -75,6 +75,9 @@ program main
! initialize parameters
call initialize(files, config)
+ ! read orography
+ call read_orog(files)
+
! read land-sea mask
call read_lsm(files)
diff --git a/source/makefile b/source/makefile
index 103ee50..9e9e29f 100644
--- a/source/makefile
+++ b/source/makefile
@@ -32,17 +32,13 @@ SRCS = mod_var.f90 \
correl.f90 \
error_cov.f90 \
read_lsm.f90 \
+ read_orog.f90 \
mod_regions.f90 \
interpz.f90 \
convertgrid.f90 \
get_initconc.f90 \
- get_initfpctm.f90 \
- get_initfpctm_hour.f90 \
- get_initcams.f90 \
init_cini.f90 \
- init_fpctm.f90 \
- init_fpctm_hour.f90 \
- init_cams.f90 \
+ init_cini_month.f90 \
init_co2.f90 \
init_ghg.f90 \
calc_conc.f90 \
diff --git a/source/mod_perturb.f90 b/source/mod_perturb.f90
index 23beee3..f454b27 100644
--- a/source/mod_perturb.f90
+++ b/source/mod_perturb.f90
@@ -98,7 +98,6 @@ module mod_perturb
! add perturbation to state vector
if ( config%lognormal ) then
! lognormal
-! states%px0(:) = states%px0(:) * exp(real(perturb(:),kind=4))
states%px0(:) = exp(real(perturb(:),kind=4))
else
! normal
diff --git a/source/mod_save.f90 b/source/mod_save.f90
index f905a85..f2298c1 100644
--- a/source/mod_save.f90
+++ b/source/mod_save.f90
@@ -465,56 +465,44 @@ module mod_save
endif
! if inc_ocean is true state vector includes all fluxes, if false need to add ocean fluxes
if ( nbox_xy(ix,jy).gt.0 ) then
- ! regrid flux increments according to prior flux distribution and add prior flux
- if ( fluxes%flx_box(nbox_xy(ix,jy),n,1).gt.smallnum ) then
- if ( config%nested ) then
- if ( config%lognormal ) then
-! ! lognormal distribution: optimize fluxes not offsets
-! fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * &
-! fluxes%flx_nest(ix,jy,n)/fluxes%flx_box(nbox_xy(ix,jy),n,1)/numscale
-! ! optimize scalars of fluxes
- fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx_nest(ix,jy,n)/numscale
- else
- ! normal distribution: optimize offsets
+ if ( config%nested ) then
+ ! nested domain
+ if ( config%lognormal ) then
+ ! lognormal: optimize scalars of fluxes
+ fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx_nest(ix,jy,n)/numscale
+ else
+ ! normal distribution: optimize offsets
+ if ( fluxes%flx_box(nbox_xy(ix,jy),n,1).gt.smallnum ) then
+ ! regrid flux increments according to prior flux distribution and add prior flux
fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy)) * &
fluxes%flx_nest(ix,jy,n)/fluxes%flx_box(nbox_xy(ix,jy),n,1)/numscale + &
fpri(ix,jy,n)
+ else
+ ! prior in box very small so no redistribution
+ fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy))/numscale + fpri(ix,jy,n)
endif
+ endif ! normal/lognormal
+ else
+ ! no nested domain
+ if ( config%lognormal ) then
+ ! lognormal: optimize scalars of fluxes
+ fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
else
- if ( config%lognormal ) then
- ! lognormal distribution: optimize fluxes not offsets
-! fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * &
-! fluxes%flx(ix0+ix-1,jy0+jy-1,n)/fluxes%flx_box(nbox_xy(ix,jy),n,1)/numscale
-! ! optimize scalars of fluxes
- fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
- else
- ! normal distribution: optimize offsets
+ ! normal distribution: optimize offsets
+ if ( fluxes%flx_box(nbox_xy(ix,jy),n,1).gt.smallnum ) then
+ ! regrid flux increments according to prior flux distribution and add prior flux
fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy)) * &
fluxes%flx(ix0+ix-1,jy0+jy-1,n)/fluxes%flx_box(nbox_xy(ix,jy),n,1)/numscale + &
fpri(ix,jy,n)
- endif
- endif
- else
- if ( config%lognormal ) then
-! ! lognormal distribution: optimize fluxes not offsets
-! fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy)))/numscale
-! ! optimize scalars of fluxes
- if ( config%nested ) then
- fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx_nest(ix,jy,n)/numscale
else
- fpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
+ ! prior in box very small so no redistribution
+ fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy))/numscale + fpri(ix,jy,n)
endif
- else
- ! normal distribution: optimize offsets
- fpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy))/numscale + fpri(ix,jy,n)
- endif
- endif
+ endif ! normal/lognormal
+ endif ! nested
! regrid without adding prior fluxes
if ( config%lognormal ) then
! lognormal distribution
-! epri(ix,jy,n) = exp(states%pxerr0((n-1)*nbox+nbox_xy(ix,jy)))/numscale
-! epos(ix,jy,n) = exp(states%pxerr((n-1)*nbox+nbox_xy(ix,jy)))/numscale
-! xpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy)))/numscale
! errors are for z = ln(x/xb) so errors in x = xb*(exp(err) - 1)
if ( config%nested ) then
epri(ix,jy,n) = fluxes%flx_nest(ix,jy,n)*abs(exp(states%pxerr0((n-1)*nbox+nbox_xy(ix,jy))) - 1.)/numscale
@@ -531,13 +519,13 @@ module mod_save
xpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy))/numscale
endif
else
- ! add ocean fluxes if not optimized
+ ! include ocean fluxes if not optimized
if ( config%nested ) then
fpos(ix,jy,n) = fluxes%flx_nest(ix,jy,n)/numscale
else
fpos(ix,jy,n) = fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
endif
- endif
+ endif ! nbox_xy
end do
end do
end do
diff --git a/source/mod_settings.f90 b/source/mod_settings.f90
index a7942ec..32a996c 100644
--- a/source/mod_settings.f90
+++ b/source/mod_settings.f90
@@ -44,6 +44,10 @@ module mod_settings
character(len=max_path_len) :: file_bg ! prior best guess file
character(len=max_name_len) :: suffix ! observation file suffix
character(len=max_name_len) :: file_log ! log file name
+ character(len=max_path_len) :: file_orog ! orography file name
+ character(len=max_name_len) :: varname_orog ! orography variable name
+ character(len=max_name_len) :: lonname_orog ! longitude variable name
+ character(len=max_name_len) :: latname_orog ! latitude variable name
character(len=max_path_len) :: file_regions ! region definitions file
character(len=max_name_len) :: varname_regs ! regions variable name
character(len=max_name_len) :: lonname_regs ! longitude variable name
@@ -95,6 +99,7 @@ module mod_settings
character(len=max_path_len) :: path_initconc ! path to initial concentrations
character(len=max_name_len) :: file_initconc ! initial concentrations generic file name
character(len=max_name_len) :: varname_init ! initial concentration variable name
+ integer :: init_type ! initial concentration product type
character(len=max_path_len) :: file_recept ! receptor list file
end type files_t
@@ -384,6 +389,20 @@ module mod_settings
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) files%file_bg = cc
+ ! read orography file settings
+ identifier = "file_orog:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) files%file_orog = cc
+ identifier = "varname_orog:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) files%varname_orog = cc
+ identifier = "lonname_orog:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) files%lonname_orog = cc
+ identifier = "latname_orog:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) files%latname_orog = cc
+
! read region mask settings
identifier = "file_regions:"
call read_content (line, identifier, cc, cn, cl, match)
@@ -558,6 +577,9 @@ module mod_settings
identifier = "varname_init:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) files%varname_init = cc
+ identifier = "init_type:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) files%init_type = cn
! read receptor file settings
identifier = "file_recept:"
@@ -597,10 +619,11 @@ module mod_settings
character(len=200), intent(in) :: filename
type(config_t), intent(in out) :: config
- integer :: ierr
+ integer :: ierr, n
character(len=200) :: line, identifier, cc
real(kind=8) :: cn
logical :: match, cl
+ character(len=max_name_len), dimension(20) :: temp ! temp string vector (arbitrary max of 20)
! default values for logical variables
config%average_fp = .false.
@@ -771,6 +794,31 @@ module mod_settings
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%globerr = real(cn,kind=4)
+ ! read initial mixing ratio optimization settings
+ identifier = "cini_lat:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) then
+ call parse_string (cc, ",", temp(:), n)
+ nycini = n
+ allocate ( cini_lat(nycini) )
+ do n = 1, nycini
+ read(temp(n),*) cini_lat(n)
+ enddo
+ endif
+ identifier = "cini_alt:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) then
+ call parse_string (cc, ",", temp(:), n)
+ nzcini = n
+ allocate ( cini_alt(nzcini) )
+ do n = 1, nzcini
+ read(temp(n),*) cini_alt(n)
+ enddo
+ endif
+ identifier = "cinires:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) cinires = real(cn,kind=4)
+
endif
end do read_loop
diff --git a/source/mod_var.f90 b/source/mod_var.f90
index cbe5f66..62773d8 100644
--- a/source/mod_var.f90
+++ b/source/mod_var.f90
@@ -33,7 +33,7 @@ module mod_var
integer, parameter :: max_name_len=100 ! max character length of variable names
integer, parameter :: recname_len=7 ! length of receptor names
real, parameter :: numscale=1.e12 ! numeric scaling factor
- integer, parameter :: maxpoint=1000 ! max number of releases per receptor per month
+ integer, parameter :: maxpoint=10000 ! max number of releases per receptor per month
integer, parameter :: maxspec=10 ! max number of species in a flexpart run
integer, parameter :: maxlev=50 ! max number of vertical levels in flexpart output
real, parameter :: trunc=1.e-9 ! truncation threshold for eigenvalues of error covariance matrix
@@ -63,6 +63,7 @@ module mod_var
real :: dy ! latitudinal resolution of flexpart grid
real, dimension(:), allocatable :: gbl_lon ! longitudes of flexpart domain
real, dimension(:), allocatable :: gbl_lat ! latitudes of flexpart domain
+ real, dimension(:,:), allocatable :: orog ! orography on flexpart grid
! nested domain
integer :: nnxgrid ! number of longitudinal grid cells
@@ -145,10 +146,13 @@ module mod_var
integer :: npvar ! number of state variables related to fluxes
integer :: ndt ! number of state vector time steps per day (for CO2)
integer :: ndvar ! number of variables per flux time step (for CO2 2*nbox else nbox)
- integer, parameter :: ncini=4 ! number of initial mixing ratio variables to optimize
- real, parameter, dimension(4) :: cini_lat=(/-30,0,30,90/) ! vector of northern edges of latitude bands for initial mixing ratio optimization
+ integer :: ncini ! number of initial mixing ratio variables to optimize
+ integer :: nycini ! number of initial mixing ratio variables to optimize in latitudinal direction
+ integer :: nzcini ! number of initial mixing ratio variables to optimize in vertical direction
+ real, dimension(:), allocatable :: cini_lat ! vector of northern edges of latitude bands for initial mixing ratio optimization
+ real, dimension(:), allocatable :: cini_alt ! vector of upper vertical levels for initial mixing ratio optimization
integer :: ntcini ! number of initial mixing ratio time steps
- integer, parameter :: cinires=28 ! time resolution for initial mixing ratio scalars (days)
+ integer :: cinires ! time resolution for initial mixing ratio scalars (days)
integer :: neig ! number of retained eigenvalues
integer :: lastiter ! number of last iteration (if restart is true)
integer :: maxiter ! max number of iterations for congrad and m1qn3
diff --git a/source/read_orog.f90 b/source/read_orog.f90
new file mode 100644
index 0000000..0a6df66
--- /dev/null
+++ b/source/read_orog.f90
@@ -0,0 +1,109 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT: read_orog
+!---------------------------------------------------------------------------------------
+! 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 2017, Rona Thompson
+!---------------------------------------------------------------------------------------
+!
+!> read_orog
+!!
+!! Purpose: Reads the orography data from file.
+!!
+!---------------------------------------------------------------------------------------
+
+subroutine read_orog(files)
+
+ use mod_var
+ use mod_settings
+ use mod_tools, only : sort
+ use mod_save, only : check
+ use netcdf
+
+ implicit none
+
+ type (files_t), intent(in) :: files
+
+ character(len=max_path_len) :: filename
+ character(len=max_name_len) :: varname, lonname, latname
+ logical :: lexist
+ integer :: nlon, nlat
+ integer :: ncid, xid, yid, varid
+ integer :: ierr
+ real, dimension(:), allocatable :: lon, lat
+ real(kind=8), dimension(:), allocatable :: lono, lato
+ integer, dimension(:), allocatable :: indx, indy
+
+ filename = files%file_orog
+ varname = files%varname_orog
+ lonname = files%lonname_orog
+ latname = files%latname_orog
+
+ inquire(file=trim(filename),exist=lexist)
+ if (.not.lexist) then
+ write(logid,*) 'ERROR: cannot find '//trim(filename)
+ stop
+ endif
+
+ ! read file
+ ! ---------
+
+ write(logid,*) 'Reading file: '//trim(filename)
+ print*, 'lonname, latname, varname = ',lonname,latname,varname
+
+ call check( nf90_open(trim(filename),nf90_NOWRITE,ncid) )
+ call check( nf90_inq_dimid(ncid,trim(lonname),xid) )
+ call check( nf90_inq_dimid(ncid,trim(latname),yid) )
+ call check( nf90_inquire_dimension(ncid,xid,len=nlon) )
+ call check( nf90_inquire_dimension(ncid,yid,len=nlat) )
+
+ ! check file resolution matches flexpart
+ if ( nlon.ne.nxgrid.or.nlat.ne.nygrid ) then
+ write(logid,*) 'ERROR read_orog: orography file resolution does not match flexpart'
+ stop
+ endif
+
+ allocate ( lon(nlon) )
+ allocate ( lat(nlat) )
+ allocate ( indx(nlon) )
+ allocate ( indy(nlat) )
+ allocate ( orog(nxgrid,nygrid), stat=ierr )
+ if ( ierr.ne.0 ) stop 'ERROR read_orog: not enough memory'
+
+ call check( nf90_inq_varid(ncid,trim(lonname),varid) )
+ call check( nf90_get_var(ncid,varid,lon) )
+ call check( nf90_inq_varid(ncid,trim(latname),varid) )
+ call check( nf90_get_var(ncid,varid,lat) )
+ call check( nf90_inq_varid(ncid,trim(varname),varid) )
+ call check( nf90_get_var(ncid,varid,orog,start=(/1,1,1/),count=(/nxgrid,nygrid,1/)) )
+ call check( nf90_close(ncid) )
+
+ allocate ( lono(nlon) )
+ allocate ( lato(nlat) )
+ lono = dble(lon)
+ lato = dble(lat)
+ where ( lon.gt.180. ) lon = lon - 360.
+ call sort(nlon, lono, indx)
+ call sort(nlat, lato, indy)
+ orog = orog(indx,indy)
+
+ deallocate(lon)
+ deallocate(lat)
+ deallocate(lono)
+ deallocate(lato)
+ deallocate(indx)
+ deallocate(indy)
+
+end subroutine read_orog
+
diff --git a/source/simulate.f90 b/source/simulate.f90
index b31c016..d6f3ed2 100644
--- a/source/simulate.f90
+++ b/source/simulate.f90
@@ -421,11 +421,6 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
hx((ns-1)*nbox+1:ns*nbox) = hx((ns-1)*nbox+1:ns*nbox) + hbox((n-1)*nbox+1:n*nbox)
end do
endif
-! if ( config%lognormal ) then
-! ! lognormal
-! ! variable is log transformed
-! px(:) = exp(px)
-! endif
endif
call system_clock(time2, countrate)
@@ -486,6 +481,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
if ( config%method.eq.'congrad'.or.config%method.eq.'m1qn3' ) then
if ( config%lognormal ) then
! lognormal
+ ! x = xb*exp(z) where xb = 1
obs%cinipos(i) = dot_product(obs%cini(i,:), exp(states%px(npvar+(ni-1)*ncini+1:npvar+ni*ncini)))
else
! normal
@@ -533,10 +529,6 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
if ( istateuni(n).gt.0 ) then
if ( config%lognormal ) then
! lognormal
-! ! Jo'(p) = diag(exp(p))H^TR^-1(H(exp(p)) - y)
-! ! calculate as: Jo'(p) = sum( exp(p)*H_i^T*ydel_i*R_i ) (note here px = exp(p))
-! grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) = grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) + &
-! px((n-1)*ndvar+1:n*ndvar) * hx((n-1)*ndvar+1:n*ndvar) * obs%delta(i) / obs%err(i)**2
! Jo'(p) = diag(p0*exp(p))H^TR^-1(H(xb*exp(p)) - y)
! calculate as: Jo'(p) = sum( xb*exp(p)*H_i^T*ydel_i*R_i ) (note here px = xb*exp(p))
grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) = grad_o((istateuni(n)-1)*ndvar+1:istateuni(n)*ndvar) + &
--
GitLab