<p><b>duda</b> 2011-02-07 15:23:18 -0700 (Mon, 07 Feb 2011)</p><p>BRANCH COMMIT<br>
<br>
- Add code from WRF height core initialization to compute<br>
hydrostatically balanced density from interpolated temperature and<br>
height fields.<br>
<br>
- Add namelist options to control whether: (1) static fields are<br>
interpolated and grid lengths are scaled by the radius of the sphere;<br>
(2) a vertical mesh is generated; (3) and meteorological fields are<br>
interpolated in three dimensions.<br>
<br>
<br>
M src/core_init_nhyd_atmos/module_test_cases.F<br>
M src/core_init_nhyd_atmos/Registry<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-02-07 21:48:34 UTC (rev 714)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-02-07 22:23:18 UTC (rev 715)
@@ -35,6 +35,9 @@
namelist character data_sources config_geog_data_path /data3/mp/wrfhelp/WPS_GEOG/
namelist character data_sources config_met_prefix GFS
namelist character data_sources config_init_date 2010-05-26_12
+namelist logical preproc_stages config_static_interp true
+namelist logical preproc_stages config_vertical_grid true
+namelist logical preproc_stages config_met_interp true
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
@@ -150,15 +153,15 @@
var persistent real dss ( nVertLevels nCells ) 0 io dss mesh - -
# Horizontally interpolated from first-guess data
-var persistent real u_fg ( nFGLevels nEdges ) 0 o u fg - -
-var persistent real v_fg ( nFGLevels nEdges ) 0 o v fg - -
-var persistent real t_fg ( nFGLevels nCells ) 0 o t fg - -
-var persistent real p_fg ( nFGLevels nCells ) 0 o p fg - -
-var persistent real z_fg ( nFGLevels nCells ) 0 o z fg - -
-var persistent real rh_fg ( nFGLevels nCells ) 0 o rh fg - -
-var persistent real soilz_fg ( nCells ) 0 o soilz fg - -
-var persistent real psfc_fg ( nCells ) 0 o psfc fg - -
-var persistent real pmsl_fg ( nCells ) 0 o pmsl fg - -
+var persistent real u_fg ( nFGLevels nEdges ) 0 - u fg - -
+var persistent real v_fg ( nFGLevels nEdges ) 0 - v fg - -
+var persistent real t_fg ( nFGLevels nCells ) 0 - t fg - -
+var persistent real p_fg ( nFGLevels nCells ) 0 - p fg - -
+var persistent real z_fg ( nFGLevels nCells ) 0 - z fg - -
+var persistent real rh_fg ( nFGLevels nCells ) 0 - rh fg - -
+var persistent real soilz_fg ( nCells ) 0 - soilz fg - -
+var persistent real psfc_fg ( nCells ) 0 - psfc fg - -
+var persistent real pmsl_fg ( nCells ) 0 - pmsl fg - -
# Horizontally interpolated from first-guess data
# and should be read in by model
@@ -173,8 +176,8 @@
var persistent real skintemp ( nCells ) 0 o skintemp fg - -
var persistent real snow ( nCells ) 0 o snow fg - -
var persistent real seaice ( nCells ) 0 o seaice fg - -
-var persistent real gfs_z ( nVertLevels nCells ) 0 o gfs_z fg - -
-var persistent real gfs_p ( nVertLevelsP1 nCells ) 0 o gfs_p fg - -
+var persistent real gfs_z ( nVertLevels nCells ) 0 - gfs_z fg - -
+#var persistent real gfs_p ( nVertLevelsP1 nCells ) 0 - gfs_p fg - -
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 io u state - -
@@ -203,6 +206,7 @@
var persistent real exner ( nVertLevels nCells Time ) 1 - exner diag - -
var persistent real exner_base ( nVertLevels nCells Time ) 1 io exner_base diag - -
var persistent real rtheta_base ( nVertLevels nCells Time ) 1 - rtheta_base diag - -
+var persistent real pressure ( nVertLevels nCells Time ) 1 io pressure diag - -
var persistent real pressure_base ( nVertLevels nCells Time ) 1 io pressure_base diag - -
var persistent real rho_base ( nVertLevels nCells Time ) 1 io rho_base diag - -
var persistent real theta_base ( nVertLevels nCells Time ) 1 io theta_base diag - -
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-02-07 21:48:34 UTC (rev 714)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-02-07 22:23:18 UTC (rev 715)
@@ -2020,7 +2020,8 @@
integer :: nsmterrain, kz, sfc_k
logical :: hybrid, smooth
- logical :: config_do_static
+ integer :: it
+ real (kind=RKIND) :: p_check
! For interpolating terrain and land use
integer :: nx, ny, nzz, iPoint, subx, suby
@@ -2125,13 +2126,11 @@
p0 = 1.e+05
- config_do_static = .true.
-
!
! Scale all distances and areas from a unit sphere to one with radius a
!
- if (config_do_static) then
+ if (config_static_interp) then
grid % xCell % array = grid % xCell % array * a
grid % yCell % array = grid % yCell % array * a
@@ -2412,12 +2411,13 @@
deallocate(rarray)
deallocate(ncat)
- end if ! config_do_static
+ end if ! config_static_interp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN ADOPT GFS TERRAIN HEIGHT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+#if 0
call read_met_init(trim(config_met_prefix), .false., trim(config_init_date), istatus)
if (istatus /= 0) then
@@ -2474,12 +2474,15 @@
end do
call read_met_close()
+#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END ADOPT GFS TERRAIN HEIGHT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ if (config_vertical_grid) then
+
!
! Vertical grid setup
!
@@ -2708,9 +2711,67 @@
! write(0,*) ' k, zx(k,1) ',k,zx(k,1)
! enddo
+
+ ! For z-metric term in omega equation
+ do iEdge = 1,grid % nEdges
+ cell1 = CellsOnEdge(1,iEdge)
+ cell2 = CellsOnEdge(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+
+ do k = 1, grid%nVertLevels
+
+ if (config_theta_adv_order == 2) then
+
+ z_edge = (zgrid(k,cell1)+zgrid(k,cell2))/2.
+
+ else !theta_adv_order == 3 or 4
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * zgrid(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * zgrid(k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ z_edge = 0.5*(zgrid(k,cell1) + zgrid(k,cell2)) &
+ - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
+
+ if (config_theta_adv_order == 3) then
+ z_edge3 = - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ else
+ z_edge3 = 0.
+ end if
+
+ end if
+
+ zb(k,1,iEdge) = (z_edge-zgrid(k,cell1))*dvEdge(iEdge)/AreaCell(cell1)
+ zb(k,2,iEdge) = (z_edge-zgrid(k,cell2))*dvEdge(iEdge)/AreaCell(cell2)
+ zb3(k,1,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell1)
+ zb3(k,2,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell2)
+
+ if (k /= 1) then
+ zf(k,1,iEdge) = ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb(k,1,iEdge)
+ zf(k,2,iEdge) = ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb(k,2,iEdge)
+ zf3(k,1,iEdge)= ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb3(k,1,iEdge)
+ zf3(k,2,iEdge)= ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb3(k,2,iEdge)
+ end if
+
+ end do
+
+ end if
+ end do
+
write(0,*) ' grid metrics setup complete '
+ end if ! config_vertical_grid
+
+ if (config_met_interp) then
+
!
! Horizontally interpolate meteorological data
!
@@ -3044,26 +3105,26 @@
call quicksort(config_nfglevels, sorted_arr)
do k=1,grid%nVertLevels
target_z = 0.5 * (grid % zgrid % array(k,iCell) + grid % zgrid % array(k+1,iCell))
- diag % pressure_base % array(k,iCell) = exp(vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=1))
+ diag % pressure % array(k,iCell) = exp(vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=1))
end do
! PRESSURE
- sorted_arr(:,:) = -999.0
- do k=1,config_nfglevels
- sorted_arr(1,k) = fg % z % array(k,iCell)
- if (vert_level(k) == 200100.0) then
-!NOSFC sorted_arr(1,k) = fg % soilz % array(iCell)
- sorted_arr(1,k) = 99999.0
- sfc_k = k
- end if
- sorted_arr(2,k) = log(fg % p % array(k,iCell))
- end do
- call quicksort(config_nfglevels, sorted_arr)
- do k=1,grid%nVertLevels+1
- target_z = grid % zgrid % array(k,iCell)
- fg % gfs_p % array(k,iCell) = exp(vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=1))
- end do
+! sorted_arr(:,:) = -999.0
+! do k=1,config_nfglevels
+! sorted_arr(1,k) = fg % z % array(k,iCell)
+! if (vert_level(k) == 200100.0) then
+!!NOSFC sorted_arr(1,k) = fg % soilz % array(iCell)
+! sorted_arr(1,k) = 99999.0
+! sfc_k = k
+! end if
+! sorted_arr(2,k) = log(fg % p % array(k,iCell))
+! end do
+! call quicksort(config_nfglevels, sorted_arr)
+! do k=1,grid%nVertLevels+1
+! target_z = grid % zgrid % array(k,iCell)
+! fg % gfs_p % array(k,iCell) = exp(vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=1))
+! end do
end do
@@ -3120,52 +3181,26 @@
!
do iCell=1,grid%nCells
do k=1,grid%nVertLevels
+
+ ! QV
es = 6.112 * exp((17.27*(state % theta % array(k,iCell) - 273.16))/(state % theta % array(k,iCell) - 35.86))
- rs = 0.622 * es / (diag % pressure_base % array(k,iCell) - es)
- state % scalars % array(state % index_qv,k,iCell) = rs * state % scalars % array(state % index_qv,k,iCell)
+ rs = 0.622 * es / (diag % pressure % array(k,iCell) - es)
+ scalars(state % index_qv,k,iCell) = rs * scalars(state % index_qv,k,iCell)
- rgas_moist = (1.0 + 0.61 * state % scalars % array(state % index_qv,k,iCell)) * rgas
- state % rho % array(k,iCell) = diag % pressure_base % array(k,iCell) / (rgas_moist * state % theta % array(k,iCell))
- state % rho % array(k,iCell) = state % rho % array(k,iCell) / (1.0 + state % scalars % array(state % index_qv,k,iCell))
- state % rho % array(k,iCell) = state % rho % array(k,iCell) / zz(k,iCell)
+ ! PI
+ p(k,iCell) = (diag % pressure % array(k,iCell) / p0) ** (rgas / cp)
- state % theta % array(k,iCell) = state % theta % array(k,iCell) * (p0 / diag % pressure_base % array(k,iCell)) ** (rgas / cp)
- state % theta % array(k,iCell) = state % theta % array(k,iCell) * (1.0 + 1.61*state % scalars % array(state % index_qv,k,iCell))
- end do
- end do
+ ! THETA - can compute this using PI instead
+! t(k,iCell) = t(k,iCell) / p(k,iCell)
+ t(k,iCell) = t(k,iCell) * (p0 / diag % pressure % array(k,iCell)) ** (rgas / cp)
-#if 0
- do iCell=1,grid%nCells
-
- fsum(:) = 0.0
- do k=2,grid%nVertLevels+1
- rgas_moist = (1.0 + 0.61 * state % scalars % array(state % index_qv,k-1,iCell)) * rgas
-!RGAS rgas_moist = rgas
- ptmp = fg % gfs_p % array(k,iCell)
- state % theta % array(k-1,iCell) = -gravity * (zgrid(k,iCell) - zgrid(k-1,iCell)) / (rgas_moist * (log(ptmp / fg % gfs_p % array(1,iCell)) + fsum(k-1)))
- fsum(k) = fsum(k-1) + gravity * (zgrid(k,iCell) - zgrid(k-1,iCell)) / (rgas_moist * state % theta % array(k-1,iCell))
-
- state % theta % array(k-1,iCell) = state % theta % array(k-1,iCell) * (p0 / diag % pressure_base % array(k-1,iCell)) ** (rgas / cp)
- state % theta % array(k-1,iCell) = state % theta % array(k-1,iCell) * (1.0 + 1.61*state % scalars % array(state % index_qv,k-1,iCell))
+ ! RHO
+ rho(k,iCell) = diag % pressure % array(k,iCell) / rgas / (p(k,iCell) * t(k,iCell))
+ rho(k,iCell) = rho(k,iCell) / (1.0 + scalars(state % index_qv,k,iCell))
end do
end do
-#endif
-#if 0
- do iCell=1,grid%nCells
- do k=2,grid%nVertLevels+1
- rgas_moist = (1.0 + 0.61 * state % scalars % array(state % index_qv,k-1,iCell)) * rgas
- state % theta % array(k-1,iCell) = -gravity * (zgrid(k,iCell) - zgrid(k-1,iCell)) / (rgas_moist * log(fg % gfs_p % array(k,iCell) / fg % gfs_p % array(k-1,iCell)))
-
- state % theta % array(k-1,iCell) = state % theta % array(k-1,iCell) * (p0 / diag % pressure_base % array(k-1,iCell)) ** (rgas / cp)
- state % theta % array(k-1,iCell) = state % theta % array(k-1,iCell) * (1.0 + 1.61*state % scalars % array(state % index_qv,k-1,iCell))
- end do
- end do
-#endif
-
-
-
!
! Reference state based on a dry isothermal atmosphere
!
@@ -3174,7 +3209,8 @@
ztemp = 0.5*(zgrid(k+1,iCell)+zgrid(k,iCell))
ppb(k,iCell) = p0*exp(-gravity*ztemp/(rgas*t0b)) ! pressure_base
pb (k,iCell) = (ppb(k,iCell)/p0)**(rgas/cp) ! exner_base
- rb (k,iCell) = ppb(k,iCell)/(rgas*t0b*zz(k,iCell)) ! rho_base
+! rb (k,iCell) = ppb(k,iCell)/(rgas*t0b*zz(k,iCell)) ! rho_base
+ rb (k,iCell) = ppb(k,iCell)/(rgas*t0b) ! rho_base
tb (k,iCell) = t0b/pb(k,iCell) ! theta_base
rtb(k,iCell) = rb(k,iCell)*tb(k,iCell) ! rtheta_base
p (k,iCell) = pb(k,iCell) ! exner
@@ -3183,66 +3219,55 @@
end do
end do
- do iEdge=1,grid%nEdges
+ do iCell=1,grid%nCells
do k=1,grid%nVertLevels
- diag % ru % array(k,iEdge) = state % u % array(k,iEdge) * 0.5*(state % rho % array(k,cellsOnEdge(1,iEdge)) + state % rho % array(k,cellsOnEdge(2,iEdge)))
+ pp(k,iCell) = diag % pressure % array(k,iCell) - ppb(k,iCell)
+ rr(k,iCell) = rho(k,iCell) - rb(k,iCell)
end do
end do
+ do iCell=1,grid%nCells
+ k = 1
+ rho(k,iCell) = ((diag % pressure % array(k,iCell) / p0)**(cv / cp)) * (p0 / rgas) / (t(k,iCell)*(1.0 + 1.61*scalars(state % index_qv,k,iCell)))
+ rr(k,iCell) = rho(k,iCell) - rb(k,iCell)
- ! For z-metric term in omega equation
- do iEdge = 1,grid % nEdges
- cell1 = CellsOnEdge(1,iEdge)
- cell2 = CellsOnEdge(2,iEdge)
- if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+ do k=2,grid % nVertLevels
+ it = 0
+ p_check = 2.0 * 0.0001
+ do while ( (it < 30) .and. (p_check > 0.0001) )
- do k = 1, grid%nVertLevels
+ p_check = pp(k,iCell)
+ dz = (zgrid(k,iCell) - zgrid(k-1,iCell))
+ pp(k,iCell) = pp(k-1,iCell) - 0.5 * (rr(k,iCell) + rr(k-1,iCell))*gravity*dz &
+ - 0.5 * (rho(k,iCell)*scalars(state % index_qv,k,iCell) + rho(k-1)*scalars(state % index_qv,k-1,iCell))*gravity*dz
+ diag % pressure % array(k,iCell) = pp(k,iCell) + ppb(k,iCell)
+ p(k,iCell) = (diag % pressure % array(k,iCell) / p0) ** (rgas / cp)
+ rho(k,iCell) = diag % pressure % array(k,iCell) / rgas / (p(k,iCell)*t(k,iCell)*(1.0 + 1.61*scalars(state % index_qv,k,iCell)))
+ rr(k,iCell) = rho(k,iCell) - rb(k,iCell)
- if (config_theta_adv_order == 2) then
+ p_check = abs(p_check - pp(k,iCell))
+
+ it = it + 1
+ end do
+ end do
+ end do
- z_edge = (zgrid(k,cell1)+zgrid(k,cell2))/2.
+ ! Compute theta_m and rho-tilde
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ t(k,iCell) = t(k,iCell) * (1.0 + 1.61*scalars(state % index_qv,k,iCell))
+ rho(k,iCell) = rho(k,iCell) / zz(k,iCell)
+ rb(k,iCell) = rb(k,iCell) / zz(k,iCell)
+ end do
+ end do
- else !theta_adv_order == 3 or 4
+ do iEdge=1,grid%nEdges
+ do k=1,grid%nVertLevels
+ diag % ru % array(k,iEdge) = state % u % array(k,iEdge) * 0.5*(state % rho % array(k,cellsOnEdge(1,iEdge)) + state % rho % array(k,cellsOnEdge(2,iEdge)))
+ end do
+ end do
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * zgrid(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * zgrid(k,cell2)
- do i=1, grid % nEdgesOnCell % array (cell1)
- if ( grid % CellsOnCell % array (i,cell1) > 0) &
- d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
- end do
- do i=1, grid % nEdgesOnCell % array (cell2)
- if ( grid % CellsOnCell % array (i,cell2) > 0) &
- d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- z_edge = 0.5*(zgrid(k,cell1) + zgrid(k,cell2)) &
- - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
- if (config_theta_adv_order == 3) then
- z_edge3 = - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 - d2fdx2_cell2) / 12.
- else
- z_edge3 = 0.
- end if
-
- end if
-
- zb(k,1,iEdge) = (z_edge-zgrid(k,cell1))*dvEdge(iEdge)/AreaCell(cell1)
- zb(k,2,iEdge) = (z_edge-zgrid(k,cell2))*dvEdge(iEdge)/AreaCell(cell2)
- zb3(k,1,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell1)
- zb3(k,2,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell2)
-
- if (k /= 1) then
- zf(k,1,iEdge) = ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb(k,1,iEdge)
- zf(k,2,iEdge) = ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb(k,2,iEdge)
- zf3(k,1,iEdge)= ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb3(k,1,iEdge)
- zf3(k,2,iEdge)= ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb3(k,2,iEdge)
- end if
-
- end do
-
- end if
- end do
-
diag % rw % array = 0.
state % w % array = 0.
do iEdge = 1,grid % nEdges
@@ -3278,7 +3303,9 @@
deallocate(vert_level)
+ end if ! config_met_interp
+
end subroutine nhyd_test_case_gfs
</font>
</pre>