Skip to content

GitLab

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

First commit of CTBTO vtables

parent f4a43397
Expand all Hide whitespace changes
Inline Side-by-side
src/calcpar.f90
View file @ e6111cbc
...@@ -32,7 +32,7 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format) ...@@ -32,7 +32,7 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! 21 May 1995 * ! 21 May 1995 *
! * ! *
! ------------------------------------------------------------------ * ! ------------------------------------------------------------------ *
! Petra Seibert, Feb 2000: * ! Petra Seibert, Feb 2000: *
! convection scheme: * ! convection scheme: *
! new variables in call to richardson * ! new variables in call to richardson *
! * ! *
...@@ -45,6 +45,12 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format) ...@@ -45,6 +45,12 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! Marian Harustak, 12.5.2017 * ! Marian Harustak, 12.5.2017 *
! - Merged calcpar and calcpar_gfs into one routine using if-then * ! - Merged calcpar and calcpar_gfs into one routine using if-then *
! for meteo-type dependent code * ! 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) ...@@ -75,12 +81,19 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
integer :: n,ix,jy,i,kz,lz,kzmin,llev,loop_start 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),hmixdummy,akzdummy 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)
real,parameter :: const=r_air/ga 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' !write(*,*) 'in calcpar writting snowheight'
!*********************************** !***********************************
! for test: write out snow depths ! for test: write out snow depths
...@@ -124,6 +137,36 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format) ...@@ -124,6 +137,36 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! 2) Calculation of inverse Obukhov length scale ! 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 if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
! NCEP version: find first level above ground ! NCEP version: find first level above ground
llev = 0 llev = 0
...@@ -136,11 +179,13 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format) ...@@ -136,11 +179,13 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
! calculate inverse Obukhov length scale with tth(llev) ! 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,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 else
llev=0 llev=0
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,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 end if
if (ol.ne.0.) then if (ol.ne.0.) then
...@@ -162,8 +207,8 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format) ...@@ -162,8 +207,8 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
! NCEP version hmix has been read in in readwind.f, is therefore not calculated here ! 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),hmixdummy,wstar(ix,jy,1,n),hmixplus,metdata_format) td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus,metdata_format)
else else
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, &
......
src/gridcheck_ecmwf.f90
View file @ e6111cbc
...@@ -41,6 +41,10 @@ subroutine gridcheck_ecmwf ...@@ -41,6 +41,10 @@ subroutine gridcheck_ecmwf
! Marian Harustak, 12.5.2017 * ! Marian Harustak, 12.5.2017 *
! - Renamed from gridcheck to gridcheck_ecmwf * ! - 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: * ! DESCRIPTION: *
...@@ -77,6 +81,7 @@ subroutine gridcheck_ecmwf ...@@ -77,6 +81,7 @@ subroutine gridcheck_ecmwf
use com_mod use com_mod
use conv_mod use conv_mod
use cmapf_mod, only: stlmbr,stcm2p use cmapf_mod, only: stlmbr,stcm2p
use class_vtable
implicit none implicit none
...@@ -100,7 +105,8 @@ subroutine gridcheck_ecmwf ...@@ -100,7 +105,8 @@ subroutine gridcheck_ecmwf
! dimension of zsec2 at least (10+nn), where nn is the number of vertical ! dimension of zsec2 at least (10+nn), where nn is the number of vertical
! coordinate parameters ! 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) real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
character(len=1) :: opt character(len=1) :: opt
...@@ -108,6 +114,19 @@ subroutine gridcheck_ecmwf ...@@ -108,6 +114,19 @@ subroutine gridcheck_ecmwf
character(len=24) :: gribErrorMsg = 'Error reading grib file' character(len=24) :: gribErrorMsg = 'Error reading grib file'
character(len=20) :: gribFunction = 'gridcheck' 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 iumax=0
iwmax=0 iwmax=0
...@@ -117,11 +136,36 @@ subroutine gridcheck_ecmwf ...@@ -117,11 +136,36 @@ subroutine gridcheck_ecmwf
else else
ifn=numbwf ifn=numbwf
endif 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) ! OPENING OF DATA FILE (GRIB CODE)
! !
5 call grib_open_file(ifile,path(3)(1:length(3)) & 5 call grib_open_file(ifile,path(3)(1:length(3)) &
//trim(wfname(ifn)),'r',iret) //trim(wfname(ifn)),'r',iret)
if (iret.ne.GRIB_SUCCESS) then if (iret.ne.GRIB_SUCCESS) then
goto 999 ! ERROR DETECTED goto 999 ! ERROR DETECTED
endif endif
...@@ -130,7 +174,7 @@ subroutine gridcheck_ecmwf ...@@ -130,7 +174,7 @@ subroutine gridcheck_ecmwf
gotGrid=0 gotGrid=0
ifield=0 ifield=0
10 ifield=ifield+1 10 ifield=ifield+1
! !
! GET NEXT FIELDS ! GET NEXT FIELDS
...@@ -142,138 +186,41 @@ subroutine gridcheck_ecmwf ...@@ -142,138 +186,41 @@ subroutine gridcheck_ecmwf
goto 999 ! ERROR DETECTED goto 999 ! ERROR DETECTED
endif 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 !!!!!!!!!!!!!!!!!!! VTABLE code
if (isec1(6).eq.77) then !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
isec1(6)=135
endif
!print*,isec1(6),isec1(8)
else
!print*,'GRiB Edition 2' !first see if we read GRIB1 or GRIB2
!read the grib2 identifiers call grib_get_int(igrib,'editionNumber',gribVer,iret)
call grib_get_int(igrib,'discipline',discipl,iret) call grib_check(iret,gribFunction,gribErrorMsg)
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)
!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 !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_get_real4_array(igrib,'values',zsec4,iret)
call grib_check(iret,gribFunction,gribErrorMsg) call grib_check(iret,gribFunction,gribErrorMsg)
endif ENDIF
if (ifield.eq.1) then 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', & call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
isec2(2),iret) isec2(2),iret)
call grib_check(iret,gribFunction,gribErrorMsg) call grib_check(iret,gribFunction,gribErrorMsg)
...@@ -287,17 +234,19 @@ subroutine gridcheck_ecmwf ...@@ -287,17 +234,19 @@ subroutine gridcheck_ecmwf
isec2(12),iret) isec2(12),iret)
call grib_check(iret,gribFunction,gribErrorMsg) 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_get_real4_array(igrib,'pv',zsec2,iret)
call grib_check(iret,gribFunction,gribErrorMsg) call grib_check(iret,gribFunction,gribErrorMsg)
nxfield=isec2(2) nxfield=isec2(2)
ny=isec2(3) ny=isec2(3)
nlev_ec=isec2(12)/2-1 nlev_ec=isec2(12)/2-1
endif endif
!HSO get the second part of the grid dimensions only from GRiB1 messages !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', & call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
xaux2in,iret) xaux2in,iret)
call grib_check(iret,gribFunction,gribErrorMsg) call grib_check(iret,gribFunction,gribErrorMsg)
...@@ -311,6 +260,8 @@ subroutine gridcheck_ecmwf ...@@ -311,6 +260,8 @@ subroutine gridcheck_ecmwf
xaux2=xaux2in xaux2=xaux2in
yaux1=yaux1in yaux1=yaux1in
yaux2=yaux2in yaux2=yaux2in
if (xaux1.gt.180.) xaux1=xaux1-360.0 if (xaux1.gt.180.) xaux1=xaux1-360.0
if (xaux2.gt.180.) xaux2=xaux2-360.0 if (xaux2.gt.180.) xaux2=xaux2-360.0
if (xaux1.lt.-180.) xaux1=xaux1+360.0 if (xaux1.lt.-180.) xaux1=xaux1+360.0
...@@ -323,10 +274,10 @@ subroutine gridcheck_ecmwf ...@@ -323,10 +274,10 @@ subroutine gridcheck_ecmwf
dxconst=180./(dx*r_earth*pi) dxconst=180./(dx*r_earth*pi)
dyconst=180./(dy*r_earth*pi) dyconst=180./(dy*r_earth*pi)
gotGrid=1 gotGrid=1
! Check whether fields are global ! Check whether fields are global
! If they contain the poles, specify polar stereographic map ! If they contain the poles, specify polar stereographic map
! projections using the stlmbr- and stcm2p-calls ! projections using the stlmbr- and stcm2p-calls
!*********************************************************** !***********************************************************
xauxa=abs(xaux2+dx-360.-xaux1) xauxa=abs(xaux2+dx-360.-xaux1)
if (xauxa.lt.0.001) then if (xauxa.lt.0.001) then
...@@ -347,8 +298,8 @@ subroutine gridcheck_ecmwf ...@@ -347,8 +298,8 @@ subroutine gridcheck_ecmwf
xauxa=abs(yaux1+90.) xauxa=abs(yaux1+90.)
if (xglobal.and.xauxa.lt.0.001) then if (xglobal.and.xauxa.lt.0.001) then
sglobal=.true. ! field contains south pole sglobal=.true. ! field contains south pole
! Enhance the map scale by factor 3 (*2=6) compared to north-south ! Enhance the map scale by factor 3 (*2=6) compared to north-south
! map scale ! map scale
sizesouth=6.*(switchsouth+90.)/dy sizesouth=6.*(switchsouth+90.)/dy
call stlmbr(southpolemap,-90.,0.) call stlmbr(southpolemap,-90.,0.)
call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, & call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
...@@ -361,8 +312,8 @@ subroutine gridcheck_ecmwf ...@@ -361,8 +312,8 @@ subroutine gridcheck_ecmwf
xauxa=abs(yaux2-90.) xauxa=abs(yaux2-90.)
if (xglobal.and.xauxa.lt.0.001) then if (xglobal.and.xauxa.lt.0.001) then
nglobal=.true. ! field contains north pole nglobal=.true. ! field contains north pole
! Enhance the map scale by factor 3 (*2=6) compared to north-south ! Enhance the map scale by factor 3 (*2=6) compared to north-south
! map scale ! map scale
sizenorth=6.*(90.-switchnorth)/dy sizenorth=6.*(90.-switchnorth)/dy
call stlmbr(northpolemap,90.,0.) call stlmbr(northpolemap,90.,0.)
call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, & call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
...@@ -378,7 +329,7 @@ subroutine gridcheck_ecmwf ...@@ -378,7 +329,7 @@ subroutine gridcheck_ecmwf
endif ! gotGrid endif ! gotGrid
if (nx.gt.nxmax) then 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(*,*) 'Reduce resolution of wind fields.'
write(*,*) 'Or change parameter settings in file par_mod.' write(*,*) 'Or change parameter settings in file par_mod.'
write(*,*) nx,nxmax write(*,*) nx,nxmax
...@@ -386,32 +337,45 @@ subroutine gridcheck_ecmwf ...@@ -386,32 +337,45 @@ subroutine gridcheck_ecmwf
endif endif
if (ny.gt.nymax) then 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(*,*) 'Reduce resolution of wind fields.'
write(*,*) 'Or change parameter settings in file par_mod.' write(*,*) 'Or change parameter settings in file par_mod.'
write(*,*) ny,nymax write(*,*) ny,nymax
stop stop
endif 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 jy=0,ny-1
do ix=0,nxfield-1 do ix=0,nxfield-1
oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga
end do end do
end do end do