Skip to content
Snippets Groups Projects
Commit 3e6ede79 authored by Christine Groot Zwaaftink's avatar Christine Groot Zwaaftink
Browse files

Adding nudging_mod.f90

parent eb39b180
Branches lcm_nudge
No related tags found
No related merge requests found
! module for handling nudging of global domain filling simulations to observed
! atmospheric concentrations
!
! Author: stephan.henne@empa.ch, cgz@nilu.no
!
! Version: 2016-12-24 initial version for FLEXPART-CTM
! 2024-06-12 implementation for FLEXPARTv11 (LCM module)
!
module nudging_mod
! Variables for handling i/o unit numbers
integer, parameter :: &
iunit_start = 21, &
iunit_end = 100
integer :: iunit_table(21:100) = 0
! ! module variables
! ! ----------------------------------------------------------------------
integer :: nobssite ! number of observations sites
integer :: nobstime ! number of observation times
integer :: nr_obs_spec ! number of observed species
integer :: maxnr_nudge_spec ! number of FLEXPART species contributing to observed species
character(len=16), allocatable, dimension(:) :: name_obs_spec ! names of observed species
! ! index of FLEXPART species (in memory not in SPECIES definitions)
! ! contributing to observed species
integer, allocatable, dimension(:,:) :: nudgeidx
real, allocatable, dimension(:) :: r_molarweight
character(len=255) :: obs_filename
integer, allocatable, dimension(:) :: obs_time ! time of observations in
! seconds since simulation start
real, allocatable, dimension(:,:,:) :: obs_mf ! observed atmospheric mole fractions
real, allocatable, dimension(:) :: xobs ! longitude of observation
real, allocatable, dimension(:) :: yobs ! latitude of observation
real, allocatable, dimension(:) :: zobs ! altitude above model ground of observation
real, allocatable, dimension(:) :: hxobs ! kernel width of observation (longitude)
real, allocatable, dimension(:) :: hyobs ! kernel width (latitude)
real, allocatable, dimension(:) :: hzobs ! kernel width (altitude)
character*16, allocatable, dimension(:) :: nameobs ! name of observation site
character(4) ::obs_test_name
character(4) ::rec_test_name
real, allocatable, dimension(:) :: ff_obs ! kernel factor
real, allocatable, dimension(:) :: tau_nudge ! time constant for nudging
real, allocatable, dimension(:) :: width_temp ! temporal nudging width (seconds)
integer, allocatable, dimension(:,:,:) :: obs_idx ! temporal index for closest valid observationo ! dimension: (obsspec,obssite)
real, allocatable, dimension(:) :: d_nudge ! summed nudging tendency for each species
real, parameter :: weightair=28.97 ! CGZ ADDED FOR DEBUG, MAYBE REMOVE IF MASS CALCULATION CHANGED OR BETTER PLACED IN COM_MOD?
contains
subroutine nudging_init()
use com_mod, only: nspec, path, weightmolar
implicit none
integer :: astat, ii, ks
! ! Read observation location and nudging definitions.
call read_nudging_options(trim(path(1))//'NUDGING')
! !Copy nudging file to output directory
call system( 'cp '//trim(path(1))//'NUDGING'//' '//trim(path(2))//'NUDGING'//'')
write(*,*) ""
write(*,*) " INFO: Nudging activated for observed species"
do ii=1,nr_obs_spec
write(*,*) " INFO: ", trim(name_obs_spec(ii))
end do
write(*,*) " INFO: Observed at", nobssite, "locations"
!Calculate ratios molar weight used in calculation of mass from observed mixing ratios
do ks=1,nr_obs_spec
print*, nudgeidx(ks,1)
print*, 'full var weightmolar:', weightmolar
print*, weightmolar(nudgeidx(ks,1))
r_molarweight(nudgeidx(ks,1))=weightmolar(nudgeidx(ks,1)) / weightmolar(1)!weightmolar(1) is molar weight of species 1 = air
print*, 'Ratio:', r_molarweight(nudgeidx(ks,1))
end do
allocate(obs_idx(nr_obs_spec, nobssite,2), stat=astat)
if (astat.ne.0) write(*,*) 'ERROR: could not allocate obs_idx'
obs_idx(:,:,:) = 0
allocate(d_nudge(nspec), stat=astat)
if (astat.ne.0) write(*,*) 'ERROR: could not allocate d_nudge'
call read_observations(obs_filename)
end subroutine nudging_init
! arguments
! fn character: input file name and path
! ext_frmt_in logical, optional: is input file in extended format?
subroutine read_nudging_options(fn, ext_frmt_in)
!use io_mod, only: unit_get_free, unit_release
use com_mod, only : lsynctime, path
use point_mod, only: xlon0, ylat0, dx, dy
use par_mod, only : pi, r_earth
use readoptions_mod, only:skplin
implicit none
! arguments
character(len=*), intent(in) :: fn ! name and path of input file
logical, intent(in), optional :: ext_frmt_in ! is input file in extended format?
! locals
integer :: funit ! unit number for file input
logical :: ext_frmt ! is input file in extended format?
character(len=16) :: obsname, nudge_receptor ! name of observation location
character(len=16) :: cdummy ! dummy for reading line of input file
real :: x, y, z, hx, hy, hz ! location and kernel widths (CTM names)
real :: lon, lat, alt, h_lon, h_lat, h_alt ! location and kernel widths (preferred in LCM?)
real :: wtemp, tau
real :: xm, ym ! helper vars for kernel calculations
integer :: ii, jj ! run index
integer :: astat, ios
real, parameter :: factor=.596831 ! factor for calculation of kernel weights (15/(8*pi))
logical :: fexist, NUDGING_old ! used for checking file status
integer :: outnmlunit
! declare namelists
!nudge ctrl includes general information on observation species and which FLEXPART species should be nudged.
!seperate from species names & numbers to allow allocation
namelist /nudge_ctrl_nr/ &
obs_filename, nr_obs_spec,maxnr_nudge_spec, nobssite
namelist /nudge_ctrl_spec/ &
name_obs_spec, nudgeidx
!nudge_receptors includes coordinates of each observation site and kernel sizes
namelist /nudge_receptors/ &
nudge_receptor, lon, lat, alt, h_lon, h_lat, h_alt, wtemp, tau
! --------------------- end of variable definition -------------------------------
! setting optional parameter
ext_frmt = .true.
if (present(ext_frmt_in)) ext_frmt = ext_frmt_in
! inquire if observation file exists
inquire(file=trim(fn), exist=fexist)
if (.not. fexist) then
write(*,*) ""
write(*,*) " INFO: File NUDGING not in input folder."
write(*,*) " INFO: Add file or switch off nudging."
stop
end if
! first reading of file to determine number of observation points
call unit_get_free(funit)
NUDGING_old=.false.
open (funit,file=trim(fn),form='formatted', &
status='old',err=999)
nr_obs_spec=-3
! try namelist input
read(funit,nudge_ctrl_nr,iostat=ios)
!First read nudge_ctrl
if ((nr_obs_spec.lt.0).or.(ios.ne.0)) then
close (funit)
write(*,*) "NUDGING file in old format"
NUDGING_old=.true.
open(funit, file=trim(fn), status='old', err=999)
! skip header
call skplin(7,funit)
read(funit,'(a255)',end=88) obs_filename
write(*,*) '"'//trim(obs_filename)//'"'
call skplin(1,funit)
read(funit,'(i10)',end=88) nr_obs_spec
print*, 'nr_obs_spec', nr_obs_spec
call skplin(1,funit)
allocate(name_obs_spec(nr_obs_spec), stat=astat)
if (astat.ne.0) write(*,*) 'ERROR: could not allocate obsspec'
do jj=1,nr_obs_spec
read(funit,'(a16)',end=88) name_obs_spec(jj)
print*, 'read name_obs_spec: ', name_obs_spec(jj)
end do
call skplin(1,funit)
read(funit,'(i10)',end=88) maxnr_nudge_spec
print*, 'read maxnr_nudge_spec: ', maxnr_nudge_spec
call skplin(1,funit)
allocate(nudgeidx(nr_obs_spec,maxnr_nudge_spec), stat=astat)
allocate(r_molarweight(nr_obs_spec),stat=astat)
if (astat.ne.0) write(*,*) 'ERROR: could not allocate nr_nudge_spec'
nudgeidx(:,:) = -1
do jj=1,nr_obs_spec
ii = 1
do
read(funit,'(a16)',end=88) cdummy
if ( cdummy(1:1).eq.'#') exit
if (ii .gt. maxnr_nudge_spec) then
write(*,*) " ERROR: Number of listed nuging species exceeds maxnr_nudge_spec", maxnr_nudge_spec
write(*,*) " ERROR: Check intput file OBS_LOC"
stop !cgz test ignoring this
end if
read(cdummy,'(i16)') nudgeidx(jj,ii)
ii = ii + 1
end do
end do
print*, 'final nudgeidx:'
do jj=1,nr_obs_spec
print*, nudgeidx(jj,:)
end do
! get the number of locations
!*******************************************
jj=0
do
jj=jj+1
read(funit,*,end=88)
if (ext_frmt) read(funit,*,end=88)
read(funit,'(4x,a16)',end=88) obsname
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=88) x
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=88) y
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=88) z
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=88) hx
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=88) hy
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=88) hz
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=88) hz
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=88) hz
if (ext_frmt) read(funit,*,end=88)
! skip obs locations without name
if ((x.eq.0.) .and. (y.eq.0.) .and. obsname.eq.' ') then
jj=jj-1
cycle
endif
end do
!********************************************
else ! onamelist format
write(*,*) "NUDGING file in namelist format"
allocate(name_obs_spec(nr_obs_spec), stat=astat)
allocate(nudgeidx(nr_obs_spec,maxnr_nudge_spec), stat=astat)
read(funit,nudge_ctrl_spec,iostat=ios)
close (funit)
print*, 'final nudgeidx:'
do jj=1,nr_obs_spec
print*, nudgeidx(jj,:)
end do
endif
88 nobssite=jj-1
close(funit)
write(*,*) "Number of observation locations:", nobssite
if (NUDGING_old) then
!prepare namelist output of nudging settings
open(121,file=trim(path(2))// &
'NUDGING.namelist',status='replace',err=1000)
write(121,nml=nudge_ctrl_nr)
write(121,nml=nudge_ctrl_spec)
close (121)
endif
! allocating module variables
call init_obsloc()
if (NUDGING_old) then
call unit_get_free(outnmlunit)
open(outnmlunit,file=trim(path(2))// &
'NUDGING.namelist',status='old', access='append')
! read file again, this time assigning values to module variables
open(funit, file=trim(fn), status='old', err=999)
! skip header
call skplin(7,funit)
! find first line with '=' separator
do
read(funit,'(a16)',end=88) cdummy
if ( cdummy(1:1).eq.'=') exit
end do
jj=0
do
jj=jj+1
if (ext_frmt) read(funit,*,end=99)
read(funit,'(4x,a16)',end=99) obsname
write(*,*) obsname
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=99) lon
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=99) lat
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=99) alt
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=99) h_lon
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=99) h_lat
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=99) h_alt
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=99) wtemp
if (ext_frmt) call skplin(3,funit)
read(funit,'(4x,f11.4)',end=99) tau
if (ext_frmt) read(funit,*,end=99)
read(funit,*,end=99)
! skip obs locations without name
if ((lon.eq.0.) .and. (lat.eq.0.) .and. obsname.eq.' ') then
jj=jj-1
cycle
endif
! check parameters
if ( tau .lt. real(2*lsynctime) ) then
write(*,*) " ERROR: Nudging relaxation time (TAU) needs to be two times"
write(*,*) " ERROR: than sync time (LSYNCTIME). Check NUDGING and/or"
write(*,*) " ERROR: COMMAND files."
write(*,*) "tau: ", tau, "lsynctime: ", lsynctime
stop
end if
if ( wtemp .le. 0.0 ) then
write(*,*) " ERROR: Nudging time window (WIDTH_TEMP) needs to be larger"
write(*,*) " ERROR: zero. Check NUDGING file."
stop
end if
obsname(len_trim(obsname)+1:len(obsname)) = &
repeat(char(0), len(obsname)-len_trim(obsname)+1)
nameobs(jj) = obsname
nudge_receptor=trim(obsname)
write(*,*) "Observations: ", jj, nameobs(jj), lon, lat, alt, h_lon, h_lat, h_alt
!WRITE TO NAMELIST
!***************************************
write(outnmlunit,nml=nudge_receptors)
!****************************************
xobs(jj) = (lon-xlon0)/dx
yobs(jj) = (lat-ylat0)/dy
zobs(jj) = alt
hxobs(jj) = h_lon
hyobs(jj) = h_lat
hzobs(jj) = h_alt
width_temp(jj) = wtemp
tau_nudge(jj) = tau
! estimate nudging area
xm=r_earth*cos(y*pi/180.)*dx/180.*pi
ym=r_earth*dy/180.*pi
!print*, 'nudging area',xm, ym
if (zobs(jj).lt.hzobs(jj)) then
ff_obs(jj) = 0.5 + 0.9375*(zobs(jj)/hzobs(jj)) &
- 0.625 * (zobs(jj)/hzobs(jj))**3 &
+ 0.1875 * (zobs(jj)/hzobs(jj))**5
ff_obs(jj) = factor/ff_obs(jj)
else
ff_obs(jj) = factor
endif
write(*,*) "ff_obs", ff_obs(jj)
end do
else
write(*,*) 'Reading NUDGING in namelist format not implemented yet. Use format from FLEXPART-CTM. '
stop
endif
99 close(funit)
close(outnmlunit)
call unit_release(funit)
call unit_release(outnmlunit)
return
999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE #### '
write(*,*) ' ', trim(fn), ' ####'
write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### '
stop
1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE #### '
write(*,*) ' NUDGING.namelist #### '
write(*,*) ' #### CANNOT BE OPENED IN OUTPUT DIRECTORY #### '
stop
end subroutine read_nudging_options
! subroutine for reading observation data within model period
!
subroutine read_observations(fn)
use netcdf
use com_mod, only: bdate, edate, nspec, path
use date_mod, only: juldate
use netcdf_output_mod, only: nf90_err
implicit none
character(len=*), intent(in) :: fn
integer :: nc_id, status, astat
! double precision :: juldate
double precision, allocatable, dimension(:) :: time_in
character(len=16), allocatable, dimension(:) :: name_in
integer :: ntime_in, time_id, spec_id, site_idx
integer :: nname_in, name_id
integer, dimension(2) :: start, count
integer :: sidx, eidx, val_obs
integer :: yyyymmdd, hhmmss, ii, jj, kk
! open netcdf file
print*, 'CGZ Start reading netcdf:', trim(fn)
call nf90_err( nf90_open(trim(fn), 0, nc_id))
call system( 'cp '//trim(fn)//' '//trim(path(2))//'/')
! prepare to read time variable
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
call nf90_err(nf90_inq_dimid(nc_id, "time", time_id))
call nf90_err(nf90_inquire_dimension(nc_id, time_id, len = ntime_in))
! allocate memory for local time variable
allocate(time_in(ntime_in), stat=astat)
if (astat.ne.0) write(*,*) 'ERROR: Could not allocate time_in'
print*, 'CGZ, nc info, size time:', ntime_in
! read all times in file
call nf90_err(nf90_get_var(nc_id, time_id, time_in))
! convert time to juldate format; assuming that time is given in days since
! 1970-01-01
time_in(:) = juldate(19700101, 0) + time_in(:)
print*, 'CGZ, nc info, read time:', time_in(1)
! prepare to read observation site names
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
call nf90_err(nf90_inq_dimid(nc_id, "nameobs", name_id))
call nf90_err(nf90_inquire_dimension(nc_id, name_id, len = nname_in))
! allocate memory for local names variable
allocate(name_in(nname_in), stat=astat)
if (astat.ne.0) write(*,*) 'ERROR: Could not allocate name_in'
print*, 'CGZ, nc info, allocate size nameobs:', nname_in
call nf90_err(nf90_inq_varid(nc_id, "nameobs", name_id))
call nf90_err(nf90_get_var(nc_id, name_id, name_in) )
print*, 'CGZ, nc info, nameobs:', name_in
! TODO: the time selected period should extend beyond start/end in order to
! use directly adjacent observations
sidx = 0
eidx = 0
do ii = 1,ntime_in
if (sidx .lt. 1 .and. time_in(ii).ge.(bdate-10.D0)) sidx = ii
if (time_in(ii) .le. (edate+10.D0)) eidx = ii
end do
! number of observations during simulation period
nobstime = eidx - sidx + 1
! allocate time variable for observations and copy from time_in
allocate(obs_time(nobstime), stat=astat)
if (astat.ne.0) write(*,*) 'ERROR: Could not allocate obs_time'
obs_time(1:nobstime) = idint((time_in(sidx:eidx) - bdate) * 86400.D0)
! allocate observed mole fractions
allocate(obs_mf(nr_obs_spec, nobstime, nobssite), stat=astat)
if (astat.ne.0) write(*,*) 'ERROR: Could not allocate obs_mf'
print*, 'CGZ, nc info:', nr_obs_spec, nobstime, nobssite
! read values from file
do ii = 1,nr_obs_spec
call nf90_err(nf90_inq_varid(nc_id, trim(name_obs_spec(ii)), spec_id))
if (spec_id.lt.0) then
write(*,*) 'ERROR: Observation file: ', trim(fn)
write(*,*) 'ERROR: Does not contain requested species/variable!'
write(*,*) 'ERROR: Check for variable: ', name_obs_spec(ii)
write(*,*) 'ERROR: or change nudging species!'
stop
end if
do jj = 1,nobssite
site_idx = -1
! find match of observation by name
do kk = 1,nname_in
obs_test_name=trim(nameobs(jj))
rec_test_name=trim(adjustl(name_in(kk))) !ori
!rec_test_name=trim(name_in(kk)) !n2o
if ( rec_test_name == obs_test_name ) then
site_idx = kk
exit
end if
end do
if (site_idx .gt. 0) then
start(1) = sidx
start(2) = site_idx
count(1) = nobstime
count(2) = 1
call nf90_err(nf90_get_var(nc_id, spec_id, obs_mf(ii,:,jj), start, count))
print*, 'CGZ: site=', site_idx, rec_test_name
print*, maxval(obs_mf(ii,:,jj))
else
write(*,*) "WARNING: No observations for site ", nameobs(jj)
obs_mf(ii,:,jj) = -9999.9
end if
val_obs = 0
do kk = 1,nobstime
if (obs_mf(ii,kk,jj).gt.0.) val_obs = val_obs + 1
end do
write(*,*) " INFO: Number of valid ", trim(name_obs_spec(ii)), " observations at ", &
nameobs(jj), ":", val_obs
end do
end do
! close netcdf file
call nf90_err(nf90_close(nc_id))
deallocate(time_in)
deallocate(name_in)
end subroutine read_observations
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! nudging_individual
!
! Nudging routine to be called from timemanager
! Nudging is done for individual particles, summing all nudging tendencies
! from different observations onto each particle in range.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine nudging_individual(itime)
use com_mod, only: numpart, weightmolar, lsynctime
use particle_mod
use, intrinsic :: ieee_arithmetic
implicit none
integer, intent(in) :: itime ! time of simulation
integer :: ii, jj, ks, kf, kt, it
real :: xd, yd, zd, r2, ww_spat
real, dimension (2) :: dtn
real, allocatable, dimension(:) :: mairtr ! holding air tracer mass at observation
real, allocatable, dimension(:,:) :: sim_mf ! model mole fraction (nr_obs_spec, nobssite)
real, allocatable, dimension(:,:) :: sim_mass ! model mass at observation
! site (nr_obs_spec, nobssite)
real, allocatable, dimension(:,:,:) :: obs_mass ! observed mass at observation
! site (nr_obs_spec, nobssite)
real, allocatable, dimension(:,:) :: ten_nudge ! nudging tendency (kg/s)
integer, allocatable, dimension(:,:) :: count_nudge ! number of observations contributing to nudging tendency
real, allocatable, dimension(:,:,:) :: ww_temp ! weight due to temporal
! distance to observation
integer :: stat
!CGZ
real :: dif_mass_after, dif_mass_before, dif_ratio_after, dif_ratio_before
real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
real :: ddx,ddy
real :: rho_moist_p(2), rho_moist_i
real :: rho_dry_p(2), rho_dry_i
integer :: il, ind, indz, ix, iy, indzp, ixp, jy, jyp
integer :: ind_a, indz_a, indzp_a
! allocate local variables
allocate(mairtr(nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: Could not allocate mairtr'
allocate(sim_mf(nr_obs_spec, nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: Could not allocate sim_mf'
allocate(sim_mass(nr_obs_spec, nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: Could not allocate sim_mass'
allocate(obs_mass(nr_obs_spec, nobssite,2), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: Could not allocate obs_mass'
allocate(ten_nudge(nr_obs_spec, maxnr_nudge_spec), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: Could not allocate ten_nudge'
allocate(count_nudge(nr_obs_spec, maxnr_nudge_spec), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: Could not allocate ten_nudge'
allocate(ww_temp(nr_obs_spec, nobssite,2), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: Could not allocate ww_temp'
! 1) GET INDEX OF CLOSEST OBSERVATION
!------------------------------------------------------------------------
! The following method only works if the time step of the simulation
! (lsynctime) is smaller than the observation time step.
do jj=1, nobssite
do ks=1, nr_obs_spec
!print*, 'obsidx at start:', obs_idx(ks,jj,1), obs_idx(ks,jj,2)
! initially we search for index with smallest difference to itime
if (obs_idx(ks,jj,1) .lt. 1) then
obs_idx(ks,jj,1) = minloc(abs(obs_time(:)-itime),1)
! check if this contains a vaild observation; if not move down until
! one is discovered
if (.not.obs_mf(ks, obs_idx(ks,jj,1), jj).gt.0) then
do kt=(obs_idx(ks,jj,1)-1),1,-1
if (obs_mf(ks, kt, jj) .gt. 0) then
obs_idx(ks,jj,1) = kt
exit
end if
end do
! if no valid obs is found before search for next
if (kt .eq. 0) then
do kt=obs_idx(ks,jj,1),nobstime
if (obs_mf(ks, kt, jj) .gt. 0) then
obs_idx(ks,jj,1) = kt
exit
end if
end do
if (kt .gt. nobstime) obs_idx(ks,jj,1) = nobstime
end if
end if
end if
!write(*,*) "Current:", obs_idx(ks,jj,1), obs_mf(ks, obs_idx(ks,jj,1), jj)
kt = obs_idx(ks,jj,1) + 1
do
if (kt .gt. nobstime) then
kt = -1
exit
end if
!write(*,*) kt, obs_mf(ks, kt, jj)
if (obs_mf(ks, kt, jj) .gt. 0.) then
exit
end if
kt = kt + 1
end do
! kt now points to next valid observation while obs_idx points to
! current observation
! Now check which one is closer to current time and adjust obs_idx if
! necessary.
!CGZ, Feb 2024: In case the simulation step is between observations, and has passed 20% of the time difference between observations,
! also add influence from the next observation. This is to be able to maintain rather large kernels for remote stations but avoid
! the step-wise change in modelled concentrations.
if (kt .gt. 0) then
! print*, 'check 1:',itime, obs_idx(ks,jj,1), abs(obs_time(obs_idx(ks,jj,1))-itime), obs_mf(ks, obs_idx(ks,jj,1), jj)
! print*, 'check 2:',itime, kt, abs(obs_time(kt)-itime), obs_mf(ks, kt, jj)
if ((abs(obs_time(obs_idx(ks,jj,1))-itime) .gt. abs(obs_time(kt)-itime)) .or. &
obs_mf(ks, obs_idx(ks,jj,1), jj) .lt. 0 ) then
obs_idx(ks,jj,1) = kt
obs_idx(ks,jj,2) = kt
else
if((obs_time(obs_idx(ks,jj,1))-itime).lt.0 .and. &
(abs(obs_time(kt)-itime))/abs(obs_time(obs_idx(ks,jj,1))-itime).lt.4) then
obs_idx(ks,jj,2) = kt
else
obs_idx(ks,jj,2) = obs_idx(ks,jj,1)
endif
end if
else
obs_idx(ks,jj,2) = obs_idx(ks,jj,1)
end if
! calculate weight for temporal influence of observation
! Using tricubic weight function
! Additional check on positive observation is used to set zero weight
! which later on switches off nudging for that site/species
do it=1,2
dtn(it) = abs(real(obs_time(obs_idx(ks,jj,it))-itime))/width_temp(jj)
if (dtn(it) .lt. 1. .and. obs_mf(ks, obs_idx(ks,jj,it), jj).gt.0) then
ww_temp(ks,jj,it) = (1-dtn(it)**3)**3
else
ww_temp(ks,jj,it) = 0.
end if
end do
write(*,*) '1:', obs_idx(ks,jj,1), obs_mf(ks,obs_idx(ks,jj,1), jj),itime, obs_time(obs_idx(ks,jj,1))
write(*,*) '2:', obs_idx(ks,jj,2), obs_mf(ks,obs_idx(ks,jj,2), jj),itime, obs_time(obs_idx(ks,jj,2))
end do
end do
write(*,*) "Observation index:" !, obs_idx
do jj=1, nobssite
do ks=1, nr_obs_spec
write(*,*) obs_idx(ks,jj,:), obs_mf(ks, obs_idx(ks,jj,:), jj), ww_temp(ks,jj,:)
end do
end do
! 2) CALCULATE OBSERVED MASS AND NUDGING TENDENCY FOR EACH PARTICLE
!------------------------------------------------------------------------
! reseting mole fractions at obs
d_nudge(:) = 0.
!$OMP PARALLEL PRIVATE(ii, jj, ks, kf, xd, yd, zd, r2, ww_spat, sim_mass, mairtr, sim_mf, &
!$OMP ten_nudge, count_nudge,obs_mass, &
!$OMP rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz, &
!$OMP ddx,ddy,rho_moist_p, rho_moist_i, &
!$OMP rho_dry_p, rho_dry_i,il, ind, indz, ix, iy, indzp, ixp, jy, jyp, &
!$OMP ind_a, indz_a, indzp_a) &
!$OMP REDUCTION(+: d_nudge)
!$OMP DO
! loop over particles
do ii=1,count%alive
ten_nudge(:,:) = 0.
count_nudge(:,:)=0.
sim_mass(:,:) = 0.
! for each observation site
do jj=1, nobssite
! west-east kernel width
xd=(part(ii)%xlon-xobs(jj))/hxobs(jj)
if (xd*xd.gt.1.) cycle ! save computing time, leave loop
! vertical distance of particle
zd=(part(ii)%z-zobs(jj))/hzobs(jj)
if (zd*zd.gt.1.) cycle ! save computing time, leave loop
! north-south kernel width
yd=(part(ii)%ylat-yobs(jj))/hyobs(jj)
if (yd*yd.gt.1.) cycle ! save computing time, leave loop
! 3d distance
r2=xd*xd+yd*yd+zd*zd
! skip receptor if distance larger 1
if (r2.ge.1.) cycle
! final kernel weight
! Kernel factor not used here in
! order to give weight 1 to particles directly at receptor
ww_spat =(1.-r2) ! *ff_obs(jj)
do ks=1,nr_obs_spec
! sum up masses of nudging species
do kf=1,maxnr_nudge_spec
! leave loop if negative index is encountered
if (nudgeidx(ks,kf)<0) exit
sim_mass(ks, jj) = sim_mass(ks, jj) + mass(ii,nudgeidx(ks,kf))
end do
end do
! sum up mass of air tracer (index 1)
mairtr(jj) = mass(ii,1)
! calculate observed mass
do ks=1,nr_obs_spec
! sim_mf(ks,jj) = sim_mass(ks,jj) / mairtr(jj) * 1E9 * &
! weightair / weightmolar(nudgeidx(ks,1))
do it=1,2
obs_mass(ks,jj,it) = obs_mf(ks, obs_idx(ks,jj,it), jj) * 1E-9 * mairtr(jj) * &
r_molarweight(nudgeidx(ks,1))
end do
end do
!write(*,*) ii, obs_mf(1, obs_idx(1,jj,1), jj), "(ppb)"
!write(*,*) 'In nudge; ', ii, ks, ww_spat, mairtr(jj), 'sim: ', sim_mass(1,jj), ', obs: ', obs_mass(1, jj,1), "(kg)"
do ks=1,nr_obs_spec
do it=1,2
if (ww_temp(ks,jj,it) .gt. 0.) then
! sum up masses of nudging species
do kf=1,maxnr_nudge_spec
if (nudgeidx(ks,kf)<0) exit
if (.not.(ieee_is_finite(sim_mass(ks,jj)))) cycle
if (isnan(sim_mass(ks,jj))) cycle
ten_nudge(ks, kf) = ten_nudge(ks,kf) + &
(obs_mass(ks,jj,it) - sim_mass(ks,jj)) * &
ww_spat * ww_temp(ks,jj,it) * &
!xmass1(ii,nudgeidx(ks,kf))/(sim_mass(ks,jj)) * &
mass(ii,nudgeidx(ks,kf))/(sim_mass(ks,jj)+1.0E-30) * &
real(lsynctime)/tau_nudge(jj)
!CGZ; keep track of how many nudging
count_nudge(ks,kf)=count_nudge(ks,kf)+1
!print*, ks, kf, ten_nudge(ks,kf), count_nudge(ks,kf)
if(abs(ten_nudge(ks,kf))>1e10 .or. isnan(ten_nudge(ks,kf)))then
write(*,*) 'Large nudging tendency:', ks, ii, jj, ten_nudge(ks,kf),&
obs_mass(ks,jj,it), sim_mass(ks,jj), obs_mf(ks, obs_idx(ks,jj,it), jj)
endif
if (obs_idx(ks,jj,2).eq.obs_idx(ks,jj,1)) cycle !don't nudge twice to same observation
end do
end if
end do!loop nearest 2 observations
end do !loop nudging species
end do ! end of loop over sites
! if (ten_nudge(1,1) .ne. 0.) then
! write(*,*) ii, sum(ten_nudge(1,1:nnudgespec)), "kg"
! end if
do ks=1,nr_obs_spec
do kf=1,maxnr_nudge_spec
if (nudgeidx(ks,kf)<0) exit
if (count_nudge(ks,kf)>1) then
ten_nudge(ks,kf)=ten_nudge(ks,kf)/count_nudge(ks,kf)
endif
!cgz debugging
if ( ten_nudge(ks,kf).gt.0.) then
write(*,*) 'Nudging ; particle: ', ii, ', obs site: ', jj, ', tend: ',ten_nudge(ks,kf),', cnt:',count_nudge(ks,kf)
endif
!CGZ 20240617 test if this causes munmap_chunk(): invalid pointer
mass(ii,nudgeidx(ks,kf)) = mass(ii,nudgeidx(ks,kf)) + ten_nudge(ks,kf)
! total nudging tendency
d_nudge(nudgeidx(ks,kf)) = d_nudge(nudgeidx(ks,kf)) + ten_nudge(ks,kf)
!CGZ!!!!!!!!!!CGZ printing debug info !!!!!!!!!!!!!!!
if(ten_nudge(ks,kf).ne.0.000)then
print*, 'Current total nudge:', ks, kf, ten_nudge(ks,kf),&
d_nudge(nudgeidx(ks,kf))- ten_nudge(ks,kf),d_nudge(nudgeidx(ks,kf)),ii
endif
!!!!!!!!!CGZ printing debug info !!!!!!!!!!!!!!!
end do
end do
end do ! end of particle loop
!$OMP END DO
!$OMP END PARALLEL
write(*,*) 'End of particle loop; total nudge:', d_nudge
deallocate(mairtr)
deallocate(sim_mf)
deallocate(sim_mass)
deallocate(obs_mass)
deallocate(ten_nudge)
deallocate(ww_temp)
end subroutine nudging_individual
! subroutine for allocation of observation location variables
subroutine init_obsloc()
implicit none
integer :: stat ! allocation status
! allocate observation positions
allocate(xobs(nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: could not allocate xobs'
allocate(yobs(nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: could not allocate yobs'
allocate(zobs(nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: could not allocate zobs'
allocate(hxobs(nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: could not allocate hxobs'
allocate(hyobs(nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: could not allocate hyobs'
allocate(hzobs(nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: could not allocate hzobs'
allocate(nameobs(nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: could not allocate nameobs'
allocate(ff_obs(nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: could not allocate ff_obs'
allocate(tau_nudge(nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: could not allocate tau_nudge'
allocate(width_temp(nobssite), stat=stat)
if (stat.ne.0) write(*,*) 'ERROR: could not allocate width_temp'
end subroutine init_obsloc
!CGZ June 2024; temporarily copied some general IO functions from FLEXPART-CTM
!>>> check if there are similar subroutines with new names or if these should be placed elsewhere
subroutine unit_get_free(iunit)
implicit none
integer, intent(OUT) :: iunit
! Variables for handling i/o unit numbers
! integer, parameter :: &
! iunit_start = 21, &
! iunit_end = 100
! integer :: iunit_table(21:100) = 0
integer :: ii
iunit = -1
do ii = iunit_start, iunit_end
if (iunit_table(ii) == 0) then
iunit_table(ii) = 1
iunit = ii
exit
endif
end do
if (iunit == -1) then
write(*,*) ' *** WARNING: No more free unit numbers available *** '
endif
end subroutine unit_get_free
subroutine unit_release(iunit)
implicit none
integer, intent(in) :: iunit
! integer, parameter :: &
! iunit_start = 21, &
! iunit_end = 100
! integer :: iunit_table(21:100) = 0
if (iunit_table(iunit) == 1) then
iunit_table(iunit) = 0
else
write (*,*) ' *** WARNING: The unit ', iunit, ' is already free *** '
endif
end subroutine unit_release
!subroutine ncdf_handle_error(status)
! use netcdf
!
! integer status
!
! write(*,*) "******* Fatal NETCDF error *******"
! write(*,*) trim(NF90_strerror(status))
! write(*,*) "*************** RIP **************"
! stop
!
!end
end module nudging_mod
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment