Commit c0884a86 authored by pesei's avatar pesei
Browse files

replace CTBTO code for checking type of GRIB

parent f251e576
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
File mode changed from 100755 to 100644
......@@ -39,6 +39,9 @@ program flexpart
! - Distinguished calls to ecmwf/gfs gridcheck versions based on *
! detected metdata format *
! - Passed metdata format down to timemanager *
! *
! Petra Seibert, 2018-06-26: simplified version met data format detection *
! *
!*****************************************************************************
! *
! Variables: *
......@@ -53,7 +56,7 @@ program flexpart
use conv_mod
use random_mod, only: gasdev1
use class_gribfile
use check_gribfile_mod
#ifdef USE_NCF
use netcdf_output_mod, only: writeheader_netcdf
......@@ -61,12 +64,9 @@ program flexpart
implicit none
integer :: i,j,ix,jy,inest
integer :: i,j,ix,jy,inest,id_centre
integer :: idummy = -320
character(len=256) :: inline_options !pathfile, flexversion, arg2
integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
integer :: detectformat
! Initialize arrays in com_mod
......@@ -189,15 +189,14 @@ program flexpart
! 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'
call cg_get_centre(path(3)(1:length(3)) // trim(wfname(1)), id_centre)
if (id_centre.eq.icg_id_ecmwf) then
print *,'ECMWF met data detected'
elseif (id_centre.eq.icg_id_ncep) then
print *,'NCEP met data detected'
else
print *,'Unknown metdata format'
print *,'Unknown met data format'
stop
endif
......@@ -219,7 +218,7 @@ program flexpart
write(*,*) 'call gridcheck'
endif
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
if (id_centre.eq.icg_id_ecmwf) then
call gridcheck_ecmwf
else
call gridcheck_gfs
......@@ -448,7 +447,7 @@ program flexpart
print*,'call timemanager'
endif
call timemanager(metdata_format)
call timemanager(id_centre)
! NIK 16.02.2005
do i=1,nspec
......
......@@ -39,6 +39,9 @@ program flexpart
! - Distinguished calls to ecmwf/gfs gridcheck versions based on *
! detected metdata format *
! - Passed metdata format down to timemanager *
! *
! Petra Seibert, 2018-06-26: simplified version met data format detection *
! *
!*****************************************************************************
! *
! Variables: *
......@@ -53,7 +56,7 @@ program flexpart
use conv_mod
use mpi_mod
use random_mod, only: gasdev1
use class_gribfile
use check_gribfile_mod
#ifdef USE_NCF
use netcdf_output_mod, only: writeheader_netcdf
......@@ -61,11 +64,9 @@ program flexpart
implicit none
integer :: i,j,ix,jy,inest
integer :: i,j,ix,jy,inest,id_centre
integer :: idummy = -320
character(len=256) :: inline_options !pathfile, flexversion, arg2
integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
integer :: detectformat
integer(selected_int_kind(16)), dimension(maxspec) :: tot_b=0, &
& tot_i=0
......@@ -202,15 +203,13 @@ program flexpart
! Detect metdata format
!**********************
metdata_format = detectformat()
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
call cg_get_centre(path(3)(1:length(3)) // trim(wfname(1)), id_centre)
if (id_centre.eq.icg_id_ecmwf) then
if (lroot) print *,'ECMWF metdata detected'
elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
if (lroot) print *,'NCEP metdata detected'
elseif (id_centre.eq.icg_id_ncep) then
if (lroot) print *,'NCEP met data detected'
else
if (lroot) print *,'Unknown metdata format'
if (lroot) print *,'Unknown met data format'
stop
endif
......@@ -232,7 +231,7 @@ program flexpart
write(*,*) 'call gridcheck'
endif
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
if (id_centre.eq.icg_id_ecmwf) then
call gridcheck_ecmwf
else
call gridcheck_gfs
......@@ -479,7 +478,7 @@ program flexpart
endif
call timemanager(metdata_format)
call timemanager(id_centre)
! NIK 16.02.2005
......
......@@ -19,7 +19,7 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcmatrix(lconv,delt,cbmf,metdata_format)
subroutine calcmatrix(lconv,delt,cbmf,id_centre)
! o i o
!*****************************************************************************
! *
......@@ -39,24 +39,27 @@ subroutine calcmatrix(lconv,delt,cbmf,metdata_format)
! Marian Harustak, 12.5.2017 *
! - Merged calcmatrix and calcmatrix_gfs into one routine using if-then *
! for meteo-type dependent code *
! *
! Petra Seibert, 2018-06-26: simplified version met data format detection *
! *
!*****************************************************************************
! *
! 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) *
! id_centre format of metdata (ecmwf/gfs) *
! *
!*****************************************************************************
use par_mod
use com_mod
use conv_mod
use class_gribfile
use check_gribfile_mod
implicit none
real :: rlevmass,summe
integer :: metdata_format
integer :: id_centre
integer :: iflag, k, kk, kuvz
......@@ -87,7 +90,7 @@ subroutine calcmatrix(lconv,delt,cbmf,metdata_format)
! 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
if (id_centre.eq.icg_id_ecmwf) then
pconv(k) = (akz(kuvz) + bkz(kuvz)*psconv)
phconv(kuvz) = (akm(kuvz) + bkm(kuvz)*psconv)
else
......
......@@ -19,7 +19,7 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
subroutine calcpar(n,uuh,vvh,pvh,id_centre)
! i i i o
!*****************************************************************************
! *
......@@ -45,8 +45,9 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! Marian Harustak, 12.5.2017 *
! - Merged calcpar and calcpar_gfs into one routine using if-then *
! for meteo-type dependent code *
!*****************************************************************************
! *
! Petra Seibert, 2018-06-26: simplified version met data format detection *
! *
!*****************************************************************************
! *
! Variables: *
......@@ -54,7 +55,7 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! uuh *
! vvh *
! pvh *
! metdata_format format of metdata (ecmwf/gfs) *
! id_centre format of metdata (ecmwf/gfs) *
! *
! Constants: *
! *
......@@ -67,11 +68,11 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
use par_mod
use com_mod
use class_gribfile
use check_gribfile_mod
implicit none
integer :: metdata_format
integer :: id_centre
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
......@@ -124,7 +125,7 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! 2) Calculation of inverse Obukhov length scale
!***********************************************
if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
if (id_centre.eq.icg_id_ncep) then
! NCEP version: find first level above ground
llev = 0
do i=1,nuvz
......@@ -136,11 +137,13 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! 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),akm,bkm,akz(llev),metdata_format)
tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n), &
akm,bkm,akz(llev),id_centre)
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,akzdummy,metdata_format)
tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n), &
akm,bkm,akzdummy,id_centre)
end if
if (ol.ne.0.) then
......@@ -160,15 +163,15 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
qvlev(i)=qvh(ix,jy,i,n)
end do
if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
if (id_centre.eq.icg_id_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),hmixdummy,wstar(ix,jy,1,n),hmixplus,metdata_format)
td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus,id_centre)
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)
td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus,id_centre)
end if
if(lsubgrid.eq.1) then
......@@ -217,7 +220,7 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
ps(ix,jy,1,n))
pold=ps(ix,jy,1,n)
zold=0.
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
if (id_centre.eq.icg_id_ecmwf) then
loop_start=2
else
loop_start=llev
......@@ -241,7 +244,7 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! to be identified as the tropopause
!************************************************************************
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
if (id_centre.eq.icg_id_ecmwf) then
loop_start=1
else
loop_start=llev
......
......@@ -19,7 +19,7 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calcpar_nests(n,uuhn,vvhn,pvhn,metdata_format)
subroutine calcpar_nests(n,uuhn,vvhn,pvhn,id_centre)
! i i i o
!*****************************************************************************
! *
......@@ -43,12 +43,15 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn,metdata_format)
! *
! Unified ECMWF and GFS builds *
! Marian Harustak, 12.5.2017 *
! - Added passing of metdata_format as it was needed by called routines *
! - Added passing of id_centre as it was needed by called routines *
! *
! Petra Seibert, 2018-06-26: simplified version met data format detection *
! *
!*****************************************************************************
! *
! Variables: *
! n temporal index for meteorological fields (1 to 3) *
! metdata_format format of metdata (ecmwf/gfs) *
! id_centre format of metdata (ecmwf/gfs) *
! *
! Constants: *
! *
......@@ -64,7 +67,7 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn,metdata_format)
implicit none
integer :: metdata_format
integer :: id_centre
integer :: n,ix,jy,i,l,kz,lz,kzmin
real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus,dummyakzllev
real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
......@@ -115,7 +118,7 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn,metdata_format)
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,dummyakzllev,metdata_format)
sshfn(ix,jy,1,n,l),akm,bkm,dummyakzllev,id_centre)
if (ol.ne.0.) then
olin(ix,jy,1,n,l)=1./ol
else
......@@ -136,7 +139,7 @@ subroutine calcpar_nests(n,uuhn,vvhn,pvhn,metdata_format)
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,metdata_format)
wstarn(ix,jy,1,n,l),hmixplus,id_centre)
if(lsubgrid.eq.1) then
subsceff=min(excessoron(ix,jy,l),hmixplus)
......
......@@ -72,3 +72,26 @@ Some minor changes
for the vertical grid. Also correct a minor bug in the calculation
of height (array assignment where it was not intended)
3) Add back the SAVE attribute to INIT, just to be sure
* checking type (EC/GFS) of grib files PS 2018-06-26
===================
replaced the CTBTO code by shorter & simper code with proper license
Affected files:
calcmatrix.f90
calcpar.f90
calcpar_nests.f90
check_gribfile_mod.f90
convmix.f90
detectformat.f90
FLEXPART.f90
FLEXPART_MPI.f90
getfields.f90
getfields_mpi.f90
obukhov.f90
richardson.f90
timemanager.f90
timemanager_mpi.f90
module check_gribfile_mod
!*
! Valid-License-Identifier: GPL-3.0-or-later
! Copyright (c) 2018 Petra Seibert (petra seibert at univie ac at)
!
! Prepared for use in FLEXPART (see flexpart.eu) version 10+
! 1. provide centre ids for ECMWF and NCEP
! 2. obtain centre id from a given grib file
! intended to replace class_gribfile from ctbto project
! requires grib_api_fortran90 and grib_api
! either from eccodes or grib_api
! I am using the convention to put an abbreviated module name (here: cg)
! in front of public entities. If they are integer, then icg.
! subroutines / functions from external libs with upper first letter
implicit none
integer, parameter :: icg_id_ncep = 7, icg_id_ecmwf = 98
!! official centre codes
integer :: icentre ! centre ID found
contains
subroutine cg_get_centre(pfname, icentre)
! this subroutine returns a code for the centre which produced the gribfile
use grib_api
integer, intent(out) :: icentre
integer :: ifile ! grib file handle
integer :: iret ! status code
integer :: igrib ! message handle
character (len=*), intent(in) :: pfname ! path+filenmae
call Grib_open_file(ifile,pfname,'r',iret)
if (iret == 0) then
call Grib_new_from_file(ifile,igrib) ! load first message
call Grib_get(igrib,'centre',icentre) ! read centre ID
call Grib_close_file(ifile) ! done
if (icentre .ne. icg_id_ecmwf .and. &
icentre .ne. icg_id_ncep) then
!! centre not foreseen in Fp
write (*,*) ' #### FLEXPART MODEL ERROR! Met input file'
write (*,*) ' #### '//trim(pfname)
write (*,*) ' #### is neither ECMWF nor NCEP grib file'
stop 'error in check_gribfile'
endif
else
!! reading gribfile failed
write (*,*) ' #### FLEXPART MODEL ERROR! Met input file'
write (*,*) ' #### '//trim(pfname)
write (*,*) ' #### cannot be opened with grib_api'
stop 'error in check_gribfile'
endif
end subroutine cg_get_centre
end module check_gribfile_mod
! what needs to be done to replace the current code with this one:
! see email
......@@ -19,7 +19,7 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine convmix(itime,metdata_format)
subroutine convmix(itime,id_centre)
! i
!**************************************************************
!handles all the calculations related to convective mixing
......@@ -38,15 +38,18 @@ subroutine convmix(itime,metdata_format)
!
! Unified ECMWF and GFS builds
! Marian Harustak, 12.5.2017
! - Merged convmix and convmix_gfs into one routine using if-then
! - Merged convmix and convmix_gfs into one routine using if-then
! for meteo-type dependent code
!
! Petra Seibert, 2018-06-26: simplified version met data format detection
!
!**************************************************************
use flux_mod
use par_mod
use com_mod
use conv_mod
use class_gribfile
use check_gribfile_mod
implicit none
......@@ -54,7 +57,7 @@ subroutine convmix(itime,metdata_format)
integer :: ipconv
integer :: jy, kpart, ktop, ngrid,kz
integer :: igrid(maxpart), ipoint(maxpart), igridn(maxpart,maxnests)
integer :: metdata_format
integer :: id_centre
! itime [s] current time
! igrid(maxpart) horizontal grid position of each particle
......@@ -115,7 +118,7 @@ subroutine convmix(itime,metdata_format)
!**********************************************************
ngrid=0
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
if (id_centre.eq.icg_id_ecmwf) then
do j=numbnests,1,-1
if ( x.gt.xln(j)+eps .and. x.lt.xrn(j)-eps .and. &
y.gt.yln(j)+eps .and. y.lt.yrn(j)-eps ) then
......@@ -188,7 +191,7 @@ subroutine convmix(itime,metdata_format)
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
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
if (id_centre.eq.icg_id_ecmwf) then
do kz=1,nuvz-1 !bugfix
tconv(kz)=(tth(ix,jy,kz+1,mind1)*dt2+ &
tth(ix,jy,kz+1,mind2)*dt1)*dtt
......@@ -207,7 +210,7 @@ subroutine convmix(itime,metdata_format)
end if
! Calculate translocation matrix
call calcmatrix(lconv,delt,cbaseflux(ix,jy),metdata_format)
call calcmatrix(lconv,delt,cbaseflux(ix,jy),id_centre)
igrold = igr
ktop = 0
endif
......@@ -284,7 +287,7 @@ subroutine convmix(itime,metdata_format)
! calculate translocation matrix
!*******************************
call calcmatrix(lconv,delt,cbasefluxn(ix,jy,inest),metdata_format)
call calcmatrix(lconv,delt,cbasefluxn(ix,jy,inest),id_centre)
igrold = igr
ktop = 0
endif
......
!**********************************************************************
! 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/>. *
!**********************************************************************
integer function detectformat()
!*****************************************************************************
! *
! This routine reads the 1st file with windfields to determine *
! the format. *
! *
! Authors: M. Harustak *
! *
! 6 May 2015 *
! *
! Unified ECMWF and GFS builds *
! Marian Harustak, 12.5.2017 *
! - Added routine to FP10 Flexpart distribution *
!*****************************************************************************
! *
! Variables: *
! fname file name of file to check *
! *
!*****************************************************************************
use par_mod
use com_mod
use class_gribfile
implicit none
character(len=255) :: filename
character(len=255) :: wfname1(maxwf)
integer :: metdata_format
! If no file is available
if ( maxwf.le.0 ) then
print*,'No wind file available'
detectformat = GRIBFILE_CENTRE_UNKNOWN
return
endif
! construct filename
filename = path(3)(1:length(3)) // trim(wfname(1))
! get format
detectformat = gribfile_centre(TRIM(filename))
end
......@@ -19,7 +19,7 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine getfields(itime,nstop,metdata_format)
subroutine getfields(itime,nstop,id_centre)
! i o
!*****************************************************************************
! *
......@@ -43,7 +43,10 @@ subroutine getfields(itime,nstop,metdata_format)
! *
! Unified ECMWF and GFS builds *
! Marian Harustak, 12.5.2017 *
! - Added passing of metdata_format as it was needed by called routines *
! - Added passing of id_centre as it was needed by called routines *
! *
! Petra Seibert, 2018-06-26: simplified version met data format detection *
! *
!*****************************************************************************
! *
! Variables: *
......@@ -61,7 +64,7 @@ subroutine getfields(itime,nstop,metdata_format)
! ww(0:nxmax,0:nymax,nwzmax,2) wind components in z-direction [deltaeta/s]*