Changeset 24 for trunk/src/calcpv_nests.f90

Timestamp:

May 23, 2014, 11:48:41 AM (10 years ago)

Author:

igpis

Message:

version 9.2 beta. Changes from HH, AST, MC, NIK, IP. Changes in vert transform. New SPECIES input includes scavenging coefficients

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,ppmk
-  real :: ppml(nuvzmax)
+  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 :: 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.
@@ -88 +97 @@
       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
@@ -97 +102 @@
 
       do kl=1,nuvz
-        ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
-        theta=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa
+        theta=tthn(ix,jy,kl,n,l)*ppmk(ix,jy,kl)
         klvrp=kl+1
         klvrm=kl-1
@@ -107 +111 @@
         if (klvrp.gt.nuvz) klvrp=nuvz
         if (klvrm.lt.1) klvrm=1
-        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))
+        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))
 
   ! Compute vertical position at pot. temperature surface on subgrid
@@ -134 +136 @@
           kch=kch+1
           k=kup
-          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
+          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
+          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
 
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
@@ -158 +158 @@
           kch=kch+1
           k=kdn
-          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
+          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
+          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
            ((thdn.le.theta).and.(thup.ge.theta))) then
@@ -212 +210 @@
           kch=kch+1
           k=kup
-          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
+          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
+          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
            ((thdn.le.theta).and.(thup.ge.theta))) then
@@ -235 +231 @@
           kch=kch+1
           k=kdn
-          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
+          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
+          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
       if (((thdn.ge.theta).and.(thup.le.theta)).or. &
            ((thdn.le.theta).and.(thup.ge.theta))) then