[Dart-dev] DART/branches Revision: 12368

[dart at ucar.edu]

hendric at ucar.edu
2018-01-19 10:34:23 -0700 (Fri, 19 Jan 2018)
46

adding oxygen ion density forward operator.




Modified: DART/branches/recam/models/cam-fv/model_mod.f90
===================================================================
--- DART/branches/recam/models/cam-fv/model_mod.f90	2018-01-18 23:06:28 UTC (rev 12367)
+++ DART/branches/recam/models/cam-fv/model_mod.f90	2018-01-19 17:34:23 UTC (rev 12368)
@@ -51,9 +51,11 @@
                                    file_exist, to_upper, E_ERR, E_MSG
 use          obs_kind_mod,  only : QTY_SURFACE_ELEVATION, QTY_PRESSURE, &
                                    QTY_GEOMETRIC_HEIGHT, QTY_VERTLEVEL, &
-                                   QTY_SURFACE_PRESSURE, QTY_PRECIPITABLE_WATER, &
+                                   QTY_SURFACE_PRESSURE, &
                                    QTY_TEMPERATURE, QTY_SPECIFIC_HUMIDITY, &
-                                   QTY_GEOPOTENTIAL_HEIGHT,  &
+                                   QTY_MOLEC_OXYGEN_MIXING_RATIO, &
+                                   QTY_ION_O_MIXING_RATIO, QTY_ATOMIC_H_MIXING_RATIO, &
+                                   QTY_ATOMIC_OXYGEN_MIXING_RATIO, QTY_NITROGEN, &
                                    get_index_for_quantity, get_num_quantities, &
                                    get_name_for_quantity, get_quantity_for_type_of_obs
 use     mpi_utilities_mod,  only : my_task_id
@@ -137,6 +139,7 @@
 character(len=32)  :: fields_to_perturb(MAX_PERT)     = "QTY_TEMPERATURE"
 real(r8)           :: perturbation_amplitude(MAX_PERT)= 0.00001_r8
 logical            :: using_chemistry                 = .false.
+logical            :: use_variable_mean_mass          = .false.
 
 ! state_variables defines the contents of the state vector.
 ! each line of this input should have the form:
@@ -166,6 +169,8 @@
    custom_routine_to_generate_ensemble, &
    fields_to_perturb,               &
    perturbation_amplitude,          &
+   use_variable_mean_mass,          &
+   using_chemistry,                   &
    debug_level
 
 ! global variables
@@ -837,8 +842,8 @@
 istatus(:)     = 99
 
 interp_handle = get_interp_handle(obs_qty)
-lon_lat_vert = get_location(location)
-which_vert = nint(query_location(location)) 
+lon_lat_vert  = get_location(location)
+which_vert    = nint(query_location(location)) 
 
 call quad_lon_lat_locate(interp_handle, lon_lat_vert(1), lon_lat_vert(2), &
                          four_lons, four_lats, lon_fract, lat_fract, istatus(1))
@@ -965,7 +970,6 @@
                                 varid, quad_vals, my_status)
 
    else ! get 3d special variables in another ways ( like QTY_PRESSURE )
-
       call get_four_nonstate_values(state_handle, ens_size, four_lons, four_lats, &
                                    four_bot_levs, four_top_levs, four_vert_fracts, &
                                    obs_qty, quad_vals, my_status)
@@ -1412,15 +1416,9 @@
 
 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) :: surface_pressure(ens_size), mbar(ens_size, nlevels)
 real(r8) :: tv(nlevels, ens_size)  ! Virtual temperature, top to bottom
 
-!>@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
 
@@ -1436,37 +1434,34 @@
 ! DEBUG
 !print*, 'lon lat surf elev ', lon_index, lat_index, surface_elevation
 
-! construct a virtual temperature column, one for each ensemble member
-do k = 1, nlevels
-   ! temperature
-   call get_staggered_values_from_qty(ens_handle, ens_size, QTY_TEMPERATURE, &
-                                     lon_index, lat_index, k, qty, temperature, status1)
+call compute_virtual_temperature(ens_handle, ens_size, lon_index, lat_index, nlevels, qty, tv, status1)
 
-   ! specific humidity
-   call get_staggered_values_from_qty(ens_handle, ens_size, QTY_SPECIFIC_HUMIDITY, &
-                                     lon_index, lat_index, k, qty, specific_humidity, status1)
-   !>tv == virtual temperature.
-   tv(k,:) = temperature(:)*(1.0_r8 + rr_factor*specific_humidity(:))
-enddo
+if (status1 /= 0) then
+   my_status = status1
+   return
+endif
 
+if (use_variable_mean_mass) then
+  call compute_mean_mass(ens_handle, ens_size, lon_index, lat_index, nlevels, qty, mbar, status1)
+
+  if (status1 /= 0) then
+     my_status = status1
+     return
+  endif
+else