<p><b>duda</b> 2012-07-27 16:25:41 -0600 (Fri, 27 Jul 2012)</p><p>BRANCH COMMIT<br>
<br>
Add potential temperature perturbation for DCMIP case 3-1.<br>
<br>
<br>
M    src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F        2012-07-27 20:41:35 UTC (rev 2071)
+++ branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F        2012-07-27 22:25:41 UTC (rev 2072)
@@ -4635,6 +4635,14 @@
       real (kind=RKIND), parameter :: t0=300., hm=250., alpha=0.
 !      real (kind=RKIND), parameter :: t0=288., hm=0., alpha=0.
 
+      ! Parameters for test case 3-1
+      real (kind=RKIND), parameter :: widthParm = 5000.0, &amp;
+                                      dTheta = 1.0,       &amp;
+                                      L_z = 20000.0,      &amp;
+                                      theta_c = 0.0,      &amp;
+                                      lambda_c = 2.0 * pii / 3.0
+
+
       real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
       real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx, cqw
       real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru 
@@ -4676,6 +4684,7 @@
 
       real (kind=RKIND) :: xi, yi, ri, xa, xc, yc, xla, zinv, xn2, xn2m, xn2l, sm, dzh, dzht, dzmin, dzmina, dzminf, &amp;
                            dzmina_global, z_edge, z_edge3, sm0
+      real (kind=RKIND) :: theta_pert, s
 
       integer, dimension(grid % nCells, 2) :: next_cell
       real (kind=RKIND),  dimension(grid % nCells) :: hxzt, pitop, ptopb
@@ -5008,7 +5017,7 @@
 !
 ! mountain wave initialization
 !
-!MGD BEGIN 3-X
+!MGD BEGIN 3-1
 !        Coefficients used to initialize 2 layer sounding based on stability
          if (trim(config_dcmip_case) == '3-1') then
             zinv = 3000.     ! Height of lower layer
@@ -5016,7 +5025,7 @@
             xn2m = 0.0001    ! N^2 for reference sounding
             xn2l = 0.0001    ! N^@ for lower layer
          end if
-!MGD END 3-X
+!MGD END 3-1
 
          if (trim(config_dcmip_case) == '2-1' .or. &amp;
              trim(config_dcmip_case) == '2-1a' .or. &amp;
@@ -5051,7 +5060,7 @@
 
                end if
 
-!MGD FOR 3-X
+!MGD FOR 3-1
 !              Initialization based on stability
                if (trim(config_dcmip_case) == '3-1') then
                   if(ztemp .le. zinv) then
@@ -5065,8 +5074,6 @@
                rh(k,i) = 0. 
             end do
 
-! MGD ADD CODE HERE FOR 3-X THERMAL PERT
-
          end do
 
       !
@@ -5365,10 +5372,20 @@
          do k=1,grid%nVertLevels
             diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
             diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(index_qv,k,iCell))
+
+! MGD ADD CODE HERE FOR 3-1 THERMAL PERT
+            if (trim(config_dcmip_case) == '3-1') then
+               s = widthParm**2.0 / (widthParm**2.0 + sphere_distance(theta_c,                   lambda_c,              &amp;
+                                                                      grid%latCell%array(iCell), grid%lonCell%array(iCell), &amp;
+                                                                      grid%sphere_radius)**2.0)
+               theta_pert = dTheta * s * sin((2.0_RKIND * pii * zgrid(k,iCell)) / L_z)
+               diag % theta % array(k,iCell) = diag % theta % array(k,iCell) + theta_pert
+            end if
+
          end do
       end do
 
-! MGD FOR 3-X:
+! MGD FOR 3-1:
 !     zt = 10000.0
 !     nVertLevels = 10
 !     X = 125

</font>
</pre>