diff --git a/src/nudging_mod.f90 b/src/nudging_mod.f90 new file mode 100755 index 0000000000000000000000000000000000000000..4e9a1765a08bd235ffd93093bc1954d36f64500c --- /dev/null +++ b/src/nudging_mod.f90 @@ -0,0 +1,971 @@ +! 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