[Dart-dev] DART/branches Revision: 12173

[dart at ucar.edu]

nancy at ucar.edu
2017-12-04 17:07:06 -0700 (Mon, 04 Dec 2017)
96
first stab at converting pressures to heights
for interpolating obs with a vertical of height.




Modified: DART/branches/recam/models/cam-fv/model_mod.f90
===================================================================
--- DART/branches/recam/models/cam-fv/model_mod.f90	2017-12-04 18:32:21 UTC (rev 12172)
+++ DART/branches/recam/models/cam-fv/model_mod.f90	2017-12-05 00:07:06 UTC (rev 12173)
@@ -1106,8 +1106,8 @@
 integer  :: k, level_one, imember, status1
 real(r8) :: surface_elevation(1)
 real(r8) :: temperature(ens_size), specific_humidity(ens_size), surface_pressure(ens_size)
-real(r8) :: phi(ens_size, nlevels)
-real(r8) :: tv(ens_size, nlevels+1)  !>@todo FIXME  ???? ! Virtual temperature, top to bottom
+real(r8) :: tv(nlevels, ens_size)  ! Virtual temperature, top to bottom
+real(r8) :: height_interf(nlevels+1, ens_size)
 
 !>@todo this should come from a model specific constant module.
 !> the forward operators and model_mod should use it.
@@ -1115,9 +1115,6 @@
 real(r8), parameter :: rv = 461.51_r8 ! wet air gas constant
 real(r8), parameter :: rr_factor = (rv/rd) - 1.0_r8
 
-!>@todo FIXME: these need to be replaced by hyam, hymb, hyai, hybi -- VERY VERY CAREFULLY
-real(r8) ::hybrid_As(nlevels+1,2), hybrid_Bs(nlevels+1,2)
-
 ! this is for surface obs
 level_one = 1
 
@@ -1140,24 +1137,22 @@
    call get_staggered_values_from_qty(ens_handle, ens_size, QTY_SPECIFIC_HUMIDITY, &
                                      lon_index, lat_index, k, qty, specific_humidity, status1)
    
-   !>@todo rename tv to something that means something to users
    !>tv == virtual temperature.
-   tv(:,k) = temperature(:)*(1.0_r8 + rr_factor*specific_humidity(:))
+   tv(k,:) = temperature(:)*(1.0_r8 + rr_factor*specific_humidity(:))
 
-   print *, 'member 1, level, t, q, tv: ', k, temperature(1), specific_humidity(1), tv(1, k)
+   print *, 'member 1, level, t, q, tv: ', k, temperature(1), specific_humidity(1), tv(k, 1)
 
 enddo
 
-! need to convert to geopotential height
+! compute the height columns for each ensemble member
 do imember = 1, ens_size
-   !>@todo refactor to just put out geometric height
-   call dcz2(nlevels, surface_pressure(imember), surface_elevation(1), tv(imember,:), &
-             grid_data%P0%vals(1), hybrid_As, hybrid_Bs, phi(imember,:))
-   do k = 1,nlevels
-      height_array(k, imember) = gph2gmh(phi(imember,k), grid_data%lat%vals(lat_index))
-   enddo
+   call build_heights(nlevels, surface_pressure(imember), surface_elevation(1), tv(:, imember), &
+                      height_array(:, imember), height_interf(:, imember))  ! can pass in variable_r
 enddo
 
+! convert entire array to geometric height (from potential height)
+call gph2gmh(height_array, grid_data%lat%vals(lat_index))
+
 do k = 1,nlevels
    print *, "member 1, level, height: ", k, height_array(k, 1)
 enddo
@@ -1192,8 +1187,8 @@
 
 
 !-----------------------------------------------------------------------
-!> Compute columns of pressures at the layer midpoints for the given number 
-!> of surface pressures. 
+!> Compute column of pressures at the layer midpoints for the given 
+!> surface pressure. 
 !>
 !>@todo FIXME unlike some other things - you could pass in an upper and lower
 !>level number and compute the pressure only at the levels between those.
@@ -1200,23 +1195,20 @@
 !>this isn't a column that has to be built from the surface up or top down.
 !>each individual pressure can be computed independently given the surface pressure.
 
-subroutine cam_p_col_midpts(num_cols, surface_pressure, n_levels, pressure_array)
+subroutine cam_p_col_midpts(surface_pressure, n_levels, pressure_array)
 
-integer,            intent(in)  :: num_cols
-real(r8),           intent(in)  :: surface_pressure(num_cols)   ! in pascals
+real(r8),           intent(in)  :: surface_pressure   ! in pascals
 integer,            intent(in)  :: n_levels
-real(r8),           intent(out) :: pressure_array(n_levels, num_cols)
+real(r8),           intent(out) :: pressure_array(n_levels)
 
-integer :: j, k
+integer :: k
 
 ! Set midpoint pressures.  This array mirrors the order of the
 ! cam model levels: 1 is the model top, N is the bottom.
 
-do j=1, num_cols
-   do k=1, n_levels
-      pressure_array(k, j) = grid_data%hyam%vals(k) * grid_data%P0%vals(1) + &
-                             grid_data%hybm%vals(k) * surface_pressure(j)
-   enddo
+do k=1, n_levels
+   pressure_array(k) = grid_data%hyam%vals(k) * grid_data%P0%vals(1) + &
+                       grid_data%hybm%vals(k) * surface_pressure
 enddo
 
 end subroutine cam_p_col_midpts
@@ -1227,23 +1219,20 @@