Commit 4d456391 authored by Espen Sollum's avatar Espen Sollum
Browse files

Merged in optimized verttransform.f90 from author Leopold Haimberger

parent adf46aef
......@@ -20,7 +20,7 @@
!**********************************************************************
subroutine verttransform(n,uuh,vvh,wwh,pvh)
! i i i i i
! i i i i i
!*****************************************************************************
! *
! This subroutine transforms temperature, dew point temperature and *
......@@ -66,27 +66,32 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
use par_mod
use com_mod
use cmapf_mod, only: cc2gll
! eso TODO:
! only used for timing of CPU measurement. Remove this (and calls to mpif_mtime below)
! as this routine should not be dependent on MPI
! use mpi_mod
! :TODO
implicit none
integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
integer :: rain_cloud_above,kz_inv
! integer :: icloudtop !PS
real :: f_qvsat,pressure
! real :: rh,lsp,convp
real :: rh,lsp,convp,prec,rhmin,precmin
real :: rhoh(nuvzmax),pinmconv(nzmax)
real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
integer :: rain_cloud_above(0:nxmax-1,0:nymax-1),kz_inv,idx(0:nxmax-1,0:nymax-1)
real :: f_qvsat,pressure,cosf(0:nymax-1)
real :: rh,lsp,convp,tim,tim2,rhmin,precmin,prec
real :: uvzlev(0:nxmax-1,0:nymax-1,nuvzmax),rhoh(0:nxmax-1,0:nymax-1,nuvzmax)
real :: pinmconv(0:nxmax-1,0:nymax-1,nzmax)
real :: ew,pint(0:nxmax-1,0:nymax-1),tv(0:nxmax-1,0:nymax-1)
real :: tvold(0:nxmax-1,0:nymax-1),pold(0:nxmax-1,0:nymax-1),dz1,dz2,dz,ui,vi
real :: xlon,ylat,xlonr,dzdx,dzdy
real :: dzdx1,dzdx2,dzdy1,dzdy2,cosf
real :: dzdx1,dzdx2,dzdy1,dzdy2
real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
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 :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
logical lconvectprec
real :: wzlev(0:nxmax-1,0:nymax-1,nwzmax)
real,parameter :: const=r_air/ga
! integer:: ihr, imin, isec, i100th,ihr2, imin2, isec2, i100th2
parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
logical :: init = .true.
......@@ -102,14 +107,22 @@ subroutine verttransform(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. *
! 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 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
! 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
!**************************************************************************************
!*****************************************************************************
do jy=0,nymin1
do ix=0,nxmin1
......@@ -123,24 +136,26 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
3 continue
tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
tvold(ixm,jym)=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)
pold(ixm,jym)=ps(ixm,jym,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-tvold).gt.0.2) then
height(kz)=height(kz-1)+const*log(pold/pint)* &
(tv-tvold)/log(tv/tvold)
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))
else
height(kz)=height(kz-1)+const*log(pold/pint)*tv
height(kz)=height(kz-1)+ &
const*log(pold(ixm,jym)/pint(ixm,jym))*tv(ixm,jym)
endif
tvold=tv
pold=pint
tvold(ixm,jym)=tv(ixm,jym)
pold(ixm,jym)=pint(ixm,jym)
end do
......@@ -166,153 +181,182 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
! Loop over the whole grid
!*************************
do jy=0,nymin1
do ix=0,nxmin1
tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ps(ix,jy,1,n))
pold=ps(ix,jy,1,n)
uvwzlev(ix,jy,1)=0.
wzlev(1)=0.
rhoh(1)=pold/(r_air*tvold)
tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
ps(ix,jy,1,n))
enddo
enddo
pold=ps(:,:,1,n)
uvzlev(:,:,1)=0.
wzlev(:,:,1)=0.
rhoh(:,:,1)=pold/(r_air*tvold)
! Compute heights of eta levels
!******************************
do kz=2,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(kz)=pint/(r_air*tv)
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)
if (abs(tv-tvold).gt.0.2) then
uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* &
(tv-tvold)/log(tv/tvold)
else
uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv
endif
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
endwhere
tvold=tv
pold=pint
end do
tvold=tv
pold=pint
end do
do kz=2,nwz-1
wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2.
end do
wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1)
do kz=2,nwz-1
wzlev(:,:,kz)=(uvzlev(:,:,kz+1)+uvzlev(:,:,kz))/2.
end do
wzlev(:,:,nwz)=wzlev(:,:,nwz-1)+ &
uvzlev(:,:,nuvz)-uvzlev(:,:,nuvz-1)
! pinmconv=(h2-h1)/(p2-p1)
pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
(aknew(1)+bknew(1)*ps(ix,jy,1,n)))
do kz=2,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
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)))
pinmconv(:,:,1)=(uvzlev(:,:,2))/ &
((aknew(2)+bknew(2)*ps(:,:,1,n))- &
(aknew(1)+bknew(1)*ps(:,:,1,n)))
do kz=2,nz-1
pinmconv(:,:,kz)=(uvzlev(:,:,kz+1)-uvzlev(:,:,kz-1))/ &
((aknew(kz+1)+bknew(kz+1)*ps(:,:,1,n))- &
(aknew(kz-1)+bknew(kz-1)*ps(:,:,1,n)))
end do
pinmconv(:,:,nz)=(uvzlev(:,:,nz)-uvzlev(:,:,nz-1))/ &
((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
!************************************
uu(ix,jy,1,n)=uuh(ix,jy,1)
vv(ix,jy,1,n)=vvh(ix,jy,1)
tt(ix,jy,1,n)=tth(ix,jy,1,n)
qv(ix,jy,1,n)=qvh(ix,jy,1,n)
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
clwc(ix,jy,1,n)=clwch(ix,jy,1,n)
ciwc(ix,jy,1,n)=ciwch(ix,jy,1,n)
clwc(:,:,1,n)=clwch(:,:,1,n)
ciwc(:,:,1,n)=ciwch(:,:,1,n)
!hg
pv(ix,jy,1,n)=pvh(ix,jy,1)
rho(ix,jy,1,n)=rhoh(1)
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(:,:,1,n)=pvh(:,:,1)
rho(:,:,1,n)=rhoh(:,:,1)
uu(:,:,nz,n)=uuh(:,:,nuvz)
vv(:,:,nz,n)=vvh(:,:,nuvz)
tt(:,:,nz,n)=tth(:,:,nuvz,n)
qv(:,:,nz,n)=qvh(:,:,nuvz,n)
!hg adding the cloud water
clwc(ix,jy,nz,n)=clwch(ix,jy,nuvz,n)
ciwc(ix,jy,nz,n)=ciwch(ix,jy,nuvz,n)
clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
!hg
pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
rho(ix,jy,nz,n)=rhoh(nuvz)
kmin=2
do iz=2,nz-1
do kz=kmin,nuvz
if(height(iz).gt.uvwzlev(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)
pv(:,:,nz,n)=pvh(:,:,nuvz)
rho(:,:,nz,n)=rhoh(:,:,nuvz)
kmin=2
idx=kmin
do iz=2,nz-1
do jy=0,nymin1
do ix=0,nxmin1
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
clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
!hg
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.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
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)=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
idx(ix,jy)=kz
exit innuvz
endif
enddo innuvz
endif
enddo
enddo
do jy=0,nymin1
do ix=0,nxmin1
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)
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
!hg adding the cloud water
clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
!hg
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
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
! Levels, where w is given
!*************************
ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
kmin=2
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
kmin=kz
goto 40
ww(:,:,1,n)=wwh(:,:,1)*pinmconv(:,:,1)
ww(:,:,nz,n)=wwh(:,:,nwz)*pinmconv(:,:,nz)
kmin=2
idx=kmin
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
idx(ix,jy)=kz
exit inn
endif
end do
40 continue
end do
enddo inn
enddo
enddo
do jy=0,nymin1
do ix=0,nxmin1
kz=idx(ix,jy)
dz1=height(iz)-wzlev(ix,jy,kz-1)
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
enddo
enddo
! Compute density gradients at intermediate levels
!*************************************************
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
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
!****************************************************************
......@@ -321,61 +365,64 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!****************************************************************
do jy=1,ny-2
cosf=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/cosf
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
kmin=kz
goto 47
cosf(jy)=1./cos((real(jy)*dy+ylat0)*pi180)
enddo
kmin=2
idx=kmin
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
idx(ix,jy)=kz
exit inneta
endif
end do
47 ix1=ix-1
enddo inneta
enddo
enddo
do jy=1,ny-2
do ix=1,nx-2
kz=idx(ix,jy)
dz1=height(iz)-uvzlev(ix,jy,kz-1)
dz2=uvzlev(ix,jy,kz)-height(iz)
dz=dz1+dz2
ix1=ix-1
jy1=jy-1
ixp=ix+1
jyp=jy+1
dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/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*ui+dzdy*vi)
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
end do
end do
! If north pole is in the domain, calculate wind velocities in polar
! stereographic coordinates
!*******************************************************************
if (nglobal) then
do jy=int(switchnorthg)-2,nymin1
ylat=ylat0+real(jy)*dy
do ix=0,nxmin1
xlon=xlon0+real(ix)*dx
do iz=1,nz
do iz=1,nz
do jy=int(switchnorthg)-2,nymin1
ylat=ylat0+real(jy)*dy
do ix=0,nxmin1
xlon=xlon0+real(ix)*dx
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
......@@ -389,9 +436,11 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!
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
else if (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
......@@ -407,7 +456,9 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
ylat=90.0
uuaux=-ffpol*sin(xlonr+ddpol)
vvaux=-ffpol*cos(xlonr+ddpol)
call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
vvpolaux)
jy=nymin1
do ix=0,nxmin1
uupol(ix,jy,iz,n)=uupolaux
......@@ -440,13 +491,14 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!*******************************************************************
if (sglobal) then
do jy=0,int(switchsouthg)+3
ylat=ylat0+real(jy)*dy
do ix=0,nxmin1
xlon=xlon0+real(ix)*dx
do iz=1,nz
do iz=1,nz
do jy=0,int(switchsouthg)+3
ylat=ylat0+real(jy)*dy
do ix=0,nxmin1
xlon=xlon0+real(ix)*dx
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
......@@ -459,11 +511,14 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
!
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)