Skip to content

GitLab

Commit 0eda0087 authored by Espen Sollum's avatar Espen Sollum
Browse files

Merge branch 'ctbto_vtables' into dev

parents 20963b15 33cfc2a3
Expand all Hide whitespace changes
Inline Side-by-side
src/calcpar.f90
View file @ 0eda0087
......@@ -32,7 +32,7 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! 21 May 1995 *
! *
! ------------------------------------------------------------------ *
! Petra Seibert, Feb 2000: *
! Petra Seibert, Feb 2000: *
! convection scheme: *
! new variables in call to richardson *
! *
......@@ -45,6 +45,12 @@ 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 *
! *
! *
! Don Morton, 13.10.2017 *
! - Repairing problems from merger and documenting the merger of *
! Harustak *
! *
!*****************************************************************************
!*****************************************************************************
......@@ -75,12 +81,19 @@ subroutine calcpar(n,uuh,vvh,pvh,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),hmixdummy,akzdummy
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
!! DJM - using these as meaningless arguments to the obukhov function call
!! For the GFS version, gfs_dummy_arg(nwzmax) is used in place of the
!! akm(nwzmax) and bkm(nwzmax) used in the call to ECMWF version
!! For the ECMWF version, ecmwf_dummy_arg is used in place of the
!! akz(llev) used in the call to the GFS version.
REAL :: ecmwf_dummy_arg, gfs_dummy_arg(nwzmax)
!write(*,*) 'in calcpar writting snowheight'
!***********************************
! for test: write out snow depths
......@@ -124,6 +137,36 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! 2) Calculation of inverse Obukhov length scale
!***********************************************
!! ..... Documentation by Don Morton, 13 Oct 2017 .....
!
!
! This subroutine is a result of merging an ECMWF and a GFS version.
! In the case of the call to the obukhov() function, originally the
! call for ECMWF looked like:
!
! 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)
!
!
! and the call for GFS looked like:
!
! 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))
!
! Harustek had also merged the ECMWF and GFS obukhov functions, and the
! new "merged" parameter list looked something like
!
! 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)
!
! For the ECMWF call, the akz(llev) argument was problematic, and the
! logic behind the argument lists was confusing and not documented. I've
! tried to resolve this by creating two new variables, gfs_dummy_arg and
! ecmwf_dummy_arg, and using those where appropriate in the call to the
! obukhov function
!
if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
! NCEP version: find first level above ground
llev = 0
......@@ -136,11 +179,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), &
& gfs_dummy_arg, gfs_dummy_arg, 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,akzdummy,metdata_format)
tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm, &
& ecmwf_dummy_arg, metdata_format)
end if
if (ol.ne.0.) then
......@@ -162,8 +207,8 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
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), &
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)
else
call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
......
src/gridcheck_ecmwf.f90
View file @ 0eda0087
......@@ -41,6 +41,10 @@ subroutine gridcheck_ecmwf
! Marian Harustak, 12.5.2017 *
! - Renamed from gridcheck to gridcheck_ecmwf *
! *
! Implementation of the Vtables approach *
! D. Morton, D. Arnold 17.11.2017 *
! - Inclusion of specific code and usage of class_vtable *
! *
!**********************************************************************
! *
! DESCRIPTION: *
......@@ -77,6 +81,7 @@ subroutine gridcheck_ecmwf
use com_mod
use conv_mod
use cmapf_mod, only: stlmbr,stcm2p
use class_vtable
implicit none
......@@ -100,7 +105,8 @@ subroutine gridcheck_ecmwf
! dimension of zsec2 at least (10+nn), where nn is the number of vertical
! coordinate parameters
integer :: isec1(56),isec2(22+nxmax+nymax)
!!! DJM -- isec1() no longer needed -- integer :: isec1(56),isec2(22+nxmax+nymax)
integer :: isec2(22+nxmax+nymax)
real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
character(len=1) :: opt
......@@ -108,6 +114,19 @@ subroutine gridcheck_ecmwf
character(len=24) :: gribErrorMsg = 'Error reading grib file'
character(len=20) :: gribFunction = 'gridcheck'
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! Vtable related variables
!
! Path to Vtable - current implementation assumes it's in cwd, named
! "Vtable"
CHARACTER(LEN=255), PARAMETER :: VTABLE_PATH = "Vtable"
CHARACTER(LEN=15) :: fpname ! stores FLEXPART name for curr grib mesg.
TYPE(Vtable) :: my_vtable ! unallocated
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! DJM
INTEGER current_grib_level ! this was isec1(8) in previous versions
iumax=0
iwmax=0
......@@ -117,11 +136,36 @@ subroutine gridcheck_ecmwf
else
ifn=numbwf
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! Vtable code
PRINT *, 'Loading Vtable: ', VTABLE_PATH
call vtable_load_by_name(VTABLE_PATH, my_vtable)
!! Debugging tool
!PRINT *, 'Dump of Vtable...'
!call vtable_dump_records(my_vtable)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!! VTABLE code
! This is diagnostic/debugging code, and will normally be commented out.
! It's purpose is to look at the provided grib file and produce an
! inventory of the FP-related messages, relative to the Vtable that's
! already been open.
!CALL vtable_gribfile_inventory(path(3)(1:length(3)) // trim(wfname(ifn)), &
!& my_vtable)
!
!!!!!!!!!!!!!!!!!!! VTABLE code
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! OPENING OF DATA FILE (GRIB CODE)
!
5 call grib_open_file(ifile,path(3)(1:length(3)) &
//trim(wfname(ifn)),'r',iret)
5 call grib_open_file(ifile,path(3)(1:length(3)) &
//trim(wfname(ifn)),'r',iret)
if (iret.ne.GRIB_SUCCESS) then
goto 999 ! ERROR DETECTED
endif
......@@ -130,7 +174,7 @@ subroutine gridcheck_ecmwf
gotGrid=0
ifield=0
10 ifield=ifield+1
10 ifield=ifield+1
!
! GET NEXT FIELDS
......@@ -142,138 +186,41 @@ subroutine gridcheck_ecmwf
goto 999 ! ERROR DETECTED
endif
!first see if we read GRIB1 or GRIB2
call grib_get_int(igrib,'editionNumber',gribVer,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
if (gribVer.eq.1) then ! GRIB Edition 1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!! VTABLE code
! Get the fpname
fpname = vtable_get_fpname(igrib, my_vtable)
!print *, 'fpname: ', trim(fpname)
!print*,'GRiB Edition 1'
!read the grib2 identifiers
call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'level',isec1(8),iret)
call grib_check(iret,gribFunction,gribErrorMsg)
!change code for etadot to code for omega
if (isec1(6).eq.77) then
isec1(6)=135
endif
!!!!!!!!!!!!!!!!!!! VTABLE code
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!print*,isec1(6),isec1(8)
else
!print*,'GRiB Edition 2'
!read the grib2 identifiers
call grib_get_int(igrib,'discipline',discipl,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'parameterCategory',parCat,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'parameterNumber',parNum,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'level',valSurf,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
call grib_get_int(igrib,'paramId',parId,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
!first see if we read GRIB1 or GRIB2
call grib_get_int(igrib,'editionNumber',gribVer,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
!print*,discipl,parCat,parNum,typSurf,valSurf
!convert to grib1 identifiers
isec1(6)=-1
isec1(7)=-1
isec1(8)=-1
isec1(8)=valSurf ! level
if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
isec1(6)=130 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
isec1(6)=131 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
isec1(6)=132 ! indicatorOfParameter
elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
isec1(6)=133 ! indicatorOfParameter
!ZHG FOR CLOUDS FROM GRIB
elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
isec1(6)=246 ! indicatorOfParameter
elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
isec1(6)=247 ! indicatorOfParameter
!ZHG end
! ESO qc(=clwc+ciwc)
elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
isec1(6)=201031 ! indicatorOfParameter
elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
isec1(6)=134 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
isec1(6)=135 ! indicatorOfParameter
elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
isec1(6)=135 ! indicatorOfParameter
elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
isec1(6)=151 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
isec1(6)=165 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
isec1(6)=166 ! indicatorOfParameter
elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
isec1(6)=167 ! indicatorOfParameter
elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
isec1(6)=168 ! indicatorOfParameter
elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
isec1(6)=141 ! indicatorOfParameter
elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
isec1(6)=164 ! indicatorOfParameter
elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
isec1(6)=142 ! indicatorOfParameter
elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
isec1(6)=143 ! indicatorOfParameter
elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
isec1(6)=146 ! indicatorOfParameter
elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
isec1(6)=176 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
isec1(6)=180 ! indicatorOfParameter
elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
isec1(6)=181 ! indicatorOfParameter
elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
isec1(6)=129 ! indicatorOfParameter
elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
isec1(6)=160 ! indicatorOfParameter
elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
(typSurf.eq.1)) then ! LSM
isec1(6)=172 ! indicatorOfParameter
else
print*,'***ERROR: undefined GRiB2 message found!',discipl, &
parCat,parNum,typSurf
endif
if(parId .ne. isec1(6) .and. parId .ne. 77) then
write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
! stop
endif
endif
!! DJM - get the current_grib_level (in previous code it was isec1(8))
!! It's the same in both GRIB1 and GRIB2
call grib_get_int(igrib,'level',current_grib_level,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
CALL grib_get_int(igrib,'numberOfPointsAlongAParallel', &
isec2(2),iret)
! nx=isec2(2)
! WRITE(*,*) nx,nxmax
IF (isec2(2).GT.nxmax) THEN
WRITE(*,*) 'FLEXPART error: Too many grid points in x direction.'
WRITE(*,*) 'Reduce resolution of wind fields.'
WRITE(*,*) 'Or change parameter settings in file ecmwf_mod.'
WRITE(*,*) nx,nxmax
! STOP
ENDIF
!get the size and data of the values array
if (isec1(6).ne.-1) then
!!!! -- original statement -- if (isec1(6).ne.-1) then
! 'NOFP' is the fpname for a grib message not recognized by FLEXPART
IF (TRIM(fpname) .NE. 'NOFP') THEN
call grib_get_real4_array(igrib,'values',zsec4,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
endif
ENDIF
if (ifield.eq.1) then
!HSO get the required fields from section 2 in a gribex compatible manner
!HSO get the required fields from section 2 in a gribex compatible manner
call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
isec2(2),iret)
call grib_check(iret,gribFunction,gribErrorMsg)
......@@ -287,17 +234,19 @@ subroutine gridcheck_ecmwf
isec2(12),iret)
call grib_check(iret,gribFunction,gribErrorMsg)
! get the size and data of the vertical coordinate array
! get the size and data of the vertical coordinate array
call grib_get_real4_array(igrib,'pv',zsec2,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
nxfield=isec2(2)
ny=isec2(3)
nlev_ec=isec2(12)/2-1
endif
endif
!HSO get the second part of the grid dimensions only from GRiB1 messages
if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
if (TRIM(fpname) .EQ. 'T2'.and. (gotGrid.eq.0)) then
!!!!! DJM --- if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
xaux2in,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
......@@ -311,6 +260,8 @@ subroutine gridcheck_ecmwf
xaux2=xaux2in
yaux1=yaux1in
yaux2=yaux2in
if (xaux1.gt.180.) xaux1=xaux1-360.0
if (xaux2.gt.180.) xaux2=xaux2-360.0
if (xaux1.lt.-180.) xaux1=xaux1+360.0
......@@ -323,10 +274,10 @@ subroutine gridcheck_ecmwf
dxconst=180./(dx*r_earth*pi)
dyconst=180./(dy*r_earth*pi)
gotGrid=1
! Check whether fields are global
! If they contain the poles, specify polar stereographic map
! projections using the stlmbr- and stcm2p-calls
!***********************************************************
! Check whether fields are global
! If they contain the poles, specify polar stereographic map
! projections using the stlmbr- and stcm2p-calls
!***********************************************************
xauxa=abs(xaux2+dx-360.-xaux1)
if (xauxa.lt.0.001) then
......@@ -347,8 +298,8 @@ subroutine gridcheck_ecmwf
xauxa=abs(yaux1+90.)
if (xglobal.and.xauxa.lt.0.001) then
sglobal=.true. ! field contains south pole
! Enhance the map scale by factor 3 (*2=6) compared to north-south
! map scale
! Enhance the map scale by factor 3 (*2=6) compared to north-south
! map scale
sizesouth=6.*(switchsouth+90.)/dy
call stlmbr(southpolemap,-90.,0.)
call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
......@@ -361,8 +312,8 @@ subroutine gridcheck_ecmwf
xauxa=abs(yaux2-90.)
if (xglobal.and.xauxa.lt.0.001) then
nglobal=.true. ! field contains north pole
! Enhance the map scale by factor 3 (*2=6) compared to north-south
! map scale
! Enhance the map scale by factor 3 (*2=6) compared to north-south
! map scale
sizenorth=6.*(90.-switchnorth)/dy
call stlmbr(northpolemap,90.,0.)
call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
......@@ -378,7 +329,7 @@ subroutine gridcheck_ecmwf
endif ! gotGrid
if (nx.gt.nxmax) then
write(*,*) 'FLEXPART error: Too many grid points in x direction.'
write(*,*) 'FLEXPART error: Too many grid points in x direction.'
write(*,*) 'Reduce resolution of wind fields.'
write(*,*) 'Or change parameter settings in file par_mod.'
write(*,*) nx,nxmax
......@@ -386,32 +337,45 @@ subroutine gridcheck_ecmwf
endif
if (ny.gt.nymax) then
write(*,*) 'FLEXPART error: Too many grid points in y direction.'
write(*,*) 'FLEXPART error: Too many grid points in y direction.'
write(*,*) 'Reduce resolution of wind fields.'
write(*,*) 'Or change parameter settings in file par_mod.'
write(*,*) ny,nymax
stop
endif
k=isec1(8)
if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
if(isec1(6).eq.129) then
!!!!!!!! DJM - orig - k=isec1(8)
k = current_grib_level
!!!!!!! DJM - orig - Iif(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
IF (TRIM(fpname) .EQ. 'UU') iumax=max(iumax,nlev_ec-k+1)
!!!!!!! DJM - orig - if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
IF (TRIM(fpname) .EQ. 'ETADOT') iwmax=max(iwmax,nlev_ec-k+1)
!!!!!!! DJM - orig - if(isec1(6).eq.129) then
IF (TRIM(fpname) .EQ. 'ORO') THEN
do jy=0,ny-1
do ix=0,nxfield-1
oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga
end do
end do
endif
if(isec1(6).eq.172) then
!!!!!!! DJM - orig - if(isec1(6).eq.172) then
IF (TRIM(fpname) .EQ. 'LSM') THEN
do jy=0,ny-1
do ix=0,nxfield-1
lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
end do
end do
endif
if(isec1(6).eq.160) then
!!!!!! DJM - orig - if(isec1(6).eq.160) then
IF (TRIM(fpname) .EQ. 'EXCESSORO') THEN
do jy=0,ny-1
do ix=0,nxfield-1
excessoro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
......@@ -425,7 +389,7 @@ subroutine gridcheck_ecmwf
! CLOSING OF INPUT DATA FILE
!
30 call grib_close_file(ifile)
30 call grib_close_file(ifile)
!error message if no fields found with correct first longitude in it
if (gotGrid.eq.0) then
......@@ -484,7 +448,7 @@ subroutine gridcheck_ecmwf
! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
numskip=nlev_ec-nuvz ! number of ecmwf model layers not used
! by trajectory model
! by trajectory model
!do 8940 i=1,244
! write (*,*) 'zsec2:',i,ifield,zsec2(i),numskip
!940 continue
......@@ -497,7 +461,7 @@ subroutine gridcheck_ecmwf
j=numskip+i
k=nlev_ec+1+numskip+i
akm(nwz-i+1)=zsec2(j)
! write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
! write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
bkm(nwz-i+1)=zsec2(k)
end do
......@@ -513,8 +477,8 @@ subroutine gridcheck_ecmwf
akz(1)=0.
bkz(1)=1.0
do i=1,nuvz
akz(i+1)=0.5*(akm(i+1)+akm(i))
bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
akz(i+1)=0.5*(akm(i+1)+akm(i))
bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
end do
nuvz=nuvz+1
......@@ -554,7 +518,7 @@ subroutine gridcheck_ecmwf
pint=akz(i)+bkz(i)*101325.
if (pint.lt.5000.) goto 96
end do
96 nconvlev=i
96 nconvlev=i
if (nconvlev.gt.nconvlevmax-1) then
nconvlev=nconvlevmax-1
write(*,*) 'Attention, convection only calculated up to ', &
......@@ -563,7 +527,7 @@ subroutine gridcheck_ecmwf
return
999 write(*,*)
999 write(*,*)
write(*,*) ' ###########################################'// &
'###### '
write(*,*) ' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
......
src/gridcheck_gfs.f90
View file @ 0eda0087
......@@ -42,6 +42,11 @@ subroutine gridcheck_gfs
! Marian Harustak, 12.5.2017 *
! - Renamed routine from gridcheck to gridcheck_gfs *
! *
! Implementation of the Vtables approach *
! D. Morton, D. Arnold 17.11.2017 *
! - Inclusion of specific code and usage of class_vtable *
! *
! *
!**********************************************************************
! *
! DESCRIPTION: *
......@@ -78,6 +83,7 @@ subroutine gridcheck_gfs
use com_mod
use conv_mod
use cmapf_mod, only: stlmbr,stcm2p
use class_vtable
implicit none
......@@ -101,7 +107,10 @@ subroutine gridcheck_gfs
! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
integer :: isec1(8),isec2(3)
!!!! DJM - isec1 no longer user - integer :: isec1(8),isec2(3)
integer :: isec2(3)
real(kind=4) :: zsec4(jpunp)
character(len=1) :: opt
......@@ -109,6 +118,23 @@ subroutine gridcheck_gfs
character(len=24) :: gribErrorMsg = 'Error reading grib file'
character(len=20) :: gribFunction = 'gridcheckwind_gfs'
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! Vtable related variables
!
! Path to Vtable - current implementation assumes it's in cwd, named
! "Vtable"
CHARACTER(LEN=255), PARAMETER :: VTABLE_PATH = "Vtable"
CHARACTER(LEN=15) :: fpname ! stores FLEXPART name for curr grib mesg.
TYPE(Vtable) :: my_vtable ! unallocated
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! DJM
INTEGER current_grib_level ! this was isec1(8) in previous versions
if (numbnests.ge.1) then
write(*,*) ' ###########################################'
write(*,*) ' FLEXPART ERROR SUBROUTINE GRIDCHECK:'
......@@ -125,6 +151,37 @@ subroutine gridcheck_gfs
else
ifn=numbwf
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! Vtable code
PRINT *, 'Loading Vtable: ', VTABLE_PATH