Skip to content

GitLab

Commit 1a8fbee0 authored by pesei's avatar pesei
Browse files

vertransform: fix #140 (ref position for grid), verttransform_nest: fix bug... 

vertransform: fix #140 (ref position for grid), verttransform_nest: fix bug with dimensions, both: smaller changes, com_mod: correct comment to nmixz
parent 77778f8d
Show whitespace changes
Inline Side-by-side
src/changelog.txt
View file @ 1a8fbee0
This is a changelog for the dev branch of FLEXPART (see flexpart.eu)
Created by Petra Seibert, 8 June 2018 to document and explain in detail why I changed what.
makefile, 2018-06-08
========
* makefile, 2018-06-08
========
1) Add a GPL3+ License statement, author statement, version information. Prudent to do that, not harmful.
......@@ -23,15 +23,42 @@ Notes and questions:
2. Would be nice to have "as fast as possible" option (thought it might be machine-dependent ...)
3. What is the purpose of -Warray-bounds ? Is it relevant for Fortran? I think if there is a compile-time array-bound violation this would be a hard error?
FLEXPART.f90
========
* FLEXPART.f90
========
Update version string!
readcommand.f90
===========
* readcommand.f90
===========
Correct misleading error message (replace "open" by "write to")
Improve error messages: combine pathname+filename, remove unaligned closing ###
Let STOP say in which subroutine we are stopping
Some minor changes
* verttransform_ecmwf.f90 PS 2018-06-13
===================
1) Remove some of commented out testing stuff
2) Fix indenting in the init if block
3) Code layout improvement
4) Change name of loops to represent the index
5) Fix ticket:140 Introduce a new way of establishing the reference position
for the vertical grid. Also correct a minor bug in the calculation
of height (array assignment where it was not intended)
6) Add back the SAVE attribute to INIT, just to be sure
7) Bring some more code under the IF (INIT) block
8) make some annotations
* com_mod.f90 PS 2018-06-13
=======
correct the comment for nmixz
* verttransform_nest.f90 PS 2018-06-17
===================
1) insert proper boundaries for implied loops in array expressions
(fixes bug annoted by ESO in 2016)
2) Code layout improvement
3) Change name of loops to represent the index
4) make some annotations
src/com_mod.f90
View file @ 1a8fbee0
......@@ -301,7 +301,7 @@ module com_mod
! nuvz,nwz vertical dimension of original ECMWF data
! nxfield same as nx for limited area fields,
! but for global fields nx=nxfield+1
! nmixz number of levels up to maximum PBL height (3500 m)
! nmixz number of levels up to maximum PBL height (hmixmax)
! nuvz is used for u,v components
! nwz is used for w components (staggered grid)
......
src/verttransform_ecmwf.f90
View file @ 1a8fbee0
!**********************************************************************
! **********************************************************************
! 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 *
......@@ -17,7 +17,7 @@
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
! **********************************************************************
subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
! i i i i i
......@@ -35,6 +35,9 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
! 12 August 1996 *
! Update: 16 January 1998 *
! *
! *
!*****************************************************************************
! CHANGES *
! Major update: 17 February 1999 *
! by G. Wotawa *
! *
......@@ -42,18 +45,30 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
! - Slope correction for vertical velocity: Modification of calculation *
! procedure *
! *
!*****************************************************************************
! Changes, Bernd C. Krueger, Feb. 2001:
! Bernd C. Krueger, Feb. 2001:
! Variables tth and qvh (on eta coordinates) from common block
!
! Sabine Eckhardt, March 2007
! *
! Sabine Eckhardt, March 2007:
! added the variable cloud for use with scavenging - descr. in com_mod
!
! *
! Who? When? *
! Unified ECMWF and GFS builds
! Marian Harustak, 12.5.2017
! - Renamed from verttransform to verttransform_ecmwf
!*****************************************************************************
! Date: 2017-05-30 modification of a bug in ew. Don Morton (CTBTO project) *
! *
! Don Morton, 2017-05-30:
! modification of a bug in ew. Don Morton (CTBTO project) *
! *
! undocumented modifications by NILU for v10 *
! *
! Petra Seibert, 2018-06-13: *
! - fix bug of ticket:140 (find reference position for vertical grid) *
! - put back SAVE attribute for INIT, just to be safe *
! - minor changes, most of them just cosmetics *
! for details see changelog.txt *
! *
! ****************************************************************************
!*****************************************************************************
! *
! Variables: *
......@@ -72,7 +87,6 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
use par_mod
use com_mod
use cmapf_mod, only: cc2gll
! use mpi_mod
implicit none
......@@ -81,12 +95,20 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
real,dimension(0:nxmax-1,0:nymax-1) :: tvold,pold,pint,tv
!> for reference profile (PS)
real :: tvoldref, poldref, pintref, psmean, psstd
integer :: ixyref(2)
integer, parameter :: ioldref = 1 ! PS 2018-06-13: set to 0 if you
!! want old method of searching reference location for the vertical grid
!! 1 for new method (options for other methods 2,... in the future)
real,dimension(0:nymax-1) :: cosf
real,dimension(0:nxmax-1,0:nymax-1) :: tvold,pold,pint,tv
!! automatic arrays introduced in v10 by ?? to achieve better loop order (PS)
integer,dimension(0:nxmax-1,0:nymax-1) :: rain_cloud_above,idx
integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixm,jym,kz_inv
integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixref,jyref,kz_inv
real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
real :: ew,dz1,dz2,dz
real :: xlon,ylat,xlonr,dzdx,dzdy
......@@ -95,9 +117,7 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
real,parameter :: const=r_air/ga
real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
logical :: init = .true.
logical :: init_w = .false.
logical :: init_r = .false.
logical, save :: init = .true. ! PS 2018-06-13: add back save attribute
!ZHG SEP 2014 tests
......@@ -119,79 +139,84 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
! If verttransform is called the first time, initialize heights of the *
! z levels in meter. The heights are the heights of model levels, where *
! u,v,T and qv are given, and of the interfaces, where w is given. So, *
! the vertical resolution in the z system is doubled. As reference point,*
! the vertical rESOlution in the z system is doubled. As reference point,*
! the lower left corner of the grid is used. *
! Unlike in the eta system, no difference between heights for u,v and *
! heights for w exists. *
!*************************************************************************
!eso measure CPU time
!ESO measure CPU time
! call mpif_mtime('verttransform',0)
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
! *****************************************************************************
if (init_r) then
open(333,file='heights.txt', &
form='formatted')
do kz=1,nuvz
read(333,*) height(kz)
end do
close(333)
write(*,*) 'height read'
else
! 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 (ioldref .eq. 0) then ! old reference grid point
do jy=0,nymin1
do ix=0,nxmin1
if (ps(ix,jy,1,n).gt.100000.) then
ixm=ix
jym=jy
if (ps(ix,jy,1,n).gt.1000.e2) then
ixref=ix
jyref=jy
goto 3
endif
end do
end do
3 continue
print*,'oldheights at' ,ixref,jyref,ps(ixref,jyref,1,n)
else ! new reference grid point
! PS: the old version fails if the pressure is <=1000 hPa in the whole
! domain. Let us find a good replacement, not just a quick fix.
! Search point near to mean pressure after excluding mountains
psmean = sum( ps(:,:,1,n) ) / (nx*ny)
print*,'height: fg psmean',psmean
psstd = sqrt(sum( (ps(:,:,1,n) - psmean)**2 ) / (nx*ny))
!> new psmean using only points within $\plusminus\sigma$
! psmean = sum( ps(:,:,1,n), abs(ps(:,:,1,n)-psmean) < psstd ) / &
! count(abs(ps(:,:,1,n)-psmean) < psstd)
!> new psmean using only points with $p\gt \overbar{p}+\sigma_p$
!> (reject mountains, accept valleys)
psmean = sum( ps(:,:,1,n), ps(:,:,1,n) > psmean - psstd ) / &
count(ps(:,:,1,n) > psmean - psstd)
print*,'height: std, new psmean',psstd,psmean
ixyref = minloc( abs( ps(:,:,1,n) - psmean ) )
ixref = ixyref(1)
jyref = ixyref(2)
print*,'newheights at' ,ixref,jyref,ps(ixref,jyref,1,n)
endif
tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
ps(ixm,jym,1,n))
pold(ixm,jym)=ps(ixm,jym,1,n)
tvoldref=tt2(ixref,jyref,1,n)* &
( 1. + 0.378*ew(td2(ixref,jyref,1,n) ) / ps(ixref,jyref,1,n))
poldref=ps(ixref,jyref,1,n)
height(1)=0.
do kz=2,nuvz
pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
if (abs(tv(ixm,jym)-tvold(ixm,jym)).gt.0.2) then
height(kz)= &
height(kz-1)+const*log(pold(ixm,jym)/pint(ixm,jym))* &
(tv(ixm,jym)-tvold(ixm,jym))/log(tv(ixm,jym)/tvold(ixm,jym))
pintref = akz(kz)+bkz(kz)*ps(ixref,jyref,1,n)
tv = tth(ixref,jyref,kz,n)*(1.+0.608*qvh(ixref,jyref,kz,n))
if (abs(tv(ixref,jyref)-tvoldref) .gt. 0.2) then
height(kz) = height(kz-1) + &
const*log( poldref/pintref ) * &
( tv(ixref,jyref) - tvoldref ) / log( tv(ixref,jyref) / tvoldref )
else
height(kz)=height(kz-1)+ &
const*log(pold(ixm,jym)/pint(ixm,jym))*tv(ixm,jym)
height(kz) = height(kz-1) + &
const*log( poldref/pintref ) * tv(ixref,jyref)
endif
tvold(ixm,jym)=tv(ixm,jym)
pold(ixm,jym)=pint(ixm,jym)
end do
tvoldref = tv(ixref,jyref)
poldref = pintref
print*,'height=',kz,height(kz),tvoldref,poldref
if (init_w) then
open(333,file='heights.txt', &
form='formatted')
do kz=1,nuvz
write(333,*) height(kz)
end do
close(333)
endif
endif ! init
! Determine highest levels that can be within PBL
!************************************************
......@@ -206,24 +231,15 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
! Do not repeat initialization of the Cartesian z grid
!*****************************************************
init=.false.
! dbg_height = height
endif
endif ! init block (vertical grid construction)
! Loop over the whole grid
!*************************
do jy=0,nymin1
do ix=0,nxmin1
tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
ps(ix,jy,1,n))
enddo
enddo
tvold(:,:)=tt2(:,:,1,n) * (1.+0.378*ew( td2(:,:,1,n) ) / ps(:,:,1,n))
pold=ps(:,:,1,n)
uvzlev(:,:,1)=0.
wzlev(:,:,1)=0.
......@@ -234,19 +250,19 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
!******************************
do kz=2,nuvz
pint=akz(kz)+bkz(kz)*ps(:,:,1,n)
tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
rhoh(:,:,kz)=pint(:,:)/(r_air*tv)
pint(:,:)=akz(kz)+bkz(kz)*ps(:,:,1,n)
tv(:,:)=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
rhoh(:,:,kz)=pint(:,:)/(r_air*tv(:,:))
where (abs(tv-tvold).gt.0.2)
uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)* &
(tv-tvold)/log(tv/tvold)
where (abs(tv(:,:)-tvold(:,:)).gt.0.2)
uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold(:,:)/pint(:,:))* &
(tv(:,:)-tvold(:,:))/log(tv(:,:)/tvold(:,:))
elsewhere
uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)*tv
uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold(:,:)/pint(:,:))*tv(:,:)
endwhere
tvold=tv
pold=pint
tvold(:,:)=tv(:,:)
pold(:,:)=pint(:,:)
end do
......@@ -270,20 +286,19 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
((aknew(nz)+bknew(nz)*ps(:,:,1,n))- &
(aknew(nz-1)+bknew(nz-1)*ps(:,:,1,n)))
! Levels, where u,v,t and q are given
!************************************
! Levels where u,v,t and q are given
!***********************************
uu(:,:,1,n)=uuh(:,:,1)
vv(:,:,1,n)=vvh(:,:,1)
tt(:,:,1,n)=tth(:,:,1,n)
qv(:,:,1,n)=qvh(:,:,1,n)
!hg adding the cloud water
!HG adding the cloud water
if (readclouds) then
clwc(:,:,1,n)=clwch(:,:,1,n)
if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)
end if
!hg
!HG
pv(:,:,1,n)=pvh(:,:,1)
rho(:,:,1,n)=rhoh(:,:,1)
......@@ -291,48 +306,52 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
vv(:,:,nz,n)=vvh(:,:,nuvz)
tt(:,:,nz,n)=tth(:,:,nuvz,n)
qv(:,:,nz,n)=qvh(:,:,nuvz,n)
!hg adding the cloud water
!HG adding the cloud water
if (readclouds) then
clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
if (.not. sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
end if
!hg
!HG
pv(:,:,nz,n)=pvh(:,:,nuvz)
rho(:,:,nz,n)=rhoh(:,:,nuvz)
kmin=2
idx=kmin
do iz=2,nz-1
iz_loop: do iz=2,nz-1
do jy=0,nymin1
do ix=0,nxmin1
if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
if (height(iz).gt.uvzlev(ix,jy,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)
!hg adding the cloud water
!HG adding the cloud water
if (readclouds) then
clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
if (.not. sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
end if
!hg
!HG
pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
else
innuvz: do kz=idx(ix,jy),nuvz
if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
(height(iz).le.uvzlev(ix,jy,kz))) then
kz_loop: do kz=idx(ix,jy),nuvz
if (idx(ix,jy) .le. kz .and. height(iz).gt.uvzlev(ix,jy,kz-1) &
.and. height(iz).le.uvzlev(ix,jy,kz)) then
idx(ix,jy)=kz
exit innuvz
exit kz_loop
endif
enddo innuvz
enddo kz_loop
endif
enddo
enddo
do jy=0,nymin1
do ix=0,nxmin1
if(height(iz).le.uvzlev(ix,jy,nuvz)) then
if (height(iz).le.uvzlev(ix,jy,nuvz)) then
kz=idx(ix,jy)
dz1=height(iz)-uvzlev(ix,jy,kz-1)
dz2=uvzlev(ix,jy,kz)-height(iz)
......@@ -343,38 +362,39 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
+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
!hg adding the cloud water
!HG adding the cloud water
if (readclouds) then
clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
if (.not.sumclouds) &
&ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
if (.not. sumclouds) ciwc(ix,jy,iz,n)= &
(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
end if
!hg
!HG
pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
endif
enddo
enddo
enddo
enddo iz_loop
! Levels, where w is given
! Levels where w is given
!*************************
ww(:,:,1,n)=wwh(:,:,1)*pinmconv(:,:,1)
ww(:,:,nz,n)=wwh(:,:,nwz)*pinmconv(:,:,nz)
kmin=2
idx=kmin
do iz=2,nz
iz_loop2: do iz=2,nz
do jy=0,nymin1
do ix=0,nxmin1
inn: do kz=idx(ix,jy),nwz
if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. &
height(iz).le.wzlev(ix,jy,kz)) then
kz_loop2: do kz=idx(ix,jy),nwz
if (idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1) &
.and. height(iz).le.wzlev(ix,jy,kz)) then
idx(ix,jy)=kz
exit inn
exit kz_loop2
endif
enddo inn
enddo kz_loop2
enddo
enddo
do jy=0,nymin1
......@@ -384,25 +404,21 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
dz2=wzlev(ix,jy,kz)-height(iz)
dz=dz1+dz2
ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(ix,jy,kz-1)*dz2 &
+wwh(ix,jy,kz)*pinmconv(ix,jy,kz)*dz1)/dz
enddo
+ wwh(ix,jy,kz) *pinmconv(ix,jy,kz)*dz1)/dz
enddo
enddo
enddo iz_loop2
! Compute density gradients at intermediate levels
!*************************************************
drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/ &
(height(2)-height(1))
drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/(height(2)-height(1))
do kz=2,nz-1
drhodz(:,:,kz,n)=(rho(:,:,kz+1,n)-rho(:,:,kz-1,n))/ &
(height(kz+1)-height(kz-1))
end do
drhodz(:,:,nz,n)=drhodz(:,:,nz-1,n)
! end do
! end do
!****************************************************************
! Compute slope of eta levels in windward direction and resulting
......@@ -410,22 +426,22 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
!****************************************************************
do jy=1,ny-2
cosf(jy)=1./cos((real(jy)*dy+ylat0)*pi180)
cosf(jy) = 1. / cos( ( real(jy)*dy + ylat0 )*pi180 )
enddo
kmin=2
idx=kmin
do iz=2,nz-1
iz_loop3: do iz=2,nz-1
do jy=1,ny-2
do ix=1,nx-2
inneta: do kz=idx(ix,jy),nz
if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
(height(iz).le.uvzlev(ix,jy,kz))) then
kz_loop3: do kz=idx(ix,jy),nz
if (idx(ix,jy) .le. kz .and. height(iz).gt.uvzlev(ix,jy,kz-1) &
.and. height(iz).le.uvzlev(ix,jy,kz)) then
idx(ix,jy)=kz
exit inneta
exit kz_loop3
endif
enddo inneta
enddo kz_loop3
enddo
enddo
......@@ -441,59 +457,57 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
jyp=jy+1
dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
dzdx2=(uvzlev(ixp,jy,kz) -uvzlev(ix1,jy,kz) )/2.
dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
dzdy2=(uvzlev(ix,jyp,kz) -uvzlev(ix,jy1,kz) )/2.
dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*uu(ix,jy,iz,n)*dxconst*cosf(jy)+dzdy*vv(ix,jy,iz,n)*dyconst)
ww(ix,jy,iz,n)=ww(ix,jy,iz,n) +( dzdx*uu(ix,jy,iz,n)*dxconst*cosf(jy) &
+ dzdy*vv(ix,jy,iz,n)*dyconst)
end do
enddo
enddo
end do
end do
enddo iz_loop3
! If north pole is in the domain, calculate wind velocities in polar
! stereographic coordinates
!*******************************************************************
if (nglobal) then
do iz=1,nz