[Dart-dev] DART/branches Revision: 12069

[dart at ucar.edu]

hendric at ucar.edu
2017-11-08 16:43:07 -0700 (Wed, 08 Nov 2017)
52
halfway through converting heights to model level.




Modified: DART/branches/recam/models/cam-fv/model_mod.f90
===================================================================
--- DART/branches/recam/models/cam-fv/model_mod.f90	2017-11-08 23:20:07 UTC (rev 12068)
+++ DART/branches/recam/models/cam-fv/model_mod.f90	2017-11-08 23:43:07 UTC (rev 12069)
@@ -634,6 +634,11 @@
    case(VERTISHEIGHT)
       ! construct a height column here and find the model levels
       ! that enclose this value
+      !@>todo put in arguments and write height_to_level
+      call cam_height_levels()
+      call height_to_level()
+      
+
       write(string1, *) 'we have not written the code yet for vertical type: ', which_vert
       call error_handler(E_ERR,'find_vertical_levels', &
                          string1,source,revision,revdate)
@@ -671,6 +676,70 @@
 end subroutine find_vertical_levels
 
 !-----------------------------------------------------------------------
+!> Compute the heights at pressure midpoints
+
+subroutine cam_height_levels(state_handle, ens_size, lon_index, lat_index, nlevels, height_array, my_status) 
+type(ensemble_type), intent(in)  :: state_handle
+integer,             intent(in)  :: ens_size
+real(r8),            intent(in)  :: lon_index
+real(r8),            intent(in)  :: lat_index
+integer,             intent(in)  :: nlevels
+real(r8),            intent(out) :: height_array(nlevels, ens_size)
+integer,             intent(out) :: my_status
+
+integer :: k, varid, level_one
+real(r8) :: temperature(ens_size), specific_humidity(ens_size), surface_pressure(ens_size)
+real(r8) :: phi(ens_size, nlevels)
+real(r8) :: surface_elevation
+
+!>@todo this should come from a model specific constant module.
+!> the forward operators and model_mod should use it.
+real(r8), parameter :: rd = 287.05_r8 ! dry air gas constant
+real(r8), parameter :: rv = 461.51_r8 ! wet air gas constant
+real(r8), parameter :: rr_factor = (rv/rd) - 1.0_r8
+
+! this is for surface obs
+level_one = 1
+
+!@>todo make into a subroutine
+! get the surface pressure from the state_handle
+varid = get_varid_from_kind(domain_id, QTY_SURFACE_PRESSURE)
+state_indx = get_dart_vector_index(lon_index, lat_index, level_one, domain_id, varid)
+surface_pressure(:) = get_state(state_indx, state_handle)
+
+! get the surface elevation from the phis
+surface_elevation = phis(lon_index, lat_index)
+
+do k = 1, nlevels
+   ! temperature
+   varid = get_varid_from_kind(domain_id, QTY_TEMPERATURE)
+   state_indx = get_dart_vector_index(lon_index, lat_index, k, domain_id, varid)
+   temperature(:) = get_state(state_indx, state_handle)
+
+   ! specific humidity
+   varid = get_varid_from_kind(domain_id, QTY_SPECIFIC_HUMIDITY)
+   state_indx = get_dart_vector_index(lon_index, lat_index, k, domain_id, varid)
+   specific_humidity(:) = get_state(state_indx, state_handle)
+   
+   !>@todo rename tv to something that mens something to users
+   tv(:,k) = temperature(:)*(1.0_r8 + rr_factor*specific_humidity(:))
+enddo
+
+! need to convert to geopotential height
+do imem = 1, ens_size
+   !>@todo refacfor to just put out geometric height
+   call dcz2(nlevels, surface_pressure(imem), surface_elevation, tv(imem,:), &
+             grid_data%P0%vals(1), hybrid_As, hybrid_Bs, phi(imem,:))
+   do k = 1,nlevels
+      height_array(k, imem) = gph2gmh(phi(imem,k), grid_data%lat%vals(lat_index))
+   enddo
+enddo
+
+istatus = 0
+
+end subroutine cam_height_levels
+
+!-----------------------------------------------------------------------
 !> Compute the pressures at the layer midpoints
 
 subroutine cam_pressure_levels(surface_pressure, n_levels, hyam, hybm, P0, pressure_array)
@@ -1509,6 +1578,168 @@
 
 !-----------------------------------------------------------------------
 
+subroutine dcz2(kmax,p_surf,h_surf,tv,hprb,hybrid_As,hybrid_Bs,z2)
+
+! Compute geopotential height for a CESM hybrid coordinate column.
+! All arrays except hybrid_As, hybrid_Bs are oriented top to bottom.
+! hybrid_[AB]s first subscript:
+!  = 1 for layer interfaces
+!  = 2 for layer midpoints
+! hybrid_As coord coeffs for P0 reference pressure term in plevs_cam
+! hybrid_Bs coord coeffs for surf pressure term in plevs_cam (in same format as hybrid_As)
+