Changes in trunk/src/calcpv_nests.f90 [24:4]

File:

• trunk/src/calcpv_nests.f90 (modified) (10 diffs)

trunk/src/calcpv_nests.f90

@@ -21 +21 @@
 
 subroutine calcpv_nests(l,n,uuhn,vvhn,pvhn)
-  !                     i i  i    i    o
+  !                        i i  i    i    o
   !*****************************************************************************
   !                                                                            *
@@ -48 +48 @@
   integer :: nlck,l
   real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
-  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt
-  real :: ppml(0:nxmaxn-1,0:nymaxn-1,nuvzmax),ppmk(0:nxmaxn-1,0:nymaxn-1,nuvzmax)
+  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
+  real :: ppml(nuvzmax)
   real :: thup,thdn
   real,parameter :: eps=1.e-5,p0=101325
@@ -61 +61 @@
   ! Loop over entire grid
   !**********************
-  do kl=1,nuvz
-    do jy=0,nyn(l)-1
-      do ix=0,nxn(l)-1
-         ppml(ix,jy,kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
-      enddo
-    enddo
-  enddo
-  ppmk=(100000./ppml)**kappa
-
   do jy=0,nyn(l)-1
     phi = (ylat0n(l) + jy * dyn(l)) * pi / 180.
@@ -97 +88 @@
       if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1
       jux=jumpx
+  ! Precalculate pressure values for efficiency
+      do kl=1,nuvz
+        ppml(kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
+      end do
   !
   ! Loop over the vertical
@@ -102 +97 @@
 
       do kl=1,nuvz
-        theta=tthn(ix,jy,kl,n,l)*ppmk(ix,jy,kl)
+        ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
+        theta=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa
         klvrp=kl+1
         klvrm=kl-1
@@ -111 +107 @@
         if (klvrp.gt.nuvz) klvrp=nuvz
         if (klvrm.lt.1) klvrm=1
-        thetap=tthn(ix,jy,klvrp,n,l)*ppmk(ix,jy,klvrp)
-        thetam=tthn(ix,jy,klvrm,n,l)*ppmk(ix,jy,klvrm)
-        dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm))
+        ppmk=akz(klvrp)+bkz(klvrp)*psn(ix,jy,1,n,l)
+        thetap=tthn(ix,jy,klvrp,n,l)*(100000./ppmk)**kappa
+        ppmk=akz(klvrm)+bkz(klvrm)*psn(ix,jy,1,n,l)
+        thetam=tthn(ix,jy,klvrm,n,l)*(100000./ppmk)**kappa
+        dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
 
   ! Compute vertical position at pot. temperature surface on subgrid
@@ -136 +134 @@
           kch=kch+1
           k=kup
-          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
-          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
+          ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
+          thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
+          thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa
 
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
@@ -158 +158 @@
           kch=kch+1
           k=kdn
-          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
-          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
+          ppmk=akz(k)+bkz(k)*psn(ivr,jy,1,n,l)
+          thdn=tthn(ivr,jy,k,n,l)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*psn(ivr,jy,1,n,l)
+          thup=tthn(ivr,jy,k+1,n,l)*(100000./ppmk)**kappa
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
            ((thdn.le.theta).and.(thup.ge.theta))) then
@@ -210 +212 @@
           kch=kch+1
           k=kup
-          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
-          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
+          ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
+          thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
+          thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
            ((thdn.le.theta).and.(thup.ge.theta))) then
@@ -231 +235 @@
           kch=kch+1
           k=kdn
-          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
-          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
+          ppmk=akz(k)+bkz(k)*psn(ix,j,1,n,l)
+          thdn=tthn(ix,j,k,n,l)*(100000./ppmk)**kappa
+          ppmk=akz(k+1)+bkz(k+1)*psn(ix,j,1,n,l)
+          thup=tthn(ix,j,k+1,n,l)*(100000./ppmk)**kappa
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
            ((thdn.le.theta).and.(thup.ge.theta))) then