Skip to content
Snippets Groups Projects
Commit 67b13e92 authored by Sabine Eckhardt's avatar Sabine Eckhardt
Browse files

update to wetdepo.f90

git-svn-id: http://flexpart.flexpart.eu:8088/svn/FlexPart90/branches/FLEXPART_9.1.3@13 ef8cc7e1-21b7-489e-abab-c1baa636049d
parents
No related branches found
No related tags found
No related merge requests found
Showing with 2970 additions and 0 deletions
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
program flexpart
!*****************************************************************************
! *
! This is the Lagrangian Particle Dispersion Model FLEXPART. *
! The main program manages the reading of model run specifications, etc. *
! All actual computing is done within subroutine timemanager. *
! *
! Author: A. Stohl *
! *
! 18 May 1996 *
! *
!*****************************************************************************
! *
! Variables: *
! *
! Constants: *
! *
!*****************************************************************************
use point_mod
use par_mod
use com_mod
use conv_mod
implicit none
integer :: i,j,ix,jy,inest
integer :: idummy = -320
character(len=256) :: pathfile
! Generate a large number of random numbers
!******************************************
do i=1,maxrand-1,2
call gasdev1(idummy,rannumb(i),rannumb(i+1))
end do
call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
! Print the GPL License statement
!*******************************************************
! print*,'Welcome to FLEXPART Version 9.1 (Build 20121029)'
print*,'Welcome to FLEXPART Version 9.1.3 (Build 2013158)'
print*,'FLEXPART is free software released under the GNU Genera'// &
'l Public License.'
! Read the pathnames where input/output files are stored
!*******************************************************
select case (iargc())
case (1)
call getarg(1,pathfile)
if (trim(pathfile).eq.'-v') then
!write(*,*) 'FLEXPART Version 9.1.3 (Build 2013158)'
stop
endif
case (0)
write(pathfile,'(a11)') './pathnames'
end select
call readpaths(pathfile)
print*,length(4)
! Read the user specifications for the current model run
!*******************************************************
call readcommand
! Read the age classes to be used
!********************************
call readageclasses
! Read, which wind fields are available within the modelling period
!******************************************************************
call readavailable
! Read the model grid specifications,
! both for the mother domain and eventual nests
!**********************************************
call gridcheck
call gridcheck_nests
! Read the output grid specifications
!************************************
call readoutgrid
if (nested_output.eq.1) call readoutgrid_nest
! Read the receptor points for which extra concentrations are to be calculated
!*****************************************************************************
call readreceptors
! Read the physico-chemical species property table
!*************************************************
!SEC: now only needed SPECIES are read in readreleases.f
!call readspecies
! Read the landuse inventory
!***************************
call readlanduse
! Assign fractional cover of landuse classes to each ECMWF grid point
!********************************************************************
call assignland
! Read the coordinates of the release locations
!**********************************************
call readreleases
! Read and compute surface resistances to dry deposition of gases
!****************************************************************
call readdepo
! Convert the release point coordinates from geografical to grid coordinates
!***************************************************************************
call coordtrafo
! Initialize all particles to non-existent
!*****************************************
do j=1,maxpart
itra1(j)=-999999999
end do
! For continuation of previous run, read in particle positions
!*************************************************************
if (ipin.eq.1) then
call readpartpositions
else
numpart=0
numparticlecount=0
endif
! Calculate volume, surface area, etc., of all output grid cells
! Allocate fluxes and OHfield if necessary
!***************************************************************
call outgrid_init
if (nested_output.eq.1) call outgrid_init_nest
! Read the OH field
!******************
if (OHREA.eqv..TRUE.) &
call readOHfield
! Write basic information on the simulation to a file "header"
! and open files that are to be kept open throughout the simulation
!******************************************************************
call writeheader
if (nested_output.eq.1) call writeheader_nest
open(unitdates,file=path(2)(1:length(2))//'dates')
call openreceptors
if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
! Releases can only start and end at discrete times (multiples of lsynctime)
!***************************************************************************
do i=1,numpoint
ireleasestart(i)=nint(real(ireleasestart(i))/ &
real(lsynctime))*lsynctime
ireleaseend(i)=nint(real(ireleaseend(i))/ &
real(lsynctime))*lsynctime
end do
! Initialize cloud-base mass fluxes for the convection scheme
!************************************************************
do jy=0,nymin1
do ix=0,nxmin1
cbaseflux(ix,jy)=0.
end do
end do
do inest=1,numbnests
do jy=0,nyn(inest)-1
do ix=0,nxn(inest)-1
cbasefluxn(ix,jy,inest)=0.
end do
end do
end do
! Calculate particle trajectories
!********************************
call timemanager
write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
&XPART MODEL RUN!'
end program flexpart
SHELL = /bin/bash
MAIN = FLEXPART_GFORTRAN
#
FC = gfortran
#INCPATH = /usr/local/ecmwf_tools/reloc/include
#LIBPATH1 = /usr/local/ecmwf_tools/reloc/lib
#LIBPATH2 = /usr/local/lib
#FFLAGS = -O3 -g -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
FFLAGS = -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
#LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -lopenjpeg
# options for nordre
INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include
LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
LIBPATH2 = /usr/lib/x86_64-linux-gnu/
#FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
#options for Lion
libs_dir=/Users/ignacio_in_EBCRPA/flexpart/libs/
INCPATH = $(libs_dir)/grib_api-1.9.9_dir/include
LIBPATH1 = $(libs_dir)/grib_api-1.9.9_dir/lib
LIBPATH2 = $(libs_dir)/jasper_dir/lib
MODOBJS = \
par_mod.o com_mod.o \
conv_mod.o hanna_mod.o \
interpol_mod.o cmapf_mod.o \
unc_mod.o oh_mod.o \
xmass_mod.o flux_mod.o \
point_mod.o outg_mod.o
OBJECTS = \
writeheader.o assignland.o\
calcpar.o part0.o \
caldate.o partdep.o \
coordtrafo.o psih.o \
raerod.o \
drydepokernel.o random.o \
readavailable.o \
ew.o readcommand.o \
advance.o readdepo.o \
releaseparticles.o psim.o \
FLEXPART.o readlanduse.o \
getfields.o init_domainfill.o\
interpol_wind.o readoutgrid.o \
interpol_all.o readpaths.o \
getrb.o readreceptors.o \
getrc.o readreleases.o \
getvdep.o readspecies.o \
interpol_misslev.o readwind.o \
conccalc.o richardson.o \
concoutput.o scalev.o \
pbl_profile.o readOHfield.o\
juldate.o timemanager.o \
interpol_vdep.o interpol_rain.o \
verttransform.o partoutput.o \
hanna.o wetdepokernel.o \
mean.o wetdepo.o \
hanna_short.o windalign.o \
obukhov.o gridcheck.o \
hanna1.o initialize.o \
gridcheck_nests.o \
readwind_nests.o calcpar_nests.o \
verttransform_nests.o interpol_all_nests.o \
interpol_wind_nests.o interpol_misslev_nests.o \
interpol_vdep_nests.o interpol_rain_nests.o \
getvdep_nests.o \
readageclasses.o readpartpositions.o \
calcfluxes.o fluxoutput.o \
qvsat.o skplin.o \
convmix.o calcmatrix.o \
convect43c.o redist.o \
sort2.o distance.o \
centerofmass.o plumetraj.o \
openouttraj.o calcpv.o \
calcpv_nests.o distance2.o \
clustering.o interpol_wind_short.o \
interpol_wind_short_nests.o shift_field_0.o \
shift_field.o outgrid_init.o \
openreceptors.o boundcond_domainfill.o\
partoutput_short.o readoutgrid_nest.o \
outgrid_init_nest.o writeheader_nest.o \
concoutput_nest.o wetdepokernel_nest.o \
drydepokernel_nest.o zenithangle.o \
ohreaction.o getvdep_nests.o \
initial_cond_calc.o initial_cond_output.o \
dynamic_viscosity.o get_settling.o
$(MAIN): $(MODOBJS) $(OBJECTS)
$(FC) *.o -o $(MAIN) $(LDFLAGS)
$(OBJECTS): $(MODOBJS)
%.o: %.f90
$(FC) -c $(FFLAGS) $<
clean:
rm *.o *.mod
SHELL = /bin/bash
MAIN = FLEXPART_GFORTRAN
#
FC = gfortran
#INCPATH = /usr/local/ecmwf_tools/reloc/include
#LIBPATH1 = /usr/local/ecmwf_tools/reloc/lib
#LIBPATH2 = /usr/local/lib
#FFLAGS = -O3 -g -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
FFLAGS = -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
#LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -lopenjpeg
# options for nordre
INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include
LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
LIBPATH2 = /usr/lib/x86_64-linux-gnu/
#FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
#options for Lion
libs_dir=/Users/ignacio_in_EBCRPA/flexpart/libs/
INCPATH = $(libs_dir)/grib_api-1.9.9_dir/include
LIBPATH1 = $(libs_dir)/grib_api-1.9.9_dir/lib
LIBPATH2 = $(libs_dir)/jasper_dir/lib
MODOBJS = \
par_mod.o com_mod.o \
conv_mod.o hanna_mod.o \
interpol_mod.o cmapf_mod.o \
unc_mod.o oh_mod.o \
xmass_mod.o flux_mod.o \
point_mod.o outg_mod.o
OBJECTS = \
writeheader.o assignland.o\
calcpar.o part0.o \
caldate.o partdep.o \
coordtrafo.o psih.o \
raerod.o \
drydepokernel.o random.o \
readavailable.o \
ew.o readcommand.o \
advance.o readdepo.o \
releaseparticles.o psim.o \
FLEXPART.o readlanduse.o \
getfields.o init_domainfill.o\
interpol_wind.o readoutgrid.o \
interpol_all.o readpaths.o \
getrb.o readreceptors.o \
getrc.o readreleases.o \
getvdep.o readspecies.o \
interpol_misslev.o readwind.o \
conccalc.o richardson.o \
concoutput.o scalev.o \
pbl_profile.o readOHfield.o\
juldate.o timemanager.o \
interpol_vdep.o interpol_rain.o \
verttransform.o partoutput.o \
hanna.o wetdepokernel.o \
mean.o wetdepo.o \
hanna_short.o windalign.o \
obukhov.o gridcheck.o \
hanna1.o initialize.o \
gridcheck_nests.o \
readwind_nests.o calcpar_nests.o \
verttransform_nests.o interpol_all_nests.o \
interpol_wind_nests.o interpol_misslev_nests.o \
interpol_vdep_nests.o interpol_rain_nests.o \
getvdep_nests.o \
readageclasses.o readpartpositions.o \
calcfluxes.o fluxoutput.o \
qvsat.o skplin.o \
convmix.o calcmatrix.o \
convect43c.o redist.o \
sort2.o distance.o \
centerofmass.o plumetraj.o \
openouttraj.o calcpv.o \
calcpv_nests.o distance2.o \
clustering.o interpol_wind_short.o \
interpol_wind_short_nests.o shift_field_0.o \
shift_field.o outgrid_init.o \
openreceptors.o boundcond_domainfill.o\
partoutput_short.o readoutgrid_nest.o \
outgrid_init_nest.o writeheader_nest.o \
concoutput_nest.o wetdepokernel_nest.o \
drydepokernel_nest.o zenithangle.o \
ohreaction.o getvdep_nests.o \
initial_cond_calc.o initial_cond_output.o \
dynamic_viscosity.o get_settling.o
$(MAIN): $(MODOBJS) $(OBJECTS)
$(FC) *.o -o $(MAIN) $(LDFLAGS)
$(OBJECTS): $(MODOBJS)
%.o: %.f90
$(FC) -c $(FFLAGS) $<
clean:
rm *.o *.mod
SHELL = /bin/bash
MAIN = FLEXPART_GFS_GFORTRAN
#
FC = gfortran
INCPATH = /usr/local/ecmwf_tools/reloc/include
LIBPATH1 = /usr/local/ecmwf_tools/reloc/lib
LIBPATH2 = /usr/local/lib
#FFLAGS = -O3 -g -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
FFLAGS = -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -lopenjpeg
#
MODOBJS = \
par_mod.o com_mod.o \
conv_mod.o hanna_mod.o \
interpol_mod.o cmapf_mod.o \
unc_mod.o oh_mod.o \
xmass_mod.o flux_mod.o \
point_mod.o outg_mod.o
OBJECTS = \
writeheader.o assignland.o\
calcpar_gfs.o part0.o \
caldate.o partdep.o \
coordtrafo.o psih.o \
raerod.o \
drydepokernel.o random.o \
erf.o readavailable.o \
ew.o readcommand.o \
advance.o readdepo.o \
releaseparticles.o psim.o \
FLEXPART.o readlanduse.o \
getfields.o init_domainfill.o\
interpol_wind.o readoutgrid.o \
interpol_all.o readpaths.o \
getrb.o readreceptors.o \
getrc.o readreleases.o \
getvdep.o readspecies.o \
interpol_misslev.o readwind_gfs.o \
conccalc.o richardson_gfs.o \
concoutput.o scalev.o \
pbl_profile.o readOHfield.o\
juldate.o timemanager.o \
interpol_vdep.o interpol_rain.o \
verttransform_gfs.o partoutput.o \
hanna.o wetdepokernel.o \
mean.o wetdepo.o \
hanna_short.o windalign.o \
obukhov_gfs.o gridcheck_gfs.o \
hanna1.o initialize.o \
gridcheck_nests.o \
readwind_nests.o calcpar_nests.o \
verttransform_nests.o interpol_all_nests.o \
interpol_wind_nests.o interpol_misslev_nests.o \
interpol_vdep_nests.o interpol_rain_nests.o \
getvdep_nests.o \
readageclasses.o readpartpositions.o \
calcfluxes.o fluxoutput.o \
qvsat.o skplin.o \
convmix_gfs.o calcmatrix_gfs.o \
convect43c.o redist.o \
sort2.o distance.o \
centerofmass.o plumetraj.o \
openouttraj.o calcpv.o \
calcpv_nests.o distance2.o \
clustering.o interpol_wind_short.o \
interpol_wind_short_nests.o shift_field_0.o \
shift_field.o outgrid_init.o \
openreceptors.o boundcond_domainfill.o\
partoutput_short.o readoutgrid_nest.o \
outgrid_init_nest.o writeheader_nest.o \
concoutput_nest.o wetdepokernel_nest.o \
drydepokernel_nest.o zenithangle.o \
ohreaction.o getvdep_nests.o \
initial_cond_calc.o initial_cond_output.o \
dynamic_viscosity.o get_settling.o
$(MAIN): $(MODOBJS) $(OBJECTS)
$(FC) *.o -o $(MAIN) $(LDFLAGS)
$(OBJECTS): $(MODOBJS)
%.o: %.f90
$(FC) -c $(FFLAGS) $<
clean:
rm *.o *.mod
This diff is collapsed.
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine assignland
!*****************************************************************************
! *
! This routine assigns fractions of the 13 landuse classes to each ECMWF *
! grid point. *
! The landuse inventory of *
! *
! Belward, A.S., Estes, J.E., and Kline, K.D., 1999, *
! The IGBP-DIS 1-Km Land-Cover Data Set DISCover: *
! A Project Overview: Photogrammetric Engineering and Remote Sensing , *
! v. 65, no. 9, p. 1013-1020 *
! *
! if there are no data in the inventory *
! the ECMWF land/sea mask is used to distinguish *
! between sea (-> ocean) and land (-> grasslands). *
! *
! Author: A. Stohl *
! *
! 5 December 1996 *
! 8 February 1999 Additional use of nests, A. Stohl *
! 29 December 2006 new landuse inventory, S. Eckhardt *
!*****************************************************************************
! *
! Variables: *
! xlanduse fractions of numclass landuses for each model grid point *
! landinvent landuse inventory (0.3 deg resolution) *
! *
!*****************************************************************************
use par_mod
use com_mod
implicit none
integer :: ix,jy,k,l,li,nrefine,iix,jjy
integer,parameter :: lumaxx=1200,lumaxy=600
integer,parameter :: xlon0lu=-180,ylat0lu=-90
real,parameter :: dxlu=0.3
real :: xlon,ylat,sumperc,p,xi,yj
real :: xlandusep(lumaxx,lumaxy,numclass)
! character*2 ck
do ix=1,lumaxx
do jy=1,lumaxy
do k=1,numclass
xlandusep(ix,jy,k)=0.
end do
sumperc=0.
do li=1,3
sumperc=sumperc+landinvent(ix,jy,li+3)
end do
do li=1,3
k=landinvent(ix,jy,li)
if (sumperc.gt.0) then
p=landinvent(ix,jy,li+3)/sumperc
else
p=0
endif
! p has values between 0 and 1
xlandusep(ix,jy,k)=p
end do
end do
end do
! do 13 k=1,11
! write (ck,'(i2.2)') k
! open(4,file='xlandusetest'//ck,form='formatted')
! do 11 ix=1,lumaxx
!11 write (4,*) (xlandusep(ix,jy,k),jy=1,lumaxy)
!11 write (4,*) (landinvent(ix,jy,k),jy=1,lumaxy)
!13 close(4)
! write (*,*) xlon0,ylat0,xlon0n(1),ylat0n(1),nxmin1,nymin1
! write (*,*) dx, dy, dxout, dyout, ylat0, xlon0
nrefine=10
do ix=0,nxmin1
do jy=0,nymin1
do k=1,numclass
sumperc=0.
xlanduse(ix,jy,k)=0.
end do
do iix=1, nrefine
xlon=(ix+(iix-1)/real(nrefine))*dx+xlon0 ! longitude, should be between -180 and 179
if (xlon.ge.(xlon0lu+lumaxx*dxlu)) then
xlon=xlon-lumaxx*dxlu
endif
do jjy=1, nrefine
ylat=(jy+(jjy-1)/real(nrefine))*dy+ylat0 ! and lat. of each gridpoint
xi=int((xlon-xlon0lu)/dxlu)+1
yj=int((ylat-ylat0lu)/dxlu)+1
if (xi.gt.lumaxx) xi=xi-lumaxx
if (yj.gt.lumaxy) yj=yj-lumaxy
if (xi.lt.0) then
write (*,*) 'problem with landuseinv sampling: ', &
xlon,xlon0lu,ix,iix,xlon0,dx,nxmax
stop
endif
do k=1,numclass
xlanduse(ix,jy,k)= &
xlanduse(ix,jy,k)+xlandusep(int(xi),int(yj),k)
sumperc=sumperc+xlanduse(ix,jy,k) ! just for the check if landuseinv. is available
end do
end do
end do
if (sumperc.gt.0) then ! detailed landuse available
sumperc=0.
do k=1,numclass
xlanduse(ix,jy,k)= &
xlanduse(ix,jy,k)/real(nrefine*nrefine)
sumperc=sumperc+xlanduse(ix,jy,k)
end do
!cc the sum of all categories should be 1 ... 100 percent ... in order to get vdep right!
if (sumperc.lt.1-1E-5) then
do k=1,numclass
xlanduse(ix,jy,k)= &
xlanduse(ix,jy,k)/sumperc
end do
endif
else
if (lsm(ix,jy).lt.0.1) then ! over sea -> ocean
xlanduse(ix,jy,3)=1.
else ! over land -> rangeland
xlanduse(ix,jy,7)=1.
endif
endif
end do
end do
!***********************************
! for test: write out xlanduse
! open(4,file='landusetest',form='formatted')
! do 56 k=1,13
! do 55 ix=0,nxmin1
!55 write (4,*) (xlanduse(ix,jy,k),jy=0,nymin1)
!56 continue
! close(4)
! write (*,*) 'landuse written'
!stop
! open(4,file='landseatest'//ck,form='formatted')
! do 57 ix=0,nxmin1
!57 write (4,*) (lsm(ix,jy),jy=0,nymin1)
! write (*,*) 'landseamask written'
!****************************************
! Same as above, but for the nested grids
!****************************************
!************** TEST ********************
! dyn(1)=dyn(1)/40
! dxn(1)=dxn(1)/40
! xlon0n(1)=1
! ylat0n(1)=50
!************** TEST ********************
do l=1,numbnests
do ix=0,nxn(l)-1
do jy=0,nyn(l)-1
do k=1,numclass
sumperc=0.
xlandusen(ix,jy,k,l)=0.
end do
do iix=1, nrefine
xlon=(ix+(iix-1)/real(nrefine))*dxn(l)+xlon0n(l)
do jjy=1, nrefine
ylat=(jy+(jjy-1)/real(nrefine))*dyn(l)+ylat0n(l)
xi=int((xlon-xlon0lu)/dxlu)+1
yj=int((ylat-ylat0lu)/dxlu)+1
if (xi.gt.lumaxx) xi=xi-lumaxx
if (yj.gt.lumaxy) yj=yj-lumaxy
do k=1,numclass
xlandusen(ix,jy,k,l)=xlandusen(ix,jy,k,l)+ &
xlandusep(int(xi),int(yj),k)
sumperc=sumperc+xlandusen(ix,jy,k,l)
end do
end do
end do
if (sumperc.gt.0) then ! detailed landuse available
sumperc=0.
do k=1,numclass
xlandusen(ix,jy,k,l)= &
xlandusen(ix,jy,k,l)/real(nrefine*nrefine)
sumperc=sumperc+xlandusen(ix,jy,k,l)
end do
!cc the sum of all categories should be 1 ... 100 percent ... in order to get vdep right!
if (sumperc.lt.1-1E-5) then
do k=1,numclass
xlandusen(ix,jy,k,l)=xlandusen(ix,jy,k,l)/sumperc
end do
endif
else ! check land/sea mask
if (lsmn(ix,jy,l).lt.0.1) then ! over sea -> ocean
xlandusen(ix,jy,3,l)=1.
else ! over land -> grasslands
xlandusen(ix,jy,7,l)=1.
endif
endif
end do
end do
end do
!***********************************
! for test: write out xlanduse
! do 66 k=1,11
! write (ck,'(i2.2)') k
! open(4,file='nlandusetest'//ck,form='formatted')
! do 65 ix=0,nxn(1)-1
!65 write (4,*) (xlandusen(ix,jy,k,1),jy=0,nyn(1)-1)
!66 close(4)
! write (*,*) 'landuse nested written'
end subroutine assignland
This diff is collapsed.
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcfluxes(nage,jpart,xold,yold,zold)
! i i i i i
!*****************************************************************************
! *
! Calculation of the gross fluxes across horizontal, eastward and *
! northward facing surfaces. The routine calculates the mass flux *
! due to the motion of only one particle. The fluxes of subsequent calls *
! to this subroutine are accumulated until the next output is due. *
! Upon output, flux fields are re-set to zero in subroutine fluxoutput.f.*
! *
! Author: A. Stohl *
! *
! 04 April 2000 *
! *
!*****************************************************************************
! *
! Variables: *
! *
! nage Age class of the particle considered *
! jpart Index of the particle considered *
! xold,yold,zold "Memorized" old positions of the particle *
! *
!*****************************************************************************
use flux_mod
use outg_mod
use par_mod
use com_mod
implicit none
integer :: jpart,nage,ixave,jyave,kz,kzave,kp
integer :: k,k1,k2,ix,ix1,ix2,ixs,jy,jy1,jy2
real :: xold,yold,zold,xmean,ymean
! Determine average positions
!****************************
if ((ioutputforeachrelease.eq.1).and.(mdomainfill.eq.0)) then
kp=npoint(jpart)
else
kp=1
endif
xmean=(xold+xtra1(jpart))/2.
ymean=(yold+ytra1(jpart))/2.
ixave=int((xmean*dx+xoutshift)/dxout)
jyave=int((ymean*dy+youtshift)/dyout)
do kz=1,numzgrid ! determine height of cell
if (outheight(kz).gt.ztra1(jpart)) goto 16
end do
16 kzave=kz
! Determine vertical fluxes
!**************************
if ((ixave.ge.0).and.(jyave.ge.0).and.(ixave.le.numxgrid-1).and. &
(jyave.le.numygrid-1)) then
do kz=1,numzgrid ! determine height of cell
if (outheighthalf(kz).gt.zold) goto 11
end do
11 k1=min(numzgrid,kz)
do kz=1,numzgrid ! determine height of cell
if (outheighthalf(kz).gt.ztra1(jpart)) goto 21
end do
21 k2=min(numzgrid,kz)
do k=1,nspec
do kz=k1,k2-1
flux(5,ixave,jyave,kz,k,kp,nage)= &
flux(5,ixave,jyave,kz,k,kp,nage)+ &
xmass1(jpart,k)
end do
do kz=k2,k1-1
flux(6,ixave,jyave,kz,k,kp,nage)= &
flux(6,ixave,jyave,kz,k,kp,nage)+ &
xmass1(jpart,k)
end do
end do
endif
! Determine west-east fluxes (fluxw) and east-west fluxes (fluxe)
!****************************************************************
if ((kzave.le.numzgrid).and.(jyave.ge.0).and. &
(jyave.le.numygrid-1)) then
! 1) Particle does not cross domain boundary
if (abs(xold-xtra1(jpart)).lt.real(nx)/2.) then
ix1=int((xold*dx+xoutshift)/dxout+0.5)
ix2=int((xtra1(jpart)*dx+xoutshift)/dxout+0.5)
do k=1,nspec
do ix=ix1,ix2-1
if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
flux(1,ix,jyave,kzave,k,kp,nage)= &
flux(1,ix,jyave,kzave,k,kp,nage) &
+xmass1(jpart,k)
endif
end do
do ix=ix2,ix1-1
if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
flux(2,ix,jyave,kzave,k,kp,nage)= &
flux(2,ix,jyave,kzave,k,kp,nage) &
+xmass1(jpart,k)
endif
end do
end do
! 2) Particle crosses domain boundary: use cyclic boundary condition
! and attribute flux to easternmost grid row only (approximation valid
! for relatively slow motions compared to output grid cell size)
else
ixs=int(((real(nxmin1)-1.e5)*dx+xoutshift)/dxout)
if ((ixs.ge.0).and.(ixs.le.numxgrid-1)) then
if (xold.gt.xtra1(jpart)) then ! west-east flux
do k=1,nspec
flux(1,ixs,jyave,kzave,k,kp,nage)= &
flux(1,ixs,jyave,kzave,k,kp,nage) &
+xmass1(jpart,k)
end do
else ! east-west flux
do k=1,nspec
flux(2,ixs,jyave,kzave,k,kp,nage)= &
flux(2,ixs,jyave,kzave,k,kp,nage) &
+xmass1(jpart,k)
end do
endif
endif
endif
endif
! Determine south-north fluxes (fluxs) and north-south fluxes (fluxn)
!********************************************************************
if ((kzave.le.numzgrid).and.(ixave.ge.0).and. &
(ixave.le.numxgrid-1)) then
jy1=int((yold*dy+youtshift)/dyout+0.5)
jy2=int((ytra1(jpart)*dy+youtshift)/dyout+0.5)
do k=1,nspec
do jy=jy1,jy2-1
if ((jy.ge.0).and.(jy.le.numygrid-1)) then
flux(3,ixave,jy,kzave,k,kp,nage)= &
flux(3,ixave,jy,kzave,k,kp,nage) &
+xmass1(jpart,k)
endif
end do
do jy=jy2,jy1-1
if ((jy.ge.0).and.(jy.le.numygrid-1)) then
flux(4,ixave,jy,kzave,k,kp,nage)= &
flux(4,ixave,jy,kzave,k,kp,nage) &
+xmass1(jpart,k)
endif
end do
end do
endif
end subroutine calcfluxes
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcmatrix(lconv,delt,cbmf)
! o i o
!*****************************************************************************
! *
! This subroutine calculates the matrix describing convective *
! redistribution of mass in a grid column, using the subroutine *
! convect43c.f provided by Kerry Emanuel. *
! *
! Petra Seibert, Bernd C. Krueger, 2000-2001 *
! *
! changed by C. Forster, November 2003 - February 2004 *
! array fmassfrac(nconvlevmax,nconvlevmax) represents *
! the convective redistribution matrix for the particles *
! *
! lconv indicates whether there is convection in this cell, or not *
! delt time step for convection [s] *
! cbmf cloud base mass flux *
! *
!*****************************************************************************
use par_mod
use com_mod
use conv_mod
implicit none
real :: rlevmass,summe
integer :: iflag, k, kk, kuvz
!1-d variables for convection
!variables for redistribution matrix
real :: cbmfold, precip, qprime
real :: tprime, wd, f_qvsat
real :: delt,cbmf
logical :: lconv
lconv = .false.
! calculate pressure at eta levels for use in convect
! and assign temp & spec. hum. to 1D workspace
! -------------------------------------------------------
! pconv(1) is the pressure at the first level above ground
! phconv(k) is the pressure between levels k-1 and k
! dpr(k) is the pressure difference "around" tconv(k)
! phconv(kmax) must also be defined 1/2 level above pconv(kmax)
! Therefore, we define k = kuvz-1 and let kuvz start from 2
! top layer cannot be used for convection because p at top of this layer is
! not given
phconv(1) = psconv
! Emanuel subroutine needs pressure in hPa, therefore convert all pressures
do kuvz = 2,nuvz
k = kuvz-1
pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
dpr(k) = phconv(k) - phconv(kuvz)
qsconv(k) = f_qvsat( pconv(k), tconv(k) )
! initialize mass fractions
do kk=1,nconvlev
fmassfrac(k,kk)=0.
enddo
enddo
!note that Emanuel says it is important
!a. to set this =0. every grid point
!b. to keep this value in the calling programme in the iteration
! CALL CONVECTION
!******************
cbmfold = cbmf
! Convert pressures to hPa, as required by Emanuel scheme
!********************************************************
!!$ do k=1,nconvlev !old
do k=1,nconvlev+1 !bugfix
pconv_hpa(k)=pconv(k)/100.
phconv_hpa(k)=phconv(k)/100.
end do
phconv_hpa(nconvlev+1)=phconv(nconvlev+1)/100.
call convect(nconvlevmax, nconvlev, delt, iflag, &
precip, wd, tprime, qprime, cbmf)
! do not update fmassfrac and cloudbase massflux
! if no convection takes place or
! if a CFL criterion is violated in convect43c.f
if (iflag .ne. 1 .and. iflag .ne. 4) then
cbmf=cbmfold
goto 200
endif
! do not update fmassfrac and cloudbase massflux
! if the old and the new cloud base mass
! fluxes are zero
if (cbmf.le.0..and.cbmfold.le.0.) then
cbmf=cbmfold
goto 200
endif
! Update fmassfrac
! account for mass displaced from level k to level k
lconv = .true.
do k=1,nconvtop
rlevmass = dpr(k)/ga
summe = 0.
do kk=1,nconvtop
fmassfrac(k,kk) = delt*fmass(k,kk)
summe = summe + fmassfrac(k,kk)
end do
fmassfrac(k,k)=fmassfrac(k,k) + rlevmass - summe
end do
200 continue
end subroutine calcmatrix
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcmatrix(lconv,delt,cbmf)
! o i o
!*****************************************************************************
! *
! This subroutine calculates the matrix describing convective *
! redistribution of mass in a grid column, using the subroutine *
! convect43c.f provided by Kerry Emanuel. *
! *
! Petra Seibert, Bernd C. Krueger, 2000-2001 *
! *
! changed by C. Forster, November 2003 - February 2004 *
! array fmassfrac(nconvlevmax,nconvlevmax) represents *
! the convective redistribution matrix for the particles *
! *
! Changes by C. Forster, November 2005, NCEP GFS version *
! *
! lconv indicates whether there is convection in this cell, or not *
! delt time step for convection [s] *
! cbmf cloud base mass flux *
! *
!*****************************************************************************
use par_mod
use com_mod
use conv_mod
implicit none
real :: rlevmass,summe
integer :: iflag, k, kk, kuvz
!1-d variables for convection
!variables for redistribution matrix
real :: cbmfold, precip, qprime
real :: tprime, wd, f_qvsat
real :: delt,cbmf
logical :: lconv
lconv = .false.
! calculate pressure at eta levels for use in convect
! and assign temp & spec. hum. to 1D workspace
! -------------------------------------------------------
! pconv(1) is the pressure at the first level above ground
! phconv(k) is the pressure between levels k-1 and k
! dpr(k) is the pressure difference "around" tconv(k)
! phconv(kmax) must also be defined 1/2 level above pconv(kmax)
! Therefore, we define k = kuvz-1 and let kuvz start from 2
! top layer cannot be used for convection because p at top of this layer is
! not given
phconv(1) = psconv
do kuvz = 2,nuvz
k = kuvz-1
phconv(kuvz) = 0.5*(pconv(kuvz)+pconv(k))
dpr(k) = phconv(k) - phconv(kuvz)
qsconv(k) = f_qvsat( pconv(k), tconv(k) )
! initialize mass fractions
do kk=1,nconvlev
fmassfrac(k,kk)=0.
enddo
end do
!note that Emanuel says it is important
!a. to set this =0. every grid point
!b. to keep this value in the calling programme in the iteration
! CALL CONVECTION
!******************
cbmfold = cbmf
! Convert pressures to hPa, as required by Emanuel scheme
!********************************************************
!!$ do k=1,nconvlev !old
do k=1,nconvlev+1 !bugfix
pconv_hpa(k)=pconv(k)/100.
phconv_hpa(k)=phconv(k)/100.
end do
phconv_hpa(nconvlev+1)=phconv(nconvlev+1)/100.
call convect(nconvlevmax, nconvlev, delt, iflag, &
precip, wd, tprime, qprime, cbmf)
! do not update fmassfrac and cloudbase massflux
! if no convection takes place or
! if a CFL criterion is violated in convect43c.f
if (iflag .ne. 1 .and. iflag .ne. 4) then
cbmf=cbmfold
goto 200
endif
! do not update fmassfrac and cloudbase massflux
! if the old and the new cloud base mass
! fluxes are zero
if (cbmf.le.0..and.cbmfold.le.0.) then
cbmf=cbmfold
goto 200
endif
! Update fmassfrac
! account for mass displaced from level k to level k
lconv = .true.
do k=1,nconvtop
rlevmass = dpr(k)/ga
summe = 0.
do kk=1,nconvtop
fmassfrac(k,kk) = delt*fmass(k,kk)
summe = summe + fmassfrac(k,kk)
end do
fmassfrac(k,k)=fmassfrac(k,k) + rlevmass - summe
end do
200 continue
end subroutine calcmatrix
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcpar(n,uuh,vvh,pvh)
! i i i o
!*****************************************************************************
! *
! Computation of several boundary layer parameters needed for the *
! dispersion calculation and calculation of dry deposition velocities. *
! All parameters are calculated over the entire grid. *
! *
! Author: A. Stohl *
! *
! 21 May 1995 *
! *
! ------------------------------------------------------------------ *
! Petra Seibert, Feb 2000: *
! convection scheme: *
! new variables in call to richardson *
! *
!*****************************************************************************
! Changes, Bernd C. Krueger, Feb. 2001:
! Variables tth and qvh (on eta coordinates) in common block
!*****************************************************************************
! *
! Variables: *
! n temporal index for meteorological fields (1 to 3) *
! *
! Constants: *
! *
! *
! Functions: *
! scalev computation of ustar *
! obukhov computatio of Obukhov length *
! *
!*****************************************************************************
use par_mod
use com_mod
implicit none
integer :: n,ix,jy,i,kz,lz,kzmin
real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
real,parameter :: const=r_air/ga
!write(*,*) 'in calcpar writting snowheight'
!***********************************
! for test: write out snow depths
! open(4,file='slandusetest',form='formatted')
! do 5 ix=0,nxmin1
!5 write (4,*) (sd(ix,jy,1,n),jy=0,nymin1)
! close(4)
! Loop over entire grid
!**********************
do jy=0,nymin1
! Set minimum height for tropopause
!**********************************
ylat=ylat0+real(jy)*dy
if ((ylat.ge.-20.).and.(ylat.le.20.)) then
altmin = 5000.
else
if ((ylat.gt.20.).and.(ylat.lt.40.)) then
altmin=2500.+(40.-ylat)*125.
else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
altmin=2500.+(40.+ylat)*125.
else
altmin=2500.
endif
endif
do ix=0,nxmin1
! 1) Calculation of friction velocity
!************************************
ustar(ix,jy,1,n)=scalev(ps(ix,jy,1,n),tt2(ix,jy,1,n), &
td2(ix,jy,1,n),surfstr(ix,jy,1,n))
if (ustar(ix,jy,1,n).le.1.e-8) ustar(ix,jy,1,n)=1.e-8
! 2) Calculation of inverse Obukhov length scale
!***********************************************
ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm)
if (ol.ne.0.) then
oli(ix,jy,1,n)=1./ol
else
oli(ix,jy,1,n)=99999.
endif
! 3) Calculation of convective velocity scale and mixing height
!**************************************************************
do i=1,nuvz
ulev(i)=uuh(ix,jy,i)
vlev(i)=vvh(ix,jy,i)
ttlev(i)=tth(ix,jy,i,n)
qvlev(i)=qvh(ix,jy,i,n)
end do
call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus)
if(lsubgrid.eq.1) then
subsceff=min(excessoro(ix,jy),hmixplus)
else
subsceff=0.0
endif
!
! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
!
hmix(ix,jy,1,n)=hmix(ix,jy,1,n)+subsceff
hmix(ix,jy,1,n)=max(hmixmin,hmix(ix,jy,1,n)) ! set minimum PBL height
hmix(ix,jy,1,n)=min(hmixmax,hmix(ix,jy,1,n)) ! set maximum PBL height
! 4) Calculation of dry deposition velocities
!********************************************
if (DRYDEP) then
! Sabine Eckhardt, Dec 06: use new index for z0 for water depending on
! windspeed
z0(7)=0.016*ustar(ix,jy,1,n)*ustar(ix,jy,1,n)/ga
! Calculate relative humidity at surface
!***************************************
rh=ew(td2(ix,jy,1,n))/ew(tt2(ix,jy,1,n))
call getvdep(n,ix,jy,ustar(ix,jy,1,n), &
tt2(ix,jy,1,n),ps(ix,jy,1,n),1./oli(ix,jy,1,n), &
ssr(ix,jy,1,n),rh,lsprec(ix,jy,1,n)+convprec(ix,jy,1,n), &
sd(ix,jy,1,n),vd)
do i=1,nspec
vdep(ix,jy,i,n)=vd(i)
end do
endif
!******************************************************
! Calculate height of thermal tropopause (Hoinka, 1997)
!******************************************************
! 1) Calculate altitudes of ECMWF model levels
!*********************************************
tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
ps(ix,jy,1,n))
pold=ps(ix,jy,1,n)
zold=0.
do kz=2,nuvz
pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n) ! pressure on model layers
tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
if (abs(tv-tvold).gt.0.2) then
zlev(kz)=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
else
zlev(kz)=zold+const*log(pold/pint)*tv
endif
tvold=tv
pold=pint
zold=zlev(kz)
end do
! 2) Define a minimum level kzmin, from which upward the tropopause is
! searched for. This is to avoid inversions in the lower troposphere
! to be identified as the tropopause
!************************************************************************
do kz=1,nuvz
if (zlev(kz).ge.altmin) then
kzmin=kz
goto 45
endif
end do
45 continue
! 3) Search for first stable layer above minimum height that fulfills the
! thermal tropopause criterion
!************************************************************************
do kz=kzmin,nuvz
do lz=kz+1,nuvz
if ((zlev(lz)-zlev(kz)).gt.2000.) then
if (((tth(ix,jy,kz,n)-tth(ix,jy,lz,n))/ &
(zlev(lz)-zlev(kz))).lt.0.002) then
tropopause(ix,jy,1,n)=zlev(kz)
goto 51
endif
goto 50
endif
end do
50 continue
end do
51 continue
end do
end do
! Calculation of potential vorticity on 3-d grid
!***********************************************
call calcpv(n,uuh,vvh,pvh)
end subroutine calcpar
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcpar(n,uuh,vvh,pvh)
! i i i o
!*****************************************************************************
! *
! Computation of several boundary layer parameters needed for the *
! dispersion calculation and calculation of dry deposition velocities. *
! All parameters are calculated over the entire grid. *
! *
! Author: A. Stohl *
! *
! 21 May 1995 *
! *
! ------------------------------------------------------------------ *
! Petra Seibert, Feb 2000: *
! convection scheme: *
! new variables in call to richardson *
! *
!*****************************************************************************
! Changes, Bernd C. Krueger, Feb. 2001:
! Variables tth and qvh (on eta coordinates) in common block
!*****************************************************************************
! *
! CHANGE 17/11/2005 Caroline Forster NCEP GFS version *
! *
! Variables: *
! n temporal index for meteorological fields (1 to 3) *
! *
! Constants: *
! *
! *
! Functions: *
! scalev computation of ustar *
! obukhov computatio of Obukhov length *
! *
!*****************************************************************************
use par_mod
use com_mod
implicit none
integer :: n,ix,jy,i,kz,lz,kzmin,llev
real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax),hmixdummy
real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
real,parameter :: const=r_air/ga
! Loop over entire grid
!**********************
do jy=0,nymin1
! Set minimum height for tropopause
!**********************************
ylat=ylat0+real(jy)*dy
if ((ylat.ge.-20.).and.(ylat.le.20.)) then
altmin = 5000.
else
if ((ylat.gt.20.).and.(ylat.lt.40.)) then
altmin=2500.+(40.-ylat)*125.
else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
altmin=2500.+(40.+ylat)*125.
else
altmin=2500.
endif
endif
do ix=0,nxmin1
! 1) Calculation of friction velocity
!************************************
ustar(ix,jy,1,n)=scalev(ps(ix,jy,1,n),tt2(ix,jy,1,n), &
td2(ix,jy,1,n),surfstr(ix,jy,1,n))
if (ustar(ix,jy,1,n).le.1.e-8) ustar(ix,jy,1,n)=1.e-8
! 2) Calculation of inverse Obukhov length scale
!***********************************************
! NCEP version: find first level above ground
llev = 0
do i=1,nuvz
if (ps(ix,jy,1,n).lt.akz(i)) llev=i
end do
llev = llev+1
if (llev.gt.nuvz) llev = nuvz-1
! NCEP version
! calculate inverse Obukhov length scale with tth(llev)
ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akz(llev))
if (ol.ne.0.) then
oli(ix,jy,1,n)=1./ol
else
oli(ix,jy,1,n)=99999.
endif
! 3) Calculation of convective velocity scale and mixing height
!**************************************************************
do i=1,nuvz
ulev(i)=uuh(ix,jy,i)
vlev(i)=vvh(ix,jy,i)
ttlev(i)=tth(ix,jy,i,n)
qvlev(i)=qvh(ix,jy,i,n)
end do
! NCEP version hmix has been read in in readwind.f, is therefore not calculated here
call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus)
if(lsubgrid.eq.1) then
subsceff=min(excessoro(ix,jy),hmixplus)
else
subsceff=0
endif
!
! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
!
hmix(ix,jy,1,n)=hmix(ix,jy,1,n)+subsceff
hmix(ix,jy,1,n)=max(hmixmin,hmix(ix,jy,1,n)) ! set minimum PBL height
hmix(ix,jy,1,n)=min(hmixmax,hmix(ix,jy,1,n)) ! set maximum PBL height
! 4) Calculation of dry deposition velocities
!********************************************
if (DRYDEP) then
! Sabine Eckhardt, Dec 06: use new index for z0 for water depending on
! windspeed
z0(7)=0.016*ustar(ix,jy,1,n)*ustar(ix,jy,1,n)/ga
! Calculate relative humidity at surface
!***************************************
rh=ew(td2(ix,jy,1,n))/ew(tt2(ix,jy,1,n))
call getvdep(n,ix,jy,ustar(ix,jy,1,n), &
tt2(ix,jy,1,n),ps(ix,jy,1,n),1./oli(ix,jy,1,n), &
ssr(ix,jy,1,n),rh,lsprec(ix,jy,1,n)+convprec(ix,jy,1,n), &
sd(ix,jy,1,n),vd)
do i=1,nspec
vdep(ix,jy,i,n)=vd(i)
end do
endif
!******************************************************
! Calculate height of thermal tropopause (Hoinka, 1997)
!******************************************************
! 1) Calculate altitudes of NCEP model levels
!*********************************************
tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
ps(ix,jy,1,n))
pold=ps(ix,jy,1,n)
zold=0.
do kz=llev,nuvz
pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n) ! pressure on model layers
tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
if (abs(tv-tvold).gt.0.2) then
zlev(kz)=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
else
zlev(kz)=zold+const*log(pold/pint)*tv
endif
tvold=tv
pold=pint
zold=zlev(kz)
end do
! 2) Define a minimum level kzmin, from which upward the tropopause is
! searched for. This is to avoid inversions in the lower troposphere
! to be identified as the tropopause
!************************************************************************
do kz=llev,nuvz
if (zlev(kz).ge.altmin) then
kzmin=kz
goto 45
endif
end do
45 continue
! 3) Search for first stable layer above minimum height that fulfills the
! thermal tropopause criterion
!************************************************************************
do kz=kzmin,nuvz
do lz=kz+1,nuvz
if ((zlev(lz)-zlev(kz)).gt.2000.) then
if (((tth(ix,jy,kz,n)-tth(ix,jy,lz,n))/ &
(zlev(lz)-zlev(kz))).lt.0.002) then
tropopause(ix,jy,1,n)=zlev(kz)
goto 51
endif
goto 50
endif
end do
50 continue
end do
51 continue
end do
end do
! Calculation of potential vorticity on 3-d grid
!***********************************************
call calcpv(n,uuh,vvh,pvh)
end subroutine calcpar
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcpar_nests(n,uuhn,vvhn,pvhn)
! i i i o
!*****************************************************************************
! *
! Computation of several boundary layer parameters needed for the *
! dispersion calculation and calculation of dry deposition velocities. *
! All parameters are calculated over the entire grid. *
! This routine is similar to calcpar, but is used for the nested grids. *
! *
! Author: A. Stohl *
! *
! 8 February 1999 *
! *
! ------------------------------------------------------------------ *
! Petra Seibert, Feb 2000: *
! convection scheme: *
! new variables in call to richardson *
! *
!*****************************************************************************
! Changes, Bernd C. Krueger, Feb. 2001:
! Variables tth and qvh (on eta coordinates) in common block
!*****************************************************************************
! *
! Variables: *
! n temporal index for meteorological fields (1 to 3) *
! *
! Constants: *
! *
! *
! Functions: *
! scalev computation of ustar *
! obukhov computatio of Obukhov length *
! *
!*****************************************************************************
use par_mod
use com_mod
implicit none
integer :: n,ix,jy,i,l,kz,lz,kzmin
real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real,parameter :: const=r_air/ga
! Loop over all nests
!********************
do l=1,numbnests
! Loop over entire grid
!**********************
do jy=0,nyn(l)-1
! Set minimum height for tropopause
!**********************************
ylat=ylat0n(l)+real(jy)*dyn(l)
if ((ylat.ge.-20.).and.(ylat.le.20.)) then
altmin = 5000.
else
if ((ylat.gt.20.).and.(ylat.lt.40.)) then
altmin=2500.+(40.-ylat)*125.
else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
altmin=2500.+(40.+ylat)*125.
else
altmin=2500.
endif
endif
do ix=0,nxn(l)-1
! 1) Calculation of friction velocity
!************************************
ustarn(ix,jy,1,n,l)=scalev(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
td2n(ix,jy,1,n,l),surfstrn(ix,jy,1,n,l))
! 2) Calculation of inverse Obukhov length scale
!***********************************************
ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), &
sshfn(ix,jy,1,n,l),akm,bkm)
if (ol.ne.0.) then
olin(ix,jy,1,n,l)=1./ol
else
olin(ix,jy,1,n,l)=99999.
endif
! 3) Calculation of convective velocity scale and mixing height
!**************************************************************
do i=1,nuvz
ulev(i)=uuhn(ix,jy,i,l)
vlev(i)=vvhn(ix,jy,i,l)
ttlev(i)=tthn(ix,jy,i,n,l)
qvlev(i)=qvhn(ix,jy,i,n,l)
end do
call richardson(psn(ix,jy,1,n,l),ustarn(ix,jy,1,n,l),ttlev, &
qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l), &
tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), &
wstarn(ix,jy,1,n,l),hmixplus)
if(lsubgrid.eq.1) then
subsceff=min(excessoron(ix,jy,l),hmixplus)
else
subsceff=0
endif
!
! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
!
hmixn(ix,jy,1,n,l)=hmixn(ix,jy,1,n,l)+subsceff
hmixn(ix,jy,1,n,l)=max(hmixmin,hmixn(ix,jy,1,n,l)) ! minim PBL height
hmixn(ix,jy,1,n,l)=min(hmixmax,hmixn(ix,jy,1,n,l)) ! maxim PBL height
! 4) Calculation of dry deposition velocities
!********************************************
if (DRYDEP) then
z0(4)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
z0(9)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
! Calculate relative humidity at surface
!***************************************
rh=ew(td2n(ix,jy,1,n,l))/ew(tt2n(ix,jy,1,n,l))
call getvdep_nests(n,ix,jy,ustarn(ix,jy,1,n,l), &
tt2n(ix,jy,1,n,l),psn(ix,jy,1,n,l),1./olin(ix,jy,1,n,l), &
ssrn(ix,jy,1,n,l),rh,lsprecn(ix,jy,1,n,l)+ &
convprecn(ix,jy,1,n,l),sdn(ix,jy,1,n,l),vd,l)
do i=1,nspec
vdepn(ix,jy,i,n,l)=vd(i)
end do
endif
!******************************************************
! Calculate height of thermal tropopause (Hoinka, 1997)
!******************************************************
! 1) Calculate altitudes of ECMWF model levels
!*********************************************
tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
psn(ix,jy,1,n,l))
pold=psn(ix,jy,1,n,l)
zold=0.
do kz=2,nuvz
pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l) ! pressure on model layers
tv=tthn(ix,jy,kz,n,l)*(1.+0.608*qvhn(ix,jy,kz,n,l))
if (abs(tv-tvold).gt.0.2) then
zlev(kz)=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
else
zlev(kz)=zold+const*log(pold/pint)*tv
endif
tvold=tv
pold=pint
zold=zlev(kz)
end do
! 2) Define a minimum level kzmin, from which upward the tropopause is
! searched for. This is to avoid inversions in the lower troposphere
! to be identified as the tropopause
!************************************************************************
do kz=1,nuvz
if (zlev(kz).ge.altmin) then
kzmin=kz
goto 45
endif
end do
45 continue
! 3) Search for first stable layer above minimum height that fulfills the
! thermal tropopause criterion
!************************************************************************
do kz=kzmin,nuvz
do lz=kz+1,nuvz
if ((zlev(lz)-zlev(kz)).gt.2000.) then
if (((tthn(ix,jy,kz,n,l)-tthn(ix,jy,lz,n,l))/ &
(zlev(lz)-zlev(kz))).lt.0.002) then
tropopausen(ix,jy,1,n,l)=zlev(kz)
goto 51
endif
goto 50
endif
end do
50 continue
end do
51 continue
end do
end do
call calcpv_nests(l,n,uuhn,vvhn,pvhn)
end do
end subroutine calcpar_nests
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcpv(n,uuh,vvh,pvh)
! i i i o
!*****************************************************************************
! *
! Calculation of potential vorticity on 3-d grid. *
! *
! Author: P. James *
! 3 February 2000 *
! *
! Adaptation to FLEXPART, A. Stohl, 1 May 2000 *
! *
!*****************************************************************************
! *
! Variables: *
! n temporal index for meteorological fields (1 to 2) *
! *
! Constants: *
! *
!*****************************************************************************
use par_mod
use com_mod
implicit none
integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
integer :: nlck
real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
real :: pvavr,ppml(nuvzmax)
real :: thup,thdn
real,parameter :: eps=1.e-5, p0=101325
real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
! Set number of levels to check for adjacent theta
nlck=nuvz/3
!
! Loop over entire grid
!**********************
do jy=0,nymin1
if (sglobal.and.jy.eq.0) goto 10
if (nglobal.and.jy.eq.nymin1) goto 10
phi = (ylat0 + jy * dy) * pi / 180.
f = 0.00014585 * sin(phi)
tanphi = tan(phi)
cosphi = cos(phi)
! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
jyvp=jy+1
jyvm=jy-1
if (jy.eq.0) jyvm=0
if (jy.eq.nymin1) jyvp=nymin1
! Define absolute gap length
jumpy=2
if (jy.eq.0.or.jy.eq.nymin1) jumpy=1
if (sglobal.and.jy.eq.1) then
jyvm=1
jumpy=1
end if
if (nglobal.and.jy.eq.ny-2) then
jyvp=ny-2
jumpy=1
end if
juy=jumpy
!
do ix=0,nxmin1
! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
ixvp=ix+1
ixvm=ix-1
jumpx=2
if (xglobal) then
ivrp=ixvp
ivrm=ixvm
if (ixvm.lt.0) ivrm=ixvm+nxmin1
if (ixvp.ge.nx) ivrp=ixvp-nx+1
else
if (ix.eq.0) ixvm=0
if (ix.eq.nxmin1) ixvp=nxmin1
ivrp=ixvp
ivrm=ixvm
! Define absolute gap length
if (ix.eq.0.or.ix.eq.nxmin1) jumpx=1
end if
jux=jumpx
! Precalculate pressure values for efficiency
do kl=1,nuvz
ppml(kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
end do
!
! Loop over the vertical
!***********************
do kl=1,nuvz
ppmk=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
theta=tth(ix,jy,kl,n)*(100000./ppmk)**kappa
klvrp=kl+1
klvrm=kl-1
klpt=kl
! If top or bottom level, dthetadp is evaluated between the current
! level and the level inside, otherwise between level+1 and level-1
!
if (klvrp.gt.nuvz) klvrp=nuvz
if (klvrm.lt.1) klvrm=1
ppmk=akz(klvrp)+bkz(klvrp)*ps(ix,jy,1,n)
thetap=tth(ix,jy,klvrp,n)*(100000./ppmk)**kappa
ppmk=akz(klvrm)+bkz(klvrm)*ps(ix,jy,1,n)
thetam=tth(ix,jy,klvrm,n)*(100000./ppmk)**kappa
dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
! Compute vertical position at pot. temperature surface on subgrid
! and the wind at that position
!*****************************************************************
! a) in x direction
ii=0
do i=ixvm,ixvp,jumpx
ivr=i
if (xglobal) then
if (i.lt.0) ivr=ivr+nxmin1
if (i.ge.nx) ivr=ivr-nx+1
end if
ii=ii+1
! Search adjacent levels for current theta value
! Spiral out from current level for efficiency
kup=klpt-1
kdn=klpt
kch=0
40 continue
! Upward branch
kup=kup+1
if (kch.ge.nlck) goto 21 ! No more levels to check,
! ! and no values found
if (kup.ge.nuvz) goto 41
kch=kch+1
k=kup
ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
if (((thdn.ge.theta).and.(thup.le.theta)).or. &
((thdn.le.theta).and.(thup.ge.theta))) then
dt1=abs(theta-thdn)
dt2=abs(theta-thup)
dt=dt1+dt2
if (dt.lt.eps) then ! Avoid division by zero error
dt1=0.5 ! G.W., 10.4.1996
dt2=0.5
dt=1.0
endif
vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
goto 20
endif
41 continue
! Downward branch
kdn=kdn-1
if (kdn.lt.1) goto 40
kch=kch+1
k=kdn
ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
if (((thdn.ge.theta).and.(thup.le.theta)).or. &
((thdn.le.theta).and.(thup.ge.theta))) then
dt1=abs(theta-thdn)
dt2=abs(theta-thup)
dt=dt1+dt2
if (dt.lt.eps) then ! Avoid division by zero error
dt1=0.5 ! G.W., 10.4.1996
dt2=0.5
dt=1.0
endif
vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
goto 20
endif
goto 40
! This section used when no values were found
21 continue
! Must use vv at current level and long. jux becomes smaller by 1
vx(ii)=vvh(ix,jy,kl)
jux=jux-1
! Otherwise OK
20 continue
end do
if (jux.gt.0) then
dvdx=(vx(2)-vx(1))/real(jux)/(dx*pi/180.)
else
dvdx=vvh(ivrp,jy,kl)-vvh(ivrm,jy,kl)
dvdx=dvdx/real(jumpx)/(dx*pi/180.)
! Only happens if no equivalent theta value
! can be found on either side, hence must use values
! from either side, same pressure level.
end if
! b) in y direction
jj=0
do j=jyvm,jyvp,jumpy
jj=jj+1
! Search adjacent levels for current theta value
! Spiral out from current level for efficiency
kup=klpt-1
kdn=klpt
kch=0
70 continue
! Upward branch
kup=kup+1
if (kch.ge.nlck) goto 51 ! No more levels to check,
! ! and no values found
if (kup.ge.nuvz) goto 71
kch=kch+1
k=kup
ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
if (((thdn.ge.theta).and.(thup.le.theta)).or. &
((thdn.le.theta).and.(thup.ge.theta))) then
dt1=abs(theta-thdn)
dt2=abs(theta-thup)
dt=dt1+dt2
if (dt.lt.eps) then ! Avoid division by zero error
dt1=0.5 ! G.W., 10.4.1996
dt2=0.5
dt=1.0
endif
uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
goto 50
endif
71 continue
! Downward branch
kdn=kdn-1
if (kdn.lt.1) goto 70
kch=kch+1
k=kdn
ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
if (((thdn.ge.theta).and.(thup.le.theta)).or. &
((thdn.le.theta).and.(thup.ge.theta))) then
dt1=abs(theta-thdn)
dt2=abs(theta-thup)
dt=dt1+dt2
if (dt.lt.eps) then ! Avoid division by zero error
dt1=0.5 ! G.W., 10.4.1996
dt2=0.5
dt=1.0
endif
uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
goto 50
endif
goto 70
! This section used when no values were found
51 continue
! Must use uu at current level and lat. juy becomes smaller by 1
uy(jj)=uuh(ix,jy,kl)
juy=juy-1
! Otherwise OK
50 continue
end do
if (juy.gt.0) then
dudy=(uy(2)-uy(1))/real(juy)/(dy*pi/180.)
else
dudy=uuh(ix,jyvp,kl)-uuh(ix,jyvm,kl)
dudy=dudy/real(jumpy)/(dy*pi/180.)
end if
!
pvh(ix,jy,kl)=dthetadp*(f+(dvdx/cosphi-dudy &
+uuh(ix,jy,kl)*tanphi)/r_earth)*(-1.e6)*9.81
!
! Resest jux and juy
jux=jumpx
juy=jumpy
end do
end do
10 continue
end do
!
! Fill in missing PV values on poles, if present
! Use mean PV of surrounding latitude ring
!
if (sglobal) then
do kl=1,nuvz
pvavr=0.
do ix=0,nxmin1
pvavr=pvavr+pvh(ix,1,kl)
end do
pvavr=pvavr/real(nx)
jy=0
do ix=0,nxmin1
pvh(ix,jy,kl)=pvavr
end do
end do
end if
if (nglobal) then
do kl=1,nuvz
pvavr=0.
do ix=0,nxmin1
pvavr=pvavr+pvh(ix,ny-2,kl)
end do
pvavr=pvavr/real(nx)
jy=nymin1
do ix=0,nxmin1
pvh(ix,jy,kl)=pvavr
end do
end do
end if
end subroutine calcpv
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcpv_nests(l,n,uuhn,vvhn,pvhn)
! i i i i o
!*****************************************************************************
! *
! Calculation of potential vorticity on 3-d nested grid *
! *
! Author: P. James *
! 22 February 2000 *
! *
!*****************************************************************************
! *
! Variables: *
! n temporal index for meteorological fields (1 to 2) *
! l index of current nest *
! *
! Constants: *
! *
!*****************************************************************************
use par_mod
use com_mod
implicit none
integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
integer :: nlck,l
real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
real :: ppml(nuvzmax)
real :: thup,thdn
real,parameter :: eps=1.e-5,p0=101325
real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
! Set number of levels to check for adjacent theta
nlck=nuvz/3
!
! Loop over entire grid
!**********************
do jy=0,nyn(l)-1
phi = (ylat0n(l) + jy * dyn(l)) * pi / 180.
f = 0.00014585 * sin(phi)
tanphi = tan(phi)
cosphi = cos(phi)
! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
jyvp=jy+1
jyvm=jy-1
if (jy.eq.0) jyvm=0
if (jy.eq.nyn(l)-1) jyvp=nyn(l)-1
! Define absolute gap length
jumpy=2
if (jy.eq.0.or.jy.eq.nyn(l)-1) jumpy=1
juy=jumpy
!
do ix=0,nxn(l)-1
! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
ixvp=ix+1
ixvm=ix-1
jumpx=2
if (ix.eq.0) ixvm=0
if (ix.eq.nxn(l)-1) ixvp=nxn(l)-1
ivrp=ixvp
ivrm=ixvm
! Define absolute gap length
if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1
jux=jumpx
! Precalculate pressure values for efficiency
do kl=1,nuvz
ppml(kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
end do
!
! Loop over the vertical
!***********************
do kl=1,nuvz
ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
theta=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa
klvrp=kl+1
klvrm=kl-1
klpt=kl
! If top or bottom level, dthetadp is evaluated between the current
! level and the level inside, otherwise between level+1 and level-1
!
if (klvrp.gt.nuvz) klvrp=nuvz
if (klvrm.lt.1) klvrm=1
ppmk=akz(klvrp)+bkz(klvrp)*psn(ix,jy,1,n,l)
thetap=tthn(ix,jy,klvrp,n,l)*(100000./ppmk)**kappa
ppmk=akz(klvrm)+bkz(klvrm)*psn(ix,jy,1,n,l)
thetam=tthn(ix,jy,klvrm,n,l)*(100000./ppmk)**kappa
dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
! Compute vertical position at pot. temperature surface on subgrid
! and the wind at that position
!*****************************************************************
! a) in x direction
ii=0
do i=ixvm,ixvp,jumpx
ivr=i
ii=ii+1
! Search adjacent levels for current theta value
! Spiral out from current level for efficiency
kup=klpt-1
kdn=klpt
kch=0
40 continue
! Upward branch
kup=kup+1
if (kch.ge.nlck) goto 21 ! No more levels to check,
! ! and no values found
if (kup.ge.nuvz) goto 41
kch=kch+1
k=kup
ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa
if (((thdn.ge.theta).and.(thup.le.theta)).or. &
((thdn.le.theta).and.(thup.ge.theta))) then
dt1=abs(theta-thdn)
dt2=abs(theta-thup)
dt=dt1+dt2
if (dt.lt.eps) then ! Avoid division by zero error
dt1=0.5 ! G.W., 10.4.1996
dt2=0.5
dt=1.0
endif
vx(ii)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt
goto 20
endif
41 continue
! Downward branch
kdn=kdn-1
if (kdn.lt.1) goto 40
kch=kch+1
k=kdn
ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa
if (((thdn.ge.theta).and.(thup.le.theta)).or. &
((thdn.le.theta).and.(thup.ge.theta))) then
dt1=abs(theta-thdn)
dt2=abs(theta-thup)
dt=dt1+dt2
if (dt.lt.eps) then ! Avoid division by zero error
dt1=0.5 ! G.W., 10.4.1996
dt2=0.5
dt=1.0
endif
vx(ii)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt
goto 20
endif
goto 40
! This section used when no values were found
21 continue
! Must use vv at current level and long. jux becomes smaller by 1
vx(ii)=vvhn(ix,jy,kl,l)
jux=jux-1
! Otherwise OK
20 continue
end do
if (jux.gt.0) then
dvdx=(vx(2)-vx(1))/real(jux)/(dxn(l)*pi/180.)
else
dvdx=vvhn(ivrp,jy,kl,l)-vvhn(ivrm,jy,kl,l)
dvdx=dvdx/real(jumpx)/(dxn(l)*pi/180.)
! Only happens if no equivalent theta value
! can be found on either side, hence must use values
! from either side, same pressure level.
end if
! b) in y direction
jj=0
do j=jyvm,jyvp,jumpy
jj=jj+1
! Search adjacent levels for current theta value
! Spiral out from current level for efficiency
kup=klpt-1
kdn=klpt
kch=0
70 continue
! Upward branch
kup=kup+1
if (kch.ge.nlck) goto 51 ! No more levels to check,
! ! and no values found
if (kup.ge.nuvz) goto 71
kch=kch+1
k=kup
ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
if (((thdn.ge.theta).and.(thup.le.theta)).or. &
((thdn.le.theta).and.(thup.ge.theta))) then
dt1=abs(theta-thdn)
dt2=abs(theta-thup)
dt=dt1+dt2
if (dt.lt.eps) then ! Avoid division by zero error
dt1=0.5 ! G.W., 10.4.1996
dt2=0.5
dt=1.0
endif
uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt
goto 50
endif
71 continue
! Downward branch
kdn=kdn-1
if (kdn.lt.1) goto 70
kch=kch+1
k=kdn
ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
if (((thdn.ge.theta).and.(thup.le.theta)).or. &
((thdn.le.theta).and.(thup.ge.theta))) then
dt1=abs(theta-thdn)
dt2=abs(theta-thup)
dt=dt1+dt2
if (dt.lt.eps) then ! Avoid division by zero error
dt1=0.5 ! G.W., 10.4.1996
dt2=0.5
dt=1.0
endif
uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt
goto 50
endif
goto 70
! This section used when no values were found
51 continue
! Must use uu at current level and lat. juy becomes smaller by 1
uy(jj)=uuhn(ix,jy,kl,l)
juy=juy-1
! Otherwise OK
50 continue
end do
if (juy.gt.0) then
dudy=(uy(2)-uy(1))/real(juy)/(dyn(l)*pi/180.)
else
dudy=uuhn(ix,jyvp,kl,l)-uuhn(ix,jyvm,kl,l)
dudy=dudy/real(jumpy)/(dyn(l)*pi/180.)
end if
!
pvhn(ix,jy,kl,l)=dthetadp*(f+(dvdx/cosphi-dudy &
+uuhn(ix,jy,kl,l)*tanphi)/r_earth)*(-1.e6)*9.81
!
! Resest jux and juy
jux=jumpx
juy=jumpy
end do
end do
end do
!
end subroutine calcpv_nests
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine caldate(juldate,yyyymmdd,hhmiss)
! i o o
!*****************************************************************************
! *
! Calculates the Gregorian date from the Julian date *
! *
! AUTHOR: Andreas Stohl (21 January 1994), adapted from Numerical Recipes*
! *
! Variables: *
! dd Day *
! hh Hour *
! hhmiss Hour, Minute, Second *
! ja,jb,jc,jd,je help variables *
! jalpha help variable *
! juldate Julian Date *
! julday help variable *
! mi Minute *
! mm Month *
! ss Seconds *
! yyyy Year *
! yyyymmdd Year, Month, Day *
! *
! Constants: *
! igreg help constant *
! *
!*****************************************************************************
use par_mod, only: dp
implicit none
integer :: yyyymmdd,yyyy,mm,dd,hhmiss,hh,mi,ss
integer :: julday,ja,jb,jc,jd,je,jalpha
real(kind=dp) :: juldate
integer,parameter :: igreg=2299161
julday=int(juldate)
if(julday.ge.igreg)then
jalpha=int(((julday-1867216)-0.25)/36524.25)
ja=julday+1+jalpha-int(0.25*jalpha)
else
ja=julday
endif
jb=ja+1524
jc=int(6680.+((jb-2439870)-122.1)/365.25)
jd=365*jc+int(0.25*jc)
je=int((jb-jd)/30.6001)
dd=jb-jd-int(30.6001*je)
mm=je-1
if (mm.gt.12) mm=mm-12
yyyy=jc-4715
if (mm.gt.2) yyyy=yyyy-1
if (yyyy.le.0) yyyy=yyyy-1
yyyymmdd=10000*yyyy+100*mm+dd
hh=int(24._dp*(juldate-real(julday,kind=dp)))
mi=int(1440._dp*(juldate-real(julday,kind=dp))-60._dp*real(hh,kind=dp))
ss=nint(86400._dp*(juldate-real(julday,kind=dp))-3600._dp*real(hh,kind=dp)- &
60._dp*real(mi,kind=dp))
if (ss.eq.60) then ! 60 seconds = 1 minute
ss=0
mi=mi+1
endif
if (mi.eq.60) then
mi=0
hh=hh+1
endif
hhmiss=10000*hh+100*mi+ss
end subroutine caldate
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine centerofmass(xl,yl,n,xcenter,ycenter)
! i i i o o
!*****************************************************************************
! *
! This routine calculates the center of mass of n points on the Earth. *
! Input are the longitudes (xl) and latitudes (yl) of the individual *
! points, output is the longitude and latitude of the centre of mass. *
! *
! Author: A. Stohl *
! *
! 24 January 2002 *
! *
!*****************************************************************************
use par_mod
implicit none
integer :: n,l
real :: xl(n),yl(n),xll,yll,xav,yav,zav,x,y,z,xcenter,ycenter
xav=0.
yav=0.
zav=0.
do l=1,n
! Convert longitude and latitude from degrees to radians
!*******************************************************
xll=xl(l)*pi180
yll=yl(l)*pi180
! Calculate 3D coordinates from longitude and latitude
!*****************************************************
x = cos(yll)*sin(xll)
y = -1.*cos(yll)*cos(xll)
z = sin(yll)
! Find the mean location in Cartesian coordinates
!************************************************
xav=xav+x
yav=yav+y
zav=zav+z
end do
xav=xav/real(n)
yav=yav/real(n)
zav=zav/real(n)
! Project the point back onto Earth's surface
!********************************************
xcenter=atan2(xav,-1.*yav)
ycenter=atan2(zav,sqrt(xav*xav+yav*yav))
! Convert back to degrees
!************************
xcenter=xcenter/pi180
ycenter=ycenter/pi180
end subroutine centerofmass
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART 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. *
! *
! FLEXPART 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 FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine clustering(xl,yl,zl,n,xclust,yclust,zclust,fclust,rms, &
rmsclust,zrms)
! i i i i o o o o o
! o o
!*****************************************************************************
! *
! This routine clusters the particle position into ncluster custers. *
! Input are the longitudes (xl) and latitudes (yl) of the individual *
! points, output are the cluster mean positions (xclust,yclust). *
! Vertical positions are not directly used for the clustering. *
! *
! For clustering, the procedure described in Dorling et al. (1992) is used.*
! *
! Dorling, S.R., Davies, T.D. and Pierce, C.E. (1992): *
! Cluster analysis: a technique for estimating the synoptic meteorological *
! controls on air and precipitation chemistry - method and applications. *
! Atmospheric Environment 26A, 2575-2581. *
! *
! *
! Author: A. Stohl *
! *
! 1 February 2002 *
! *
! Variables: *
! fclust fraction of particles belonging to each cluster *
! ncluster number of clusters to be used *
! rms total horizontal rms distance after clustering *
! rmsclust horizontal rms distance for each individual cluster *
! zrms total vertical rms distance after clustering *
! xclust,yclust, Cluster centroid positions *
! zclust *
! xl,yl,zl particle positions *
! *
!*****************************************************************************
use par_mod
implicit none
integer :: n,i,j,l,nclust(maxpart),numb(ncluster),ncl
real :: xl(n),yl(n),zl(n),xclust(ncluster),yclust(ncluster),x,y,z
real :: zclust(ncluster),distance2,distances,distancemin,rms,rmsold
real :: xav(ncluster),yav(ncluster),zav(ncluster),fclust(ncluster)
real :: rmsclust(ncluster)
real :: zdist,zrms
if (n.lt.ncluster) return
rmsold=-5.
! Convert longitude and latitude from degrees to radians
!*******************************************************
do i=1,n
nclust(i)=i
xl(i)=xl(i)*pi180
yl(i)=yl(i)*pi180
end do
! Generate a seed for each cluster
!*********************************
do j=1,ncluster
zclust(j)=0.
xclust(j)=xl(j*n/ncluster)
yclust(j)=yl(j*n/ncluster)
end do
! Iterative loop to compute the cluster means
!********************************************
do l=1,100
! Assign each particle to a cluster: criterion minimum distance to the
! cluster mean position
!*********************************************************************
do i=1,n
distancemin=10.**10.
do j=1,ncluster
distances=distance2(yl(i),xl(i),yclust(j),xclust(j))
if (distances.lt.distancemin) then
distancemin=distances
ncl=j
endif
end do
nclust(i)=ncl
end do
! Recalculate the cluster centroid position: convert to 3D Cartesian coordinates,
! calculate mean position, and re-project this point onto the Earth's surface
!*****************************************************************************
do j=1,ncluster
xav(j)=0.
yav(j)=0.
zav(j)=0.
rmsclust(j)=0.
numb(j)=0
end do
rms=0.
do i=1,n
numb(nclust(i))=numb(nclust(i))+1
distances=distance2(yl(i),xl(i), &
yclust(nclust(i)),xclust(nclust(i)))
! rms is the total rms of all particles
! rmsclust is the rms for a particular cluster
!*********************************************
rms=rms+distances*distances
rmsclust(nclust(i))=rmsclust(nclust(i))+distances*distances
! Calculate Cartesian 3D coordinates from longitude and latitude
!***************************************************************
x = cos(yl(i))*sin(xl(i))
y = -1.*cos(yl(i))*cos(xl(i))
z = sin(yl(i))
xav(nclust(i))=xav(nclust(i))+x
yav(nclust(i))=yav(nclust(i))+y
zav(nclust(i))=zav(nclust(i))+z
end do
rms=sqrt(rms/real(n))
! Find the mean location in Cartesian coordinates
!************************************************
do j=1,ncluster
if (numb(j).gt.0) then
rmsclust(j)=sqrt(rmsclust(j)/real(numb(j)))
xav(j)=xav(j)/real(numb(j))
yav(j)=yav(j)/real(numb(j))
zav(j)=zav(j)/real(numb(j))
! Project the point back onto Earth's surface
!********************************************
xclust(j)=atan2(xav(j),-1.*yav(j))
yclust(j)=atan2(zav(j),sqrt(xav(j)*xav(j)+yav(j)*yav(j)))
endif
end do
! Leave the loop if the RMS distance decreases only slightly between 2 iterations
!*****************************************************************************
if ((l.gt.1).and.(abs(rms-rmsold)/rmsold.lt.0.005)) goto 99
rmsold=rms
end do
99 continue
! Convert longitude and latitude from radians to degrees
!*******************************************************
do i=1,n
xl(i)=xl(i)/pi180
yl(i)=yl(i)/pi180
zclust(nclust(i))=zclust(nclust(i))+zl(i)
end do
do j=1,ncluster
xclust(j)=xclust(j)/pi180
yclust(j)=yclust(j)/pi180
if (numb(j).gt.0) zclust(j)=zclust(j)/real(numb(j))
fclust(j)=100.*real(numb(j))/real(n)
end do
! Determine total vertical RMS deviation
!***************************************
zrms=0.
do i=1,n
zdist=zl(i)-zclust(nclust(i))
zrms=zrms+zdist*zdist
end do
if (zrms.gt.0.) zrms=sqrt(zrms/real(n))
end subroutine clustering
This diff is collapsed.
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment