Skip to content
GitLab
Commit 5b34509e authored by Espen Sollum's avatar Espen Sollum
Browse files

Deleted unused files

parent 6ecb30a4
Expand all Hide whitespace changes
Inline Side-by-side
src/calcmatrix_gfs.f90 deleted 100644 → 0
View file @ 6ecb30a4
!**********************************************************************
! 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
src/calcpar_gfs.f90 deleted 100644 → 0
View file @ 6ecb30a4
!**********************************************************************
! 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
src/convmix_gfs.f90 deleted 100644 → 0
View file @ 6ecb30a4
!**********************************************************************
! 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 convmix(itime)
! i
!**************************************************************
!handles all the calculations related to convective mixing
!Petra Seibert, Bernd C. Krueger, Feb 2001
!nested grids included, Bernd C. Krueger, May 2001
!
!Changes by Caroline Forster, April 2004 - February 2005:
! convmix called every lsynctime seconds
!CHANGES by A. Stohl:
! various run-time optimizations - February 2005
!CHANGES by C. Forster, November 2005, NCEP GFS version
! in the ECMWF version convection is calculated on the
! original eta-levels
! in the GFS version convection is calculated on the
! FLEXPART levels
!**************************************************************
use par_mod
use com_mod
use conv_mod
implicit none
integer :: igr,igrold, ipart, itime, ix, j, inest
integer :: ipconv
integer :: jy, kpart, ktop, ngrid,kz
integer :: igrid(maxpart), ipoint(maxpart), igridn(maxpart,maxnests)
! itime [s] current time
! igrid(maxpart) horizontal grid position of each particle
! igridn(maxpart,maxnests) dto. for nested grids
! ipoint(maxpart) pointer to access particles according to grid position
logical :: lconv
real :: x, y, xtn,ytn, ztold, delt
real :: dt1,dt2,dtt
integer :: mind1,mind2
! dt1,dt2,dtt,mind1,mind2 variables used for time interpolation
integer :: itage,nage
!monitoring variables
!real sumconv,sumall
! Calculate auxiliary variables for time interpolation
!*****************************************************
dt1=real(itime-memtime(1))
dt2=real(memtime(2)-itime)
dtt=1./(dt1+dt2)
mind1=memind(1)
mind2=memind(2)
delt=real(abs(lsynctime))
lconv = .false.
! if no particles are present return after initialization
!********************************************************
if (numpart.le.0) return
! Assign igrid and igridn, which are pseudo grid numbers indicating particles
! that are outside the part of the grid under consideration
! (e.g. particles near the poles or particles in other nests).
! Do this for all nests but use only the innermost nest; for all others
! igrid shall be -1
! Also, initialize index vector ipoint
!************************************************************************
do ipart=1,numpart
igrid(ipart)=-1
do j=numbnests,1,-1
igridn(ipart,j)=-1
end do
ipoint(ipart)=ipart
! do not consider particles that are (yet) not part of simulation
if (itra1(ipart).ne.itime) goto 20
x = xtra1(ipart)
y = ytra1(ipart)
! Determine which nesting level to be used
!**********************************************************
ngrid=0
do j=numbnests,1,-1
if ( x.gt.xln(j) .and. x.lt.xrn(j) .and. &
y.gt.yln(j) .and. y.lt.yrn(j) ) then
ngrid=j
goto 23
endif
end do
23 continue
! Determine nested grid coordinates
!**********************************
if (ngrid.gt.0) then
! nested grids
xtn=(x-xln(ngrid))*xresoln(ngrid)
ytn=(y-yln(ngrid))*yresoln(ngrid)
ix=nint(xtn)
jy=nint(ytn)
igridn(ipart,ngrid) = 1 + jy*nxn(ngrid) + ix
else if(ngrid.eq.0) then
! mother grid
ix=nint(x)
jy=nint(y)
igrid(ipart) = 1 + jy*nx + ix
endif
20 continue
end do
!sumall = 0.
!sumconv = 0.
!*****************************************************************************
! 1. Now, do everything for the mother domain and, later, for all of the nested domains
! While all particles have to be considered for redistribution, the Emanuel convection
! scheme only needs to be called once for every grid column where particles are present.
! Therefore, particles are sorted according to their grid position. Whenever a new grid
! cell is encountered by looping through the sorted particles, the convection scheme is called.
!*****************************************************************************
! sort particles according to horizontal position and calculate index vector IPOINT
call sort2(numpart,igrid,ipoint)
! Now visit all grid columns where particles are present
! by going through the sorted particles
igrold = -1
do kpart=1,numpart
igr = igrid(kpart)
if (igr .eq. -1) goto 50
ipart = ipoint(kpart)
! sumall = sumall + 1
if (igr .ne. igrold) then
! we are in a new grid column
jy = (igr-1)/nx
ix = igr - jy*nx - 1
! Interpolate all meteorological data needed for the convection scheme
psconv=(ps(ix,jy,1,mind1)*dt2+ps(ix,jy,1,mind2)*dt1)*dtt
tt2conv=(tt2(ix,jy,1,mind1)*dt2+tt2(ix,jy,1,mind2)*dt1)*dtt
td2conv=(td2(ix,jy,1,mind1)*dt2+td2(ix,jy,1,mind2)*dt1)*dtt
!!$ do kz=1,nconvlev+1 !old
do kz=1,nuvz-1 !bugfix
pconv(kz)=(pplev(ix,jy,kz,mind1)*dt2+ &
pplev(ix,jy,kz,mind2)*dt1)*dtt
tconv(kz)=(tt(ix,jy,kz,mind1)*dt2+ &
tt(ix,jy,kz,mind2)*dt1)*dtt
qconv(kz)=(qv(ix,jy,kz,mind1)*dt2+ &
qv(ix,jy,kz,mind2)*dt1)*dtt
end do
! Calculate translocation matrix
call calcmatrix(lconv,delt,cbaseflux(ix,jy))
igrold = igr
ktop = 0
endif
! treat particle only if column has convection
if (lconv .eqv. .true.) then
! assign new vertical position to particle
ztold=ztra1(ipart)
call redist(ipart,ktop,ipconv)
! if (ipconv.le.0) sumconv = sumconv+1
! Calculate the gross fluxes across layer interfaces
!***************************************************
if (iflux.eq.1) then
itage=abs(itra1(ipart)-itramem(ipart))
do nage=1,nageclass
if (itage.lt.lage(nage)) goto 37
end do
37 continue
if (nage.le.nageclass) &
call calcfluxes(nage,ipart,real(xtra1(ipart)), &
real(ytra1(ipart)),ztold)
endif
endif !(lconv .eqv. .true)
50 continue
end do
!*****************************************************************************
! 2. Nested domains
!*****************************************************************************
! sort particles according to horizontal position and calculate index vector IPOINT
do inest=1,numbnests
do ipart=1,numpart
ipoint(ipart)=ipart