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 1247)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_real.F	(revision 1248)
@@ -145,5 +145,7 @@
 !****MARS
       INTEGER :: hypsometric_opt = 1 ! classic
+      INTEGER :: interp_theta = .true. ! classic
       !INTEGER :: hypsometric_opt = 2 ! Wee et al. 2012 correction
+      !INTEGER :: interp_theta = .false. ! Wee et al. 2012 correction
       REAL :: pfu, pfd, phm
 
@@ -1173,5 +1175,6 @@
                                its , ite , jts , jte , kts , kte )
          END IF
-   
+  
+IF ( interp_theta ) THEN 
          !  The input field is temperature, we want potential temp.
 !****MARS: here em_p_gc is really needed ! 
@@ -1180,5 +1183,6 @@
                            ims , ime , jms , jme , 1   , num_metgrid_levels , &
                            its , ite , jts , jte , 1   , num_metgrid_levels )
-   
+ENDIF   
+
 
 
@@ -1261,5 +1265,10 @@
                             ims , ime , jms , jme , kms , kme , &
                             its , ite , jts , jte , kts , kte )
- 
+
+         !  Depending on the setting of interp_theta = T/F, t_gc is is either theta Xor 
+         !  temperature, and that means that the t_2 field is also the associated field.
+         !  It is better to interpolate temperature and potential temperature in LOG(p),
+         !  regardless of requested default.
+
          !!****MARS: normalement c'est vert_interp
          CALL vert_interp_old ( grid%em_t_gc , grid%em_pd_gc , grid%em_t_2               , grid%em_pb , &
@@ -1270,4 +1279,31 @@
                             ims , ime , jms , jme , kms , kme , &
                             its , ite , jts , jte , kts , kte )
+
+         IF ( .NOT. interp_theta ) THEN
+
+            !! correction Wee et al. 2012 
+            !! first interpolate temperature (see above)
+            !! then interpolate pressure
+            !! and in the end compute potential temperature
+
+            !! scalar just an intermediate thing for interpolated pressure 
+            !! -- it is reinitialized afterwards
+            !! It is better to interpolate pressure in p regardless default options
+            CALL vert_interp_old ( grid%em_p_gc , grid%em_pd_gc , scalar(:,:,:,1)         , grid%em_pb , &
+                            num_metgrid_levels , 'T' , &
+                            1, lagrange_order , lowest_lev_from_sfc , &
+                            zap_close_levels , force_sfc_in_vinterp , &
+                            ids , ide , jds , jde , kds , kde , &
+                            ims , ime , jms , jme , kms , kme , &
+                            its , ite , jts , jte , kts , kte )
+
+            CALL t_to_theta ( grid%em_t_2 , scalar(:,:,:,1) , p00 , &
+                              ids , ide , jds , jde , kds , kde , &
+                              ims , ime , jms , jme , kms , kme , &
+                              its , ite , jts , jte , kts , kte )
+      
+            scalar(:,:,:,1) = 0.
+         END IF
+