Commit 1278c7be authored by ronesy's avatar ronesy
Browse files

Added sorting of observations chronologically in init_cini to increase efficiency

parent 005382b2
......@@ -77,13 +77,14 @@ subroutine convertgrid(midpoints, nx_in, ny_in, lon_in, lat_in, work_in, &
real :: dy_out
real :: llx_out
real :: lly_out
real, dimension(nx_in) :: lono
real, dimension(ny_in) :: lato
real(kind=8), dimension(nx_in) :: lono
real(kind=8), dimension(ny_in) :: lato
real, dimension(nx_in) :: wlon_in
real, dimension(ny_in) :: slat_in
real, dimension(nx_out) :: wlon_out
real, dimension(ny_out) :: slat_out
real, dimension(nx_in,ny_in) :: area_xy
integer, dimension(:), allocatable :: ind
ntime_in = 1
ntime_out = 1
......@@ -113,21 +114,25 @@ subroutine convertgrid(midpoints, nx_in, ny_in, lon_in, lat_in, work_in, &
endif
! check lon and lat dimensions
lato = lat_in
call sort(ny_in, lato)
lato = dble(lat_in)
allocate( ind(ny_in) )
call sort(ny_in, lato, ind)
deallocate(ind)
! indices to lat
do jy = 1, ny_in
indy(jy) = minloc(lato,dim=1,mask=lato.eq.lat_in(jy))
indy(jy) = minloc(sngl(lato),dim=1,mask=sngl(lato).eq.lat_in(jy))
end do
lat_in = lato
lat_in = sngl(lato)
where (lon_in.gt.180.) lon_in = lon_in - 360.
lono = lon_in
call sort(nx_in, lono)
lono = dble(lon_in)
allocate( ind(nx_in) )
call sort(nx_in, lono, ind)
deallocate(ind)
! indices to lon
do ix = 1, nx_in
indx(ix) = minloc(lono,dim=1,mask=lono.eq.lon_in(ix))
indx(ix) = minloc(sngl(lono),dim=1,mask=sngl(lono).eq.lon_in(ix))
end do
lon_in = lono
lon_in = sngl(lono)
! reorder flx
work_in(:,:,:) = work_in(indx,indy,:)
......
......@@ -43,6 +43,7 @@ subroutine init_cini(files, obs)
use mod_var
use mod_dates
use mod_obs
use mod_tools
implicit none
......@@ -58,6 +59,9 @@ subroutine init_cini(files, obs)
integer :: ix, jy, kz, i
real :: conc
real(kind=8) :: jdmid, jdpre, jdlast
real(kind=8), dimension(nobs) :: jdates
real, dimension(nobs) :: cini
integer, dimension(nobs) :: ind
real, dimension(nxgrid,nygrid,nzgrid) :: gridini
real, dimension(nxgrid,nygrid,nzgrid) :: concini, preconcini, proconcini
......@@ -69,16 +73,25 @@ subroutine init_cini(files, obs)
concini(:,:,:) = 0.
preconcini(:,:,:) = 0.
proconcini(:,:,:) = 0.
cini(:) = 0.
! sort obs chronologically for efficiency
jdates = obs%jdate(:)
call sort(nobs,jdates,ind)
! print*, 'init_cini: ind = ',ind
! print*, 'init_cini: jdates = ',jdates
! print*, 'init_cini: obs%jdate = ',obs%jdate
do i = 1, nobs
! define file name
call caldate(obs%jdate(i), jjjjmmdd, hhmiss)
call caldate(jdates(i), jjjjmmdd, hhmiss)
write(adate,'(I6.6)') jjjjmmdd/100
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(i))//'/'//trim(adate)//'/'
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 = ',filename
! read flexpart initial concentrations weighting
inquire(file=trim(filename),exist=lexist)
......@@ -92,7 +105,7 @@ subroutine init_cini(files, obs)
endif
! calculate termination time of trajectory
call caldate(obs%jdate(i) - trajdays, jjjjmmdd, hhmiss)
call caldate(jdates(i) - trajdays, jjjjmmdd, hhmiss)
jjjjmm = jjjjmmdd/100
! preceding and proceding months
......@@ -176,14 +189,14 @@ subroutine init_cini(files, obs)
do jy = 1, nygrid
do ix = 1, nxgrid
! interpolate between concentrations
if( (obs%jdate(i) - trajdays).lt.jdmid ) then
if( (jdates(i) - trajdays).lt.jdmid ) then
conc = preconcini(ix,jy,kz) + (concini(ix,jy,kz) - preconcini(ix,jy,kz)) * &
real((obs%jdate(i) - trajdays) - jdpre, kind=4)/30.
real((jdates(i) - trajdays) - jdpre, kind=4)/30.
else
conc = concini(ix,jy,kz) + (proconcini(ix,jy,kz) - concini(ix,jy,kz)) * &
real((obs%jdate(i) - trajdays) - jdmid, kind=4)/30.
real((jdates(i) - trajdays) - jdmid, kind=4)/30.
endif
obs%cini(i) = obs%cini(i) + gridini(ix,jy,kz)*conc
cini(i) = cini(i) + gridini(ix,jy,kz)*conc
end do
end do
end do
......@@ -194,6 +207,14 @@ subroutine init_cini(files, obs)
end do
! reorder cini for observations
! print*, 'init_cini: ind = ',ind
! print*, 'init_cini: cini = ',cini
do i = 1, nobs
obs%cini(ind(i)) = cini(i)
end do
! print*, 'init_cini: obs%cini = ',obs%cini
end subroutine init_cini
......@@ -62,17 +62,18 @@ module mod_regions
implicit none
type (files_t), intent(in) :: files
type (config_t), intent(in) :: config
integer :: ix, jy
integer, dimension(nxregrid,nyregrid) :: hloc_xy
real, dimension(nxregrid*nyregrid) :: work
integer :: hh
integer, dimension(:), allocatable :: cnt
character(len=max_path_len) :: rowfmt
integer :: n, i, ierr
real :: area
type (files_t), intent(in) :: files
type (config_t), intent(in) :: config
integer :: ix, jy
integer, dimension(nxregrid,nyregrid) :: hloc_xy
real(kind=8), dimension(nxregrid*nyregrid) :: work
integer, dimension(nxregrid*nyregrid) :: ind
integer :: hh
integer, dimension(:), allocatable :: cnt
character(len=max_path_len) :: rowfmt
integer :: n, i, ierr
real :: area
! define region boxes
! -------------------
......@@ -134,19 +135,19 @@ module mod_regions
cnt(:) = 0
! local time of region is median time of grid cells in box
work(:) = 0.
work(:) = 0d0
do n = 1, nbox
i = 0
do jy = 1, nyregrid
do ix = 1, nxregrid
if( nbox_xy(ix,jy).eq.n ) then
i = i + 1
work(i) = real(hloc_xy(ix,jy))
work(i) = dble(hloc_xy(ix,jy))
endif
end do
end do
call sort(i, work(1:i))
hloc(n) = work(i/2+1)
call sort(i, work(1:i), ind(1:i))
hloc(n) = sngl(work(i/2+1))
end do
deallocate(cnt)
......
......@@ -46,30 +46,38 @@ module mod_tools
!!
! --------------------------------------------------
subroutine sort(n, arr)
subroutine sort(n, arr, ind)
implicit none
integer, intent(in) :: n
real, dimension(n), intent(in out) :: arr
integer, intent(in) :: n
real(kind=8), dimension(n), intent(in out) :: arr
integer, dimension(n), intent(in out) :: ind
integer, parameter :: m=7, nstack=50
integer :: i, ir, j, k, l, jstack
integer, dimension(nstack) :: istack
real :: a, temp
integer, parameter :: m=7, nstack=50
integer :: i, ir, j, k, l, jstack
integer, dimension(nstack) :: istack
real(kind=8) :: a, temp
integer :: ia, itemp
do i=1,n
ind(i)=i
end do
jstack=0
l=1
ir=n
1 if(ir-l.lt.m) then
do j=l+1,ir
a=arr(j)
ia=ind(j)
do i=j-1,1,-1
if(arr(i).le.a) go to 2
arr(i+1)=arr(i)
ind(i+1)=ind(i)
end do
i=0
2 arr(i+1)=a
ind(i+1)=ia
end do
if(jstack.eq.0) return
ir=istack(jstack)
......@@ -78,26 +86,39 @@ module mod_tools
else
k=(l+ir)/2
temp=arr(k)
itemp=ind(k)
arr(k)=arr(l+1)
arr(l+1)=temp
ind(k)=ind(l+1)
ind(l+1)=itemp
if(arr(l+1).gt.arr(ir)) then
temp=arr(l+1)
itemp=ind(l+1)
arr(l+1)=arr(ir)
arr(ir)=temp
ind(l+1)=ind(ir)
ind(ir)=itemp
endif
if(arr(l).gt.arr(ir)) then
temp=arr(l)
itemp=ind(l)
arr(l)=arr(ir)
arr(ir)=temp
ind(l)=ind(ir)
ind(ir)=itemp
endif
if(arr(l+1).gt.arr(l)) then
temp=arr(l+1)
itemp=ind(l+1)
arr(l+1)=arr(l)
arr(l)=temp
ind(l+1)=ind(l)
ind(l)=itemp
endif
i=l+1
j=ir
a=arr(l)
ia=ind(l)
3 continue
i=i+1
if(arr(i).lt.a) go to 3
......@@ -106,11 +127,16 @@ module mod_tools
if(arr(j).gt.a) go to 4
if(j.lt.i) go to 5
temp=arr(i)
itemp=ind(i)
arr(i)=arr(j)
arr(j)=temp
ind(i)=ind(j)
ind(j)=itemp
go to 3
5 arr(l)=arr(j)
arr(j)=a
ind(l)=ind(j)
ind(j)=ia
jstack=jstack+2
if(jstack.gt.nstack) print*, 'ERROR sort: nstack too small'
if(ir-i+1.ge.j-l) then
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment