[Dart-dev] DART/branches Revision: 11727

[dart at ucar.edu]

nancy at ucar.edu
2017-06-13 17:36:49 -0600 (Tue, 13 Jun 2017)
76
down to about 2 pages of errors, mostly of the 
array shape-matching kind.




Modified: DART/branches/coamps/models/coamps_nest/coamps_interp_mod.f90
===================================================================
--- DART/branches/coamps/models/coamps_nest/coamps_interp_mod.f90	2017-06-13 23:05:35 UTC (rev 11726)
+++ DART/branches/coamps/models/coamps_nest/coamps_interp_mod.f90	2017-06-13 23:36:49 UTC (rev 11727)
@@ -157,8 +157,9 @@
     ! I can keep them straight: 
     !  "values" arrays index as (neighbor, level)
     !  availability array index as (level, variable)
-    integer, parameter :: VALUES_DIM_NEIGHBOR       = 1
-    integer, parameter :: VALUES_DIM_LEVEL          = 2
+    integer, parameter :: VALUES_DIM_ENSEMBLE       = 1
+    integer, parameter :: VALUES_DIM_NEIGHBOR       = 2
+    integer, parameter :: VALUES_DIM_LEVEL          = 3
     integer, parameter :: AVAILABILITY_DIM_LEVEL    = 1
     integer, parameter :: AVAILABILITY_DIM_VARIABLE = 2
   
@@ -728,13 +729,15 @@
         select case (obs_kind)
           case (QTY_GEOPOTENTIAL_HEIGHT)
               call calculate_heights(interpolator, AVAILABLE_INDEX_TARGET,      &
-                                     interpolator%target_values)
+                                     interpolator%target_values(1,:,:))
+              interpolator%target_values = spread(interpolator%target_values(1,:,:), &
+                                     VALUES_DIM_ENSEMBLE, interpolator%ensemble_size)
           case (QTY_VERTLEVEL)
              interpolator%target_values(1,:,:)  =                      &
                 spread(get_domain_msigma(interpolator%model_domain),   &
                        VALUES_DIM_NEIGHBOR, NUM_NEIGHBORS)
              interpolator%target_values = spread(interpolator%target_values(1,:,:), &
-                                                 1, interpolator%ensemble_size)
+                                           VALUES_DIM_ENSEMBLE, interpolator%ensemble_size)
 
              interpolator%vars_available(:, AVAILABLE_INDEX_TARGET) = .true.
 
@@ -769,7 +772,9 @@
             AVAILABLE_INDEX_HEIGHT = get_next_availability_index()
 
             call calculate_heights(interpolator, AVAILABLE_INDEX_HEIGHT,     &
-                                   interpolator%vcoord_values)
+                                   interpolator%vcoord_values(1,:,:))
+            interpolator%vcoord_values = spread(interpolator%vcoord_values(1,:,:), &
+                          VALUES_DIM_ENSEMBLE, interpolator%ensemble_size)
 
         case (INTERPOLATE_TO_SIGMA)
             AVAILABLE_INDEX_SIGMA = get_next_availability_index()
@@ -796,9 +801,10 @@
 
         character(len=*), parameter :: routine = 'calculate_surface_heights'
 
+        ! all ensemble_members have the same terrain
         call get_terrain_height_at_points(get_nest(interpolator%interp_point),                      &
                                           interpolator%neighbors_i, interpolator%neighbors_j,       &
-                                          interpolator%vinterp_values(:, SINGLE_LEVEL) )
+                                          interpolator%vinterp_values(1, :, SINGLE_LEVEL) )
     end subroutine calculate_surface_heights
 
     ! calculate_surface_pressure
@@ -811,10 +817,10 @@
         real(kind=r8), dimension(NUM_NEIGHBORS) :: zsfc
         real(kind=r8), dimension(SINGLE_POINT, SINGLE_POINT) :: sfc_pres 
 
-        real(kind=r8), dimension(:,:), pointer :: mean_exner_values
-        real(kind=r8), dimension(:,:), pointer :: mean_theta_values
-        real(kind=r8), dimension(:,:), pointer :: theta_values
-        real(kind=r8), dimension(:,:), pointer :: exner_values
+        real(kind=r8), dimension(:,:,:), pointer :: mean_exner_values
+        real(kind=r8), dimension(:,:,:), pointer :: mean_theta_values
+        real(kind=r8), dimension(:,:,:), pointer :: theta_values
+        real(kind=r8), dimension(:,:,:), pointer :: exner_values
 
         logical, parameter :: IS_NOT_MEAN  = .false.
         logical, parameter :: IS_MEAN      = .true.
@@ -821,7 +827,7 @@
         logical, parameter :: IS_M_LEVEL   = .true.
         logical, parameter :: IS_W_LEVEL   = .false.
 
-        integer :: n, k
+        integer :: n, k, e
 
         character(len=*), parameter :: routine = 'calculate_surface_pressure'
         integer                     :: alloc_status
@@ -857,19 +863,22 @@
         call get_terrain_height_at_points(get_nest(interpolator%interp_point),                      &
                                           interpolator%neighbors_i, interpolator%neighbors_j, zsfc) 
 
+        ! FIXME: can we skip the outer loop?
+        do e=1, interpolator%ensemble_size
         do n=1,NUM_NEIGHBORS
 
-          call sfcp(theta_values(n,:), exner_values(n,:),              &
-                    mean_theta_values(n,:), mean_exner_values(n,:),    &
-                    get_domain_dsigmaw(interpolator%model_domain),     &
-                    get_domain_wsigma(interpolator%model_domain),      &
-                    zsfc(n), SINGLE_POINT, SINGLE_POINT,               &
+          call sfcp(theta_values(e,n,:), exner_values(e,n,:),           &
+                    mean_theta_values(e,n,:), mean_exner_values(e,n,:), &
+                    get_domain_dsigmaw(interpolator%model_domain),      &
+                    get_domain_wsigma(interpolator%model_domain),       &
+                    zsfc(n), SINGLE_POINT, SINGLE_POINT,                &
                     num_model_levels, sfc_pres)