diff --git a/flexpart_code/calcmatrix_gfs.f90 b/flexpart_code/calcmatrix_gfs.f90
index ff08a40..31eeb1f 100644
--- a/flexpart_code/calcmatrix_gfs.f90
+++ b/flexpart_code/calcmatrix_gfs.f90
@@ -74,21 +74,16 @@ subroutine calcmatrix(lconv,delt,cbmf)
! not given
phconv(1) = psconv
-
-
- ! JMA / SH Bugfix phconv was set in loop with access to undefined pconv(nuvz)
- phconv(1) = psconv
- phconv(2:nuvz-1) = 0.5*(pconv(2:nuvz-1)+pconv(1:nuvz-2))
- phconv(nuvz) = pconv(nuvz-1)
-
- dpr(1:nuvz-1) = phconv(1:nuvz-1)-phconv(2:nuvz)
-
- do k = 1,nuvz-1
+ do kuvz = 2,nuvz
+ k = kuvz-1
+ phconv(kuvz) = 0.5*(pconv(kuvz)+pconv(k))
+ dpr(k) = phconv(k) - phconv(kuvz)
qsconv(k) = f_qvsat( pconv(k), tconv(k) )
- end do
-
! initialize mass fractions
- fmassfrac(1:nuvz-1,1:nconvlev)=0.0
+ do kk=1,nconvlev
+ fmassfrac(k,kk)=0.
+ enddo
+ end do
!note that Emanuel says it is important
!a. to set this =0. every grid point
diff --git a/flexpart_code/calcpv.f90 b/flexpart_code/calcpv.f90
index ceab700..d3f57a0 100644
--- a/flexpart_code/calcpv.f90
+++ b/flexpart_code/calcpv.f90
@@ -1,4 +1,3 @@
-
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
@@ -20,295 +19,319 @@
! along with FLEXPART. If not, see . *
!**********************************************************************
- subroutine calcpv(n,uuh,vvh,pvh)
- ! i i i o
- !*****************************************************************************
- ! *
- ! Calculation of potential vorticity on 3-d grid. *
- ! *
- ! Author: P. James *
- ! 3 February 2000 *
- ! *
- ! Adaptation to FLEXPART, A. Stohl, 1 May 2000 *
- ! *
- !*****************************************************************************
- ! *
- ! Variables: *
- ! n temporal index for meteorological fields (1 to 2) *
- ! *
- ! Constants: *
- ! *
- !*****************************************************************************
-
- use par_mod
- use com_mod
- use omp_lib
-
- implicit none
-
- integer :: n,ix,jy,i,j,k,kl,jj,klvrp,klvrm,klpt,kup,kdn,kch
- integer :: jux,juy,ivrm,ivrp,ivr
- integer :: nlck
- !SH additions/replacements
- integer :: khi,klo
- integer :: jyvp_arr(0:nymin1),jyvm_arr(0:nymin1),jumpy_arr(0:nymin1)
- integer :: jybords(2),ixvp_arr(0:nxmin1),ixvm_arr(0:nxmin1),jumpx_arr(0:nxmin1)
- !end SH
- real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
- real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
- real :: pvavr
- real :: thup,thdn
- !SH additions
- real :: th_pre(0:nxmax-1,0:nymax-1,nuvzmax)
- real :: ppmstore(0:nxmin1,0:nymin1,1:nuvz)
- !end SH
- !PARAMETERS
- real,parameter :: eps=1.e-5
- !DUMMY ARGUMENTS
- 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)
-
- write(*,*) 'enter calcpv parallelised...'
- pvh(:,:,:)=0.0_8
-
- ! Set number of levels to check for adjacent theta
- nlck=nuvz/3
-
- ! *** Precalculate all theta levels for efficiency
- do jy=0,nymin1
- do ix=0,nxmin1
- ppmstore(ix,jy,1:nuvz) = akz(1:nuvz)+bkz(1:nuvz)*ps(ix,jy,1,n)
- end do
- end do
- !write(*,*) 'calcpv kappa exponentiation with kappa=',kappa
+subroutine calcpv(n,uuh,vvh,pvh)
+ ! i i i o
+ !*****************************************************************************
+ ! *
+ ! Calculation of potential vorticity on 3-d grid. *
+ ! *
+ ! Author: P. James *
+ ! 3 February 2000 *
+ ! *
+ ! Adaptation to FLEXPART, A. Stohl, 1 May 2000 *
+ ! *
+ !*****************************************************************************
+ ! *
+ ! Variables: *
+ ! n temporal index for meteorological fields (1 to 2) *
+ ! *
+ ! Constants: *
+ ! *
+ !*****************************************************************************
-!$omp parallel default(none) shared(th_pre,tth,ppmstore, &
-!$omp nxmin1,nymin1,nuvz,n) num_threads(omp_get_max_threads())
-!$omp workshare
- th_pre(0:nxmin1,0:nymin1,1:nuvz)=tth(0:nxmin1,0:nymin1,1:nuvz,n)* &
- (100000./ppmstore(0:nxmin1,0:nymin1,1:nuvz))**kappa
-!$omp end workshare
-!$omp end parallel
- write(*,*) 'calcpv rest...'
-
- jybords=[0,nymin1]
- if (sglobal) jybords(1)=1
- if (nglobal) jybords(2)=nymin1-1
- jyvp_arr(jybords(1):jybords(2)) = [ (jy,jy=jybords(1)+1,jybords(2)+1) ]
- jyvm_arr(jybords(1):jybords(2)) = [ (jy,jy=jybords(1)-1,jybords(2)-1) ]
- jyvp_arr(jybords(2))=jybords(2)
- jyvm_arr(jybords(1))=jybords(1)
- jumpy_arr=2
- jumpy_arr(jybords(1:2))=1
+ use par_mod
+ use com_mod
- ixvp_arr(0:nxmin1) = [ (ix,ix=1,nxmin1+1) ]
- ixvm_arr(0:nxmin1) = [ (ix,ix=-1,nxmin1-1) ]
- jumpx_arr=2
- if (.not. xglobal) then
- ixvp_arr(nxmin1)=nxmin1
- ixvm_arr(0)=0
- jumpx_arr([0,nxmin1])=1
- else
- ixvp_arr(nxmin1)=1
- ixvm_arr(0)=nxmin1-1
- end if
+ implicit none
- ! Loop over entire grid
- !**********************
-!$omp parallel default(none) shared(jybords,jumpx_arr,jumpy_arr, &
-!$omp ixvp_arr,ixvm_arr,jyvp_arr,jyvm_arr,th_pre,ppmstore, &
-!$omp vvh,pvh,nuvz,nxmin1,ylat0,dx,dy,n,nlck,xglobal,uuh) &
-!$omp private(f,cosphi,tanphi,phi,juy,jux,theta,thetap,thetam, &
-!$omp thdn,thup,vx,dvdx,kup,kdn,ivrp,ivr,ivrm, &
-!$omp uy,dudy,dthetadp,jj,j,khi,klo, &
-!$omp klvrp,klvrm,k,klpt,kl,kch,i,dt1,dt2) &
-!$omp num_threads(omp_get_max_threads())
-!$omp do
- DO jy=jybords(1),jybords(2)
- phi = (ylat0 + jy * dy) * pi / 180.
- f = SIN(phi)
- cosphi = COS(phi)
- tanphi = f/cosphi
- f = 0.00014585 * f
-
- juy=jumpy_arr(jy)
-
- DO ix=0,nxmin1
- ! Define absolute gap length
- jux=jumpx_arr(ix)
-
- ! Loop over the vertical
- !***********************
- DO kl=1,nuvz
- theta=th_pre(ix,jy,kl)
- klvrp=kl+1
- klvrm=kl-1
- klpt=kl
- ! If top or bottom level, dthetadp is evaluated between the current
- ! level and the level inside, otherwise between level+1 and level-1
-
- IF (klvrp > nuvz) klvrp=nuvz
- IF (klvrm < 1) klvrm=1
- thetap=th_pre(ix,jy,klvrp)
- thetam=th_pre(ix,jy,klvrm)
+ integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
+ integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
+ integer :: nlck
+ real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
+ real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
+ real :: pvavr,ppml(nuvzmax)
+ real :: thup,thdn
+ real,parameter :: eps=1.e-5, p0=101325
+ 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)
- dthetadp=(thetap-thetam)/(ppmstore(ix,jy,klvrp)-ppmstore(ix,jy,klvrm))
-
- ! Compute vertical position at pot. temperature surface on subgrid
- ! and the wind at that position
- !*****************************************************************
- ! a) in x direction
- xloop: DO i=1,2
- if (i==1) ivr=ixvm_arr(ix)
- if (i==2) ivr=ixvp_arr(ix)
+ ! Set number of levels to check for adjacent theta
+ nlck=nuvz/3
+ !
+ ! Loop over entire grid
+ !**********************
+ do jy=0,nymin1
+ if (sglobal.and.jy.eq.0) goto 10
+ if (nglobal.and.jy.eq.nymin1) goto 10
+ phi = (ylat0 + jy * dy) * pi / 180.
+ f = 0.00014585 * sin(phi)
+ tanphi = tan(phi)
+ cosphi = cos(phi)
+ ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
+ jyvp=jy+1
+ jyvm=jy-1
+ if (jy.eq.0) jyvm=0
+ if (jy.eq.nymin1) jyvp=nymin1
+ ! Define absolute gap length
+ jumpy=2
+ if (jy.eq.0.or.jy.eq.nymin1) jumpy=1
+ if (sglobal.and.jy.eq.1) then
+ jyvm=1
+ jumpy=1
+ end if
+ if (nglobal.and.jy.eq.ny-2) then
+ jyvp=ny-2
+ jumpy=1
+ end if
+ juy=jumpy
+ !
+ do ix=0,nxmin1
+ ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
+ ixvp=ix+1
+ ixvm=ix-1
+ jumpx=2
+ if (xglobal) then
+ ivrp=ixvp
+ ivrm=ixvm
+ if (ixvm.lt.0) ivrm=ixvm+nxmin1
+ if (ixvp.ge.nx) ivrp=ixvp-nx+1
+ else
+ if (ix.eq.0) ixvm=0
+ if (ix.eq.nxmin1) ixvp=nxmin1
+ ivrp=ixvp
+ ivrm=ixvm
+ ! Define absolute gap length
+ if (ix.eq.0.or.ix.eq.nxmin1) jumpx=1
+ end if
+ jux=jumpx
+ ! Precalculate pressure values for efficiency
+ do kl=1,nuvz
+ ppml(kl)=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
+ end do
+ !
+ ! Loop over the vertical
+ !***********************
- ! Search adjacent levels for current theta value
- ! Spiral out from current level for efficiency
- kup=klpt
- kdn=klpt-1
- kch=0
+ do kl=1,nuvz
+ ppmk=akz(kl)+bkz(kl)*ps(ix,jy,1,n)
+ theta=tth(ix,jy,kl,n)*(100000./ppmk)**kappa
+ klvrp=kl+1
+ klvrm=kl-1
+ klpt=kl
+ ! If top or bottom level, dthetadp is evaluated between the current
+ ! level and the level inside, otherwise between level+1 and level-1
+ !
+ if (klvrp.gt.nuvz) klvrp=nuvz
+ if (klvrm.lt.1) klvrm=1
+ ppmk=akz(klvrp)+bkz(klvrp)*ps(ix,jy,1,n)
+ thetap=tth(ix,jy,klvrp,n)*(100000./ppmk)**kappa
+ ppmk=akz(klvrm)+bkz(klvrm)*ps(ix,jy,1,n)
+ thetam=tth(ix,jy,klvrm,n)*(100000./ppmk)**kappa
+ dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
- do while (kch=1) then
- kch=kch+1
- if (th_pre(ivr,jy,k)<=th_pre(ivr,jy,k+1)) then
- khi=k+1
- klo=k
- else
- khi=k
- klo=k+1
- end if
+ if (((thdn.ge.theta).and.(thup.le.theta)).or. &
+ ((thdn.le.theta).and.(thup.ge.theta))) then
+ dt1=abs(theta-thdn)
+ dt2=abs(theta-thup)
+ dt=dt1+dt2
+ if (dt.lt.eps) then ! Avoid division by zero error
+ dt1=0.5 ! G.W., 10.4.1996
+ dt2=0.5
+ dt=1.0
+ endif
+ vx(ii)=(vvh(ivr,jy,k)*dt2+vvh(ivr,jy,k+1)*dt1)/dt
+ goto 20
+ endif
+41 continue
+ ! Downward branch
+ kdn=kdn-1
+ if (kdn.lt.1) goto 40
+ kch=kch+1
+ k=kdn
+ ppmk=akz(k)+bkz(k)*ps(ivr,jy,1,n)
+ thdn=tth(ivr,jy,k,n)*(100000./ppmk)**kappa
+ ppmk=akz(k+1)+bkz(k+1)*ps(ivr,jy,1,n)
+ thup=tth(ivr,jy,k+1,n)*(100000./ppmk)**kappa
- dt1=th_pre(ivr,jy,khi)-theta
- dt2=theta-th_pre(ivr,jy,klo)
- if (dt1>=0 .and. dt2>=0) then
- if (dt1+dt2 0) THEN
- dvdx=(vx(2)-vx(1))/real(jux)/(dx*pi/180.)
- ELSE
- ivrp=ixvp_arr(ix)
- ivrm=ixvm_arr(ix)
- if (ix == 0) then
- if (xglobal) then
- ivrm=nxmin1-1
- end if
- else if (ix == nxmin1) then
- if (xglobal) then
- ivrp=1
- end if
- end if
- dvdx=vvh(ivrp,jy,kl)-vvh(ivrm,jy,kl)
- dvdx=dvdx/real(jumpx_arr(ix))/(dx*pi/180.)
- ! Only happens if no equivalent theta value
- ! can be found on either side, hence must use values
- ! from either side, same pressure level.
- END IF
-
- ! b) in y direction
- yloop: DO jj=1,2
- if (jj==1) j=jyvm_arr(jy)
- if (jj==2) j=jyvp_arr(jy)
- ! Search adjacent levels for current theta value
- ! Spiral out from current level for efficiency
- kup=klpt
- kdn=klpt-1
- kch=0
- DO WHILE (kch=1) then
- kch=kch+1
+ jj=0
+ do j=jyvm,jyvp,jumpy
+ jj=jj+1
+ ! Search adjacent levels for current theta value
+ ! Spiral out from current level for efficiency
+ kup=klpt-1
+ kdn=klpt
+ kch=0
+70 continue
+ ! Upward branch
+ kup=kup+1
+ if (kch.ge.nlck) goto 51 ! No more levels to check,
+ ! ! and no values found
+ if (kup.ge.nuvz) goto 71
+ kch=kch+1
+ k=kup
+ ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
+ thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
+ ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
+ thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
+ if (((thdn.ge.theta).and.(thup.le.theta)).or. &
+ ((thdn.le.theta).and.(thup.ge.theta))) then
+ dt1=abs(theta-thdn)
+ dt2=abs(theta-thup)
+ dt=dt1+dt2
+ if (dt.lt.eps) then ! Avoid division by zero error
+ dt1=0.5 ! G.W., 10.4.1996
+ dt2=0.5
+ dt=1.0
+ endif
+ uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
+ goto 50
+ endif
+71 continue
+ ! Downward branch
+ kdn=kdn-1
+ if (kdn.lt.1) goto 70
+ kch=kch+1
+ k=kdn
+ ppmk=akz(k)+bkz(k)*ps(ix,j,1,n)
+ thdn=tth(ix,j,k,n)*(100000./ppmk)**kappa
+ ppmk=akz(k+1)+bkz(k+1)*ps(ix,j,1,n)
+ thup=tth(ix,j,k+1,n)*(100000./ppmk)**kappa
+ if (((thdn.ge.theta).and.(thup.le.theta)).or. &
+ ((thdn.le.theta).and.(thup.ge.theta))) then
+ dt1=abs(theta-thdn)
+ dt2=abs(theta-thup)
+ dt=dt1+dt2
+ if (dt.lt.eps) then ! Avoid division by zero error
+ dt1=0.5 ! G.W., 10.4.1996
+ dt2=0.5
+ dt=1.0
+ endif
+ uy(jj)=(uuh(ix,j,k)*dt2+uuh(ix,j,k+1)*dt1)/dt
+ goto 50
+ endif
+ goto 70
+ ! This section used when no values were found
+51 continue
+ ! Must use uu at current level and lat. juy becomes smaller by 1
+ uy(jj)=uuh(ix,jy,kl)
+ juy=juy-1
+ ! Otherwise OK
+50 continue
+ end do
+ if (juy.gt.0) then
+ dudy=(uy(2)-uy(1))/real(juy)/(dy*pi/180.)
+ else
+ dudy=uuh(ix,jyvp,kl)-uuh(ix,jyvm,kl)
+ dudy=dudy/real(jumpy)/(dy*pi/180.)
+ end if
+ !
+ pvh(ix,jy,kl)=dthetadp*(f+(dvdx/cosphi-dudy &
+ +uuh(ix,jy,kl)*tanphi)/r_earth)*(-1.e6)*9.81
- if (th_pre(ix,j,k)<=th_pre(ix,j,k+1)) then
- khi=k+1
- klo=k
- else
- khi=k
- klo=k+1
- end if
- dt1=th_pre(ix,j,khi)-theta
- dt2=theta-th_pre(ix,j,klo)
- if (dt1>=0 .and. dt2>=0) then
- if (dt1+dt2 0) THEN
- dudy=(uy(2)-uy(1))/real(juy)/(dy*pi/180.)
- ELSE
- dudy=uuh(ix,jyvp_arr(jy),kl)-uuh(ix,jyvm_arr(jy),kl)
- dudy=dudy/real(jumpy_arr(jy))/(dy*pi/180.)
- END IF
-
- pvh(ix,jy,kl)=dthetadp*(f+(dvdx/cosphi-dudy+uuh(ix,jy,kl)*tanphi)/r_earth)*(-1.e6)*9.81
-
- !"RESET"
- jux=jumpx_arr(ix)
- juy=jumpy_arr(jy)
- END DO
- END DO
- END DO
-!$omp end do
-!$omp end parallel
-
- ! Fill in missing PV values on poles, if present
- ! Use mean PV of surrounding latitude ring
-
- if (sglobal) then
- do kl=1,nuvz
- pvh(0:nxmin1,0,kl)=sum(pvh(0:nxmin1,1,kl))/real(nxmin1+1)
- end do
- end if
- if (nglobal) then
- do kl=1,nuvz
- pvh(0:nxmin1,nymin1,kl)=sum(pvh(0:nxmin1,nymin1-1,kl))/real(nxmin1+1)
- end do
- end if
- write(*,*) 'leave calcpv'
-
- end subroutine
-
+end subroutine calcpv
diff --git a/flexpart_code/cmapf_mod.f90 b/flexpart_code/cmapf_mod.f90
index 4a914a0..576b291 100644
--- a/flexpart_code/cmapf_mod.f90
+++ b/flexpart_code/cmapf_mod.f90
@@ -26,12 +26,10 @@
module cmapf_mod
use par_mod, only: dp
- !this module is only a procedure module -> needs no "save" attribute
implicit none
private
- !routine names public declaration
public :: cc2gll, cll2xy, cgszll, cxy2ll, stlmbr, stcm2p
real,parameter :: rearth=6371.2, almst1=.9999999
@@ -49,7 +47,8 @@ subroutine cc2gll (strcmp, xlat,xlong, ue,vn, ug,vg)
implicit none
real :: strcmp(9), xlat, xlong, ue, vn, ug, vg
- real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot,tempdb,atempdb
+
+ real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot
along = cspanf( xlong - strcmp(2), -180., 180.)
if (xlat.gt.89.985) then
@@ -61,28 +60,8 @@ subroutine cc2gll (strcmp, xlat,xlong, ue,vn, ug,vg)
else
rot = - strcmp(1) * along
endif
-
- tempdb=radpdg*rot
- atempdb=abs(tempdb)
-
- if (atempdb<=pi) then
- if (atempdb<0.05E0_dp) then
- slong=tempdb*tempdb
- clong=1.0E0_dp-slong/2.0E0_dp+slong*slong/24.0E0_dp
- slong=tempdb-slong*tempdb/6.0E0_dp+slong*slong*tempdb/120.0E0_dp
- else
- clong = cos(tempdb)
- if (tempdb<0.0E0_dp) then
- slong = -1.0E0_dp * sqrt( 1.0E0_dp - clong*clong )
- else
- slong = sqrt( 1.0E0_dp - clong*clong )
- end if
- end if
- else
- clong=cos(tempdb)
- slong=sin(tempdb)
- end if
-
+ slong = sin( radpdg * rot )
+ clong = cos( radpdg * rot )
xpolg = slong * strcmp(5) + clong * strcmp(6)
ypolg = clong * strcmp(5) - slong * strcmp(6)
ug = ypolg * ue + xpolg * vn
diff --git a/flexpart_code/cmapf_mod.mod b/flexpart_code/cmapf_mod.mod
deleted file mode 100644
index f56a2a3..0000000
Binary files a/flexpart_code/cmapf_mod.mod and /dev/null differ
diff --git a/flexpart_code/com_mod.f90 b/flexpart_code/com_mod.f90
index 4b71c21..e9b79fc 100644
--- a/flexpart_code/com_mod.f90
+++ b/flexpart_code/com_mod.f90
@@ -17,7 +17,6 @@ module com_mod
maxreceptor, maxpart, maxrand, nwzmax, nuvzmax
implicit none
- save
!****************************************************************
! Variables defining where FLEXPART input/output files are stored
@@ -305,10 +304,10 @@ module com_mod
! Fixed fields, unchangeable with time
!*************************************
- real, target :: oro(0:nxmax-1,0:nymax-1)
- real, target :: excessoro(0:nxmax-1,0:nymax-1)
- real, target :: lsm(0:nxmax-1,0:nymax-1)
- real, target :: xlanduse(0:nxmax-1,0:nymax-1,numclass)
+ real :: oro(0:nxmax-1,0:nymax-1)
+ real :: excessoro(0:nxmax-1,0:nymax-1)
+ real :: lsm(0:nxmax-1,0:nymax-1)
+ real :: xlanduse(0:nxmax-1,0:nymax-1,numclass)
! oro [m] orography of the ECMWF model
! excessoro excess orography mother domain
@@ -328,8 +327,8 @@ module com_mod
real :: pv(0:nxmax-1,0:nymax-1,nzmax,2)
real :: rho(0:nxmax-1,0:nymax-1,nzmax,2)
real :: drhodz(0:nxmax-1,0:nymax-1,nzmax,2)
- real, target :: tth(0:nxmax-1,0:nymax-1,nuvzmax,2)
- real, target :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
+ real :: tth(0:nxmax-1,0:nymax-1,nuvzmax,2)
+ real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
integer :: cloudsh(0:nxmax-1,0:nymax-1,2)
@@ -351,25 +350,25 @@ module com_mod
! 2d fields
!**********
- real, target :: ps(0:nxmax-1,0:nymax-1,1,2)
- real, target :: sd(0:nxmax-1,0:nymax-1,1,2)
- real, target :: msl(0:nxmax-1,0:nymax-1,1,2)
- real, target :: tcc(0:nxmax-1,0:nymax-1,1,2)
- real, target :: u10(0:nxmax-1,0:nymax-1,1,2)
- real, target :: v10(0:nxmax-1,0:nymax-1,1,2)
- real, target :: tt2(0:nxmax-1,0:nymax-1,1,2)
- real, target :: td2(0:nxmax-1,0:nymax-1,1,2)
- real, target :: lsprec(0:nxmax-1,0:nymax-1,1,2)
- real, target :: convprec(0:nxmax-1,0:nymax-1,1,2)
- real, target :: sshf(0:nxmax-1,0:nymax-1,1,2)
- real, target :: ssr(0:nxmax-1,0:nymax-1,1,2)
- real, target :: surfstr(0:nxmax-1,0:nymax-1,1,2)
- real, target :: ustar(0:nxmax-1,0:nymax-1,1,2)
- real, target :: wstar(0:nxmax-1,0:nymax-1,1,2)
- real, target :: hmix(0:nxmax-1,0:nymax-1,1,2)
- real, target :: tropopause(0:nxmax-1,0:nymax-1,1,2)
- real, target :: oli(0:nxmax-1,0:nymax-1,1,2)
- real, target :: diffk(0:nxmax-1,0:nymax-1,1,2)
+ real :: ps(0:nxmax-1,0:nymax-1,1,2)
+ real :: sd(0:nxmax-1,0:nymax-1,1,2)
+ real :: msl(0:nxmax-1,0:nymax-1,1,2)
+ real :: tcc(0:nxmax-1,0:nymax-1,1,2)
+ real :: u10(0:nxmax-1,0:nymax-1,1,2)
+ real :: v10(0:nxmax-1,0:nymax-1,1,2)
+ real :: tt2(0:nxmax-1,0:nymax-1,1,2)
+ real :: td2(0:nxmax-1,0:nymax-1,1,2)
+ real :: lsprec(0:nxmax-1,0:nymax-1,1,2)
+ real :: convprec(0:nxmax-1,0:nymax-1,1,2)
+ real :: sshf(0:nxmax-1,0:nymax-1,1,2)
+ real :: ssr(0:nxmax-1,0:nymax-1,1,2)
+ real :: surfstr(0:nxmax-1,0:nymax-1,1,2)
+ real :: ustar(0:nxmax-1,0:nymax-1,1,2)
+ real :: wstar(0:nxmax-1,0:nymax-1,1,2)
+ real :: hmix(0:nxmax-1,0:nymax-1,1,2)
+ real :: tropopause(0:nxmax-1,0:nymax-1,1,2)
+ real :: oli(0:nxmax-1,0:nymax-1,1,2)
+ real :: diffk(0:nxmax-1,0:nymax-1,1,2)
! ps surface pressure
! sd snow depth
@@ -671,6 +670,5 @@ module com_mod
! rannumb field of normally distributed random numbers
- logical :: verttransform_inited = .false.
end module com_mod
diff --git a/flexpart_code/com_mod.mod b/flexpart_code/com_mod.mod
deleted file mode 100644
index 880a4fc..0000000
Binary files a/flexpart_code/com_mod.mod and /dev/null differ
diff --git a/flexpart_code/conv_mod.f90 b/flexpart_code/conv_mod.f90
index 16417ec..7b1bbc8 100644
--- a/flexpart_code/conv_mod.f90
+++ b/flexpart_code/conv_mod.f90
@@ -13,7 +13,6 @@ module conv_mod
use par_mod, only: nconvlevmax, na, nxmax, nymax, nxmaxn, nymaxn, maxnests
implicit none
- save
!integer,parameter :: nconvlevmax = nuvzmax-1, &
! na = nconvlevmax+1
diff --git a/flexpart_code/conv_mod.mod b/flexpart_code/conv_mod.mod
deleted file mode 100644
index 4f65e9e..0000000
Binary files a/flexpart_code/conv_mod.mod and /dev/null differ
diff --git a/flexpart_code/ew.f90 b/flexpart_code/ew.f90
index 6c61ca0..0d44e82 100644
--- a/flexpart_code/ew.f90
+++ b/flexpart_code/ew.f90
@@ -1,159 +1,47 @@
+!**********************************************************************
+! 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 . *
+!**********************************************************************
- function ew(t_in)
- implicit none
- real :: ew
- real, intent(in) :: t_in
- interface
- elemental function ew_nocheck(t_in)
- implicit none
- real, intent(in) :: t_in
- real :: ew_nocheck,w1,w2
- end function
- end interface
- if(t_in <= 0.) then
- write(*,*) 'TEMP: ',t_in
- stop 'SORRY: T NOT IN [K]'
- end if
- ew=ew_nocheck(t_in)
- end function
+real function ew(x)
- elemental function ew_nocheck(t_in)
- implicit none
- real, intent(in) :: t_in
- real :: ew_nocheck,w1,w2
- interface
- elemental function ew_gg(x)
- implicit none
- real :: ew_gg
- real, intent(in) :: x
- end function
- elemental function ew_wobus_refined(t_in_wrf) result(ewrf)
- implicit none
- real :: ewrf
- real, intent(in) :: t_in_wrf
- end function
- end interface
+ !****************************************************************
+ !SAETTIGUNGSDAMPFDRUCK UEBER WASSER IN PA. X IN KELVIN.
+ !NACH DER GOFF-GRATCH-FORMEL.
+ !****************************************************************
- if (abs(t_in-300.0)<100.0) then
- ew_nocheck=ew_wobus_refined(t_in)
- else if (t_in>399.9) then
- w1=min(t_in,410.0)
- w2=(w1-400.0)
- w1=(410.0-w1)
- ew_nocheck=(w1*ew_wobus_refined(t_in)+w2*ew_gg(t_in))/(410.0-400.0)
- else if (t_in<200.1) then
- w1=max(t_in,194.0)
- w2=(w1-194.0)
- w1=(200.0-w1)
- ew_nocheck=(w2*ew_wobus_refined(t_in)+w1*ew_gg(t_in))/(200.0-194.0)
- end if
- end function
+ implicit none
- elemental FUNCTION ew_wobus_refined(t_in_wrf) result(ewrf)
- implicit none
- real :: ewrf
- real, intent(in) :: t_in_wrf
- real(8) :: ref_pol,td
- real(8), parameter :: acoeff_ewrf(0:10)=[ &
- -6958.391726914337D0, 245.45012981061453D0, &
- -3.8799554209110156D0, 3.619958715118083D-2, &
- -2.207540481565966D-4, 9.194341350199019D-7, &
- -2.6487536137880456D-9, 5.211800006142822D-12, &
- -6.703450445729384D-15, 5.089509368882543D-18, &
- -1.7321846403884217D-21 ]
- interface
- elemental function ew_wobus(t_in) result(esw)
- real, intent(in) :: t_in
- real :: esw
- end function ew_wobus
- end interface
- td=real(t_in_wrf,8)
-
- ref_pol = acoeff_ewrf(0)+ &
- td*(acoeff_ewrf(1)+ &
- td*(acoeff_ewrf(2)+ &
- td*(acoeff_ewrf(3)+ &
- td*(acoeff_ewrf(4)+ &
- td*(acoeff_ewrf(5)+ &
- td*(acoeff_ewrf(6)+ &
- td*(acoeff_ewrf(7)+ &
- td*(acoeff_ewrf(8)+ &
- td*(acoeff_ewrf(9)+ &
- td*acoeff_ewrf(10) )))))))))
-
- ewrf=ew_wobus(t_in_wrf)*real(ref_pol)
- end function
+ real :: x, y, a, c, d
- elemental function ew_wobus(t_in)
- ! this function returns the saturation vapor pressure esw (millibars)
- ! over liquid water given the temperature t (kelvin). the polynomial
- ! approximation below is due to herman wobus, a mathematician who
- ! worked at the navy weather research facility, norfolk, virginia,
- ! but who is now retired. the coefficients of the polynomial were
- ! chosen to fit the values in table 94 on pp. 351-353 of the smith-
- ! sonian meteorological tables by roland list (6th edition). the
- ! approximation is valid for -50 < t < 100c.
- !
- ! baker, schlatter 17-may-1982 original version.
- !
- ! es0 = saturation vapor ressure over liquid water at 0c
- implicit none
- real, intent(in) :: t_in
- real :: ew_wobus
- double precision :: pol,t
- double precision :: es0
- es0=6.1078D0
- t=dble(t_in-273.15)
- pol = 0.99999683D0 + t*(-0.90826951D-02 + &
- t*(0.78736169D-04 + t*(-0.61117958D-06 + &
- t*(0.43884187D-08 + t*(-0.29883885D-10 + &
- t*(0.21874425D-12 + t*(-0.17892321D-14 + &
- t*(0.11112018D-16 + t*(-0.30994571D-19)))))))))
- ew_wobus = real(100.0D0*es0 / pol**8)
- end function ew_wobus
+ ew=0.
+ if(x.le.0.) stop 'sorry: t not in [k]'
+ y=373.16/x
+ a=-7.90298*(y-1.)
+ a=a+(5.02808*0.43429*alog(y))
+ c=(1.-(1./y))*11.344
+ c=-1.+(10.**c)
+ c=-1.3816*c/(10.**7)
+ d=(1.-y)*3.49149
+ d=-1.+(10.**d)
+ d=8.1328*d/(10.**3)
+ y=a+c+d
+ ew=101324.6*(10.**y) ! Saettigungsdampfdruck in Pa
- elemental function ew_gg(x)
- ! ****************************************************************
- ! SAETTIGUNGSDAMPFDRUCK UEBER WASSER IN PA. X IN KELVIN.
- ! NACH DER GOFF-GRATCH-FORMEL.
- ! ****************************************************************
- implicit none
- real :: ew_gg
- real, intent(in) :: x
- real :: y,a,c,d
- y = 373.16/x
- a = -7.90298 * (y-1.)
- a = a + (5.02808*0.43429*LOG(y))
- c = -1. + 10.0**( 11.344 - ( (11.344/373.16) * x ) )
- c = -1.3816E-7 * c
- d = (1.-y)*3.49149
- d = 8.1328E-3 * ( -1. + (10.**d) )
- y = a + c + d
- ew_gg=1013.246*(10.**y)*100. ! Saettigungsdampfdruck in Pa
- end function
-
- elemental FUNCTION ew_fitted(t_in_ewf)
- ! a 14th order polynomial fit to GG -- performes worse than Wobus however
- implicit none
- integer :: i
- real, intent(in) :: t_in_ewf
- real :: ew_fitted
- real(8) :: t_in_ewf_dble,ref_pol
- real(8), parameter :: a_ewfit(0:14) = [ &
- -1.2241382741591574D6, 6.2760007286859395D4, &
- -1.4467525836456739D3, 1.9735651467825693D1, &
- -0.17595965373830060D0, 0.10629184255214715D-2, &
- -4.3022048695559316D-6, 1.0568061161848794D-8, &
- -7.9252050616823303D-12, -4.5604622615056078D-14,&
- 1.9821751947758235D-16, -4.0185846042019225D-19,&
- 4.7379290439536481D-22, -3.1275571418968759D-25,&
- 8.9966308005149801D-29 ]
- t_in_ewf_dble=real(t_in_ewf,8)
- ref_pol=0.0D0
- do i=14,1,-1
- ref_pol = t_in_ewf_dble*(a_ewfit(i)+ref_pol)
- end do
- ref_pol=ref_pol+a_ewfit(0)
- ew_fitted=real(ref_pol)
- end function
-
+end function ew
diff --git a/flexpart_code/flux_mod.f90 b/flexpart_code/flux_mod.f90
index c36fee5..220dd4f 100644
--- a/flexpart_code/flux_mod.f90
+++ b/flexpart_code/flux_mod.f90
@@ -26,7 +26,6 @@ module flux_mod
! areaeast,areanorth [m2] side areas of each grid cell
implicit none
- save
real,allocatable, dimension (:,:,:,:,:,:,:) :: flux
diff --git a/flexpart_code/flux_mod.mod b/flexpart_code/flux_mod.mod
deleted file mode 100644
index d522127..0000000
Binary files a/flexpart_code/flux_mod.mod and /dev/null differ
diff --git a/flexpart_code/getfields.f90 b/flexpart_code/getfields.f90
index bee7ef1..fad4e96 100644
--- a/flexpart_code/getfields.f90
+++ b/flexpart_code/getfields.f90
@@ -71,10 +71,10 @@ subroutine getfields(itime,nstop)
integer :: indj,itime,nstop,memaux
- real, target :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
- real, target :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
- real, target :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
- real, target :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
+ 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 :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
diff --git a/flexpart_code/gridcheck_gfs.f90 b/flexpart_code/gridcheck_gfs.f90
index 41b417a..f1fce80 100644
--- a/flexpart_code/gridcheck_gfs.f90
+++ b/flexpart_code/gridcheck_gfs.f90
@@ -88,12 +88,12 @@ subroutine gridcheck
integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
real :: sizesouth,sizenorth,xauxa,pint
real :: akm_usort(nwzmax)
- real,parameter :: eps=spacing(2.0_4*360.0_4)
+ real,parameter :: eps=0.0001
! NCEP GFS
real :: pres(nwzmax), help
- integer :: i180
+ integer :: i179,i180,i181
! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
@@ -104,16 +104,14 @@ subroutine gridcheck
!HSO grib api error messages
character(len=24) :: gribErrorMsg = 'Error reading grib file'
character(len=20) :: gribFunction = 'gridcheckwind_gfs'
-
- write(*,*) 'enter gridcheck...'
-
+ !
if (numbnests.ge.1) then
write(*,*) ' ###########################################'
write(*,*) ' FLEXPART ERROR SUBROUTINE GRIDCHECK:'
- write(*,*) ' NO NESTED WINDFIELDS ALLOWED FOR GFS! '
+ write(*,*) ' NO NESTED WINDFIELDAS ALLOWED FOR GFS! '
write(*,*) ' ###########################################'
stop
- end if
+ endif
iumax=0
iwmax=0
@@ -122,255 +120,253 @@ subroutine gridcheck
ifn=1
else
ifn=numbwf
- end if
-
- ! OPENING OF DATA FILE (GRIB CODE) ... loop (acc to user choice)
- ! if file not present (e.g. CD ROM)
- iret=GRIB_SUCCESS+1
- do while (iret/=GRIB_SUCCESS)
- call grib_open_file(ifile,path(3)(1:length(3))//trim(wfname(ifn)),'r',iret)
- if (iret.ne.GRIB_SUCCESS) then ! ERROR DETECTED
- write(*,*)
- write(*,*) ' ################################################# '
- write(*,*) ' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
- write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
- write(*,*) ' ################################################# '
- write(*,*)
- write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND !!!'
- write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE... !!!'
- write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
- write(*,'(a)') '!!! THE "X" KEY... !!!'
- write(*,*)
- read(*,'(a)') opt
- if(opt.eq.'X') stop
- end if
- end do
-
+ endif
+ !
+ ! OPENING OF DATA FILE (GRIB CODE)
+ !
+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
!turn on support for multi fields messages
call grib_multi_support_on
ifield=0
- do
- ifield=ifield+1
-
- ! GET NEXT FIELDS
- call grib_new_from_file(ifile,igrib,iret)
- if (iret.eq.GRIB_END_OF_FILE ) then
- exit ! EOF DETECTED
- elseif (iret.ne.GRIB_SUCCESS) then ! ERROR DETECTED
- write(*,*)
- write(*,*) ' ################################################# '
- write(*,*) ' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
- write(*,*) ' ERROR READING NEXT FIELD FROM FILE '//wfname(ifn)
- write(*,*) ' ################################################# '
- write(*,*)
- write(*,'(a)') '!!! HANDLING OF SPLIT/PARTIAL GFS !!!'
- write(*,'(a)') '!!! FILES IS REMOVED FROM THIS VERSION !!!'
- write(*,'(a)') '!!! FLEXPART TERMINATING... !!!'
- write(*,*)
- stop
- end if
-
- !first see if we read GRIB1 or GRIB2
- call grib_get_int(igrib,'editionNumber',gribVer,iret)
+10 ifield=ifield+1
+ !
+ ! GET NEXT FIELDS
+ !
+ call grib_new_from_file(ifile,igrib,iret)
+ if (iret.eq.GRIB_END_OF_FILE ) then
+ goto 30 ! EOF DETECTED
+ elseif (iret.ne.GRIB_SUCCESS) then
+ 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
+
+ !read the grib1 identifiers
+ call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'level',isec1(8),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+
+ !get the size and data of the values array
+ call grib_get_real4_array(igrib,'values',zsec4,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+
+ else ! 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,'scaledValueOfFirstFixedSurface', &
+ valSurf,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+
+ !convert to grib1 identifiers
+ isec1(6)=-1
+ isec1(7)=-1
+ isec1(8)=-1
+ if ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
+ isec1(6)=33 ! indicatorOfParameter
+ isec1(7)=100 ! indicatorOfTypeOfLevel
+ isec1(8)=valSurf/100 ! level, convert to hPa
+ elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
+ isec1(6)=7 ! indicatorOfParameter
+ isec1(7)=1 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
+ .and.(discipl.eq.2)) then ! LSM
+ isec1(6)=81 ! indicatorOfParameter
+ isec1(7)=1 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ endif
+
+ if (isec1(6).ne.-1) then
+ ! get the size and data of the values array
+ call grib_get_real4_array(igrib,'values',zsec4,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
+ endif
+
+ endif ! gribVer
+
+ if(ifield.eq.1) then
+
+ !get the required fields from section 2
+ !store compatible to gribex input
+ call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
+ isec2(2),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
+ isec2(3),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
+ xaux1in,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
+ xaux2in,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
+ yaux1in,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
+ yaux2in,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+
+ xaux1=xaux1in
+ xaux2=xaux2in
+ yaux1=yaux1in
+ yaux2=yaux2in
+
+ nxfield=isec2(2)
+ ny=isec2(3)
+ if((abs(xaux1).lt.eps).and.(xaux2.ge.359)) then ! NCEP DATA FROM 0 TO
+ xaux1=-179.0 ! 359 DEG EAST ->
+ xaux2=-179.0+360.-360./real(nxfield) ! TRANSFORMED TO -179
+ endif ! TO 180 DEG EAST
+ if (xaux1.gt.180) xaux1=xaux1-360.0
+ if (xaux2.gt.180) xaux2=xaux2-360.0
+ if (xaux1.lt.-180) xaux1=xaux1+360.0
+ if (xaux2.lt.-180) xaux2=xaux2+360.0
+ if (xaux2.lt.xaux1) xaux2=xaux2+360.
+ xlon0=xaux1
+ ylat0=yaux1
+ dx=(xaux2-xaux1)/real(nxfield-1)
+ dy=(yaux2-yaux1)/real(ny-1)
+ dxconst=180./(dx*r_earth*pi)
+ dyconst=180./(dy*r_earth*pi)
+ !HSO end edits
+
+
+ ! 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
+ nx=nxfield+1 ! field is cyclic
+ xglobal=.true.
+ if (abs(nxshift).ge.nx) &
+ stop 'nxshift in file par_mod is too large'
+ xlon0=xlon0+real(nxshift)*dx
+ else
+ nx=nxfield
+ xglobal=.false.
+ if (nxshift.ne.0) &
+ stop 'nxshift (par_mod) must be zero for non-global domain'
+ endif
+ nxmin1=nx-1
+ nymin1=ny-1
+ if (xlon0.gt.180.) xlon0=xlon0-360.
+ 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
+ sizesouth=6.*(switchsouth+90.)/dy
+ call stlmbr(southpolemap,-90.,0.)
+ call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
+ sizesouth,switchsouth,180.)
+ switchsouthg=(switchsouth-ylat0)/dy
+ else
+ sglobal=.false.
+ switchsouthg=999999.
+ endif
+ 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
+ sizenorth=6.*(90.-switchnorth)/dy
+ call stlmbr(northpolemap,90.,0.)
+ call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
+ sizenorth,switchnorth,180.)
+ switchnorthg=(switchnorth-ylat0)/dy
+ else
+ nglobal=.false.
+ switchnorthg=999999.
+ endif
+ endif ! ifield.eq.1
+
+ if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative'
+ if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
+
+ ! NCEP ISOBARIC LEVELS
+ !*********************
+
+ if((isec1(6).eq.33).and.(isec1(7).eq.100)) then ! check for U wind
+ iumax=iumax+1
+ pres(iumax)=real(isec1(8))*100.0
+ endif
+
- if (gribVer.eq.1) then ! GRIB Edition 1
-
- !read the grib1 identifiers
- call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igrib,'level',isec1(8),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
-
- !get the size and data of the values array
- if ( (isec1(6)==33 .and. isec1(7)==100) .or. ((isec1(6)==81 .or. isec1(6)==7) .and. isec1(7)==1) ) then
- call grib_get_real4_array(igrib,'values',zsec4,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- end if
- else ! 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,'scaledValueOfFirstFixedSurface',valSurf,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
-
- !convert to grib1 identifiers
- isec1(6:8)=-1
- if ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
- isec1(6)=33 ! indicatorOfParameter
- isec1(7)=100 ! indicatorOfTypeOfLevel
- isec1(8)=valSurf/100 ! level, convert to hPa
- else if ( (parCat.eq.3 .and. parNum.eq.5 .and. typSurf.eq.1) .or. &
- (parCat.eq.0 .and. parNum.eq.0 .and. typSurf.eq.1 .and. discipl.eq.2) ) then
- isec1(6)=81 ! indicatorOfParameter
- isec1(7)=1 ! indicatorOfTypeOfLevel
- isec1(8)=0
- if (parCat.eq.3 .and. parNum.eq.5 .and. typSurf.eq.1) isec1(6)=7
- end if
-
- if (isec1(6)/=-1) then ! LSM
- !get the size and data of the values array
- call grib_get_real4_array(igrib,'values',zsec4,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- end if
-
- end if ! gribVer
-
- if(ifield.eq.1) then
-
- !get the required fields from section 2
- !store compatible to gribex input
- call grib_get_int(igrib,'numberOfPointsAlongAParallel',isec2(2),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igrib,'numberOfPointsAlongAMeridian',isec2(3),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees',xaux1in,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees',xaux2in,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees',yaux1in,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees',yaux2in,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
-
- xaux1=xaux1in
- xaux2=xaux2in
- yaux1=yaux1in
- yaux2=yaux2in
-
- nxfield=isec2(2)
- ny=isec2(3)
- if((abs(xaux1).lt.eps).and.(xaux2.ge.359.)) then ! NCEP DATA FROM 0 TO
- xaux1=-180.0 ! 359 DEG EAST ->
- xaux2=-180.0+360.-360./real(nxfield) ! TRANSFORMED TO -180
- end if ! TO 179 DEG EAST
- if (xaux1.gt.180) xaux1=xaux1-360.0
- if (xaux2.gt.180) xaux2=xaux2-360.0
- if (xaux1.lt.-180) xaux1=xaux1+360.0
- if (xaux2.lt.-180) xaux2=xaux2+360.0
- if (xaux2.lt.xaux1) xaux2=xaux2+360.
- xlon0=xaux1
- ylat0=yaux1
- dx=(xaux2-xaux1)/real(nxfield-1)
- dy=(yaux2-yaux1)/real(ny-1)
- dxconst=180./(dx*r_earth*pi)
- dyconst=180./(dy*r_earth*pi)
- !HSO end edits
- i180=nint(180./dx) ! 0.5 deg data
-
-
- ! 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
- nx=nxfield+1 ! field is cyclic
- xglobal=.true.
- if (abs(nxshift).ge.nx) stop 'nxshift in file par_mod is too large'
- xlon0=xlon0+real(nxshift)*dx
- else
- nx=nxfield
- xglobal=.false.
- if (nxshift.ne.0) stop 'nxshift (par_mod) must be zero for non-global domain'
- end if
- nxmin1=nx-1
- nymin1=ny-1
- if (xlon0.gt.180.) xlon0=xlon0-360.
- 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
- sizesouth=6.*(switchsouth+90.)/dy
- call stlmbr(southpolemap,-90.,0.)
- call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth,sizesouth,switchsouth,180.)
- switchsouthg=(switchsouth-ylat0)/dy
- else
- sglobal=.false.
- switchsouthg=999999.
- end if
- 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
- sizenorth=6.*(90.-switchnorth)/dy
- call stlmbr(northpolemap,90.,0.)
- call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth,sizenorth,switchnorth,180.)
- switchnorthg=(switchnorth-ylat0)/dy
- else
- nglobal=.false.
- switchnorthg=999999.
- end if
- end if ! ifield.eq.1
-
- if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative'
- if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
-
- ! NCEP ISOBARIC LEVELS
- !*********************
-
- if((isec1(6).eq.33).and.(isec1(7).eq.100)) then ! check for U wind
- iumax=iumax+1
- pres(iumax)=real(isec1(8))*100.0
- end if
-
-
-
-
- ! NCEP TERRAIN
- !*************
-
- if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
- do jy=0,ny-1
- do ix=0,nxfield-1
- help=zsec4(nxfield*(ny-jy-1)+ix+1)
- if(ix.lt.i180) then
- oro(i180+ix,jy)=help
- excessoro(i180+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
- else
- oro(ix-i180,jy)=help
- excessoro(ix-i180,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
- end if
- end do
+ i179=nint(179./dx)
+ 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
+ !*************
+
+ if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
+ do jy=0,ny-1
+ do ix=0,nxfield-1
+ help=zsec4(nxfield*(ny-jy-1)+ix+1)
+ if(ix.le.i180) then
+ oro(i179+ix,jy)=help
+ excessoro(i179+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
+ else
+ oro(ix-i181,jy)=help
+ excessoro(ix-i181,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
+ endif
end do
- end if
-
- ! NCEP LAND SEA MASK
- !*******************
-
- if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
- do jy=0,ny-1
- do ix=0,nxfield-1
- help=zsec4(nxfield*(ny-jy-1)+ix+1)
- if(ix.lt.i180) then
- lsm(i180+ix,jy)=help
- else
- lsm(ix-i180,jy)=help
- end if
- end do
+ end do
+ endif
+
+ ! NCEP LAND SEA MASK
+ !*******************
+
+ if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
+ do jy=0,ny-1
+ do ix=0,nxfield-1
+ help=zsec4(nxfield*(ny-jy-1)+ix+1)
+ if(ix.le.i180) then
+ lsm(i179+ix,jy)=help
+ else
+ lsm(ix-i181,jy)=help
+ endif
end do
- end if
+ end do
+ endif
- call grib_release(igrib)
+ call grib_release(igrib)
- end do !! READ NEXT LEVEL OR PARAMETER
+ goto 10 !! READ NEXT LEVEL OR PARAMETER
+ !
+ ! CLOSING OF INPUT DATA FILE
+ !
- ! CLOSING OF INPUT DATA FILE (HSO)
+ ! HSO
+30 continue
call grib_close_file(ifile)
! HSO end edits
@@ -378,19 +374,39 @@ subroutine gridcheck
nwz =iumax
nlev_ec=iumax
- if (nx>nxmax .or. ny>nymax) write(*,*) 'FLEXPART error: Too many grid points in x/y direction:'
- if (nuvz.gt.nuvzmax) write(*,*) 'FLEXPART error: Too many u,v grid points in z direction.'
- if (nwz.gt.nwzmax) write(*,*) 'FLEXPART error: Too many w grid points in z direction.'
- if (nx>nxmax) write(*,*) 'nx,nxmax:',nx,nxmax
- if (ny>nymax) write(*,*) 'ny,nymax:',ny,nymax
- if (nuvz>nuvzmax) write(*,*) 'nuvz,nuvzmax:',nuvz,nuvzmax
- if (nwz>nwzmax) write(*,*) 'nwz,nwzmax:',nwz,nwzmax
- if (nx>nxmax .or. ny>nymax .or. nuvz>nuvzmax .or. nwz>nwzmax) then
- write(*,*) 'STOPPING!!! Reduce resolution of wind fields.'
+ if (nx.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 par_mod.'
+ write(*,*) nx,nxmax
+ stop
+ endif
+
+ if (ny.gt.nymax) then
+ 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
- end if
+ endif
+ if (nuvz.gt.nuvzmax) then
+ write(*,*) 'FLEXPART error: Too many u,v grid points in z '// &
+ 'direction.'
+ write(*,*) 'Reduce resolution of wind fields.'
+ write(*,*) 'Or change parameter settings in file par_mod.'
+ write(*,*) nuvz,nuvzmax
+ stop
+ endif
+
+ if (nwz.gt.nwzmax) then
+ write(*,*) 'FLEXPART error: Too many w grid points in z '// &
+ 'direction.'
+ write(*,*) 'Reduce resolution of wind fields.'
+ write(*,*) 'Or change parameter settings in file par_mod.'
+ write(*,*) nwz,nwzmax
+ stop
+ endif
! If desired, shift all grids by nxshift grid cells
!**************************************************
@@ -399,18 +415,21 @@ subroutine gridcheck
call shift_field_0(oro,nxfield,ny)
call shift_field_0(lsm,nxfield,ny)
call shift_field_0(excessoro,nxfield,ny)
- end if
+ endif
! Output of grid info
!********************
write(*,*)
write(*,*)
- write(*,'(a,2i7)') '# of vertical levels in NCEP data: ',nuvz,nwz
+ write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
+ nuvz,nwz
write(*,*)
write(*,'(a)') 'Mother domain:'
- write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ',xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
- write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range: ',ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
+ write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', &
+ xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
+ write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range: ', &
+ ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
write(*,*)
@@ -429,14 +448,15 @@ subroutine gridcheck
!******************************
! change Sabine Eckhardt: akm should always be in descending order ... readwind adapted!
!******************************
- do i=1,nwz
- if (akm_usort(1).gt.akm_usort(2)) then
- akm(i)=akm_usort(i)
- else
- akm(i)=akm_usort(nwz-i+1)
- end if
- end do
+ do i=1,nwz
+ if (akm_usort(1).gt.akm_usort(2)) then
+ akm(i)=akm_usort(i)
+ else
+ akm(i)=akm_usort(nwz-i+1)
+ endif
+ end do
+ !
! CALCULATION OF AKZ, BKZ
! AKZ,BKZ: model discretization parameters at the center of each model
! layer
@@ -465,19 +485,54 @@ subroutine gridcheck
bknew(i)=bkz(i)
end do
+ ! Switch on following lines to use doubled vertical resolution
+ !*************************************************************
+ !nz=nuvz+nwz-1
+ !if (nz.gt.nzmax) stop 'nzmax too small'
+ !do 100 i=1,nwz
+ ! aknew(2*(i-1)+1)=akm(i)
+ !00 bknew(2*(i-1)+1)=bkm(i)
+ !do 110 i=2,nuvz
+ ! aknew(2*(i-1))=akz(i)
+ !10 bknew(2*(i-1))=bkz(i)
+ ! End doubled vertical resolution
+
+
! Determine the uppermost level for which the convection scheme shall be applied
! by assuming that there is no convection above 50 hPa (for standard SLP)
!*****************************************************************************
do i=1,nuvz-2
pint=akz(i)+bkz(i)*101325.
- if (pint.lt.5000.) exit
+ if (pint.lt.5000.) goto 96
end do
- nconvlev=i
+96 nconvlev=i
if (nconvlev.gt.nconvlevmax-1) then
nconvlev=nconvlevmax-1
write(*,*) 'Attention, convection only calculated up to ', &
akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
- end if
- write(*,*) 'exit gridcheck...'
+ endif
+
+ return
+
+999 write(*,*)
+ write(*,*) ' ###########################################'// &
+ '###### '
+ write(*,*) ' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
+ write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
+ write(*,*) ' ###########################################'// &
+ '###### '
+ write(*,*)
+ write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND !!!'
+ write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE... !!!'
+ write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
+ write(*,'(a)') '!!! THE "X" KEY... !!!'
+ write(*,*)
+ read(*,'(a)') opt
+ if(opt.eq.'X') then
+ stop
+ else
+ goto 5
+ endif
+
end subroutine gridcheck
diff --git a/flexpart_code/hanna_mod.f90 b/flexpart_code/hanna_mod.f90
index c1f5b0e..91717ec 100644
--- a/flexpart_code/hanna_mod.f90
+++ b/flexpart_code/hanna_mod.f90
@@ -1,7 +1,6 @@
module hanna_mod
implicit none
- save
real :: ust,wst,ol,h,zeta,sigu,sigv,tlu,tlv,tlw
real :: sigw,dsigwdz,dsigw2dz
diff --git a/flexpart_code/hanna_mod.mod b/flexpart_code/hanna_mod.mod
deleted file mode 100644
index a26ffa3..0000000
Binary files a/flexpart_code/hanna_mod.mod and /dev/null differ
diff --git a/flexpart_code/interpol_mod.f90 b/flexpart_code/interpol_mod.f90
index f22df86..0a9c329 100644
--- a/flexpart_code/interpol_mod.f90
+++ b/flexpart_code/interpol_mod.f90
@@ -3,7 +3,6 @@ module interpol_mod
use par_mod, only: nzmax, maxspec
implicit none
- save
real :: uprof(nzmax),vprof(nzmax),wprof(nzmax)
real :: usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
diff --git a/flexpart_code/interpol_mod.mod b/flexpart_code/interpol_mod.mod
deleted file mode 100644
index 012737e..0000000
Binary files a/flexpart_code/interpol_mod.mod and /dev/null differ
diff --git a/flexpart_code/oh_mod.f90 b/flexpart_code/oh_mod.f90
index 2ba0e5d..0cb0b01 100644
--- a/flexpart_code/oh_mod.f90
+++ b/flexpart_code/oh_mod.f90
@@ -25,7 +25,6 @@ module oh_mod
!for this field
implicit none
- save
real,allocatable, dimension (:,:,:,:) :: OH_field
real,allocatable, dimension (:) :: OH_field_height
diff --git a/flexpart_code/oh_mod.mod b/flexpart_code/oh_mod.mod
deleted file mode 100644
index 5c53bd2..0000000
Binary files a/flexpart_code/oh_mod.mod and /dev/null differ
diff --git a/flexpart_code/outg_mod.f90 b/flexpart_code/outg_mod.f90
index 4b515d7..6557c37 100644
--- a/flexpart_code/outg_mod.f90
+++ b/flexpart_code/outg_mod.f90
@@ -22,7 +22,6 @@
module outg_mod
implicit none
- save
real,allocatable, dimension (:) :: outheight
real,allocatable, dimension (:) :: outheighthalf
diff --git a/flexpart_code/outg_mod.mod b/flexpart_code/outg_mod.mod
deleted file mode 100644
index fb4d8bd..0000000
Binary files a/flexpart_code/outg_mod.mod and /dev/null differ
diff --git a/flexpart_code/par_mod.f90 b/flexpart_code/par_mod.f90
index 1a4ec0d..d06d066 100644
--- a/flexpart_code/par_mod.f90
+++ b/flexpart_code/par_mod.f90
@@ -34,7 +34,6 @@
module par_mod
implicit none
- save
!****************************************************************
! Parameter defining KIND parameter for "double precision"
@@ -123,12 +122,9 @@ module par_mod
!integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92
!integer,parameter :: nxmax=361,nymax=181,nuvzmax=61,nwzmax=61,nzmax=61
- !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
- integer,parameter :: nxmax=1441,nymax=721,nuvzmax=34,nwzmax=34,nzmax=34 !for GFS as of 2016
- !integer,parameter :: nxmax=1801,nymax=901,nuvzmax=64,nwzmax=64,nzmax=64 !for ECMWF 0.2deg
-
+ integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
!integer,parameter :: nxshift=359 ! for ECMWF
- integer,parameter :: nxshift=0 ! for GFS and pre-processed ECMWF
+ integer,parameter :: nxshift=0 ! for GFS
integer,parameter :: nconvlevmax = nuvzmax-1
integer,parameter :: na = nconvlevmax+1
diff --git a/flexpart_code/par_mod.mod b/flexpart_code/par_mod.mod
deleted file mode 100644
index 9275361..0000000
Binary files a/flexpart_code/par_mod.mod and /dev/null differ
diff --git a/flexpart_code/plumetraj.f90 b/flexpart_code/plumetraj.f90
index f6de13b..0df59c6 100644
--- a/flexpart_code/plumetraj.f90
+++ b/flexpart_code/plumetraj.f90
@@ -68,7 +68,6 @@ subroutine plumetraj(itime)
real :: topo,topocenter,hm(2),hmixi,hmixfract,hmixcenter
real :: pv1(2),pvprof(2),pvi,pvcenter,pvfract,tr(2),tri,tropofract
real :: tropocenter
- character(len=255) :: out_fmtstr
dt1=real(itime-memtime(1))
@@ -233,13 +232,12 @@ subroutine plumetraj(itime)
! Write out results in trajectory data file
!******************************************
- out_fmtstr=''
- write(out_fmtstr,'(i0)') ncluster
- out_fmtstr='(i5,i8,2f9.4,4f8.1,f8.2,4f8.1,3f6.1,'//trim(out_fmtstr)//'(2f8.3,f7.0,f6.1,f8.1))'
-
- write(unitouttraj,out_fmtstr) j,itime-(ireleasestart(j)+ireleaseend(j))/2, &
+ write(unitouttraj,'(i5,i8,2f9.4,4f8.1,f8.2,4f8.1,3f6.1,&
+ &5(2f8.3,f7.0,f6.1,f8.1))')&
+ &j,itime-(ireleasestart(j)+ireleaseend(j))/2, &
xcenter,ycenter,zcenter,topocenter,hmixcenter,tropocenter, &
- pvcenter,rmsdist,rms,zrmsdist,zrms,hmixfract,pvfract,tropofract, &
+ pvcenter,rmsdist,rms,zrmsdist,zrms,hmixfract,pvfract, &
+ tropofract, &
(xclust(k),yclust(k),zclust(k),fclust(k),rmsclust(k), &
k=1,ncluster)
endif
diff --git a/flexpart_code/point_mod.f90 b/flexpart_code/point_mod.f90
index ae25917..0e1a2a1 100644
--- a/flexpart_code/point_mod.f90
+++ b/flexpart_code/point_mod.f90
@@ -22,7 +22,6 @@
module point_mod
implicit none
- save
integer, allocatable, dimension (:) :: ireleasestart
integer, allocatable, dimension (:) :: ireleaseend
diff --git a/flexpart_code/point_mod.mod b/flexpart_code/point_mod.mod
deleted file mode 100644
index 75a712f..0000000
Binary files a/flexpart_code/point_mod.mod and /dev/null differ
diff --git a/flexpart_code/readreleases.f90 b/flexpart_code/readreleases.f90
index 113c4f9..bb317f1 100644
--- a/flexpart_code/readreleases.f90
+++ b/flexpart_code/readreleases.f90
@@ -185,7 +185,7 @@ subroutine readreleases
write (*,*) ' Releasepoints allocated: ', numpoint
do i=1,numpoint
- xmasssave(i)=0.0_dp
+ xmasssave(i)=0.
end do
!now save the information
@@ -385,9 +385,9 @@ subroutine readreleases
if (ldirect.eq.1) then
if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
write(*,*) 'FLEXPART MODEL ERROR'
- write(*,*) 'Release starts before simulation begins or ends (1)'
+ write(*,*) 'Release starts before simulation begins or ends'
write(*,*) 'after simulation stops.'
- write(*,*) 'Make files COMMAND and RELEASES consistent (1).'
+ write(*,*) 'Make files COMMAND and RELEASES consistent.'
stop
endif
ireleasestart(numpoint)=int((jul1-bdate)*86400.)
@@ -395,9 +395,9 @@ subroutine readreleases
else if (ldirect.eq.-1) then
if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
write(*,*) 'FLEXPART MODEL ERROR'
- write(*,*) 'Release starts before simulation begins or ends (2)'
+ write(*,*) 'Release starts before simulation begins or ends'
write(*,*) 'after simulation stops.'
- write(*,*) 'Make files COMMAND and RELEASES consistent (2).'
+ write(*,*) 'Make files COMMAND and RELEASES consistent.'
stop
endif
ireleasestart(numpoint)=int((jul1-bdate)*86400.)
diff --git a/flexpart_code/readwind_gfs.f90 b/flexpart_code/readwind_gfs.f90
index ce532b9..f89e0d5 100644
--- a/flexpart_code/readwind_gfs.f90
+++ b/flexpart_code/readwind_gfs.f90
@@ -63,58 +63,39 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
use grib_api
use par_mod
use com_mod
- use omp_lib
implicit none
-
!HSO new parameters for grib_api
integer :: ifile
integer :: iret
integer :: igrib
-
- !SH
- integer, dimension(1:100000) :: igribind
- integer :: inumfields,iextf,isubfield,chunksize,mynumthreads
- logical, parameter :: parallel = .true.
- !SH W fix
- !0: fix by zeroes
- !1: fix by last valid value
- !2: no fix
- integer, parameter :: wfix = 0
- integer :: numpwmax, match_zlev
- !SH compactification with pointers
- real, dimension(:,:), pointer :: rarr2_pnt
-
integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
!HSO end edits
- real, target :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
- real, target :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
- real, target :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
+ real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
+ real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
+ real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
integer :: ii,indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
! NCEP
integer :: numpt,numpu,numpv,numpw,numprh
real :: help, temp, ew
real :: elev
- real, target :: ulev1(0:nxmax-1,0:nymax-1),vlev1(0:nxmax-1,0:nymax-1)
- real, target :: tlev1(0:nxmax-1,0:nymax-1)
- real, target :: qvh2(0:nxmax-1,0:nymax-1)
+ real :: ulev1(0:nxmax-1,0:nymax-1),vlev1(0:nxmax-1,0:nymax-1)
+ real :: tlev1(0:nxmax-1,0:nymax-1)
+ real :: qvh2(0:nxmax-1,0:nymax-1)
- integer :: i180
+ integer :: i179,i180,i181
! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
!HSO kept isec1, isec2 and zsec4 for consistency with gribex GRIB input
- integer :: isec2(3)
- integer, dimension(:,:), allocatable :: isec1_arr
- real(kind=4), dimension(:), allocatable :: xsec18_arr
- real(kind=4), dimension(:,:), allocatable :: zsec4_arr
-
- real(kind=4) :: xaux,yaux,xaux0,yaux0
- real,parameter :: eps=spacing(2.0_4*360.0_4)
- real(kind=8) :: xauxin,yauxin
- real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
+ integer :: isec1(8),isec2(3)
+ real(kind=4) :: zsec4(jpunp)
+ real(kind=4) :: xaux,yaux,xaux0,yaux0
+ real(kind=8) :: xauxin,yauxin
+ real,parameter :: eps=1.e-4
+ real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
real :: plev1,hlev1,ff10m,fflev1
logical :: hflswitch,strswitch
@@ -124,44 +105,20 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
character(len=20) :: gribFunction = 'readwind_gfs'
- interface
- subroutine readwind_grib2_to_1(parcat_in,parnum_in,typsurf_in,valsurf_in,discipl_in,i6_out,i7_out,x8_out)
- implicit none
- integer, intent(in) :: parcat_in,parnum_in,typsurf_in,valsurf_in,discipl_in
- integer, intent(out) :: i6_out,i7_out
- real(kind=4), intent(out) :: x8_out
- end subroutine
- end interface
-
- write(*,*) 'enter readwind_gfs'
-
- if (parallel) then
- chunksize = 24
- mynumthreads = omp_get_max_threads()
- else
- chunksize = 24
- mynumthreads = 1
- end if
-
- numpwmax = 0
-
hflswitch=.false.
strswitch=.false.
levdiff2=nlev_ec-nwz+1
iumax=0
iwmax=0
- i180=nint(180./dx)
+
! OPENING OF DATA FILE (GRIB CODE)
!HSO
- call grib_open_file(ifile,path(3)(1:length(3)) &
+5 call grib_open_file(ifile,path(3)(1:length(3)) &
//trim(wfname(indj)),'r',iret)
if (iret.ne.GRIB_SUCCESS) then
- write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
- write(*,*) ' #### ',wfname(indj),' #### '
- write(*,*) ' #### IS NOT GRIB FORMAT !!! #### '
- stop 'Execution terminated'
+ goto 888 ! ERROR DETECTED
endif
!turn on support for multi fields messages
call grib_multi_support_on
@@ -171,326 +128,441 @@ subroutine readwind_grib2_to_1(parcat_in,parnum_in,typsurf_in,valsurf_in,discipl
numpv=0
numpw=0
numprh=0
- inumfields=0
- do
- inumfields=inumfields+1
- call grib_new_from_file(ifile,igribind(inumfields),iret)
- if (iret.eq.GRIB_END_OF_FILE) then
- inumfields=inumfields-1
- exit ! EOF DETECTED
- elseif (iret.ne.GRIB_SUCCESS) then
- write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
- write(*,*) ' #### ',wfname(indj),' #### '
- write(*,*) ' #### IS NOT GRIB FORMAT !!! #### '
- stop 'Execution terminated'
- endif
- end do
-
- allocate(isec1_arr(1:chunksize,8))
- allocate(xsec18_arr(1:chunksize))
- allocate(zsec4_arr(1:jpunp,1:chunksize))
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! numfield = 1 w/ preparations !!!!!!!!
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- write(*,*) 'readwind_gfs number of fields:',inumfields
-
- ifield=1
- isubfield=1
+ ifield=0
+10 ifield=ifield+1
+ !
+ ! GET NEXT FIELDS
+ !
+ call grib_new_from_file(ifile,igrib,iret)
+ if (iret.eq.GRIB_END_OF_FILE) then
+ goto 50 ! EOF DETECTED
+ elseif (iret.ne.GRIB_SUCCESS) then
+ goto 888 ! ERROR DETECTED
+ endif
!first see if we read GRIB1 or GRIB2
- call grib_get_int(igribind(ifield),'editionNumber',gribVer,iret)
+ call grib_get_int(igrib,'editionNumber',gribVer,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
-
+
if (gribVer.eq.1) then ! GRIB Edition 1
- call grib_get_int(igribind(ifield),'indicatorOfParameter',isec1_arr(isubfield,6),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'indicatorOfTypeOfLevel',isec1_arr(isubfield,7),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'level',isec1_arr(isubfield,8),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- xsec18_arr(ifield) = real(isec1_arr(isubfield,8))
+
+ !read the grib1 identifiers
+ call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'level',isec1(8),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+
else ! GRIB Edition 2
- call grib_get_int(igribind(ifield),'discipline',discipl,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'parameterCategory',parCat,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'parameterNumber',parNum,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'typeOfFirstFixedSurface',typSurf,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'scaledValueOfFirstFixedSurface', valSurf,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- !convert
- isec1_arr(isubfield,6:8) = -1
- xsec18_arr(isubfield) = -1.0
- call readwind_grib2_to_1(parCat,parNum,typSurf,valSurf,discipl,isec1_arr(isubfield,6),isec1_arr(isubfield,7), &
- xsec18_arr(isubfield))
- isec1_arr(isubfield,8) = nint(xsec18_arr(isubfield))
+
+ !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,'scaledValueOfFirstFixedSurface', &
+ valSurf,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+
+ !convert to grib1 identifiers
+ isec1(6)=-1
+ isec1(7)=-1
+ isec1(8)=-1
+ if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T
+ isec1(6)=11 ! indicatorOfParameter
+ isec1(7)=100 ! indicatorOfTypeOfLevel
+ isec1(8)=valSurf/100 ! level, convert to hPa
+ elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
+ isec1(6)=33 ! indicatorOfParameter
+ isec1(7)=100 ! indicatorOfTypeOfLevel
+ isec1(8)=valSurf/100 ! level, convert to hPa
+ elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.100)) then ! V
+ isec1(6)=34 ! indicatorOfParameter
+ isec1(7)=100 ! indicatorOfTypeOfLevel
+ isec1(8)=valSurf/100 ! level, convert to hPa
+ elseif ((parCat.eq.2).and.(parNum.eq.8).and.(typSurf.eq.100)) then ! W
+ isec1(6)=39 ! indicatorOfParameter
+ isec1(7)=100 ! indicatorOfTypeOfLevel
+ isec1(8)=valSurf/100 ! level, convert to hPa
+ elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.100)) then ! RH
+ isec1(6)=52 ! indicatorOfParameter
+ isec1(7)=100 ! indicatorOfTypeOfLevel
+ isec1(8)=valSurf/100 ! level, convert to hPa
+ elseif ((parCat.eq.1).and.(parNum.eq.1).and.(typSurf.eq.103)) then ! RH2
+ isec1(6)=52 ! indicatorOfParameter
+ isec1(7)=105 ! indicatorOfTypeOfLevel
+ isec1(8)=2
+ elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! T2
+ isec1(6)=11 ! indicatorOfParameter
+ isec1(7)=105 ! indicatorOfTypeOfLevel
+ isec1(8)=2
+ elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! U10
+ isec1(6)=33 ! indicatorOfParameter
+ isec1(7)=105 ! indicatorOfTypeOfLevel
+ isec1(8)=10
+ elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! V10
+ isec1(6)=34 ! indicatorOfParameter
+ isec1(7)=105 ! indicatorOfTypeOfLevel
+ isec1(8)=10
+ elseif ((parCat.eq.3).and.(parNum.eq.1).and.(typSurf.eq.101)) then ! SLP
+ isec1(6)=2 ! indicatorOfParameter
+ isec1(7)=102 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then ! SP
+ isec1(6)=1 ! indicatorOfParameter
+ isec1(7)=1 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ elseif ((parCat.eq.1).and.(parNum.eq.13).and.(typSurf.eq.1)) then ! SNOW
+ isec1(6)=66 ! indicatorOfParameter
+ isec1(7)=1 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.104)) then ! T sigma 0
+ isec1(6)=11 ! indicatorOfParameter
+ isec1(7)=107 ! indicatorOfTypeOfLevel
+ isec1(8)=0.995 ! lowest sigma level
+ elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.104)) then ! U sigma 0
+ isec1(6)=33 ! indicatorOfParameter
+ isec1(7)=107 ! indicatorOfTypeOfLevel
+ isec1(8)=0.995 ! lowest sigma level
+ elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.104)) then ! V sigma 0
+ isec1(6)=34 ! indicatorOfParameter
+ isec1(7)=107 ! indicatorOfTypeOfLevel
+ isec1(8)=0.995 ! lowest sigma level
+ elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
+ isec1(6)=7 ! indicatorOfParameter
+ isec1(7)=1 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
+ .and.(discipl.eq.2)) then ! LSM
+ isec1(6)=81 ! indicatorOfParameter
+ isec1(7)=1 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ elseif ((parCat.eq.3).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! BLH
+ isec1(6)=221 ! indicatorOfParameter
+ isec1(7)=1 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ elseif ((parCat.eq.1).and.(parNum.eq.7).and.(typSurf.eq.1)) then ! LSP/TP
+ isec1(6)=62 ! indicatorOfParameter
+ isec1(7)=1 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ elseif ((parCat.eq.1).and.(parNum.eq.196).and.(typSurf.eq.1)) then ! CP
+ isec1(6)=63 ! indicatorOfParameter
+ isec1(7)=1 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ endif
+
endif ! gribVer
- if (isec1_arr(isubfield,6).ne.-1) then
- !get the size and data of the values array
- call grib_get_real4_array(igribind(ifield),'values',zsec4_arr(1:jpunp,isubfield),iret)
+
+ if (isec1(6).ne.-1) then
+ ! get the size and data of the values array
+ call grib_get_real4_array(igrib,'values',zsec4,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
endif
- call grib_get_int(igribind(ifield),'numberOfPointsAlongAParallel',isec2(2),iret)
+ if(ifield.eq.1) then
+
+ !get the required fields from section 2
+ !store compatible to gribex input
+ call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
+ isec2(2),iret)
call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'numberOfPointsAlongAMeridian',isec2(3),iret)
+ call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
+ isec2(3),iret)
call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_real8(igribind(ifield),'longitudeOfFirstGridPointInDegrees',xauxin,iret)
+ call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
+ xauxin,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_real8(igribind(ifield),'latitudeOfLastGridPointInDegrees',yauxin,iret)
+ call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
+ yauxin,iret)
call grib_check(iret,gribFunction,gribErrorMsg)
xaux=xauxin+real(nxshift)*dx
yaux=yauxin
-
+
! CHECK GRID SPECIFICATIONS
-
- if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
- if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
- if(xaux.eq.0.) xaux=-180.0 ! NCEP DATA
- xaux0=xlon0
- yaux0=ylat0
- if(xaux.lt.0.) xaux=xaux+360.
- if(yaux.lt.0.) yaux=yaux+360.
- if(xaux0.lt.0.) xaux0=xaux0+360.
- if(yaux0.lt.0.) yaux0=yaux0+360.
- if(abs(xaux-xaux0).gt.eps .or. abs(yaux-yaux0).gt.eps) then
- write (*, *) 'xaux/xaux0/yaux/yaux0: xaux, xaux0, yaux, yaux0'
- stop 'READWIND: LOWER LEFT LONGITUDE or LATITUDE NOT CONSISTENT'
+
+ if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
+ if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
+ if(xaux.eq.0.) xaux=-179.0 ! NCEP DATA
+ xaux0=xlon0
+ yaux0=ylat0
+ if(xaux.lt.0.) xaux=xaux+360.
+ if(yaux.lt.0.) yaux=yaux+360.
+ if(xaux0.lt.0.) xaux0=xaux0+360.
+ if(yaux0.lt.0.) yaux0=yaux0+360.
+ if(abs(xaux-xaux0).gt.eps) &
+ stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
+ if(abs(yaux-yaux0).gt.eps) &
+ stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
+ endif
+ !HSO end of edits
+
+ i179=nint(179./dx)
+ 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
+
+ if (isec1(6).ne.-1) then
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! OK done now do real stuff
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- write(*,*) 'grib main readin loop: ',mynumthreads,' parallel threads process chunks of ',chunksize,' fields'
-
- do iextf=1,inumfields,chunksize
-!$omp parallel default(none) &
-!$omp shared(isec1_arr,xsec18_arr,zsec4_arr,iextf,inumfields,igribind,chunksize) &
-!$omp private(ifield,isubfield,iret,discipl,gribVer,parCat,parNum, &
-!$omp typSurf,ValSurf,gribfunction,griberrormsg) num_threads(mynumthreads)
-!$omp do schedule(guided)
- do ifield=iextf,min(inumfields,iextf+chunksize-1)
- isubfield=ifield-iextf+1
- !write(*,*) isubfield
- !next field
- !HERE, isec1, xsec18 are read
- !and zsec4
- !and (locally needed) discipl, parCat, parNum, typSurf, valSurf
-
- !write(*,*) ifield
- !write(*,*) igribind(ifield),igribind(ifield+1)
- !first see if we read GRIB1 or GRIB2
- call grib_get_int(igribind(ifield),'editionNumber',gribVer,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
-
- if (gribVer.eq.1) then ! GRIB Edition 1
- !read the grib1 identifiers
- call grib_get_int(igribind(ifield),'indicatorOfParameter',isec1_arr(isubfield,6),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'indicatorOfTypeOfLevel',isec1_arr(isubfield,7),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'level',isec1_arr(isubfield,8),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- !JMA / SH: isec1(8) not evaluated any more below
- !b/c with GRIB 2 this may be a real variable
- xsec18_arr(isubfield) = real(isec1_arr(isubfield,8))
- else ! GRIB Edition 2
- !read the grib2 identifiers
- call grib_get_int(igribind(ifield),'discipline',discipl,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'parameterCategory',parCat,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'parameterNumber',parNum,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'typeOfFirstFixedSurface',typSurf,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- call grib_get_int(igribind(ifield),'scaledValueOfFirstFixedSurface', valSurf,iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- isec1_arr(isubfield,6:8) = -1
- xsec18_arr(isubfield) = -1.0
- !convert
- call readwind_grib2_to_1(parCat,parNum,typSurf,valSurf,discipl,isec1_arr(isubfield,6),isec1_arr(isubfield,7), &
- xsec18_arr(isubfield))
- endif ! gribVer
- if (isec1_arr(isubfield,6).ne.-1) then
- if ((isec1_arr(isubfield,6).eq.011 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.033 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.034 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.052 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.001 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.039 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.066 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.002 .and. isec1_arr(isubfield,7).eq.102) .or. &
- (isec1_arr(isubfield,6).eq.071 .and. isec1_arr(isubfield,7).eq.244) .or. &
- (isec1_arr(isubfield,6).eq.033 .and. isec1_arr(isubfield,7).eq.105 .and. nint(xsec18_arr(isubfield)).eq.10) .or. &
- (isec1_arr(isubfield,6).eq.034 .and. isec1_arr(isubfield,7).eq.105 .and. nint(xsec18_arr(isubfield)).eq.10) .or. &
- (isec1_arr(isubfield,6).eq.011 .and. isec1_arr(isubfield,7).eq.105 .and. nint(xsec18_arr(isubfield)).eq.2) .or. &
- (isec1_arr(isubfield,6).eq.017 .and. isec1_arr(isubfield,7).eq.105 .and. nint(xsec18_arr(isubfield)).eq.2) .or. &
- (isec1_arr(isubfield,6).eq.062 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.063 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.007 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.081 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.221 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.052 .and. isec1_arr(isubfield,7).eq.105 .and. nint(xsec18_arr(isubfield)).eq.02) .or. &
- (isec1_arr(isubfield,6).eq.011 .and. isec1_arr(isubfield,7).eq.107) .or. &
- (isec1_arr(isubfield,6).eq.033 .and. isec1_arr(isubfield,7).eq.107) .or. &
- (isec1_arr(isubfield,6).eq.034 .and. isec1_arr(isubfield,7).eq.107)) then
- !write(*,*) 'read',isec1_arr(isubfield,6),isec1_arr(isubfield,7)
-
- !get the size and data of the values array
- call grib_get_real4_array(igribind(ifield),'values',zsec4_arr(1:jpunp,isubfield),iret)
- call grib_check(iret,gribFunction,gribErrorMsg)
- end if
+ do j=0,nymin1
+ do i=0,nxfield-1
+ if((isec1(6).eq.011).and.(isec1(7).eq.100)) then
+ ! TEMPERATURE
+ if((i.eq.0).and.(j.eq.0)) then
+ do ii=1,nuvz
+ if ((isec1(8)*100.0).eq.akz(ii)) numpt=ii
+ end do
+ endif
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ tth(i179+i,j,numpt,n)=help
+ else
+ tth(i-i181,j,numpt,n)=help
+ endif
endif
- end do
-!$omp end do
-!$omp end parallel
-
- do ifield=iextf,min(inumfields,iextf+chunksize-1)
- isubfield=ifield-iextf+1
- if (isec1_arr(isubfield,6).ne.-1) then
- if ((isec1_arr(isubfield,6).eq.011 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.033 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.034 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.052 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.039 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.001 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.066 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.002 .and. isec1_arr(isubfield,7).eq.102) .or. &
- (isec1_arr(isubfield,6).eq.071 .and. isec1_arr(isubfield,7).eq.244) .or. &
- (isec1_arr(isubfield,6).eq.033 .and. isec1_arr(isubfield,7).eq.105 .and. nint(xsec18_arr(isubfield)).eq.10) .or. &
- (isec1_arr(isubfield,6).eq.034 .and. isec1_arr(isubfield,7).eq.105 .and. nint(xsec18_arr(isubfield)).eq.10) .or. &
- (isec1_arr(isubfield,6).eq.011 .and. isec1_arr(isubfield,7).eq.105 .and. nint(xsec18_arr(isubfield)).eq.2) .or. &
- (isec1_arr(isubfield,6).eq.017 .and. isec1_arr(isubfield,7).eq.105 .and. nint(xsec18_arr(isubfield)).eq.2) .or. &
- (isec1_arr(isubfield,6).eq.062 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.063 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.007 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.081 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.221 .and. isec1_arr(isubfield,7).eq.001) .or. &
- (isec1_arr(isubfield,6).eq.052 .and. isec1_arr(isubfield,7).eq.105 .and. nint(xsec18_arr(isubfield)).eq.2) .or. &
- (isec1_arr(isubfield,6).eq.011 .and. isec1_arr(isubfield,7).eq.107) .or. &
- (isec1_arr(isubfield,6).eq.033 .and. isec1_arr(isubfield,7).eq.107) .or. &
- (isec1_arr(isubfield,6).eq.034 .and. isec1_arr(isubfield,7).eq.107)) then
-
- if ((isec1_arr(isubfield,6).eq.011 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.033 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.034 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.052 .and. isec1_arr(isubfield,7).eq.100) .or. &
- (isec1_arr(isubfield,6).eq.039 .and. isec1_arr(isubfield,7).eq.100)) then
+ if((isec1(6).eq.033).and.(isec1(7).eq.100)) then
+ ! U VELOCITY
+ if((i.eq.0).and.(j.eq.0)) then
do ii=1,nuvz
- if (abs(xsec18_arr(isubfield)*100.0-akz(ii)) < &
- 10.0*max(spacing(akz(ii)),spacing(xsec18_arr(isubfield)*100.0))) then
- match_zlev=ii
- end if
+ if ((isec1(8)*100.0).eq.akz(ii)) numpu=ii
end do
- end if
-
- ! TEMPERATURE
- if ((isec1_arr(isubfield,6).eq.011).and.(isec1_arr(isubfield,7).eq.100)) rarr2_pnt => tth(:,:,match_zlev,n)
- ! U VELOCITY
- if ((isec1_arr(isubfield,6).eq.033).and.(isec1_arr(isubfield,7).eq.100)) then
- rarr2_pnt => uuh(:,:,match_zlev)
- !NCEP ISOBARIC LEVELS
- iumax=iumax+1
- end if
- ! V VELOCITY
- if ((isec1_arr(isubfield,6).eq.034).and.(isec1_arr(isubfield,7).eq.100)) rarr2_pnt => vvh(:,:,match_zlev)
- ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
- if ((isec1_arr(isubfield,6).eq.052).and.(isec1_arr(isubfield,7).eq.100)) rarr2_pnt => qvh(:,:,match_zlev,n)
- ! W VELOCITY
- if ((isec1_arr(isubfield,6).eq.039).and.(isec1_arr(isubfield,7).eq.100)) then
- numpwmax=max(match_zlev,numpwmax)
- rarr2_pnt => wwh(:,:,match_zlev)
- end if
- ! SURFACE PRESSURE
- if ((isec1_arr(isubfield,6).eq.001).and.(isec1_arr(isubfield,7).eq.001)) rarr2_pnt => ps(:,:,1,n)
- ! SNOW DEPTH
- if ((isec1_arr(isubfield,6).eq.066).and.(isec1_arr(isubfield,7).eq.001)) rarr2_pnt => sd(:,:,1,n)
- ! MEAN SEA LEVEL PRESSURE
- if ((isec1_arr(isubfield,6).eq.002).and.(isec1_arr(isubfield,7).eq.102)) rarr2_pnt => msl(:,:,1,n)
- ! TOTAL CLOUD COVER
- if ((isec1_arr(isubfield,6).eq.071).and.(isec1_arr(isubfield,7).eq.244)) rarr2_pnt => tcc(:,:,1,n)
- ! 10M U VELOCITY
- if ((isec1_arr(isubfield,6).eq.033).and.(isec1_arr(isubfield,7).eq.105).and.(nint(xsec18_arr(isubfield)).eq.10)) &
- rarr2_pnt => u10(:,:,1,n)
- ! 10M V VELOCITY
- if ((isec1_arr(isubfield,6).eq.034).and.(isec1_arr(isubfield,7).eq.105).and.(nint(xsec18_arr(isubfield)).eq.10)) &
- rarr2_pnt => v10(:,:,1,n)
- ! 2M TEMPERATURE
- if ((isec1_arr(isubfield,6).eq.011).and.(isec1_arr(isubfield,7).eq.105).and.(nint(xsec18_arr(isubfield)).eq.2)) &
- rarr2_pnt => tt2(:,:,1,n)
- ! 2M DEW POINT TEMPERATURE
- if ((isec1_arr(isubfield,6).eq.017).and.(isec1_arr(isubfield,7).eq.105).and.(nint(xsec18_arr(isubfield)).eq.2)) &
- rarr2_pnt => td2(:,:,1,n)
- ! LARGE SCALE PREC.
- if ((isec1_arr(isubfield,6).eq.062).and.(isec1_arr(isubfield,7).eq.001)) rarr2_pnt => lsprec(:,:,1,n)
- ! CONVECTIVE PREC.
- if ((isec1_arr(isubfield,6).eq.063).and.(isec1_arr(isubfield,7).eq.001)) rarr2_pnt => convprec(:,:,1,n)
- ! TOPOGRAPHY
- if ((isec1_arr(isubfield,6).eq.007).and.(isec1_arr(isubfield,7).eq.001)) then
- rarr2_pnt => oro(:,:)
- excessoro(i180:i180+min(nxfield,i180)-1,:) = 0.0
- if (nxfield-1>=i180) excessoro(0:nxfield-1-i180,:) = 0.0
- end if
- ! LAND SEA MASK
- if ((isec1_arr(isubfield,6).eq.081).and.(isec1_arr(isubfield,7).eq.001)) rarr2_pnt => lsm(:,:)
- ! MIXING HEIGHT
- if ((isec1_arr(isubfield,6).eq.221).and.(isec1_arr(isubfield,7).eq.001)) rarr2_pnt => hmix(:,:,1,n)
- ! 2M REALTIVE HUMIDITY
- if ((isec1_arr(isubfield,6).eq.052).and.(isec1_arr(isubfield,7).eq.105).and.(nint(xsec18_arr(isubfield)).eq.2)) &
- rarr2_pnt => qvh2(:,:)
- ! TEMPERATURE LOWEST SIGMA LEVEL
- if ((isec1_arr(isubfield,6).eq.011).and.(isec1_arr(isubfield,7).eq.107)) rarr2_pnt => tlev1(:,:)
- ! U VELOCITY LOWEST SIGMA LEVEL
- if ((isec1_arr(isubfield,6).eq.033).and.(isec1_arr(isubfield,7).eq.107)) rarr2_pnt => ulev1(:,:)
- ! V VELOCITY LOWEST SIGMA LEVEL
- if ((isec1_arr(isubfield,6).eq.034).and.(isec1_arr(isubfield,7).eq.107)) rarr2_pnt => vlev1(:,:)
-
- do j=0,nymin1
- rarr2_pnt(i180+1:i180+min(nxfield,i180), j+1) = &
- zsec4_arr(nxfield*(ny-j-1)+1:nxfield*(ny-j-1)+min(nxfield,i180), isubfield)
- end do
- if (nxfield-1>=i180) then
- do j=0,nymin1
- rarr2_pnt(1:nxfield-i180, j+1) = zsec4_arr(nxfield*(ny-j-1)+i180+1:nxfield*(ny-j-1)+nxfield, isubfield)
+ endif
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ uuh(i179+i,j,numpu)=help
+ else
+ uuh(i-i181,j,numpu)=help
+ endif
+ endif
+ if((isec1(6).eq.034).and.(isec1(7).eq.100)) then
+ ! V VELOCITY
+ if((i.eq.0).and.(j.eq.0)) then
+ do ii=1,nuvz
+ if ((isec1(8)*100.0).eq.akz(ii)) numpv=ii
end do
- end if
-
- nullify(rarr2_pnt)
- end if
+ endif
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ vvh(i179+i,j,numpv)=help
+ else
+ vvh(i-i181,j,numpv)=help
+ endif
+ endif
+ if((isec1(6).eq.052).and.(isec1(7).eq.100)) then
+ ! RELATIVE HUMIDITY -> CONVERT TO SPECIFIC HUMIDITY LATER
+ if((i.eq.0).and.(j.eq.0)) then
+ do ii=1,nuvz
+ if ((isec1(8)*100.0).eq.akz(ii)) numprh=ii
+ end do
+ endif
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ qvh(i179+i,j,numprh,n)=help
+ else
+ qvh(i-i181,j,numprh,n)=help
+ endif
+ endif
+ if((isec1(6).eq.001).and.(isec1(7).eq.001)) then
+ ! SURFACE PRESSURE
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ ps(i179+i,j,1,n)=help
+ else
+ ps(i-i181,j,1,n)=help
+ endif
+ endif
+ if((isec1(6).eq.039).and.(isec1(7).eq.100)) then
+ ! W VELOCITY
+ if((i.eq.0).and.(j.eq.0)) then
+ do ii=1,nuvz
+ if ((isec1(8)*100.0).eq.akz(ii)) numpw=ii
+ end do
+ endif
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ wwh(i179+i,j,numpw)=help
+ else
+ wwh(i-i181,j,numpw)=help
+ endif
+ endif
+ if((isec1(6).eq.066).and.(isec1(7).eq.001)) then
+ ! SNOW DEPTH
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ sd(i179+i,j,1,n)=help
+ else
+ sd(i-i181,j,1,n)=help
+ endif
+ endif
+ if((isec1(6).eq.002).and.(isec1(7).eq.102)) then
+ ! MEAN SEA LEVEL PRESSURE
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ msl(i179+i,j,1,n)=help
+ else
+ msl(i-i181,j,1,n)=help
+ endif
+ endif
+ if((isec1(6).eq.071).and.(isec1(7).eq.244)) then
+ ! TOTAL CLOUD COVER
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ tcc(i179+i,j,1,n)=help
+ else
+ tcc(i-i181,j,1,n)=help
+ endif
+ endif
+ if((isec1(6).eq.033).and.(isec1(7).eq.105).and. &
+ (isec1(8).eq.10)) then
+ ! 10 M U VELOCITY
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ u10(i179+i,j,1,n)=help
+ else
+ u10(i-i181,j,1,n)=help
+ endif
+ endif
+ if((isec1(6).eq.034).and.(isec1(7).eq.105).and. &
+ (isec1(8).eq.10)) then
+ ! 10 M V VELOCITY
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ v10(i179+i,j,1,n)=help
+ else
+ v10(i-i181,j,1,n)=help
+ endif
+ endif
+ if((isec1(6).eq.011).and.(isec1(7).eq.105).and. &
+ (isec1(8).eq.02)) then
+ ! 2 M TEMPERATURE
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ tt2(i179+i,j,1,n)=help
+ else
+ tt2(i-i181,j,1,n)=help
+ endif
+ endif
+ if((isec1(6).eq.017).and.(isec1(7).eq.105).and. &
+ (isec1(8).eq.02)) then
+ ! 2 M DEW POINT TEMPERATURE
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ td2(i179+i,j,1,n)=help
+ else
+ td2(i-i181,j,1,n)=help
+ endif
+ endif
+ if((isec1(6).eq.062).and.(isec1(7).eq.001)) then
+ ! LARGE SCALE PREC.
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ lsprec(i179+i,j,1,n)=help
+ else
+ lsprec(i-i181,j,1,n)=help
+ endif
+ endif
+ if((isec1(6).eq.063).and.(isec1(7).eq.001)) then
+ ! CONVECTIVE PREC.
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ convprec(i179+i,j,1,n)=help
+ else
+ convprec(i-i181,j,1,n)=help
+ endif
+ endif
+ if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
+ ! TOPOGRAPHY
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ oro(i179+i,j)=help
+ excessoro(i179+i,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
+ else
+ oro(i-i181,j)=help
+ excessoro(i-i181,j)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
+ endif
+ endif
+ if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
+ ! LAND SEA MASK
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ lsm(i179+i,j)=help
+ else
+ lsm(i-i181,j)=help
+ endif
+ endif
+ if((isec1(6).eq.221).and.(isec1(7).eq.001)) then
+ ! MIXING HEIGHT
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ hmix(i179+i,j,1,n)=help
+ else
+ hmix(i-i181,j,1,n)=help
+ endif
+ endif
+ if((isec1(6).eq.052).and.(isec1(7).eq.105).and. &
+ (isec1(8).eq.02)) then
+ ! 2 M RELATIVE HUMIDITY
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ qvh2(i179+i,j)=help
+ else
+ qvh2(i-i181,j)=help
+ endif
+ endif
+ if((isec1(6).eq.011).and.(isec1(7).eq.107)) then
+ ! TEMPERATURE LOWEST SIGMA LEVEL
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ tlev1(i179+i,j)=help
+ else
+ tlev1(i-i181,j)=help
+ endif
+ endif
+ if((isec1(6).eq.033).and.(isec1(7).eq.107)) then
+ ! U VELOCITY LOWEST SIGMA LEVEL
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ ulev1(i179+i,j)=help
+ else
+ ulev1(i-i181,j)=help
+ endif
+ endif
+ if((isec1(6).eq.034).and.(isec1(7).eq.107)) then
+ ! V VELOCITY LOWEST SIGMA LEVEL
+ help=zsec4(nxfield*(ny-j-1)+i+1)
+ if(i.le.i180) then
+ vlev1(i179+i,j)=help
+ else
+ vlev1(i-i181,j)=help
+ endif
endif
- call grib_release(igribind(ifield))
- end do !! READ NEXT LEVEL OR PARAMETER
+ end do
end do
-
- deallocate(isec1_arr)
- deallocate(xsec18_arr)
- deallocate(zsec4_arr)
-
- !HSO close grib file
- call grib_close_file(ifile)
+ endif
+ if((isec1(6).eq.33).and.(isec1(7).eq.100)) then
+ ! NCEP ISOBARIC LEVELS
+ iumax=iumax+1
+ endif
- !fix W velocities at uppermost levles
- if (numpwmax= nuvz:
- ! JMA: In that case: k is set to nuvz
-
- if (k .gt. nuvz) k = nuvz ! JMA
zl1=zold
theta1=thetaold
diff --git a/flexpart_code/unc_mod.f90 b/flexpart_code/unc_mod.f90
index ecf2f9e..3ba2de3 100644
--- a/flexpart_code/unc_mod.f90
+++ b/flexpart_code/unc_mod.f90
@@ -22,7 +22,6 @@
module unc_mod
implicit none
- save
real,allocatable, dimension (:,:,:,:,:,:,:) :: gridunc
real,allocatable, dimension (:,:,:,:,:,:,:) :: griduncn
diff --git a/flexpart_code/unc_mod.mod b/flexpart_code/unc_mod.mod
deleted file mode 100644
index 4db589e..0000000
Binary files a/flexpart_code/unc_mod.mod and /dev/null differ
diff --git a/flexpart_code/verttransform.f90 b/flexpart_code/verttransform.f90
index b2a52db..8d7c323 100644
--- a/flexpart_code/verttransform.f90
+++ b/flexpart_code/verttransform.f90
@@ -82,9 +82,7 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
- real :: invcos, invdz
real,parameter :: const=r_air/ga
- real,parameter :: inv_one_eighty = 1. / 180.
logical :: init = .true.
@@ -228,7 +226,7 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
do kz=2,nwz-1
- wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))*0.5
+ wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
end do
wzlev(nwz)=wzlev(nwz-1)+ &
uvzlev(nuvz)-uvzlev(nuvz-1)
@@ -277,34 +275,34 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
rho(ix,jy,nz,n)=rhoh(nuvz)
kmin=2
do iz=2,nz-1
- if(height(iz).gt.uvzlev(nuvz)) then
- uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
- vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
- tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
- qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
- pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
- rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
- else
- do kz=kmin,nuvz
- if ((height(iz).gt.uvzlev(kz-1)).and. &
- (height(iz).le.uvzlev(kz))) then
- dz1=height(iz)-uvzlev(kz-1)
- dz2=uvzlev(kz)-height(iz)
- dz=dz1+dz2
- invdz = 1./dz
- uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)*invdz
- vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)*invdz
- tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
- +tth(ix,jy,kz,n)*dz1)*invdz
- qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
- +qvh(ix,jy,kz,n)*dz1)*invdz
- pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)*invdz
- rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)*invdz
- kmin=kz
- exit
- endif
- end do
- end if
+ do kz=kmin,nuvz
+ if(height(iz).gt.uvzlev(nuvz)) then
+ uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
+ vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
+ tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
+ qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
+ pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
+ rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
+ goto 30
+ endif
+ if ((height(iz).gt.uvzlev(kz-1)).and. &
+ (height(iz).le.uvzlev(kz))) then
+ dz1=height(iz)-uvzlev(kz-1)
+ dz2=uvzlev(kz)-height(iz)
+ dz=dz1+dz2
+ uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
+ vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
+ tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
+ +tth(ix,jy,kz,n)*dz1)/dz
+ qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
+ +qvh(ix,jy,kz,n)*dz1)/dz
+ pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
+ rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
+ kmin=kz
+ goto 30
+ endif
+ end do
+30 continue
end do
@@ -351,13 +349,12 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!****************************************************************
do jy=1,ny-2
- invcos = 1./cos((real(jy)*dy+ylat0)*pi180)
do ix=1,nx-2
kmin=2
do iz=2,nz-1
- ui=uu(ix,jy,iz,n)*dxconst*invcos
+ ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
vi=vv(ix,jy,iz,n)*dyconst
do kz=kmin,nz
@@ -378,12 +375,12 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
ixp=ix+1
jyp=jy+1
- dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))*0.5
- dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))*0.5
+ dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
+ dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
- dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))*0.5
- dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))*0.5
+ dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
+ dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
@@ -399,11 +396,11 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!*******************************************************************
if (nglobal) then
- do iz=1,nz
- do jy=int(switchnorthg)-2,nymin1
- ylat=ylat0+real(jy)*dy
- do ix=0,nxmin1
- xlon=xlon0+real(ix)*dx
+ do jy=int(switchnorthg)-2,nymin1
+ ylat=ylat0+real(jy)*dy
+ do ix=0,nxmin1
+ xlon=xlon0+real(ix)*dx
+ do iz=1,nz
call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
vvpol(ix,jy,iz,n))
@@ -419,7 +416,7 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
! AMSnauffer Nov 18 2004 Added check for case vv=0
!
xlon=xlon0+real(nx/2-1)*dx
- xlonr=xlon*pi*inv_one_eighty
+ xlonr=xlon*pi/180.
ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
vv(nx/2-1,nymin1,iz,n)**2)
if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
@@ -436,7 +433,7 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
xlon=180.0
- xlonr=xlon*pi*inv_one_eighty
+ xlonr=xlon*pi/180.
ylat=90.0
uuaux=-ffpol*sin(xlonr+ddpol)
vvaux=-ffpol*cos(xlonr+ddpol)
@@ -475,11 +472,11 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!*******************************************************************
if (sglobal) then
- do iz=1,nz
- do jy=0,int(switchsouthg)+3
- ylat=ylat0+real(jy)*dy
- do ix=0,nxmin1
- xlon=xlon0+real(ix)*dx
+ do jy=0,int(switchsouthg)+3
+ ylat=ylat0+real(jy)*dy
+ do ix=0,nxmin1
+ xlon=xlon0+real(ix)*dx
+ do iz=1,nz
call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
vvpol(ix,jy,iz,n))
@@ -494,7 +491,7 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
! AMSnauffer Nov 18 2004 Added check for case vv=0
!
xlon=xlon0+real(nx/2-1)*dx
- xlonr=xlon*pi*inv_one_eighty
+ xlonr=xlon*pi/180.
ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
vv(nx/2-1,0,iz,n)**2)
if (vv(nx/2-1,0,iz,n).lt.0.) then
@@ -511,16 +508,17 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
xlon=180.0
- xlonr=xlon*pi*inv_one_eighty
+ xlonr=xlon*pi/180.
ylat=-90.0
uuaux=+ffpol*sin(xlonr-ddpol)
vvaux=-ffpol*cos(xlonr-ddpol)
call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
vvpolaux)
+ jy=0
do ix=0,nxmin1
- uupol(ix,0,iz,n)=uupolaux
- vvpol(ix,0,iz,n)=vvpolaux
+ uupol(ix,jy,iz,n)=uupolaux
+ vvpol(ix,jy,iz,n)=vvpolaux
end do
end do
@@ -530,12 +528,14 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
do iz=1,nz
wdummy=0.
+ jy=1
do ix=0,nxmin1
- wdummy=wdummy+ww(ix,1,iz,n)
+ wdummy=wdummy+ww(ix,jy,iz,n)
end do
wdummy=wdummy/real(nx)
+ jy=0
do ix=0,nxmin1
- ww(ix,0,iz,n)=wdummy
+ ww(ix,jy,iz,n)=wdummy
end do
end do
endif
@@ -544,14 +544,14 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
! create a cloud and rainout/washout field, clouds occur where rh>80%
! total cloudheight is stored at level 0
- do kz_inv=1,nz-1
- kz=nz-kz_inv+1
- do jy=0,nymin1
- do ix=0,nxmin1
- rain_cloud_above=0
- lsp=lsprec(ix,jy,1,n)
- convp=convprec(ix,jy,1,n)
- cloudsh(ix,jy,n)=0
+ do jy=0,nymin1
+ do ix=0,nxmin1
+ rain_cloud_above=0
+ lsp=lsprec(ix,jy,1,n)
+ convp=convprec(ix,jy,1,n)
+ cloudsh(ix,jy,n)=0
+ do kz_inv=1,nz-1
+ kz=nz-kz_inv+1
pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
clouds(ix,jy,kz,n)=0
diff --git a/flexpart_code/verttransform_gfs.f90 b/flexpart_code/verttransform_gfs.f90
index e1905e5..4826b3e 100644
--- a/flexpart_code/verttransform_gfs.f90
+++ b/flexpart_code/verttransform_gfs.f90
@@ -63,24 +63,15 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
use par_mod
use com_mod
use cmapf_mod
- use omp_lib
implicit none
- integer :: ix,jy,kz,iz,n,ix1,jy1,ixp,jyp,ixm,jym
+ integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
integer :: rain_cloud_above,kz_inv
real :: f_qvsat,pressure
real :: rh,lsp,convp
+ real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
- !SH
- real :: tmp_real
- real :: tv_arr(1:nuvz),pint_arr(1:nuvz)
- real :: rhoh_arr(0:nxmax-1,0:nymax-1,nuvzmax)
- real :: pinmconv_arr(0:nxmax-1,0:nymax-1,nzmax)
- integer :: llev_arr(0:nxmax-1,0:nymax-1),hind_above_arr(0:nxmin1,0:nymin1)
- integer, dimension(0:nxmax-1,0:nymax-1,1:nz) :: match_height_arr
- integer :: loopcnt
-
real :: xlon,ylat,xlonr,dzdx,dzdy
real :: dzdx1,dzdx2,dzdy1,dzdy2
real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
@@ -88,15 +79,15 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
- real :: uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
+ real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
+ real,parameter :: const=r_air/ga
! NCEP version
- ! in this case, nz==nuvz is crucial for correctness and assumed below.
+ integer :: llev, i
+
+ logical :: init = .true.
+
- if (nz/=nuvz) stop 'FATAL FLEXPART-GFS, nz != nuvz'
- write(*,*) 'enter verttransform...'
- uvwzlev(:,:,:) = 0.0
-
!*************************************************************************
! If verttransform is called the first time, initialize heights of the *
! z levels in meter. The heights are the heights of model levels, where *
@@ -106,361 +97,381 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
! Unlike in the eta system, no difference between heights for u,v and *
! heights for w exists. *
!*************************************************************************
- if (.not. verttransform_inited) then
- ! Search for a point with high surface pressure (i.e. not above significant topography)
- ! Then, use this point to construct a reference z profile, to be used at all times
- !*****************************************************************************
+ if (init) then
+
+ ! Search for a point with high surface pressure (i.e. not above significant topography)
+ ! Then, use this point to construct a reference z profile, to be used at all times
+ !*****************************************************************************
- yloop_init: do jy=0,nymin1
+ do jy=0,nymin1
do ix=0,nxmin1
if (ps(ix,jy,1,n).gt.100000.) then
ixm=ix
jym=jy
- exit yloop_init
+ goto 3
endif
end do
- end do yloop_init
- write(*,*) 'verttransform start: init ixm/jym:',ixm,jym
+ end do
+3 continue
- height(1) = 0.
- tv_arr(1) = tt2(ixm,jym,1,n) * (1. + 0.378*ew(td2(ixm,jym,1,n))/ps(ixm,jym,1,n))
- tv_arr(2:nuvz) = tth(ixm,jym,2:nuvz,n) * (1. + 0.608*qvh(ixm,jym,2:nuvz,n))
-
- pint_arr(1) = ps(ixm,jym,1,n)
- pint_arr(2:nuvz) = akz(2:nuvz)+bkz(2:nuvz)*ps(ixm,jym,1,n)
+
+ tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
+ ps(ixm,jym,1,n))
+ pold=ps(ixm,jym,1,n)
+ height(1)=0.
do kz=2,nuvz
- if (abs(tv_arr(kz)-tv_arr(kz-1)).gt.0.2) then
- height(kz) = height(kz-1) + (r_air/ga)*log(pint_arr(kz-1)/pint_arr(kz)) * (tv_arr(kz)-tv_arr(kz-1)) / &
- log(tv_arr(kz)/tv_arr(kz-1))
+ pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
+ tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
+
+
+ ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
+ ! upon the transformation to z levels. In order to save computer memory, this is
+ ! not done anymore in the standard version. However, this option can still be
+ ! switched on by replacing the following lines with those below, that are
+ ! currently commented out.
+ ! Note that two more changes are necessary in this subroutine below.
+ ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
+ !*****************************************************************************
+
+ if (abs(tv-tvold).gt.0.2) then
+ height(kz)= &
+ height(kz-1)+const*log(pold/pint)* &
+ (tv-tvold)/log(tv/tvold)
else
- height(kz) = height(kz-1) + (r_air/ga)*log(pint_arr(kz-1)/pint_arr(kz)) * tv_arr(kz)
+ height(kz)=height(kz-1)+ &
+ const*log(pold/pint)*tv
endif
+
+ ! Switch on following lines to use doubled vertical resolution
+ !*************************************************************
+ ! if (abs(tv-tvold).gt.0.2) then
+ ! height((kz-1)*2)=
+ ! + height(max((kz-2)*2,1))+const*log(pold/pint)*
+ ! + (tv-tvold)/log(tv/tvold)
+ ! else
+ ! height((kz-1)*2)=height(max((kz-2)*2,1))+
+ ! + const*log(pold/pint)*tv
+ ! endif
+ ! End doubled vertical resolution
+
+ tvold=tv
+ pold=pint
end do
- do kz=nuvz+1,ubound(height,1)
- height(kz)=height(nuvz)
- end do
-
- ! Determine highest levels that can be within PBL
- !************************************************
+
+
+ ! Switch on following lines to use doubled vertical resolution
+ !*************************************************************
+ ! do 7 kz=3,nz-1,2
+ ! height(kz)=0.5*(height(kz-1)+height(kz+1))
+ ! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
+ ! End doubled vertical resolution
+
+
+ ! Determine highest levels that can be within PBL
+ !************************************************
+
do kz=1,nz
- if (height(kz).gt.hmixmax) exit
+ if (height(kz).gt.hmixmax) then
+ nmixz=kz
+ goto 9
+ endif
end do
- !FIXME CHECK NEW SH: formerly was undefined if condition never fulfilled; now is nz+1
- nmixz=kz
+9 continue
+
+ ! Do not repeat initialization of the Cartesian z grid
+ !*****************************************************
+
+ init=.false.
- ! Do not repeat initialization of the Cartesian z grid
- !*****************************************************
- verttransform_inited=.true.
endif
- write(*,*) 'verttransform_gfs.f90 mostly parallel using ',omp_get_max_threads(),' threads'
-!$omp parallel default(none) private(ix,jy,kz,loopcnt,tvold,pold,tv,pint) &
-!$omp shared(n,nz,nxmin1,nymin1,nuvz,akz,bkz,ps,llev_arr, &
-!$omp tth,qvh,uvwzlev,rhoh_arr,pinmconv_arr,aknew,bknew) &
-!$omp num_threads(omp_get_max_threads())
-!$omp do
+ ! Loop over the whole grid
+ !*************************
+
do jy=0,nymin1
do ix=0,nxmin1
- !NCEP version: find first level above ground
- !FIXME Assume AKZ ist absteigend in i
- do loopcnt=1,nuvz-2
- llev_arr(ix,jy)=loopcnt
- if (akz(loopcnt) <= ps(ix,jy,1,n)) exit
+
+ ! 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
- end do
- end do
-!$omp end do
+ llev = llev+1
+ if (llev.gt.nuvz-2) llev = nuvz-2
+ ! if (llev.eq.nuvz-2) write(*,*) 'verttransform
+ ! +WARNING: LLEV eq NUZV-2'
+ ! NCEP version
-!$omp do
- do jy=0,nymin1
- do ix=0,nxmin1
- ! compute height of pressure levels above ground
- !***********************************************
- tvold=tth(ix,jy,llev_arr(ix,jy),n)*(1.+0.608*qvh(ix,jy,llev_arr(ix,jy),n))
- pold=akz(llev_arr(ix,jy))
- uvwzlev(ix,jy,llev_arr(ix,jy))=0.
- rhoh_arr(ix,jy,llev_arr(ix,jy))=pold/(r_air*tvold)
-
- do kz=llev_arr(ix,jy)+1,nuvz
+
+ ! compute height of pressure levels above ground
+ !***********************************************
+
+ tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n))
+ pold=akz(llev)
+ uvzlev(llev)=0.
+ wzlev(llev)=0.
+ uvwzlev(ix,jy,llev)=0.
+ rhoh(llev)=pold/(r_air*tvold)
+
+ do kz=llev+1,nuvz
pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
- rhoh_arr(ix,jy,kz)=pint/(r_air*tv)
+ rhoh(kz)=pint/(r_air*tv)
if (abs(tv-tvold).gt.0.2) then
- uvwzlev(ix,jy,kz) = uvwzlev(ix,jy,kz-1) + (r_air/ga)*log(pold/pint) * (tv-tvold)/log(tv/tvold)
+ uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
+ (tv-tvold)/log(tv/tvold)
else
- uvwzlev(ix,jy,kz) = uvwzlev(ix,jy,kz-1) + (r_air/ga)*log(pold/pint) * tv
+ uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
endif
+ wzlev(kz)=uvzlev(kz)
+ uvwzlev(ix,jy,kz)=uvzlev(kz)
tvold=tv
pold=pint
end do
- pinmconv_arr(ix,jy,llev_arr(ix,jy))=(uvwzlev(ix,jy,llev_arr(ix,jy)+1)-uvwzlev(ix,jy,llev_arr(ix,jy)))/ &
- ((aknew(llev_arr(ix,jy)+1)+bknew(llev_arr(ix,jy)+1)*ps(ix,jy,1,n))- &
- (aknew(llev_arr(ix,jy))+bknew(llev_arr(ix,jy))*ps(ix,jy,1,n)))
- do kz=llev_arr(ix,jy)+1,nz-1
- pinmconv_arr(ix,jy,kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
- (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
- end do
- pinmconv_arr(ix,jy,nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
- (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
- end do
- end do
-!$omp end do
-!$omp end parallel
+ ! Switch on following lines to use doubled vertical resolution
+ ! Switch off the three lines above.
+ !*************************************************************
+ !22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
+ ! do 23 kz=2,nwz
+ !23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
+ ! End doubled vertical resolution
-! write(*,*) 'verttransform_gfs.f90 mostly parallel 2/3 using ',omp_get_max_threads(),' threads'
+ ! pinmconv=(h2-h1)/(p2-p1)
-!$omp parallel default(none) private(ix,jy,iz,kz,dz,dz1,dz2) &
-!$omp shared(nymin1,nxmin1,nz,nuvz,nwz,n, &
-!$omp height,hind_above_arr,match_height_arr,llev_arr,uvwzlev, &
-!$omp uu,uuh,vv,vvh,tt,tth,qv,qvh,pv,pvh,rho,rhoh_arr,pplev,akz, &
-!$omp ww,wwh,pinmconv_arr,drhodz) &
-!$omp num_threads(omp_get_max_threads())
-!$omp do
- do jy=0,nymin1
- do ix=0,nxmin1
- !1 is always OK (zero - zero), 2 is thus the first number which can appear
- do iz=2,nz
- if(height(iz) > uvwzlev(ix,jy,nuvz)) exit
+ pinmconv(llev)=(uvwzlev(ix,jy,llev+1)-uvwzlev(ix,jy,llev))/ &
+ ((aknew(llev+1)+bknew(llev+1)*ps(ix,jy,1,n))- &
+ (aknew(llev)+bknew(llev)*ps(ix,jy,1,n)))
+ do kz=llev+1,nz-1
+ pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
+ ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
+ (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
end do
- hind_above_arr(ix,jy)=iz
- end do
- end do
-!$omp end do
+ pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
+ ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
+ (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
-!$omp do
- !matching level z in uvwzlev to each height level in height
- !default -1 so that anything going wrong will be found by array bounds checker
- do jy=0,nymin1
- do ix=0,nxmin1
- do iz=1,nz
- match_height_arr(ix,jy,iz)=-1
- end do
- end do
- end do
-!$omp end do
-!$omp do
- do jy=0,nymin1
- do ix=0,nxmin1
- do iz=2,min(nz,hind_above_arr(ix,jy))-1
- do kz=llev_arr(ix,jy)+1,nuvz
- !hind_above_arr such that condition is *always* fulfilled at some point
- !meaning kz-based uwvzlev goes from llev_arr height (zero) to nuvz-height
- !and iz-based height goes from 2-height to some hind_above_arr-height
- !so a matching kz can always be found
- !
- !WHEN OPTIMISING MIND THE FACT THAT SOME OF THESE ARRAYS MAY HAVE ZERO PADDING BEFORE/AFTER!
- !
- if ( height(iz) > uvwzlev(ix,jy,kz-1) .and. height(iz) <= uvwzlev(ix,jy,kz) ) exit
- end do
- if (kz==nuvz+1) then
- write(*,*) 'FATAL hindex search verttransform'
- stop
- end if
- match_height_arr(ix,jy,iz)=kz
- end do
- end do
- end do
-!$omp end do
-!$omp do
! Levels, where u,v,t and q are given
!************************************
- do iz=min(minval(hind_above_arr),nz),nz
- uu( 0:nxmin1,0:nymin1,iz,n)=uuh(0:nxmin1,0:nymin1,nuvz)
- vv( 0:nxmin1,0:nymin1,iz,n)=vvh(0:nxmin1,0:nymin1,nuvz)
- tt( 0:nxmin1,0:nymin1,iz,n)=tth(0:nxmin1,0:nymin1,nuvz,n)
- qv( 0:nxmin1,0:nymin1,iz,n)=qvh(0:nxmin1,0:nymin1,nuvz,n)
- pv( 0:nxmin1,0:nymin1,iz,n)=pvh(0:nxmin1,0:nymin1,nuvz)
- rho( 0:nxmin1,0:nymin1,iz,n)=rhoh_arr(0:nxmin1,0:nymin1,nuvz)
- pplev(0:nxmin1,0:nymin1,iz,n)=akz(nuvz)
- end do
-!$omp end do
-!$omp do
- do jy=0,nymin1
- do ix=0,nxmin1
- uu(ix,jy,1,n)=uuh(ix,jy,llev_arr(ix,jy))
- vv(ix,jy,1,n)=vvh(ix,jy,llev_arr(ix,jy))
- tt(ix,jy,1,n)=tth(ix,jy,llev_arr(ix,jy),n)
- qv(ix,jy,1,n)=qvh(ix,jy,llev_arr(ix,jy),n)
- pv(ix,jy,1,n)=pvh(ix,jy,llev_arr(ix,jy))
- rho(ix,jy,1,n)=rhoh_arr(ix,jy,llev_arr(ix,jy))
- pplev(ix,jy,1,n)=akz(llev_arr(ix,jy))
- end do
- end do
-!$omp end do
-!$omp do
- do jy=0,nymin1
- do ix=0,nxmin1
- do iz=2,min(nz,hind_above_arr(ix,jy))-1
- kz=match_height_arr(ix,jy,iz)
- dz1=height(iz)-uvwzlev(ix,jy,kz-1)
- dz2=uvwzlev(ix,jy,kz)-height(iz)
- dz=dz1+dz2
- uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
- vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
- tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2+tth(ix,jy,kz,n)*dz1)/dz
- qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2+qvh(ix,jy,kz,n)*dz1)/dz
- pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
- rho(ix,jy,iz,n)=(rhoh_arr(ix,jy,kz-1)*dz2+rhoh_arr(ix,jy,kz)*dz1)/dz
- pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
+
+ uu(ix,jy,1,n)=uuh(ix,jy,llev)
+ vv(ix,jy,1,n)=vvh(ix,jy,llev)
+ tt(ix,jy,1,n)=tth(ix,jy,llev,n)
+ qv(ix,jy,1,n)=qvh(ix,jy,llev,n)
+ pv(ix,jy,1,n)=pvh(ix,jy,llev)
+ rho(ix,jy,1,n)=rhoh(llev)
+ pplev(ix,jy,1,n)=akz(llev)
+ uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
+ vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
+ tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
+ qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
+ pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
+ rho(ix,jy,nz,n)=rhoh(nuvz)
+ pplev(ix,jy,nz,n)=akz(nuvz)
+ kmin=llev+1
+ do iz=2,nz-1
+ do kz=kmin,nuvz
+ if(height(iz).gt.uvzlev(nuvz)) then
+ uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
+ vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
+ tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
+ qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
+ pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
+ rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
+ pplev(ix,jy,iz,n)=pplev(ix,jy,nz,n)
+ goto 30
+ endif
+ if ((height(iz).gt.uvzlev(kz-1)).and. &
+ (height(iz).le.uvzlev(kz))) then
+ dz1=height(iz)-uvzlev(kz-1)
+ dz2=uvzlev(kz)-height(iz)
+ dz=dz1+dz2
+ uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
+ vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
+ tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
+ +tth(ix,jy,kz,n)*dz1)/dz
+ qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
+ +qvh(ix,jy,kz,n)*dz1)/dz
+ pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
+ rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
+ pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
+ endif
+ end do
+30 continue
end do
- end do
- end do
-!$omp end do
-
+
+
! Levels, where w is given
!*************************
-!$omp workshare
- ww(0:nxmin1,0:nymin1,nz,n)=wwh(0:nxmin1,0:nymin1,nwz)*pinmconv_arr(0:nxmin1,0:nymin1,nz)
-!$omp end workshare
-!$omp do
- do iz=minval(hind_above_arr),nz-1
- !SH: FIXME shouldn't it be nwz here in the last thingy??? -> unclear
- ww(0:nxmin1,0:nymin1,iz,n)=wwh(0:nxmin1,0:nymin1,nwz)*pinmconv_arr(0:nxmin1,0:nymin1,nz)
- end do
-!$omp end do
-!$omp do
- do jy=0,nymin1
- do ix=0,nxmin1
- ww(ix,jy,1,n)=wwh(ix,jy,llev_arr(ix,jy))*pinmconv_arr(ix,jy,llev_arr(ix,jy))
- end do
- end do
-!$omp end do
-!$omp do
- do jy=0,nymin1
- do ix=0,nxmin1
- do iz=2,min(nz,hind_above_arr(ix,jy))-1
- kz=match_height_arr(ix,jy,iz)
- dz1=height(iz)-uvwzlev(ix,jy,kz-1)
- dz2=uvwzlev(ix,jy,kz)-height(iz)
- ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv_arr(ix,jy,kz-1)*dz2 + wwh(ix,jy,kz)*pinmconv_arr(ix,jy,kz)*dz1)/(dz1+dz2)
+
+ ww(ix,jy,1,n)=wwh(ix,jy,llev)*pinmconv(llev)
+ ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
+ kmin=llev+1
+ do iz=2,nz
+ do kz=kmin,nwz
+ if ((height(iz).gt.wzlev(kz-1)).and. &
+ (height(iz).le.wzlev(kz))) then
+ dz1=height(iz)-wzlev(kz-1)
+ dz2=wzlev(kz)-height(iz)
+ dz=dz1+dz2
+ ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
+ +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
+
+ endif
+ end do
end do
- end do
- end do
-!$omp end do
+
! Compute density gradients at intermediate levels
!*************************************************
-!$omp workshare
- drhodz(0:nxmin1,0:nymin1,1,n)=(rho(0:nxmin1,0:nymin1,2,n)-rho(0:nxmin1,0:nymin1,1,n))/(height(2)-height(1))
-!$omp end workshare
-!$omp do
- do kz=2,nz-1
- drhodz(0:nxmin1,0:nymin1,kz,n)=(rho(0:nxmin1,0:nymin1,kz+1,n)-rho(0:nxmin1,0:nymin1,kz-1,n))/(height(kz+1)-height(kz-1))
- end do
-!$omp end do
-!$omp workshare
- drhodz(0:nxmin1,0:nymin1,nz,n)=drhodz(0:nxmin1,0:nymin1,nz-1,n)
-!$omp end workshare
-!$omp end parallel
+ drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
+ (height(2)-height(1))
+ do kz=2,nz-1
+ drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
+ (height(kz+1)-height(kz-1))
+ end do
+ drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
+ end do
+ end do
-! write(*,*) 'verttransform_gfs.f90 mostly parallel 3/3 using ',omp_get_max_threads(),' threads'
-!$omp parallel default(none) private(ix,jy,iz,kz,dz,dz1,dz2,ix1,jy1, &
-!$omp ixp,jyp,dzdx1,dzdx2,dzdx,dzdy1,dzdy2,dzdy,ui,vi,tmp_real) &
-!$omp shared(n,nz,nx,ny,uu,vv,dxconst,dyconst,dy,ylat0,ww, &
-!$omp hind_above_arr,match_height_arr,height,uvwzlev) &
-!$omp num_threads(omp_get_max_threads())
-!$omp do
!****************************************************************
! Compute slope of eta levels in windward direction and resulting
! vertical wind correction
!****************************************************************
+
do jy=1,ny-2
- tmp_real=cos((real(jy)*dy+ylat0)*pi180)
do ix=1,nx-2
- do iz=2,min(nz,hind_above_arr(ix,jy))-1
- kz=match_height_arr(ix,jy,iz)
- dz1=height(iz)-uvwzlev(ix,jy,kz-1)
- dz2=uvwzlev(ix,jy,kz)-height(iz)
- dz=dz1+dz2
+ ! 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-2) llev = nuvz-2
+ ! if (llev.eq.nuvz-2) write(*,*) 'verttransform
+ ! +WARNING: LLEV eq NUZV-2'
+ ! NCEP version
- ix1=ix-1
+ kmin=llev+1
+ do iz=2,nz-1
+
+ ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
+ vi=vv(ix,jy,iz,n)*dyconst
+
+ do kz=kmin,nz
+ if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
+ (height(iz).le.uvwzlev(ix,jy,kz))) then
+ dz1=height(iz)-uvwzlev(ix,jy,kz-1)
+ dz2=uvwzlev(ix,jy,kz)-height(iz)
+ dz=dz1+dz2
+ kl=kz-1
+ klp=kz
+ goto 47
+ endif
+ end do
+
+47 ix1=ix-1
jy1=jy-1
ixp=ix+1
jyp=jy+1
- dzdx1=(uvwzlev(ixp,jy,kz-1)-uvwzlev(ix1,jy,kz-1))/2.
- dzdx2=(uvwzlev(ixp,jy,kz)-uvwzlev(ix1,jy,kz))/2.
+ dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
+ dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
- dzdy1=(uvwzlev(ix,jyp,kz-1)-uvwzlev(ix,jy1,kz-1))/2.
- dzdy2=(uvwzlev(ix,jyp,kz)-uvwzlev(ix,jy1,kz))/2.
+ dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
+ dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
- ui=uu(ix,jy,iz,n)*dxconst/tmp_real
- vi=vv(ix,jy,iz,n)*dyconst
ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
+
end do
+
end do
end do
-!$omp end do
-!$omp end parallel
-! write(*,*) 'verttransform_gfs.f90 mostly parallel 3a/3 using ',omp_get_max_threads(),' threads'
! If north pole is in the domain, calculate wind velocities in polar
! stereographic coordinates
!*******************************************************************
if (nglobal) then
-!$omp parallel default(none) private(jy,iz,ylat,xlon) &
-!$omp shared(n,nz,switchnorthg,nxmin1,nymin1,xlon0,ylat0,dy,dx, &
-!$omp northpolemap,uu,vv,uupol,vvpol) &
-!$omp num_threads(omp_get_max_threads())
-!$omp do
do jy=int(switchnorthg)-2,nymin1
ylat=ylat0+real(jy)*dy
do ix=0,nxmin1
xlon=xlon0+real(ix)*dx
do iz=1,nz
call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
- vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
+ vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+ vvpol(ix,jy,iz,n))
end do
end do
end do
-!$omp end do
-!$omp end parallel
+
do iz=1,nz
- ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+
+ ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
xlon=xlon0+real(nx/2-1)*dx
xlonr=xlon*pi/180.
- ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2)
+ ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
+ vv(nx/2-1,nymin1,iz,n)**2)
if(vv(nx/2-1,nymin1,iz,n).lt.0.) then
- ddpol = atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
+ ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
+ vv(nx/2-1,nymin1,iz,n))-xlonr
elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
- ddpol = pi + atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
+ ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
+ vv(nx/2-1,nymin1,iz,n))-xlonr
else
- ddpol = pi/2 - xlonr
+ ddpol=pi/2-xlonr
endif
if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
- ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
- uuaux=-ffpol*sin(pi+ddpol)
- vvaux=-ffpol*cos(pi+ddpol)
- call cc2gll(northpolemap,90.0,180.0,uuaux,vvaux,uupolaux,vvpolaux)
+ ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+ xlon=180.0
+ xlonr=xlon*pi/180.
+ ylat=90.0
+ uuaux=-ffpol*sin(xlonr+ddpol)
+ vvaux=-ffpol*cos(xlonr+ddpol)
+ call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+ vvpolaux)
- uupol(0:nxmin1,nymin1,iz,n)=uupolaux
- vvpol(0:nxmin1,nymin1,iz,n)=vvpolaux
+ jy=nymin1
+ do ix=0,nxmin1
+ uupol(ix,jy,iz,n)=uupolaux
+ vvpol(ix,jy,iz,n)=vvpolaux
+ end do
end do
- ! Fix: Set W at pole to the zonally averaged W of the next equator-
- ! ward parallel of latitude
- do iz=1,nz
- ww(0:nxmin1,nymin1,iz,n)=sum(ww(0:nxmin1,ny-2,iz,n))/real(nx)
- end do
+
+ ! Fix: Set W at pole to the zonally averaged W of the next equator-
+ ! ward parallel of latitude
+
+ do iz=1,nz
+ wdummy=0.
+ jy=ny-2
+ do ix=0,nxmin1
+ wdummy=wdummy+ww(ix,jy,iz,n)
+ end do
+ wdummy=wdummy/real(nx)
+ jy=nymin1
+ do ix=0,nxmin1
+ ww(ix,jy,iz,n)=wdummy
+ end do
+ end do
+
endif
@@ -469,56 +480,73 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!*******************************************************************
if (sglobal) then
-!$omp parallel default(none) private(jy,iz,ylat,xlon) &
-!$omp shared(n,nz,switchsouthg,nymin1,nxmin1,xlon0,ylat0,dy,dx, &
-!$omp southpolemap,uu,vv,uupol,vvpol) &
-!$omp num_threads(omp_get_max_threads())
-!$omp do
do jy=0,int(switchsouthg)+3
ylat=ylat0+real(jy)*dy
do ix=0,nxmin1
xlon=xlon0+real(ix)*dx
do iz=1,nz
call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
- vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
+ vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+ vvpol(ix,jy,iz,n))
end do
end do
end do
-!$omp end do
-!$omp end parallel
do iz=1,nz
- ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+
+ ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
xlon=xlon0+real(nx/2-1)*dx
xlonr=xlon*pi/180.
- ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2)
+ ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
+ vv(nx/2-1,0,iz,n)**2)
if(vv(nx/2-1,0,iz,n).lt.0.) then
- ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
+ ddpol=atan(uu(nx/2-1,0,iz,n)/ &
+ vv(nx/2-1,0,iz,n))+xlonr
elseif (vv(nx/2-1,0,iz,n).gt.0.) then
- ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))-xlonr
+ ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
+ vv(nx/2-1,0,iz,n))-xlonr
else
ddpol=pi/2-xlonr
endif
if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
- ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
- uuaux=+ffpol*sin(pi-ddpol)
- vvaux=-ffpol*cos(pi-ddpol)
- call cc2gll(northpolemap,-90.0,180.0,uuaux,vvaux,uupolaux,vvpolaux)
+ ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+ xlon=180.0
+ xlonr=xlon*pi/180.
+ ylat=-90.0
+ uuaux=+ffpol*sin(xlonr-ddpol)
+ vvaux=-ffpol*cos(xlonr-ddpol)
+ call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+ vvpolaux)
- uupol(0:nxmin1,0,iz,n)=uupolaux
- vvpol(0:nxmin1,0,iz,n)=vvpolaux
+ jy=0
+ do ix=0,nxmin1
+ uupol(ix,jy,iz,n)=uupolaux
+ vvpol(ix,jy,iz,n)=vvpolaux
+ end do
end do
- ! Fix: Set W at pole to the zonally averaged W of the next equator-
- ! ward parallel of latitude
+
+ ! Fix: Set W at pole to the zonally averaged W of the next equator-
+ ! ward parallel of latitude
+
do iz=1,nz
- ww(0:nxmin1,0,iz,n)=sum(ww(0:nxmin1,1,iz,n))/real(nx)
+ wdummy=0.
+ jy=1
+ do ix=0,nxmin1
+ wdummy=wdummy+ww(ix,jy,iz,n)
+ end do
+ wdummy=wdummy/real(nx)
+ jy=0
+ do ix=0,nxmin1
+ ww(ix,jy,iz,n)=wdummy
+ end do
end do
endif
+ ! write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
! create a cloud and rainout/washout field, clouds occur where rh>80%
! total cloudheight is stored at level 0
do jy=0,nymin1
@@ -533,32 +561,30 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
clouds(ix,jy,kz,n)=0
if (rh.gt.0.8) then ! in cloud
- if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
- rain_cloud_above=1
- cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
- height(kz)-height(kz-1)
- if (lsp.ge.convp) then
- clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
- else
- clouds(ix,jy,kz,n)=2 ! convp dominated rainout
- endif
- else ! no precipitation
- clouds(ix,jy,kz,n)=1 ! cloud
- endif
+ if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+ rain_cloud_above=1
+ cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
+ height(kz)-height(kz-1)
+ if (lsp.ge.convp) then
+ clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+ else
+ clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+ endif
+ else ! no precipitation
+ clouds(ix,jy,kz,n)=1 ! cloud
+ endif
else ! no cloud
- if (rain_cloud_above.eq.1) then ! scavenging
- if (lsp.ge.convp) then
- clouds(ix,jy,kz,n)=5 ! lsp dominated washout
- else
- clouds(ix,jy,kz,n)=4 ! convp dominated washout
- endif
- endif
+ if (rain_cloud_above.eq.1) then ! scavenging
+ if (lsp.ge.convp) then
+ clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+ else
+ clouds(ix,jy,kz,n)=4 ! convp dominated washout
+ endif
+ endif
endif
end do
end do
end do
-
- write(*,*) 'exit verttransform...'
-
-end subroutine verttransform
+
+end subroutine verttransform
diff --git a/flexpart_code/xmass_mod.f90 b/flexpart_code/xmass_mod.f90
index d6fa7db..d97c63b 100644
--- a/flexpart_code/xmass_mod.f90
+++ b/flexpart_code/xmass_mod.f90
@@ -21,9 +21,8 @@
module xmass_mod
- use par_mod, only : dp
implicit none
- real(kind=dp), allocatable, dimension (:) :: xmasssave
+ real,allocatable, dimension (:) :: xmasssave
end module xmass_mod
diff --git a/flexpart_code/xmass_mod.mod b/flexpart_code/xmass_mod.mod
deleted file mode 100644
index 29caa21..0000000
Binary files a/flexpart_code/xmass_mod.mod and /dev/null differ