Skip to content

GitLab

Commit bb579a9d authored by Espen Sollum's avatar Espen Sollum
Browse files

Merge branch 'ctbto' into dev

parents aa8c34a4 dd6a4aae
Hide whitespace changes
Inline Side-by-side
src/FLEXPART.f90
View file @ bb579a9d
......@@ -32,6 +32,14 @@ program flexpart
! 18 May 1996 *
! *
!*****************************************************************************
! Changes: *
! Unified ECMWF and GFS builds *
! Marian Harustak, 12.5.2017 *
! - Added detection of metdata format using gributils routines *
! - Distinguished calls to ecmwf/gfs gridcheck versions based on *
! detected metdata format *
! - Passed metdata format down to timemanager *
!*****************************************************************************
! *
! Variables: *
! *
......@@ -45,12 +53,16 @@ program flexpart
use conv_mod
use netcdf_output_mod, only: writeheader_netcdf
use random_mod, only: gasdev1
use class_gribfile
implicit none
integer :: i,j,ix,jy,inest
integer :: idummy = -320
character(len=256) :: inline_options !pathfile, flexversion, arg2
integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
integer :: detectformat
! Initialize arrays in com_mod
......@@ -69,7 +81,7 @@ program flexpart
! FLEXPART version string
flexversion_major = '10' ! Major version number, also used for species file names
flexversion='Version '//trim(flexversion_major)//'.1beta (2016-11-02)'
flexversion='Version '//trim(flexversion_major)//'.2beta (2017-08-01)'
verbosity=0
! Read the pathnames where input/output files are stored
......@@ -171,6 +183,22 @@ program flexpart
endif
call readavailable
! Detect metdata format
!**********************
metdata_format = detectformat()
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
print *,'ECMWF metdata detected'
elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
print *,'NCEP metdata detected'
else
print *,'Unknown metdata format'
return
endif
! If nested wind fields are used, allocate arrays
!************************************************
......@@ -187,7 +215,11 @@ program flexpart
write(*,*) 'call gridcheck'
endif
call gridcheck
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
call gridcheck_ecmwf
else
call gridcheck_gfs
end if
if (verbosity.gt.1) then
CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
......@@ -410,7 +442,7 @@ program flexpart
print*,'call timemanager'
endif
call timemanager
call timemanager(metdata_format)
! NIK 16.02.2005
do i=1,nspec
......
src/FLEXPART_MPI.f90
View file @ bb579a9d
......@@ -32,6 +32,14 @@ program flexpart
! 18 May 1996 *
! *
!*****************************************************************************
! Changes: *
! Unified ECMWF and GFS builds *
! Marian Harustak, 12.5.2017 *
! - Added detection of metdata format using gributils routines *
! - Distinguished calls to ecmwf/gfs gridcheck versions based on *
! detected metdata format *
! - Passed metdata format down to timemanager *
!*****************************************************************************
! *
! Variables: *
! *
......@@ -46,12 +54,16 @@ program flexpart
use mpi_mod
use netcdf_output_mod, only: writeheader_netcdf
use random_mod, only: gasdev1
use class_gribfile
implicit none
integer :: i,j,ix,jy,inest
integer :: idummy = -320
character(len=256) :: inline_options !pathfile, flexversion, arg2
integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
integer :: detectformat
! Initialize mpi
......@@ -79,7 +91,7 @@ program flexpart
! FLEXPART version string
flexversion_major = '10' ! Major version number, also used for species file names
flexversion='Ver. '//trim(flexversion_major)//'.1beta MPI (2016-11-02)'
flexversion='Ver. '//trim(flexversion_major)//'.2beta MPI (2017-08-01)'
verbosity=0
! Read the pathnames where input/output files are stored
......@@ -196,6 +208,22 @@ program flexpart
endif
call readavailable
! Detect metdata format
!**********************
metdata_format = detectformat()
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
print *,'ECMWF metdata detected'
elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
print *,'NCEP metdata detected'
else
print *,'Unknown metdata format'
return
endif
! If nested wind fields are used, allocate arrays
!************************************************
......@@ -212,7 +240,12 @@ program flexpart
write(*,*) 'call gridcheck'
endif
call gridcheck
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
call gridcheck_ecmwf
else
call gridcheck_gfs
end if
if (verbosity.gt.1 .and. lroot) then
call system_clock(count_clock, count_rate, count_max)
......@@ -455,7 +488,7 @@ program flexpart
endif
call timemanager
call timemanager(metdata_format)
! NIK 16.02.2005
......
src/calcmatrix.f90
View file @ bb579a9d
......@@ -19,7 +19,7 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcmatrix(lconv,delt,cbmf)
subroutine calcmatrix(lconv,delt,cbmf,metdata_format)
! o i o
!*****************************************************************************
! *
......@@ -29,23 +29,34 @@ subroutine calcmatrix(lconv,delt,cbmf)
! *
! Petra Seibert, Bernd C. Krueger, 2000-2001 *
! *
!*****************************************************************************
! Changes: *
! changed by C. Forster, November 2003 - February 2004 *
! array fmassfrac(nconvlevmax,nconvlevmax) represents *
! the convective redistribution matrix for the particles *
! *
! Unified ECMWF and GFS builds *
! Marian Harustak, 12.5.2017 *
! - Merged calcmatrix and calcmatrix_gfs into one routine using if-then *
! for meteo-type dependent code *
!*****************************************************************************
! *
! lconv indicates whether there is convection in this cell, or not *
! delt time step for convection [s] *
! cbmf cloud base mass flux *
! metdata_format format of metdata (ecmwf/gfs) *
! *
!*****************************************************************************
use par_mod
use com_mod
use conv_mod
use class_gribfile
implicit none
real :: rlevmass,summe
integer :: metdata_format
integer :: iflag, k, kk, kuvz
......@@ -76,16 +87,20 @@ subroutine calcmatrix(lconv,delt,cbmf)
! Emanuel subroutine needs pressure in hPa, therefore convert all pressures
do kuvz = 2,nuvz
k = kuvz-1
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
else
phconv(kuvz) = 0.5*(pconv(kuvz)+pconv(k))
endif
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
end do
end do
!note that Emanuel says it is important
......
src/calcmatrix_gfs.f90 deleted 100644 → 0
View file @ aa8c34a4
!**********************************************************************
! 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.f90
View file @ bb579a9d
......@@ -19,7 +19,7 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcpar(n,uuh,vvh,pvh)
subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! i i i o
!*****************************************************************************
! *
......@@ -36,13 +36,25 @@ subroutine calcpar(n,uuh,vvh,pvh)
! convection scheme: *
! new variables in call to richardson *
! *
! CHANGE 17/11/2005 Caroline Forster NCEP GFS version *
! *
! Changes, Bernd C. Krueger, Feb. 2001: *
! Variables tth and qvh (on eta coordinates) in common block *
! *
! Unified ECMWF and GFS builds *
! Marian Harustak, 12.5.2017 *
! - Merged calcpar and calcpar_gfs into one routine using if-then *
! for meteo-type dependent code *
!*****************************************************************************
! 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) *
! uuh *
! vvh *
! pvh *
! metdata_format format of metdata (ecmwf/gfs) *
! *
! Constants: *
! *
......@@ -55,13 +67,15 @@ subroutine calcpar(n,uuh,vvh,pvh)
use par_mod
use com_mod
use class_gribfile
implicit none
integer :: n,ix,jy,i,kz,lz,kzmin
integer :: metdata_format
integer :: n,ix,jy,i,kz,lz,kzmin,llev,loop_start
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 :: 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)
......@@ -110,8 +124,25 @@ subroutine calcpar(n,uuh,vvh,pvh)
! 2) Calculation of inverse Obukhov length scale
!***********************************************
if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
! 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,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm)
tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akz(llev),metdata_format)
else
llev=0
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,akz(llev),metdata_format)
end if
if (ol.ne.0.) then
oli(ix,jy,1,n)=1./ol
else
......@@ -129,9 +160,16 @@ subroutine calcpar(n,uuh,vvh,pvh)
qvlev(i)=qvh(ix,jy,i,n)
end do
if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
! 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),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus)
td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus,metdata_format)
else
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,metdata_format)
end if
if(lsubgrid.eq.1) then
subsceff=min(excessoro(ix,jy),hmixplus)
......@@ -172,14 +210,19 @@ subroutine calcpar(n,uuh,vvh,pvh)
! Calculate height of thermal tropopause (Hoinka, 1997)
!******************************************************
! 1) Calculate altitudes of ECMWF model levels
!*********************************************
! 1) Calculate altitudes of 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
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
loop_start=2
else
loop_start=llev
end if
do kz=loop_start,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))
......@@ -198,7 +241,13 @@ subroutine calcpar(n,uuh,vvh,pvh)
! to be identified as the tropopause
!************************************************************************
do kz=1,nuvz
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
loop_start=1
else
loop_start=llev
end if
do kz=loop_start,nuvz
if (zlev(kz).ge.altmin) then
kzmin=kz
goto 45
......
src/calcpar_gfs.f90 deleted 100644 → 0
View file @ aa8c34a4
!**********************************************************************
! 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