[Dart-dev] [5046] DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90: Use the original profile when computing the line integral

[nancy at ucar.edu]

Revision: 5046
Author:   nancy
Date:     2011-06-30 10:01:06 -0600 (Thu, 30 Jun 2011)
Log Message:
-----------
Use the original profile when computing the line integral
for the non-local operator.  This code was mistakenly using
only the requested output observation heights as the profile.

Modified Paths:
--------------
    DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90

-------------- next part --------------
Modified: DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90
===================================================================
--- DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90	2011-06-29 19:14:31 UTC (rev 5045)
+++ DART/trunk/observations/gps/convert_cosmic_gps_cdf.f90	2011-06-30 16:01:06 UTC (rev 5046)
@@ -169,7 +169,7 @@
 
 end if
 
-allocate(hghtp(nlevels))  ;  allocate(refrp(nlevels))
+
 did_obs = .false.
 
 ! main loop that does either a single file or a list of files
@@ -206,6 +206,7 @@
    call nc_check( nf90_get_att(ncid,nf90_global,'rfict',rfict),'get_att rfict')
    rfict = rfict * 1000.0_r8
    
+   allocate(hghtp(nobs)) ;  allocate(refrp(nobs))
    allocate( lat(nobs))  ;  allocate( lon(nobs))
    allocate(hght(nobs))  ;  allocate(refr(nobs))
    allocate(azim(nobs))
@@ -235,16 +236,10 @@
    
    call nc_check( nf90_close(ncid) , 'close file')
    
-   obsloop: do k = 1, nlevels
-   
-     call interp_height_wght(hght, obs_levels(k), nobs, zloc, wght)
-     if ( zloc < 1 )  cycle obsloop
-     hghtp(nlevels-k+1) = obs_levels(k) * 1000.0_r8
-     refrp(nlevels-k+1) = exp( wght * log(refr(zloc)) + & 
-                       (1.0_r8 - wght) * log(refr(zloc+1)) ) * 1.0e-6_r8
-   
-   end do obsloop
-   
+   ! convert units here.
+   hghtp(:) = hght(:) * 1000.0_r8
+   refrp(:) = refr(:) * 1.0e-6_r8
+    
    first_obs = .true.
    
    obsloop2: do k = 1, nlevels
@@ -281,7 +276,7 @@
    
        !  compute the excess phase
        call excess(refrp, hghtp, lono, lato, hghto, nx, & 
-                   ny, nz, rfict, ray_ds, ray_htop, phs, nlevels)
+                   ny, nz, rfict, ray_ds, ray_htop, phs, nobs)
    
        !  if too high, phs will return as 0.  cycle loop here.
        if (phs <= 0) cycle obsloop2
@@ -316,7 +311,7 @@
    end do obsloop2
 
   ! clean up and loop if there is another input file
-  deallocate( lat, lon, hght, refr, azim )
+  deallocate( lat, lon, hght, refr, azim, hghtp, refrp)
 
   filenum = filenum + 1
 
@@ -664,14 +659,21 @@
 bot_lev = -1
 fract = 0.0_r8
 
+! make sure it's not higher than the highest available level.
+if (height > hghtp(1)) then
+   write(msgstring, *) 'requested height is ', height, '; highest available is ', hghtp(1)
+   call error_handler(E_ERR, 'bad level, above highest in file', &
+                      source, revision, revdate, text2=msgstring)
+endif
+
 !  Search down through height levels
-do k = 2, nobs
+heights: do k = 2, nobs
   if ( height >= hghtp(k) ) then
     bot_lev = k
     fract = (hghtp(k) - height) / (hghtp(k) - hghtp(k-1))
-    exit
+    exit heights
   endif
-end do
+end do heights
 
 ! the hghtp() array is currently an interpolated list of levels
 ! and on at least 1 PGI compiler version computing the lowest value 
@@ -680,12 +682,13 @@
 ! close to the lowest level then return it as equal; otherwise it's 
 ! an internally inconsistent input file.
 if (bot_lev < 0) then
-  if (abs(height - hghtp(nobs)) < 0.00001) then
+  if (abs(height - hghtp(nobs)) < 0.00001_r8) then
      bot_lev = nobs
      fract = 0.0_r8
   else
+     write(msgstring, *) 'requested height is ', height, '; lowest available is ', hghtp(nobs)
      call error_handler(E_ERR, 'bad level, below lowest in file', &
-                       source, revision, revdate)
+                       source, revision, revdate, text2=msgstring)
   endif
 endif