Commit 467460af authored by Espen Sollum's avatar Espen Sollum
Browse files

First commit of files from Hachinger

parent 759df5f2
...@@ -70,8 +70,8 @@ subroutine calcmatrix(lconv,delt,cbmf,metdata_format) ...@@ -70,8 +70,8 @@ subroutine calcmatrix(lconv,delt,cbmf,metdata_format)
do kuvz = 2,nuvz do kuvz = 2,nuvz
k = kuvz-1 k = kuvz-1
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then 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 else
phconv(kuvz) = 0.5*(pconv(kuvz)+pconv(k)) phconv(kuvz) = 0.5*(pconv(kuvz)+pconv(k))
endif endif
......
...@@ -13,7 +13,11 @@ real function ew(x) ...@@ -13,7 +13,11 @@ real function ew(x)
real :: x, y, a, c, d real :: x, y, a, c, d
ew=0. ew=0.
if(x.le.0.) stop 'sorry: t not in [k]' if(x.le.0.) then
write(*,*) 'in ew.f90: x=', x
stop 'sorry: t not in [k]'
end if
y=373.16/x y=373.16/x
a=-7.90298*(y-1.) a=-7.90298*(y-1.)
a=a+(5.02808*0.43429*alog(y)) a=a+(5.02808*0.43429*alog(y))
......
...@@ -74,12 +74,12 @@ subroutine gridcheck_gfs ...@@ -74,12 +74,12 @@ subroutine gridcheck_gfs
integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
real :: sizesouth,sizenorth,xauxa,pint real :: sizesouth,sizenorth,xauxa,pint
real :: akm_usort(nwzmax) real :: akm_usort(nwzmax)
real,parameter :: eps=0.0001 real,parameter :: eps=spacing(2.0_4*360.0_4)
! NCEP GFS ! NCEP GFS
real :: pres(nwzmax), help real :: pres(nwzmax), help
integer :: i179,i180,i181 integer :: i180
! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
...@@ -223,9 +223,9 @@ subroutine gridcheck_gfs ...@@ -223,9 +223,9 @@ subroutine gridcheck_gfs
nxfield=isec2(2) nxfield=isec2(2)
ny=isec2(3) ny=isec2(3)
if((abs(xaux1).lt.eps).and.(xaux2.ge.359)) then ! NCEP DATA FROM 0 TO if((abs(xaux1).lt.eps).and.(xaux2.ge.359.)) then ! NCEP DATA FROM 0 TO
xaux1=-179.0 ! 359 DEG EAST -> xaux1=-180.0 ! 359 DEG EAST ->
xaux2=-179.0+360.-360./real(nxfield) ! TRANSFORMED TO -179 xaux2=-180.0+360.-360./real(nxfield) ! TRANSFORMED TO -179
endif ! TO 180 DEG EAST endif ! TO 180 DEG EAST
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
...@@ -304,13 +304,7 @@ subroutine gridcheck_gfs ...@@ -304,13 +304,7 @@ subroutine gridcheck_gfs
endif endif
i179=nint(179./dx) i180=nint(180./dx) ! 0.5 deg data
if (dx.lt.0.7) then
i180=nint(180./dx)+1 ! 0.5 deg data
else
i180=nint(179./dx)+1 ! 1 deg data
endif
i181=i180+1
! NCEP TERRAIN ! NCEP TERRAIN
...@@ -320,12 +314,12 @@ subroutine gridcheck_gfs ...@@ -320,12 +314,12 @@ subroutine gridcheck_gfs
do jy=0,ny-1 do jy=0,ny-1
do ix=0,nxfield-1 do ix=0,nxfield-1
help=zsec4(nxfield*(ny-jy-1)+ix+1) help=zsec4(nxfield*(ny-jy-1)+ix+1)
if(ix.le.i180) then if(ix.lt.i180) then
oro(i179+ix,jy)=help oro(i180+ix,jy)=help
excessoro(i179+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED excessoro(i180+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
else else
oro(ix-i181,jy)=help oro(ix-i180,jy)=help
excessoro(ix-i181,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED excessoro(ix-i180,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
endif endif
end do end do
end do end do
...@@ -338,10 +332,10 @@ subroutine gridcheck_gfs ...@@ -338,10 +332,10 @@ subroutine gridcheck_gfs
do jy=0,ny-1 do jy=0,ny-1
do ix=0,nxfield-1 do ix=0,nxfield-1
help=zsec4(nxfield*(ny-jy-1)+ix+1) help=zsec4(nxfield*(ny-jy-1)+ix+1)
if(ix.le.i180) then if(ix.lt.i180) then
lsm(i179+ix,jy)=help lsm(i180+ix,jy)=help
else else
lsm(ix-i181,jy)=help lsm(ix-i180,jy)=help
endif endif
end do end do
end do end do
...@@ -412,7 +406,7 @@ subroutine gridcheck_gfs ...@@ -412,7 +406,7 @@ subroutine gridcheck_gfs
if (lroot) then if (lroot) then
write(*,*) write(*,*)
write(*,*) write(*,*)
write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', & write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
nuvz,nwz nuvz,nwz
write(*,*) write(*,*)
write(*,'(a)') 'Mother domain:' write(*,'(a)') 'Mother domain:'
......
...@@ -145,7 +145,9 @@ module par_mod ...@@ -145,7 +145,9 @@ module par_mod
! integer,parameter :: nxmax=181,nymax=91,nuvzmax=92,nwzmax=92,nzmax=92,nxshift=0 ! CERA 2.0 degree 92 level ! integer,parameter :: nxmax=181,nymax=91,nuvzmax=92,nwzmax=92,nzmax=92,nxshift=0 ! CERA 2.0 degree 92 level
! GFS ! GFS
integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 ! integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138
! GFS 0.25
integer,parameter :: nxmax=1441,nymax=721,nuvzmax=138,nwzmax=138,nzmax=138
integer :: nxshift=0 ! shift not fixed for the executable integer :: nxshift=0 ! shift not fixed for the executable
...@@ -211,7 +213,7 @@ module par_mod ...@@ -211,7 +213,7 @@ module par_mod
! Maximum number of particles, species, and similar ! Maximum number of particles, species, and similar
!************************************************** !**************************************************
integer,parameter :: maxpart=100000 integer,parameter :: maxpart=10000000
integer,parameter :: maxspec=6 integer,parameter :: maxspec=6
real,parameter :: minmass=0.0001 real,parameter :: minmass=0.0001
...@@ -255,7 +257,7 @@ module par_mod ...@@ -255,7 +257,7 @@ module par_mod
! Dimension of random number field ! Dimension of random number field
!********************************* !*********************************
integer,parameter :: maxrand=1000000 integer,parameter :: maxrand=10000000
! maxrand number of random numbers used ! maxrand number of random numbers used
......
...@@ -545,9 +545,9 @@ subroutine readreleases ...@@ -545,9 +545,9 @@ subroutine readreleases
if (ldirect.eq.1) then if (ldirect.eq.1) then
if ((jul1.lt.bdate).or.(jul2.gt.edate)) then if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
write(*,*) 'FLEXPART MODEL ERROR' write(*,*) 'FLEXPART MODEL ERROR'
write(*,*) 'Release starts before simulation begins or ends' write(*,*) 'Release starts before simulation begins or ends (1)'
write(*,*) 'after simulation stops.' write(*,*) 'after simulation stops.'
write(*,*) 'Make files COMMAND and RELEASES consistent.' write(*,*) 'Make files COMMAND and RELEASES consistent (1).'
stop stop
endif endif
if (npart(numpoint).gt.num_min_discrete) then if (npart(numpoint).gt.num_min_discrete) then
...@@ -560,9 +560,9 @@ subroutine readreleases ...@@ -560,9 +560,9 @@ subroutine readreleases
else if (ldirect.eq.-1) then else if (ldirect.eq.-1) then
if ((jul1.lt.edate).or.(jul2.gt.bdate)) then if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
write(*,*) 'FLEXPART MODEL ERROR' write(*,*) 'FLEXPART MODEL ERROR'
write(*,*) 'Release starts before simulation begins or ends' write(*,*) 'Release starts before simulation begins or ends (2)'
write(*,*) 'after simulation stops.' write(*,*) 'after simulation stops.'
write(*,*) 'Make files COMMAND and RELEASES consistent.' write(*,*) 'Make files COMMAND and RELEASES consistent (2).'
stop stop
endif endif
if (npart(numpoint).gt.num_min_discrete) then if (npart(numpoint).gt.num_min_discrete) then
......
This diff is collapsed.
...@@ -140,11 +140,15 @@ subroutine richardson(psurf,ust,ttlev,qvlev,ulev,vlev,nuvz, & ...@@ -140,11 +140,15 @@ subroutine richardson(psurf,ust,ttlev,qvlev,ulev,vlev,nuvz, &
thetaold=theta thetaold=theta
zold=z zold=z
end do end do
k=k-1 ! ESO: make sure k <= nuvz (ticket #139) ! k=k-1 ! ESO: make sure k <= nuvz (ticket #139)
20 continue 20 continue
! Determine Richardson number between the critical levels ! Determine Richardson number between the critical levels
!******************************************************** !********************************************************
! JMA: It may happen that k >= nuvz:
! JMA: In that case: k is set to nuvz
if (k .gt. nuvz) k = nuvz ! JMA
zl1=zold zl1=zold
theta1=thetaold theta1=thetaold
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment