[Dart-dev] [5800] DART/branches/development/obs_def/obs_def_tower_mod.f90: First attempt at supporting the latent and sensible heat flux observations .
nancy at ucar.edu
nancy at ucar.edu
Wed Jul 18 14:57:02 MDT 2012
Revision: 5800
Author: thoar
Date: 2012-07-18 14:57:01 -0600 (Wed, 18 Jul 2012)
Log Message:
-----------
First attempt at supporting the latent and sensible heat flux observations.
Modified Paths:
--------------
DART/branches/development/obs_def/obs_def_tower_mod.f90
-------------- next part --------------
Modified: DART/branches/development/obs_def/obs_def_tower_mod.f90
===================================================================
--- DART/branches/development/obs_def/obs_def_tower_mod.f90 2012-07-18 20:15:46 UTC (rev 5799)
+++ DART/branches/development/obs_def/obs_def_tower_mod.f90 2012-07-18 20:57:01 UTC (rev 5800)
@@ -317,15 +317,154 @@
real(r8), intent(out) :: obs_val
integer, intent(out) :: istatus
+integer, dimension(NF90_MAX_VAR_DIMS) :: dimids
+real(r8), dimension(3) :: loc
+integer, dimension(3) :: ncstart, nccount
+integer, dimension(1) :: loninds, latinds, timeinds
+integer :: gridloni, gridlatj, timei
+integer :: varid, xtype, ndims, natts, dimlen
+integer :: io1, io2, second, day
+real(r8) :: loc_lon, loc_lat
+real(r4), dimension(1) :: hyperslab
+real(r4) :: spvalR4
+real(r8) :: scale_factor, add_offset
+real(digits12) :: otime
+character(len=20) :: strshort
+
if ( .not. module_initialized ) call initialize_module(state_time)
obs_val = MISSING_R8
istatus = 1
-call error_handler(E_ERR, 'get_expected_latent_heat_flux', &
- 'not implemented yet.', &
- source, revision, revdate)
+write(strshort,'(''ens_index '',i4)')ens_index
+if (ens_index > ens_size) then
+ write(string1,*)'believed to have ',ens_size,'ensemble members for observation operator.'
+ write(string2,*)'asking to use operator for ensemble member ',ens_index
+ call error_handler(E_ERR, 'get_expected_latent_heat_flux', &
+ string1, source, revision, revdate, text2=string2)
+endif
+
+! bombproofing ... make sure the netcdf file is open.
+
+write(*,*)'ncid(',ens_index,') is ',ncid(ens_index)
+call nc_check(nf90_inquire(ncid(ens_index)), &
+ 'get_expected_latent_heat_flux', 'inquire '//trim(strshort))
+
+! bombproofing ... make sure the variable is the shape and size we expect
+
+call nc_check(nf90_inq_varid(ncid(ens_index), 'EFLX_LH_TOT_R', varid), &
+ 'get_expected_latent_heat_flux', 'inq_varid EFLX_LH_TOT_R '//trim(strshort))
+call nc_check(nf90_inquire_variable(ncid(ens_index), varid, xtype=xtype, ndims=ndims, &
+ dimids=dimids, natts=natts),'get_expected_latent_heat_flux','inquire variable EFLX_LH_TOT_R '//trim(strshort))
+
+if (ndims /= 3) then
+ write(string1,*)'EFLX_LH_TOT_R is supposed to have 3 dimensions, it has',ndims
+ call error_handler(E_ERR, 'get_expected_latent_heat_flux', &
+ string1, source, revision, revdate)
+endif
+
+! If the variable is not a NF90_FLOAT, then the assumptions for processing
+! the missing_value, _FillValue, etc., may not be correct.
+if (xtype /= NF90_FLOAT) then
+ write(string1,*)'EFLX_LH_TOT_R supposed to be a 32 bit real. xtype = ',NF90_FLOAT,' it is ',xtype
+ call error_handler(E_ERR, 'get_expected_latent_heat_flux', &
+ string1, source, revision, revdate)
+endif
+
+! Dimension 1 is longitude
+call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(1), len=dimlen), &
+ 'get_expected_latent_heat_flux', 'inquire_dimension EFLX_LH_TOT_R 1'//trim(strshort))
+if (dimlen /= nlon) then
+ write(string1,*)'LON has length',nlon,'EFLX_LH_TOT_R has ',dimlen,'longitudes.'
+ call error_handler(E_ERR, 'get_expected_latent_heat_flux', &
+ string1, source, revision, revdate)
+endif
+
+! Dimension 2 is latitude
+call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(2), len=dimlen), &
+ 'get_expected_latent_heat_flux', 'inquire_dimension EFLX_LH_TOT_R 2'//trim(strshort))
+if (dimlen /= nlat) then
+ write(string1,*)'LAT has length',nlat,'EFLX_LH_TOT_R has ',dimlen,'latitudes.'
+ call error_handler(E_ERR, 'get_expected_latent_heat_flux', &
+ string1, source, revision, revdate)
+endif
+
+! Dimension 3 is time
+call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(3), len=dimlen), &
+ 'get_expected_latent_heat_flux', 'inquire_dimension EFLX_LH_TOT_R 3'//trim(strshort))
+if (dimlen /= ntime) then
+ write(string1,*)'TIME has length',ntime,'EFLX_LH_TOT_R has ',dimlen,'times.'
+ call error_handler(E_ERR, 'get_expected_latent_heat_flux', &
+ string1, source, revision, revdate)
+endif
+
+! Find the grid cell and timestep of interest
+! Get the individual locations values
+
+call get_time(obs_time, second, day)
+otime = real(day,digits12) + real(second,digits12)/86400.0_digits12
+loc = get_location(location) ! loc is in DEGREES
+loc_lon = loc(1)
+loc_lat = loc(2)
+
+latinds = minloc(abs(lat - loc_lat)) ! these return 'arrays' ...
+loninds = minloc(abs(lon - loc_lon)) ! these return 'arrays' ...
+timeinds = minloc(abs(rtime - otime)) ! these return 'arrays' ...
+
+gridlatj = latinds(1)
+gridloni = loninds(1)
+timei = timeinds(1)
+
+if (debug .and. do_output()) then
+ write(*,*)'get_expected_latent_heat_flux:targetlon, lon, lon index is ', &
+ loc_lon,lon(gridloni),gridloni
+ write(*,*)'get_expected_latent_heat_flux:targetlat, lat, lat index is ', &
+ loc_lat,lat(gridlatj),gridlatj
+ write(*,*)'get_expected_latent_heat_flux: targetT, T, T index is ', &
+ otime,rtime(timei),timei
+endif
+
+if ( abs(otime - rtime(timei)) > 30*60 ) then
+ if (debug .and. do_output()) then
+ write(*,*)'get_expected_latent_heat_flux: no close time ... skipping observation'
+ call print_time(obs_time,'get_expected_latent_heat_flux:observation time')
+ call print_date(obs_time,'get_expected_latent_heat_flux:observation date')
+ endif
+ istatus = 2
+ return
+endif
+
+! Grab exactly the scalar we want.
+
+ncstart = (/ gridloni, gridlatj, timei /)
+nccount = (/ 1, 1, 1 /)
+
+call nc_check(nf90_get_var(ncid(ens_index), varid, hyperslab, start=ncstart, count=nccount), &
+ 'get_expected_latent_heat_flux', 'get_var')
+
+obs_val = hyperslab(1)
+
+! Apply any netCDF attributes ...
+
+io1 = nf90_get_att(ncid(ens_index), varid, '_FillValue' , spvalR4)
+if ((io1 == NF90_NOERR) .and. (hyperslab(1) == spvalR4)) obs_val = MISSING_R8
+
+io2 = nf90_get_att(ncid(ens_index), varid, 'missing_value' , spvalR4)
+if ((io2 == NF90_NOERR) .and. (hyperslab(1) == spvalR4)) obs_val = MISSING_R8
+
+io1 = nf90_get_att(ncid(ens_index), varid, 'scale_factor', scale_factor)
+io2 = nf90_get_att(ncid(ens_index), varid, 'add_offset' , add_offset)
+
+if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then
+ if (obs_val /= MISSING_R8) obs_val = obs_val * scale_factor + add_offset
+elseif (io1 == NF90_NOERR) then
+ if (obs_val /= MISSING_R8) obs_val = obs_val * scale_factor
+elseif (io2 == NF90_NOERR) then
+ if (obs_val /= MISSING_R8) obs_val = obs_val + add_offset
+endif
+
+if (obs_val /= MISSING_R8) istatus = 0
end subroutine get_expected_latent_heat_flux
@@ -343,15 +482,155 @@
real(r8), intent(out) :: obs_val
integer, intent(out) :: istatus
+integer, dimension(NF90_MAX_VAR_DIMS) :: dimids
+real(r8), dimension(3) :: loc
+integer, dimension(3) :: ncstart, nccount
+integer, dimension(1) :: loninds, latinds, timeinds
+integer :: gridloni, gridlatj, timei
+integer :: varid, xtype, ndims, natts, dimlen
+integer :: io1, io2, second, day
+real(r8) :: loc_lon, loc_lat
+real(r4), dimension(1) :: hyperslab
+real(r4) :: spvalR4
+real(r8) :: scale_factor, add_offset
+real(digits12) :: otime
+character(len=20) :: strshort
+
if ( .not. module_initialized ) call initialize_module(state_time)
obs_val = MISSING_R8
istatus = 1
-call error_handler(E_ERR, 'get_expected_sensible_heat_flux', &
- 'not implemented yet.', &
- source, revision, revdate)
+write(strshort,'(''ens_index '',i4)')ens_index
+if (ens_index > ens_size) then
+ write(string1,*)'believed to have ',ens_size,'ensemble members for observation operator.'
+ write(string2,*)'asking to use operator for ensemble member ',ens_index
+ call error_handler(E_ERR, 'get_expected_sensible_heat_flux', &
+ string1, source, revision, revdate, text2=string2)
+endif
+
+! bombproofing ... make sure the netcdf file is open.
+
+write(*,*)'ncid(',ens_index,') is ',ncid(ens_index)
+call nc_check(nf90_inquire(ncid(ens_index)), &
+ 'get_expected_sensible_heat_flux', 'inquire '//trim(strshort))
+
+! bombproofing ... make sure the variable is the shape and size we expect
+
+call nc_check(nf90_inq_varid(ncid(ens_index), 'FSH', varid), &
+ 'get_expected_sensible_heat_flux', 'inq_varid FSH '//trim(strshort))
+call nc_check(nf90_inquire_variable(ncid(ens_index), varid, xtype=xtype, ndims=ndims, &
+ dimids=dimids, natts=natts),'get_expected_sensible_heat_flux','inquire variable FSH '//trim(strshort))
+
+if (ndims /= 3) then
+ write(string1,*)'FSH is supposed to have 3 dimensions, it has',ndims
+ call error_handler(E_ERR, 'get_expected_sensible_heat_flux', &
+ string1, source, revision, revdate)
+endif
+
+! If the variable is not a NF90_FLOAT, then the assumptions for processing
+! the missing_value, _FillValue, etc., may not be correct.
+if (xtype /= NF90_FLOAT) then
+ write(string1,*)'FSH supposed to be a 32 bit real. xtype = ',NF90_FLOAT,' it is ',xtype
+ call error_handler(E_ERR, 'get_expected_sensible_heat_flux', &
+ string1, source, revision, revdate)
+endif
+
+! Dimension 1 is longitude
+call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(1), len=dimlen), &
+ 'get_expected_sensible_heat_flux', 'inquire_dimension FSH 1'//trim(strshort))
+if (dimlen /= nlon) then
+ write(string1,*)'LON has length',nlon,'FSH has ',dimlen,'longitudes.'
+ call error_handler(E_ERR, 'get_expected_sensible_heat_flux', &
+ string1, source, revision, revdate)
+endif
+
+! Dimension 2 is latitude
+call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(2), len=dimlen), &
+ 'get_expected_sensible_heat_flux', 'inquire_dimension FSH 2'//trim(strshort))
+if (dimlen /= nlat) then
+ write(string1,*)'LAT has length',nlat,'FSH has ',dimlen,'latitudes.'
+ call error_handler(E_ERR, 'get_expected_sensible_heat_flux', &
+ string1, source, revision, revdate)
+endif
+
+! Dimension 3 is time
+call nc_check(nf90_inquire_dimension(ncid(ens_index), dimids(3), len=dimlen), &
+ 'get_expected_sensible_heat_flux', 'inquire_dimension FSH 3'//trim(strshort))
+if (dimlen /= ntime) then
+ write(string1,*)'TIME has length',ntime,'FSH has ',dimlen,'times.'
+ call error_handler(E_ERR, 'get_expected_sensible_heat_flux', &
+ string1, source, revision, revdate)
+endif
+
+! Find the grid cell and timestep of interest
+! Get the individual locations values
+
+call get_time(obs_time, second, day)
+otime = real(day,digits12) + real(second,digits12)/86400.0_digits12
+loc = get_location(location) ! loc is in DEGREES
+loc_lon = loc(1)
+loc_lat = loc(2)
+
+latinds = minloc(abs(lat - loc_lat)) ! these return 'arrays' ...
+loninds = minloc(abs(lon - loc_lon)) ! these return 'arrays' ...
+timeinds = minloc(abs(rtime - otime)) ! these return 'arrays' ...
+
+gridlatj = latinds(1)
+gridloni = loninds(1)
+timei = timeinds(1)
+
+if (debug .and. do_output()) then
+ write(*,*)'get_expected_sensible_heat_flux:targetlon, lon, lon index is ', &
+ loc_lon,lon(gridloni),gridloni
+ write(*,*)'get_expected_sensible_heat_flux:targetlat, lat, lat index is ', &
+ loc_lat,lat(gridlatj),gridlatj
+ write(*,*)'get_expected_sensible_heat_flux: targetT, T, T index is ', &
+ otime,rtime(timei),timei
+endif
+
+if ( abs(otime - rtime(timei)) > 30*60 ) then
+ if (debug .and. do_output()) then
+ write(*,*)'get_expected_sensible_heat_flux: no close time ... skipping observation'
+ call print_time(obs_time,'get_expected_sensible_heat_flux:observation time')
+ call print_date(obs_time,'get_expected_sensible_heat_flux:observation date')
+ endif
+ istatus = 2
+ return
+endif
+
+! Grab exactly the scalar we want.
+
+ncstart = (/ gridloni, gridlatj, timei /)
+nccount = (/ 1, 1, 1 /)
+
+call nc_check(nf90_get_var(ncid(ens_index), varid, hyperslab, start=ncstart, count=nccount), &
+ 'get_expected_sensible_heat_flux', 'get_var')
+
+obs_val = hyperslab(1)
+
+! Apply any netCDF attributes ...
+
+io1 = nf90_get_att(ncid(ens_index), varid, '_FillValue' , spvalR4)
+if ((io1 == NF90_NOERR) .and. (hyperslab(1) == spvalR4)) obs_val = MISSING_R8
+
+io2 = nf90_get_att(ncid(ens_index), varid, 'missing_value' , spvalR4)
+if ((io2 == NF90_NOERR) .and. (hyperslab(1) == spvalR4)) obs_val = MISSING_R8
+
+io1 = nf90_get_att(ncid(ens_index), varid, 'scale_factor', scale_factor)
+io2 = nf90_get_att(ncid(ens_index), varid, 'add_offset' , add_offset)
+
+if ( (io1 == NF90_NOERR) .and. (io2 == NF90_NOERR) ) then
+ if (obs_val /= MISSING_R8) obs_val = obs_val * scale_factor + add_offset
+elseif (io1 == NF90_NOERR) then
+ if (obs_val /= MISSING_R8) obs_val = obs_val * scale_factor
+elseif (io2 == NF90_NOERR) then
+ if (obs_val /= MISSING_R8) obs_val = obs_val + add_offset
+endif
+
+if (obs_val /= MISSING_R8) istatus = 0
+
end subroutine get_expected_sensible_heat_flux
More information about the Dart-dev
mailing list