<p><b>duda</b> 2012-07-26 11:36:53 -0600 (Thu, 26 Jul 2012)</p><p>BRANCH COMMIT<br>
<br>
Merge all DCMIP test cases under config_test_case=9, and distinguish <br>
between cases using a new namelist variable config_dcmip_case.<br>
<br>
Still more work needed in DCMIP case 3-1 to add a thermal perturbation.<br>
<br>
<br>
M namelist.input.init_nhyd_atmos<br>
M src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F<br>
M src/core_init_nhyd_atmos/Registry<br>
</p><hr noshade><pre><font color="gray">Modified: branches/dcmip/namelist.input.init_nhyd_atmos
===================================================================
--- branches/dcmip/namelist.input.init_nhyd_atmos 2012-07-26 14:34:46 UTC (rev 2058)
+++ branches/dcmip/namelist.input.init_nhyd_atmos 2012-07-26 17:36:53 UTC (rev 2059)
@@ -1,26 +1,18 @@
&nhyd_model
- config_test_case = 9
+ config_start_time = '0000-01-01_00:00:00'
config_theta_adv_order = 3
- config_start_time = '2010-10-23_00:00:00'
- config_stop_time = '2010-10-23_00:00:00'
/
&dcmip
+ config_dcmip_case = '3-1'
config_planet_scale = 500.0
/
&dimensions
config_nvertlevels = 41
- config_nsoillevels = 4
- config_nfglevels = 38
- config_nfgsoillevels = 4
/
&data_sources
- config_geog_data_path = '/mmm/users/wrfhelp/WPS_GEOG/'
- config_met_prefix = 'CFSR'
- config_sfc_prefix = 'SST'
- config_fg_interval = 21600
/
&vertical_grid
@@ -30,10 +22,6 @@
/
&preproc_stages
- config_static_interp = .false.
- config_vertical_grid = .true.
- config_met_interp = .true.
- config_input_sst = .false.
/
&io
Modified: branches/dcmip/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/dcmip/src/core_init_nhyd_atmos/Registry 2012-07-26 14:34:46 UTC (rev 2058)
+++ branches/dcmip/src/core_init_nhyd_atmos/Registry 2012-07-26 17:36:53 UTC (rev 2059)
@@ -1,12 +1,13 @@
%
% namelist type namelist_record name default_value
%
-namelist integer nhyd_model config_test_case 7
+namelist integer nhyd_model config_test_case 9
namelist character nhyd_model config_calendar_type gregorian
namelist character nhyd_model config_start_time none
namelist character nhyd_model config_stop_time none
namelist integer nhyd_model config_theta_adv_order 3
namelist real nhyd_model config_coef_3rd_order 0.25
+namelist character dcmip config_dcmip_case 2-0-0
namelist real dcmip config_planet_scale 1.0
namelist integer dimensions config_nvertlevels 26
namelist integer dimensions config_nsoillevels 4
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-26 14:34:46 UTC (rev 2058)
+++ branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-07-26 17:36:53 UTC (rev 2059)
@@ -123,42 +123,59 @@
else if (config_test_case == 9 ) then
- write(0,*) ' reduced radius mountain wave test case '
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- write(0,*) ' calling test case setup '
- call init_atm_test_case_reduced_radius(block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, config_test_case)
- write(0,*) ' returned from test case setup '
- do i=2,nTimeLevs
- call mpas_copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
+ write(0,*) ' '
+ write(0,*) ' '
+ write(0,*) ' Setting up DCMIP test case '//trim(config_dcmip_case)
+ write(0,*) ' '
+ write(0,*) ' '
+
+ if (trim(config_dcmip_case) == '2-0-0' .or. &
+ trim(config_dcmip_case) == '2-0-1') then
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call init_atm_test_case_resting_atmosphere(block_ptr % mesh, block_ptr % state % time_levs(1) % state, &
+ block_ptr % diag, config_test_case)
+ block_ptr => block_ptr % next
end do
- block_ptr => block_ptr % next
- end do
+ else 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
- else if (config_test_case == 10 ) then
-
- write(0,*) ' resting atmosphere test case '
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- write(0,*) ' calling test case setup '
- call init_atm_test_case_resting_atmosphere(block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, config_test_case)
- write(0,*) ' returned from test case setup '
- do i=2,nTimeLevs
- call mpas_copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call init_atm_test_case_reduced_radius(block_ptr % mesh, block_ptr % state % time_levs(1) % state, &
+ block_ptr % diag, config_test_case)
+ block_ptr => block_ptr % next
end do
- block_ptr => block_ptr % next
- end do
+ else
+ write(0,*) ' '
+ write(0,*) ' *************'
+ write(0,*) ' Unrecognized DCMIP case '//trim(config_dcmip_case)
+ write(0,*) ' Please choose either 2-0-0, 2-0-1, 2-1, 2-1a, 2-2, or 3-1'
+ write(0,*) ' *************'
+ write(0,*) ' '
+ call mpas_dmpar_abort(domain % dminfo)
+
+ end if
+
else
- write(0,*) ' Only test cases 1 through 10 are currently supported for the nonhydrostatic core'
- stop
+ write(0,*) ' '
+ write(0,*) ' *************'
+ write(0,*) ' Only test cases 1 through 9 are currently supported for the nonhydrostatic core'
+ write(0,*) ' *************'
+ write(0,*) ' '
+ call mpas_dmpar_abort(domain % dminfo)
end if
+ ! Copy initialized state to all time levels
block_ptr => domain % blocklist
do while (associated(block_ptr))
do i=2,nTimeLevs
@@ -3260,9 +3277,7 @@
tempField % recvList => parinfo % cellsToRecv
tempField % array => hs
- call mpas_timer_start("EXCHANGE_1D_REAL")
call mpas_dmpar_exch_halo_field(tempField)
- call mpas_timer_stop("EXCHANGE_1D_REAL")
! dzmina = minval(hs(:)-hx(k-1,:))
dzmina = minval(zw(k)+ah(k)*hs(1:grid%nCellsSolve)-zw(k-1)-ah(k-1)*hx(k-1,1:grid%nCellsSolve))
@@ -4740,7 +4755,6 @@
! metrics for hybrid coordinate and vertical stretching
str = 1.0
-! zt = 20000.
zt = 30000.
dz = zt/float(nz1)
@@ -4809,8 +4823,13 @@
! setting for terrain
- xa = 5000.
- xla = 4000.
+! MGD for both 2-1 and 2-1a (and 2-2)
+ if (trim(config_dcmip_case) == '2-1' .or. &
+ trim(config_dcmip_case) == '2-1a' .or. &
+ trim(config_dcmip_case) == '2-2') then
+ xa = 5000.
+ xla = 4000.
+ end if
do iCell=1,grid % nCells
@@ -4824,15 +4843,27 @@
ri = sphere_distance(yi, xi, 0., 0., grid % sphere_radius)
+! MGD BEGIN 2-1
! Circular mountain with Schar mtn cross section
- hx(1,iCell) = hm*exp(-(ri/xa)**2)*cos(pii*ri/xla)**2
+ if (trim(config_dcmip_case) == '2-1') then
+ hx(1,iCell) = hm*exp(-(ri/xa)**2)*cos(pii*ri/xla)**2
+ end if
+! MGD END 2-1
+! MGD BEGIN 2-1a
+! proposed to be run with x333 rather than x500
! Ridge mountain with Schar mtn cross section
-! hx(1,iCell) = hm*exp(-(xc/xac)**2)*cos(pii*xc/xlac)**2*cos(yc/grid % sphere_radius)
+ if (trim(config_dcmip_case) == '2-1a') then
+ hx(1,iCell) = hm*exp(-(xc/xac)**2)*cos(pii*xc/xlac)**2*cos(yc/grid % sphere_radius)
+ end if
+! MGD END 2-1a
hx(nz,iCell) = zt
!***** SHP -> get the temporary point information for the neighbor cell ->> should be changed!!!!!
+
+!!! MGD WE NEED TO REPLACE THIS TERRAIN SMOOTHING WITH TC9
+
if (terrain_smooth) then
do i=1,grid % nCells
!option 1
@@ -4925,18 +4956,28 @@
!
! mountain wave initialization
!
+!MGD BEGIN 3-X
! Coefficients used to initialize 2 layer sounding based on stability
- zinv = 3000. ! Height of lower layer
- xn2 = 0.0001 ! N^2 for upper layer
- xn2m = 0.0000 ! N^2 for reference sounding
- xn2l = 0.0001 ! N^@ for lower layer
+ if (trim(config_dcmip_case) == '3-1') then
+ zinv = 3000. ! Height of lower layer
+ xn2 = 0.0001 ! N^2 for upper layer
+ xn2m = 0.0001 ! N^2 for reference sounding
+ xn2l = 0.0001 ! N^@ for lower layer
+ end if
+!MGD END 3-X
-! um = 10.
- um = 20.
-! shear = 0.00025
- shear = 0.
- us = 0.
+ if (trim(config_dcmip_case) == '2-1' .or. &
+ trim(config_dcmip_case) == '2-1a' .or. &
+ trim(config_dcmip_case) == '3-1') then
+ um = 20. ! base wind for 2-1, 2-1a, and 3-1
+ end if
+ if (trim(config_dcmip_case) == '2-2') then
+ shear = 0.00025 ! MGD 2-2
+ else
+ shear = 0. ! MGD everything else, 2-1, ...
+ end if
+
do i=1,grid % nCells
! Surface temp and Exner function as function of latitude to balance wind fed
@@ -4947,21 +4988,33 @@
do k=1,nz1
ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
+!MGD FOR 2-1, 2-1a, 2-2
! Isothermal temerature initialization
- t (k,i) = tsurf/pis*exp(gravity*ztemp/(cp*tsurf))
- tb(k,i) = t(k,i)
+ if (trim(config_dcmip_case) == '2-1' .or. &
+ trim(config_dcmip_case) == '2-1a' .or. &
+ trim(config_dcmip_case) == '2-2') then
+ t (k,i) = tsurf/pis*exp(gravity*ztemp/(cp*tsurf))
+ tb(k,i) = t(k,i)
+
+ end if
+
+!MGD FOR 3-X
! Initialization based on stability
-! if(ztemp .le. zinv) then
-! t (k,i) = t0*(1.+xn2l/gravity*ztemp)
-! else
-! t (k,i) = t0*(1.+xn2l/gravity*zinv+xn2/gravity*(ztemp-zinv))
-! end if
-! tb(k,i) = t0*(1. + xn2m/gravity*ztemp)
+ if (trim(config_dcmip_case) == '3-1') then
+ if(ztemp .le. zinv) then
+ t (k,i) = t0*(1.+xn2l/gravity*ztemp)
+ else
+ t (k,i) = t0*(1.+xn2l/gravity*zinv+xn2/gravity*(ztemp-zinv))
+ end if
+ tb(k,i) = t0*(1. + xn2m/gravity*ztemp)
+ end if
+ rh(k,i) = 0.
+ end do
- rh(k,i) = 0.
- end do
+! MGD ADD CODE HERE FOR 3-X THERMAL PERT
+
end do
!
@@ -4986,22 +5039,22 @@
+zgrid(k,cell2)+zgrid(k+1,cell2))
! Top of shear layer set at 10 km
- if(ztemp.lt.10000.) then
+! if(ztemp.lt.10000.) then
u(k,iEdge) = usurf * sqrt(1.+2.*shear*ztemp)
- else
- u(k,iEdge) = usurf * sqrt(1.+2.*shear*10000.)
- end if
+! else
+! u(k,iEdge) = usurf * sqrt(1.+2.*shear*10000.)
+! end if
end do
end do
deallocate(psiVertex)
do k=1,nz1
ztemp = .5*( zw(k)+zw(k+1))
- if(ztemp.lt.10000.) then
+! if(ztemp.lt.10000.) then
grid % u_init % array(k) = um * sqrt(1.+2.*shear*ztemp)
- else
- grid % u_init % array(k) = um * sqrt(1.+2.*shear*10000.)
- end if
+! else
+! grid % u_init % array(k) = um * sqrt(1.+2.*shear*10000.)
+! end if
end do
!
@@ -5263,10 +5316,18 @@
end do
end do
+! MGD FOR 3-X:
+! zt = 10000.0
+! nVertLevels = 10
+! X = 125
+! dt = 12.
+! nso = 8
+! 2nd-order horiz mixing = 50.0
+
end subroutine init_atm_test_case_reduced_radius
-!--------------------- TEST CASE 10 -----------------------------------------------
+!--------------------- TEST CASE 9 -----------------------------------------------
subroutine init_atm_test_case_resting_atmosphere(grid, state, diag, test_case)
@@ -5281,7 +5342,8 @@
type (diag_type), intent(inout) :: diag
integer, intent(in) :: test_case
- real (kind=RKIND), parameter :: t0=300., hm=2000., alpha=0.
+ real (kind=RKIND), parameter :: t0=300., alpha=0.
+ real (kind=RKIND) :: hm
real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx, cqw
@@ -5446,7 +5508,7 @@
! ah(k) = cos(.5*pii*zw(k)/zh)**6
ah(k) = cos(.5*pii*zw(k)/zh)**2
!
- ah(k) = ah(k)*(1.-zw(k)/zt)
+! ah(k) = ah(k)*(1.-zw(k)/zt)
!
else
ah(k) = 0.
@@ -5469,45 +5531,62 @@
write(0,*) 'EARTH RADIUS = ', grid % sphere_radius
- ! for hx computation
- r1m = .75*pii
- r2m = pii/16.
+! MGD 2-0-0, not used in 2-0-1
+ if (trim(config_dcmip_case) == '2-0-0') then
+ ! for hx computation
+ r1m = .75*pii
+ r2m = pii/16.
+ end if
+! MGD 2-0-1, not used in 2-0-0
+ if (trim(config_dcmip_case) == '2-0-1') then
! setting for terrain
+! xa = pii/16. ! for specifying mtn with in degrees
+ xa = pii*grid%sphere_radius/16. ! corresponds to ~11 grid intervals across entire mtn with 2 deg res
+ end if
- xa = pii/16. ! for specifying mtn with in degrees
- xa = pii*grid%sphere_radius/16. ! corresponds to ~11 grid intervals across entire mtn with 2 deg res
+! MGD both 2-0-0 and 2-0-1
+ hm = 2000.0
do iCell=1,grid % nCells
+
+ if (trim(config_dcmip_case) == '2-0-0') then
! Comb mountain as specified for DCMIP case 2.0
+! MGD BEGIN 2-0-0
+ xi = grid % lonCell % array(iCell)
+ yi = grid % latCell % array(iCell)
-! xi = grid % lonCell % array(iCell)
-! yi = grid % latCell % array(iCell)
+ rad = acos(cos(xi)*cos(yi))
-! rad = acos(cos(xi)*cos(yi))
+ if (rad.lt.r1m) THEN
+ hx(1,iCell) = hm*cos(.5*pii*rad/r1m)**2.*cos(pii*rad/r2m)**2
+ else
+ hx(1,iCell) = 0.
+ end if
+! MGD END 2-0-0
+ end if
-! if(rad.lt.r1m) THEN
-! hx(1,iCell) = hm*cos(.5*pii*rad/r1m)**2.*cos(pii*rad/r2m)**2
-! else
-! hx(1,iCell) = 0.
-! end if
+ if (trim(config_dcmip_case) == '2-0-1') then
+! cosine**2 ridge
+! MGD BEGIN 2-0-1
-! cosine**2 ridge, tapering toward the poles
+ xi = grid % lonCell % array(iCell)
+ yi = grid % latCell % array(iCell)
+ xc = sphere_distance(yi, xi, yi, 0., grid % sphere_radius)
+ yc = sphere_distance(yi, xi, 0., xi, grid % sphere_radius)
- xi = grid % lonCell % array(iCell)
- yi = grid % latCell % array(iCell)
- xc = sphere_distance(yi, xi, yi, 0., grid % sphere_radius)
- yc = sphere_distance(yi, xi, 0., xi, grid % sphere_radius)
- if(abs(xc).ge.xa) then
-! if(abs(xi).ge.xa.and.abs(2.*pii-xi).ge.xa) then
- hx(1,iCell) = 0.
- else
-! for mtn ridge with uniform width in km
- hx(1,iCell) = hm*cos(.5*pii*xc/xa)**2*cos(yc/grid % sphere_radius)
-! for mtn ridge with uniform width in degrees
-! hx(1,iCell) = hm*cos(.5*pii*xi/xa)**2*cos(yc/grid % sphere_radius)
+ if (abs(xc).ge.xa) then ! for mtn ridge with uniform width in km
+! if (abs(xi).ge.xa.and.abs(2.*pii-xi).ge.xa) then ! for mtn ridge with uniform width in degrees
+ hx(1,iCell) = 0.
+ else
+! for mtn ridge with uniform width in km
+ hx(1,iCell) = hm*cos(.5*pii*xc/xa)**2*cos(yc/grid % sphere_radius)
+! for mtn ridge with uniform width in degrees
+! hx(1,iCell) = hm*cos(.5*pii*xi/xa)**2*cos(yc/grid % sphere_radius)
+ end if
+! MGD END 2-0-1
end if
hx(:,iCell) = hx(1,iCell)
@@ -5524,25 +5603,15 @@
allocate(hs (grid % nCells+1))
allocate(hs1(grid % nCells+1))
-!! dzmin = .2
-! dzmin = .3
dzmin = .5
-! dzmin = .6
-! dzmin = .7
-! sm0 = .1
- sm0 = .05
-! sm0 = .02
-! sm0 = .01
+ sm0 = .05
-! nsm = 50
nsm = 30
-! nsm = 10
write(6,*) 'dzmin = ',dzmin,' sm0 = ',sm0,' nsm = ',nsm
do k=2,kz-1
-! do k=3,kz-1
hx(k,:) = hx(k-1,:)
dzminf = zw(k)-zw(k-1)
@@ -5661,12 +5730,22 @@
grid % cf3 % scalar = cf3
um = 0.
- gamma = .0065
+ gamma = .0065 ! temp lapse rate in K/km
- zinb = 3000.
-! zinb = zt
+! MGD BEGIN 2-0-0
+ if (trim(config_dcmip_case) == '2-0-0') then
+! zinb = zt
+ zint = zt ! no inversion layer
+ end if
+! MGD END 2-0-0
+! MGD BEGIN 2-0-1
+ if (trim(config_dcmip_case) == '2-0-1') then
+ zinb = 3000. ! bottom of inversion layer
+ zint = 5000. ! top of inversion layer
+ end if
+! MGD END 2-0-1
- zint = 5000.
+ ! computing intermediate T and Theta used to build sounding that includes inversion layer
tinv = t0-gamma*zinb
th_inb = t0*(1.-gamma*zinb/t0)**(1.-gravity/(cp*gamma))
th_int = th_inb*exp((gravity*(zint-zinb))/(cp*tinv))
</font>
</pre>