<p><b>duda</b> 2012-08-01 14:24:12 -0600 (Wed, 01 Aug 2012)</p><p>Fixes to 2-2 and 3-1 from Bill and Joe.<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-08-01 17:20:58 UTC (rev 2076)
+++ branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-08-01 20:24:12 UTC (rev 2077)
@@ -2074,7 +2074,7 @@
write(0,*) ' *** sounding for the simulation ***'
write(0,*) ' z theta pres qv rho_m u rr'
do k=1,nz1
- write(6,'(8(f14.9,2x))') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
+ write(0,'(8(f14.9,2x))') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
t(k,1)/(1.+1.61*scalars(index_qv,k,1)), &
.01*p0*p(k,1)**(1./rcp), &
1000.*scalars(index_qv,k,1), &
@@ -4773,9 +4773,34 @@
call atm_initialize_advection_rk(grid)
call atm_initialize_deformation_weights(grid)
- xnutr = 0.1 ! Coefficient for implicit w damping in absorbing layer
- zd = 20000. ! Bottom of absorbing layer
+ if (trim(config_dcmip_case) == '2-1') then
+ zt = 30000.
+ xnutr = 0.1 ! Coefficient for implicit w damping in absorbing layer
+ zd = 20000. ! Bottom of absorbing layer
+ write(0,*) ' test case 2-1, zt, zd, xnutr ', zt,zd,xnutr
+ end if
+ if (trim(config_dcmip_case) == '2-1a') then
+ zt = 20000.
+ xnutr = 0.1 ! Coefficient for implicit w damping in absorbing layer
+ zd = 10000. ! Bottom of absorbing layer
+ write(0,*) ' test case 2-1a, zt, zd, xnutr ', zt,zd,xnutr
+ end if
+
+ if (trim(config_dcmip_case) == '2-2') then
+ zt = 30000.
+ xnutr = 0.1 ! Coefficient for implicit w damping in absorbing layer
+ zd = 20000. ! Bottom of absorbing layer
+ write(0,*) ' test case 2-2, zt, zd, xnutr ', zt,zd,xnutr
+ end if
+
+ if (trim(config_dcmip_case) == '3-1') then
+ zt = 10000.
+ xnutr = 0.0 ! Coefficient for implicit w damping in absorbing layer
+ zd = 10000. ! Bottom of absorbing layer
+ write(0,*) ' test case 3-1, zt, zd, xnutr ', zt,zd,xnutr
+ end if
+
p0 = 1.e+05
rcp = rgas/cp
rcv = rgas/(cp-rgas)
@@ -4783,7 +4808,6 @@
! metrics for hybrid coordinate and vertical stretching
str = 1.0
- zt = 30000.
dz = zt/float(nz1)
! write(0,*) ' dz = ',dz
@@ -4859,7 +4883,10 @@
xla = 4000.
end if
+ write(0,*) ' hm, xa, xla ',hm,xa,xla
+ hx = 0.
+
do iCell=1,grid % nCells
xi = grid % lonCell % array(iCell)
@@ -4878,6 +4905,13 @@
end if
! MGD END 2-1
+! MGD BEGIN 2-2
+! Circular mountain with Schar mtn cross section
+ if (trim(config_dcmip_case) == '2-2') then
+ hx(1,iCell) = hm*exp(-(ri/xa)**2)*cos(pii*ri/xla)**2
+ end if
+! MGD END 2-2
+
! MGD BEGIN 2-1a
! proposed to be run with x333 rather than x500
! Ridge mountain with Schar mtn cross section
@@ -5029,8 +5063,9 @@
if (trim(config_dcmip_case) == '2-1' .or. &
trim(config_dcmip_case) == '2-1a' .or. &
+ trim(config_dcmip_case) == '2-2' .or. &
trim(config_dcmip_case) == '3-1') then
- um = 20. ! base wind for 2-1, 2-1a, and 3-1
+ um = 20. ! base wind for 2-1, 2-1a, 2-2, and 3-1
end if
if (trim(config_dcmip_case) == '2-2') then
@@ -5056,7 +5091,9 @@
trim(config_dcmip_case) == '2-2') then
t (k,i) = tsurf/pis*exp(gravity*ztemp/(cp*tsurf))
- tb(k,i) = t(k,i)
+ tb (k,i) = t0*exp(gravity*ztemp/(cp*t0))
+!! JBK fix, 20120801
+ !! tb(k,i) = t(k,i)
end if
@@ -5074,6 +5111,21 @@
rh(k,i) = 0.
end do
+
+! MGD ADD CODE HERE FOR 3-1 THERMAL PERT
+ if (trim(config_dcmip_case) == '3-1') then
+ do k=1,nz1
+ s = widthParm**2.0 / (widthParm**2.0 + sphere_distance(theta_c, lambda_c, &
+ grid%latCell%array(i), grid%lonCell%array(i), &
+ grid%sphere_radius)**2.0)
+ theta_pert = dTheta * s * sin((2.0_RKIND * pii * 0.5*(zgrid(k,i)+zgrid(k+1,i))) / L_z)
+ ! diag % theta % array(k,i) = diag % theta % array(k,i) + theta_pert
+ t(k,i) = t(k,i) + theta_pert
+ end do
+ end if
+
+
+
end do
!
@@ -5122,9 +5174,11 @@
do i=1, grid % nCells
tsurf = t0*exp(-shear*um**2/gravity*sin(grid%latCell%array(i))**2)
- pis = exp(-um**2*sin(grid%latCell%array(i))**2/(2.*cp*tsurf))
-! pis = 1.
+!! JBK fix 20120801
+!! pis = exp(-um**2*sin(grid%latCell%array(i))**2/(2.*cp*tsurf))
+ pis = 1.
+
pitop(i) = pis-.5*dzw(1)*gravity/(cp*tb(1,1)*zz(1,1))
do k=2,nz1
pitop(i) = pitop(i)-dzu(k)*gravity/(cp*(fzm(k)*tb(k,1)+fzp(k)*tb(k-1,1)) &
@@ -5215,7 +5269,7 @@
write(0,*) ' *** sounding for the simulation ***'
write(0,*) ' z theta pres qv rho_m u rr'
do k=1,nz1
- write(6,'(8(f14.9,2x))') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
+ write(0,'(8(f14.9,2x))') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
t(k,1)/(1.+1.61*scalars(index_qv,k,1)), &
.01*p0*p(k,1)**(1./rcp), &
1000.*scalars(index_qv,k,1), &
@@ -5372,16 +5426,6 @@
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, &
- grid%latCell%array(iCell), grid%lonCell%array(iCell), &
- grid%sphere_radius)**2.0)
- theta_pert = dTheta * s * sin((2.0_RKIND * pii * 0.5*(zgrid(k,iCell)+zgrid(k+1,iCell))) / L_z)
- diag % theta % array(k,iCell) = diag % theta % array(k,iCell) + theta_pert
- end if
-
end do
end do
@@ -5976,7 +6020,7 @@
write(0,*) ' *** sounding for the simulation ***'
write(0,*) ' z temp theta pres rho_m u rr'
do k=1,nz1
- write(6,'(8(f14.9,2x))') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
+ write(0,'(8(f14.9,2x))') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
t(k,1)/(1.+1.61*scalars(index_qv,k,1))*p(k,1), &
t(k,1)/(1.+1.61*scalars(index_qv,k,1)), &
.01*p0*p(k,1)**(1./rcp), &
</font>
</pre>