Commit 6ecb30a4 authored by Espen Sollum's avatar Espen Sollum
Browse files

Merged changes from CTBTO project

parent 61e07bac
...@@ -32,6 +32,14 @@ program flexpart ...@@ -32,6 +32,14 @@ program flexpart
! 18 May 1996 * ! 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: * ! Variables: *
! * ! *
...@@ -45,12 +53,16 @@ program flexpart ...@@ -45,12 +53,16 @@ program flexpart
use conv_mod use conv_mod
use netcdf_output_mod, only: writeheader_netcdf use netcdf_output_mod, only: writeheader_netcdf
use random_mod, only: gasdev1 use random_mod, only: gasdev1
use class_gribfile
implicit none implicit none
integer :: i,j,ix,jy,inest integer :: i,j,ix,jy,inest
integer :: idummy = -320 integer :: idummy = -320
character(len=256) :: inline_options !pathfile, flexversion, arg2 character(len=256) :: inline_options !pathfile, flexversion, arg2
integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
integer :: detectformat
! Initialize arrays in com_mod ! Initialize arrays in com_mod
...@@ -171,6 +183,22 @@ program flexpart ...@@ -171,6 +183,22 @@ program flexpart
endif endif
call readavailable 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 ! If nested wind fields are used, allocate arrays
!************************************************ !************************************************
...@@ -187,7 +215,11 @@ program flexpart ...@@ -187,7 +215,11 @@ program flexpart
write(*,*) 'call gridcheck' write(*,*) 'call gridcheck'
endif 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 if (verbosity.gt.1) then
CALL SYSTEM_CLOCK(count_clock, count_rate, count_max) CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
...@@ -410,7 +442,7 @@ program flexpart ...@@ -410,7 +442,7 @@ program flexpart
print*,'call timemanager' print*,'call timemanager'
endif endif
call timemanager call timemanager(metdata_format)
! NIK 16.02.2005 ! NIK 16.02.2005
do i=1,nspec do i=1,nspec
......
...@@ -32,6 +32,14 @@ program flexpart ...@@ -32,6 +32,14 @@ program flexpart
! 18 May 1996 * ! 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: * ! Variables: *
! * ! *
...@@ -46,12 +54,16 @@ program flexpart ...@@ -46,12 +54,16 @@ program flexpart
use mpi_mod use mpi_mod
use netcdf_output_mod, only: writeheader_netcdf use netcdf_output_mod, only: writeheader_netcdf
use random_mod, only: gasdev1 use random_mod, only: gasdev1
use class_gribfile
implicit none implicit none
integer :: i,j,ix,jy,inest integer :: i,j,ix,jy,inest
integer :: idummy = -320 integer :: idummy = -320
character(len=256) :: inline_options !pathfile, flexversion, arg2 character(len=256) :: inline_options !pathfile, flexversion, arg2
integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
integer :: detectformat
! Initialize mpi ! Initialize mpi
...@@ -196,6 +208,22 @@ program flexpart ...@@ -196,6 +208,22 @@ program flexpart
endif endif
call readavailable 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 ! If nested wind fields are used, allocate arrays
!************************************************ !************************************************
...@@ -212,7 +240,12 @@ program flexpart ...@@ -212,7 +240,12 @@ program flexpart
write(*,*) 'call gridcheck' write(*,*) 'call gridcheck'
endif 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 if (verbosity.gt.1 .and. lroot) then
call system_clock(count_clock, count_rate, count_max) call system_clock(count_clock, count_rate, count_max)
...@@ -455,7 +488,7 @@ program flexpart ...@@ -455,7 +488,7 @@ program flexpart
endif endif
call timemanager call timemanager(metdata_format)
! NIK 16.02.2005 ! NIK 16.02.2005
......
...@@ -19,7 +19,7 @@ ...@@ -19,7 +19,7 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * ! 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 ! o i o
!***************************************************************************** !*****************************************************************************
! * ! *
...@@ -29,23 +29,34 @@ subroutine calcmatrix(lconv,delt,cbmf) ...@@ -29,23 +29,34 @@ subroutine calcmatrix(lconv,delt,cbmf)
! * ! *
! Petra Seibert, Bernd C. Krueger, 2000-2001 * ! Petra Seibert, Bernd C. Krueger, 2000-2001 *
! * ! *
!*****************************************************************************
! Changes: *
! changed by C. Forster, November 2003 - February 2004 * ! changed by C. Forster, November 2003 - February 2004 *
! array fmassfrac(nconvlevmax,nconvlevmax) represents * ! array fmassfrac(nconvlevmax,nconvlevmax) represents *
! the convective redistribution matrix for the particles * ! 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 * ! lconv indicates whether there is convection in this cell, or not *
! delt time step for convection [s] * ! delt time step for convection [s] *
! cbmf cloud base mass flux * ! cbmf cloud base mass flux *
! metdata_format format of metdata (ecmwf/gfs) *
! * ! *
!***************************************************************************** !*****************************************************************************
use par_mod use par_mod
use com_mod use com_mod
use conv_mod use conv_mod
use class_gribfile
implicit none implicit none
real :: rlevmass,summe real :: rlevmass,summe
integer :: metdata_format
integer :: iflag, k, kk, kuvz integer :: iflag, k, kk, kuvz
...@@ -76,8 +87,12 @@ subroutine calcmatrix(lconv,delt,cbmf) ...@@ -76,8 +87,12 @@ subroutine calcmatrix(lconv,delt,cbmf)
! Emanuel subroutine needs pressure in hPa, therefore convert all pressures ! Emanuel subroutine needs pressure in hPa, therefore convert all pressures
do kuvz = 2,nuvz do kuvz = 2,nuvz
k = kuvz-1 k = kuvz-1
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv) pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
phconv(kuvz) = (akm(kuvz) + bkm(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) dpr(k) = phconv(k) - phconv(kuvz)
qsconv(k) = f_qvsat( pconv(k), tconv(k) ) qsconv(k) = f_qvsat( pconv(k), tconv(k) )
......
...@@ -19,7 +19,7 @@ ...@@ -19,7 +19,7 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * ! 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 ! i i i o
!***************************************************************************** !*****************************************************************************
! * ! *
...@@ -36,13 +36,25 @@ subroutine calcpar(n,uuh,vvh,pvh) ...@@ -36,13 +36,25 @@ subroutine calcpar(n,uuh,vvh,pvh)
! convection scheme: * ! convection scheme: *
! new variables in call to richardson * ! 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: * ! Variables: *
! n temporal index for meteorological fields (1 to 3) * ! n temporal index for meteorological fields (1 to 3) *
! uuh *
! vvh *
! pvh *
! metdata_format format of metdata (ecmwf/gfs) *
! * ! *
! Constants: * ! Constants: *
! * ! *
...@@ -55,13 +67,15 @@ subroutine calcpar(n,uuh,vvh,pvh) ...@@ -55,13 +67,15 @@ subroutine calcpar(n,uuh,vvh,pvh)
use par_mod use par_mod
use com_mod use com_mod
use class_gribfile
implicit none 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 :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat 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 :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
real :: vvh(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 :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
...@@ -110,8 +124,25 @@ subroutine calcpar(n,uuh,vvh,pvh) ...@@ -110,8 +124,25 @@ subroutine calcpar(n,uuh,vvh,pvh)
! 2) Calculation of inverse Obukhov length scale ! 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), & 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 if (ol.ne.0.) then
oli(ix,jy,1,n)=1./ol oli(ix,jy,1,n)=1./ol
else else
...@@ -129,9 +160,16 @@ subroutine calcpar(n,uuh,vvh,pvh) ...@@ -129,9 +160,16 @@ subroutine calcpar(n,uuh,vvh,pvh)
qvlev(i)=qvh(ix,jy,i,n) qvlev(i)=qvh(ix,jy,i,n)
end do 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, & 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), & 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 if(lsubgrid.eq.1) then
subsceff=min(excessoro(ix,jy),hmixplus) subsceff=min(excessoro(ix,jy),hmixplus)
...@@ -172,14 +210,19 @@ subroutine calcpar(n,uuh,vvh,pvh) ...@@ -172,14 +210,19 @@ subroutine calcpar(n,uuh,vvh,pvh)
! Calculate height of thermal tropopause (Hoinka, 1997) ! 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))/ & tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
ps(ix,jy,1,n)) ps(ix,jy,1,n))
pold=ps(ix,jy,1,n) pold=ps(ix,jy,1,n)
zold=0. 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 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)) tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
...@@ -198,7 +241,13 @@ subroutine calcpar(n,uuh,vvh,pvh) ...@@ -198,7 +241,13 @@ subroutine calcpar(n,uuh,vvh,pvh)
! to be identified as the tropopause ! 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 if (zlev(kz).ge.altmin) then
kzmin=kz kzmin=kz
goto 45 goto 45
......
...@@ -19,7 +19,7 @@ ...@@ -19,7 +19,7 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!********************************************************************** !**********************************************************************
subroutine calcpar_nests(n,uuhn,vvhn,pvhn) subroutine calcpar_nests(n,uuhn,vvhn,pvhn,metdata_format)
! i i i o ! i i i o
!***************************************************************************** !*****************************************************************************
! * ! *
...@@ -38,12 +38,17 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn) ...@@ -38,12 +38,17 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn)
! new variables in call to richardson * ! new variables in call to richardson *
! * ! *
!***************************************************************************** !*****************************************************************************
! Changes, Bernd C. Krueger, Feb. 2001: ! Changes, Bernd C. Krueger, Feb. 2001: *
! Variables tth and qvh (on eta coordinates) in common block ! Variables tth and qvh (on eta coordinates) in common block *
! *
! Unified ECMWF and GFS builds *
! Marian Harustak, 12.5.2017 *
! - Added passing of metdata_format as it was needed by called routines *
!***************************************************************************** !*****************************************************************************
! * ! *
! Variables: * ! Variables: *
! n temporal index for meteorological fields (1 to 3) * ! n temporal index for meteorological fields (1 to 3) *
! metdata_format format of metdata (ecmwf/gfs) *
! * ! *
! Constants: * ! Constants: *
! * ! *
...@@ -59,8 +64,9 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn) ...@@ -59,8 +64,9 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn)
implicit none implicit none
integer :: metdata_format
integer :: n,ix,jy,i,l,kz,lz,kzmin integer :: n,ix,jy,i,l,kz,lz,kzmin
real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus,dummyakzllev
real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat 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)
real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
...@@ -109,7 +115,7 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn) ...@@ -109,7 +115,7 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn)
ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), & 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), & 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) sshfn(ix,jy,1,n,l),akm,bkm,dummyakzllev,metdata_format)
if (ol.ne.0.) then if (ol.ne.0.) then
olin(ix,jy,1,n,l)=1./ol olin(ix,jy,1,n,l)=1./ol
else else
...@@ -130,7 +136,7 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn) ...@@ -130,7 +136,7 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn)
call richardson(psn(ix,jy,1,n,l),ustarn(ix,jy,1,n,l),ttlev, & 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), & 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), & 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) wstarn(ix,jy,1,n,l),hmixplus,metdata_format)
if(lsubgrid.eq.1) then if(lsubgrid.eq.1) then
subsceff=min(excessoron(ix,jy,l),hmixplus) subsceff=min(excessoron(ix,jy,l),hmixplus)
......
...@@ -19,7 +19,7 @@ ...@@ -19,7 +19,7 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!********************************************************************** !**********************************************************************
subroutine convmix(itime) subroutine convmix(itime,metdata_format)
! i ! i
!************************************************************** !**************************************************************
!handles all the calculations related to convective mixing !handles all the calculations related to convective mixing
...@@ -30,12 +30,23 @@ subroutine convmix(itime) ...@@ -30,12 +30,23 @@ subroutine convmix(itime)
! convmix called every lsynctime seconds ! convmix called every lsynctime seconds
!CHANGES by A. Stohl: !CHANGES by A. Stohl:
! various run-time optimizations - February 2005 ! 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
!
! Unified ECMWF and GFS builds
! Marian Harustak, 12.5.2017
! - Merged convmix and convmix_gfs into one routine using if-then
! for meteo-type dependent code
!************************************************************** !**************************************************************
use flux_mod use flux_mod
use par_mod use par_mod
use com_mod use com_mod
use conv_mod use conv_mod
use class_gribfile