Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F	(revision 1198)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F	(revision 1199)
@@ -129,4 +129,9 @@
  INTEGER :: ierr
 !!MARS
+
+      REAL :: pfu, pfd, phm
+      INTEGER :: hypsometric_opt = 1 ! classic
+      !INTEGER :: hypsometric_opt = 2 ! Wee et al. 2012 correction
+
 
 
@@ -537,7 +542,17 @@
 !  respect to the model eqns.
 
+   IF (hypsometric_opt == 1) THEN
     DO k  = 2,kte
       grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)*grid%em_mub(i,j)*grid%em_alb(i,k-1,j)
     ENDDO
+   ELSE IF (hypsometric_opt == 2) THEN
+    DO k = 2,kte
+      pfu = grid%em_mub(i,j)*grid%em_znw(k)   + grid%p_top
+      pfd = grid%em_mub(i,j)*grid%em_znw(k-1)   + grid%p_top
+      phm = grid%em_mub(i,j)*grid%em_znu(k-1)   + grid%p_top
+      grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) + grid%em_alb(i,k-1,j)*phm*LOG(pfd/pfu)
+    END DO
+   END IF
+
 
   ENDDO
@@ -622,4 +637,5 @@
 
     grid%em_ph_1(i,1,j) = 0.
+   IF (hypsometric_opt == 1) THEN
     DO k  = 2,kte
       grid%em_ph_1(i,k,j) = grid%em_ph_1(i,k-1,j) - (1./grid%em_rdnw(k-1))*(       &
@@ -630,4 +646,26 @@
       grid%em_ph0(i,k,j) = grid%em_ph_1(i,k,j) + grid%em_phb(i,k,j)
     ENDDO
+   ELSE IF (hypsometric_opt == 2) THEN
+
+             ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
+             ! Note that al*p approximates Rd*T and dLOG(p) does z.
+             ! Here T varies mostly linear with z, the first-order integration produces better result.
+
+               grid%em_ph_2(i,1,j) = grid%em_phb(i,1,j)
+               DO k = 2,kte
+                  pfu = grid%em_mu0(i,j)*grid%em_znw(k)   + grid%p_top
+                  pfd = grid%em_mu0(i,j)*grid%em_znw(k-1) + grid%p_top
+                  phm = grid%em_mu0(i,j)*grid%em_znu(k-1) + grid%p_top
+                  grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) + grid%em_alt(i,k-1,j)*phm*LOG(pfd/pfu)
+               END DO
+
+               DO k = 1,kte
+                  grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k,j) - grid%em_phb(i,k,j)
+                  grid%em_ph_1(i,k,j) = grid%em_ph_2(i,k,j)
+               END DO
+
+   END IF
+
+
 
     IF ( wrf_dm_on_monitor() ) THEN
@@ -677,4 +715,6 @@
 !  rebalance hydrostatically
 
+   IF (hypsometric_opt == 1) THEN
+
       DO k  = 2,kte
         grid%em_ph_1(i,k,j) = grid%em_ph_1(i,k-1,j) - (1./grid%em_rdnw(k-1))*(       &
@@ -685,4 +725,26 @@
         grid%em_ph0(i,k,j) = grid%em_ph_1(i,k,j) + grid%em_phb(i,k,j)
       ENDDO
+
+   ELSE IF (hypsometric_opt == 2) THEN
+
+             ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
+             ! Note that al*p approximates Rd*T and dLOG(p) does z.
+             ! Here T varies mostly linear with z, the first-order integration produces better result.
+
+               grid%em_ph_2(i,1,j) = grid%em_phb(i,1,j)
+               DO k = 2,kte
+                  pfu = grid%em_mu0(i,j)*grid%em_znw(k)   + grid%p_top
+                  pfd = grid%em_mu0(i,j)*grid%em_znw(k-1) + grid%p_top
+                  phm = grid%em_mu0(i,j)*grid%em_znu(k-1) + grid%p_top
+                  grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) + grid%em_alt(i,k-1,j)*phm*LOG(pfd/pfu)
+               END DO
+
+               DO k = 1,kte
+                  grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k,j) - grid%em_phb(i,k,j)
+                  grid%em_ph_1(i,k,j) = grid%em_ph_2(i,k,j)
+               END DO
+
+   END IF
+
 
     ENDDO
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_real.F
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_real.F	(revision 1198)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_real.F	(revision 1199)
@@ -144,5 +144,7 @@
       REAl :: yeah, yeahc, yeah2
 !****MARS
-
+      INTEGER :: hypsometric_opt = 1 ! classic
+      !INTEGER :: hypsometric_opt = 2 ! Wee et al. 2012 correction
+      REAL :: pfu, pfd, phm
 
 
@@ -2264,7 +2266,17 @@
 
             grid%em_phb(i,1,j) = grid%ht(i,j) * g
-            DO k  = 2,kte
-               grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)*grid%em_mub(i,j)*grid%em_alb(i,k-1,j)
-            END DO
+            IF (hypsometric_opt == 1) THEN
+               DO k  = 2,kte
+                  grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)*grid%em_mub(i,j)*grid%em_alb(i,k-1,j)
+               END DO
+            ELSE IF (hypsometric_opt == 2) THEN
+               DO k = 2,kte
+                  pfu = grid%em_mub(i,j)*grid%em_znw(k)   + grid%p_top
+                  pfd = grid%em_mub(i,j)*grid%em_znw(k-1)   + grid%p_top
+                  phm = grid%em_mub(i,j)*grid%em_znu(k-1)   + grid%p_top
+                  grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) + grid%em_alb(i,k-1,j)*phm*LOG(pfd/pfu)
+               END DO
+            END IF
+
          END DO
       END DO
@@ -2382,14 +2394,33 @@
                   grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
                END DO
-         
+        
                !  This is the hydrostatic equation used in the model after the small timesteps.  In 
                !  the model, grid%em_al (inverse density) is computed from the geopotential.
          
-               DO k  = 2,kte
-                  grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) - &
-                                grid%em_dnw(k-1) * ( (grid%em_mub(i,j)+grid%em_mu_2(i,j))*grid%em_al(i,k-1,j) &
-                              + grid%em_mu_2(i,j)*grid%em_alb(i,k-1,j) )
-                  grid%em_ph0(i,k,j) = grid%em_ph_2(i,k,j) + grid%em_phb(i,k,j)
-               END DO
+               IF (hypsometric_opt == 1) THEN
+                  DO k  = 2,kte
+                     grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) - &
+                                   grid%em_dnw(k-1) * ( (grid%em_mub(i,j)+grid%em_mu_2(i,j))*grid%em_al(i,k-1,j) &
+                                 + grid%em_mu_2(i,j)*grid%em_alb(i,k-1,j) )
+                     grid%em_ph0(i,k,j) = grid%em_ph_2(i,k,j) + grid%em_phb(i,k,j)
+                  END DO
+               ELSE IF (hypsometric_opt == 2) THEN
+                  ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
+                  ! Note that al*p approximates Rd*T and dLOG(p) does z.
+                  ! Here T varies mostly linear with z, the first-order integration produces better result.
+                  PRINT*,"WEE ET AL. 2012 CORRECTION."
+                  grid%em_ph_2(i,1,j) = grid%em_phb(i,1,j)
+                  DO k = 2,kte
+                     pfu = grid%em_mu0(i,j)*grid%em_znw(k)   + grid%p_top
+                     pfd = grid%em_mu0(i,j)*grid%em_znw(k-1) + grid%p_top
+                     phm = grid%em_mu0(i,j)*grid%em_znu(k-1) + grid%p_top
+                     grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) + grid%em_alt(i,k-1,j)*phm*LOG(pfd/pfu)
+                  END DO
+
+                  DO k = 1,kte
+                     grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k,j) - grid%em_phb(i,k,j)
+                  END DO
+               END IF
+
 
 !               !  Adjust the column pressure so that the computed 500 mb height is close to the
@@ -2610,4 +2641,10 @@
                                      i, j, k
 
+      REAL :: pfu, pfd, phm
+      INTEGER :: hypsometric_opt = 1 ! classic
+      !INTEGER :: hypsometric_opt = 2 ! Wee et al. 2012 correction
+
+
+
 #ifdef DM_PARALLEL
 #    include "em_data_calls.inc"
@@ -2716,7 +2753,16 @@
 
             grid%em_phb(i,1,j) = grid%ht_fine(i,j) * g
-            DO k  = 2,kte
-               grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)*grid%em_mub(i,j)*grid%em_alb(i,k-1,j)
-            END DO
+            IF (hypsometric_opt == 1) THEN
+              DO k  = 2,kte
+                 grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)*grid%em_mub(i,j)*grid%em_alb(i,k-1,j)
+              END DO
+            ELSE IF (hypsometric_opt == 2) THEN
+              DO k = 2,kte
+                 pfu = grid%em_mub(i,j)*grid%em_znw(k)   + grid%p_top
+                 pfd = grid%em_mub(i,j)*grid%em_znw(k-1) + grid%p_top
+                 phm = grid%em_mub(i,j)*grid%em_znu(k-1) + grid%p_top
+                 grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) + grid%em_alb(i,k-1,j)*phm*LOG(pfd/pfu)
+              END DO
+            END IF
          END DO
       END DO
@@ -2771,14 +2817,34 @@
             END DO
       
-            !  This is the hydrostatic equation used in the model after the small timesteps.  In 
-            !  the model, grid%em_al (inverse density) is computed from the geopotential.
-      
-            DO k  = 2,kte
-               grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) - &
-                             grid%em_dnw(k-1) * ( (grid%em_mub(i,j)+grid%em_mu_2(i,j))*grid%em_al(i,k-1,j) &
-                           + grid%em_mu_2(i,j)*grid%em_alb(i,k-1,j) )
-               grid%em_ph0(i,k,j) = grid%em_ph_2(i,k,j) + grid%em_phb(i,k,j)
-            END DO
-       
+            !  This is the hydrostatic equation used in the model after the small timesteps.  In
+            !  the model, grid%al (inverse density) is computed from the geopotential.
+
+            IF (hypsometric_opt == 1) THEN
+               DO k  = 2,kte
+                  grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) - &
+                                grid%em_dnw(k-1) * ( (grid%em_mub(i,j)+grid%em_mu_2(i,j))*grid%em_al(i,k-1,j) &
+                              + grid%em_mu_2(i,j)*grid%em_alb(i,k-1,j) )
+                  grid%em_ph0(i,k,j) = grid%em_ph_2(i,k,j) + grid%em_phb(i,k,j)
+               END DO
+            ELSE IF (hypsometric_opt == 2) THEN
+
+             ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
+             ! Note that al*p approximates Rd*T and dLOG(p) does z.
+             ! Here T varies mostly linear with z, the first-order integration produces better result.
+
+               grid%em_ph_2(i,1,j) = grid%em_phb(i,1,j)
+               DO k = 2,kte
+                  pfu = grid%em_mu0(i,j)*grid%em_znw(k)   + grid%p_top
+                  pfd = grid%em_mu0(i,j)*grid%em_znw(k-1) + grid%p_top
+                  phm = grid%em_mu0(i,j)*grid%em_znu(k-1) + grid%p_top
+                  grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) + grid%em_alt(i,k-1,j)*phm*LOG(pfd/pfu)
+               END DO
+
+               DO k = 1,kte
+                  grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k,j) - grid%em_phb(i,k,j)
+               END DO
+
+            END IF
+
          END DO
       END DO