| blob | c52dc7dee8cd0b6626ded285a69cd998d0606d24 |
1 ! Program Name:
2 ! Author(s)/Contact(s):
3 ! Abstract:
4 ! History Log:
5 !
6 ! Usage:
7 ! Parameters: <Specify typical arguments passed>
8 ! Input Files:
9 ! <list file names and briefly describe the data they include>
10 ! Output Files:
11 ! <list file names and briefly describe the information they include>
12 !
13 ! Condition codes:
14 ! <list exit condition or error codes returned >
15 ! If appropriate, descriptive troubleshooting instructions or
16 ! likely causes for failures could be mentioned here with the
17 ! appropriate error code
18 !
19 ! User controllable options: <if applicable>
21 module module_stream_nudging
23 use config_base, only: nlst
24 use module_nudging_io, only: lastObsStructure, lastObsStructure_SoA
25 #ifdef MPP_LAND
26 use module_mpp_land
27 use module_mpp_reachls, only: ReachLS_write_io
28 #endif
29 use module_hydro_stop, only:HYDRO_stop
31 implicit none
32 !===================================================================================================
33 ! Module Variables
35 !========================
36 ! obs and obsTime data structures. Each entry in obsTime holds a timeslice file, with the
37 ! individual obs in the obsStr.
38 type obsStructure
39 character(len=15) :: usgsId ! the 15 char USGS identifier.
40 character(len=19) :: obsTime ! observation at gage dims: nGages
41 real :: obsDischarge ! observation at gage dims: nGages
42 real :: obsQC ! quality control factpr [0,1]
43 integer :: obsStaticInd ! the index to the obsStaticStr where static info is kept
44 real :: innov ! obs-modeled
45 end type obsStructure
47 type obsTimeStructure
48 character(len=19) :: time
49 character(len=19) :: updateTime
50 integer, allocatable, dimension(:) :: allCellInds ! cell indices affected at this time
51 integer, allocatable, dimension(:) :: nGageCell ! number of gages for each affected cell ind
52 type(obsStructure), allocatable, dimension(:) :: obsStr ! the obs at this file time / timeslice
53 end type obsTimeStructure
55 ! The top level structure used to solve the nudges
56 type(obsTimeStructure), allocatable, dimension(:) :: obsTimeStr ! size=nObsTimes
58 !========================
59 ! lastObs structure, corresponding to nudgingLastObs.YYYY-mm-dd_HH:MM:ss.nc
60 ! How observations from the past are carried forward.
61 ! Type defined in module_nudging_io.F
62 type(lastObsStructure), allocatable, dimension(:) :: lastObsStr
63 type(lastObsStructure_SoA) :: lastObsStr_SoA
65 !========================
66 ! This holds static information for a given gage for a given cycle (when R, G, and tau) do not
67 ! change.
68 ! The "lastObs" variables are not exactly static... they were an afterthought and this was
69 ! by far the best place for them to live.
70 type, extends(lastObsStructure) :: obsStaticStructure
71 !! Inherited components
72 !!character(len=15) :: usgsId ! the 15 char USGS identifier.
73 !!real :: lastObsDischarge(:) ! last observed discharge
74 !!character(len=19) :: lastObsTime(:) ! time of last obs discharge (.le.hydroTime)
75 !!real :: lastObsQuality(:) ! quality of the last obs discharge
76 !!real :: lastObsModelDischarge(:) ! the modeled discharge value at the time of the last obs
77 !! New components
78 integer, allocatable, dimension(:) :: cellsAffected ! indices of cells affected
79 real, allocatable, dimension(:) :: dist ! optional: dist to affected cells, optional
80 real, allocatable, dimension(:) :: ws ! spatial cressman weights at affected cells discharge.
81 end type obsStaticStructure
83 type, extends(lastObsStructure_SoA) :: obsStaticStructure_SoA
84 !! Inherited components
85 !!character(len=15) :: usgsId ! the 15 char USGS identifier.
86 !!real :: lastObsDischarge(:) ! last observed discharge
87 !!character(len=19) :: lastObsTime(:) ! time of last obs discharge (.le.hydroTime)
88 !!real :: lastObsQuality(:) ! quality of the last obs discharge
89 !!real :: lastObsModelDischarge(:) ! the modeled discharge value at the time of the last obs
90 !! New components
91 integer, allocatable :: obsCellInd(:) ! index of obs on model channel network, for distance calc
92 real,allocatable :: R(:), G(:), tau(:) ! the nudging parameters at this gage.
93 end type obsStaticStructure_SoA
95 ! The static obs/gage information - store here to perform calculations
96 ! only once per cycle.
97 ! Currently the dimensions of this are fixed to 10000 which should
98 ! work for the forseeable future. May want to consider some routine for
99 ! augmenting this size if necessary.
100 type(obsStaticStructure), allocatable, dimension(:) :: obsStaticStr
101 type(obsStaticStructure_SoA) :: obsStaticStr_Soa
103 !! The number of up/down stream links which can be collected
104 !! using R to solve which links are "neighboring".
105 !! This value applies to both directions, it's the number
106 !! of links per R.
107 integer, parameter :: maxNeighLinksDim=5000
109 !========================
110 ! Node and gage collocation - For reach based routing, corresponds to the "gage" column
111 ! of RouteLink.nc.
112 type nodeGageStructure
113 integer, allocatable, dimension(:) :: nodeId
114 character(len=15), allocatable, dimension(:) :: usgsId
115 end type nodeGageStructure
116 type(nodeGageStructure) :: nodeGageTmp, nodeGageStr
117 integer :: nGagesDomain ! the number of gages specified
119 !========================
120 ! Nudging parameters structure, corresponding to NudgeParams.nc file.
121 !! for dealloction purposes, might be better to put the dimension onthe derived type.
122 type nudgingParamStructure
123 character(len=15), allocatable, dimension(:) :: usgsId
124 real, allocatable, dimension(:) :: R
125 real, allocatable, dimension(:) :: G
126 real, allocatable, dimension(:) :: tau
127 real, allocatable, dimension(:,:,:) :: qThresh !! gage, month, nThresh
128 real, allocatable, dimension(:,:,:) :: expCoeff !! gage, month, nThresh
129 end type nudgingParamStructure
130 type(nudgingParamStructure) :: nudgingParamsTmp, nudgingParamsStr
132 !========================
133 ! Network reExpression structure, corresponding to netwkReExFile.nc
134 type netwkReExpStructure
135 integer*4, allocatable, dimension(:) :: go
136 integer*4, allocatable, dimension(:) :: start
137 integer*4, allocatable, dimension(:) :: end
138 end type netwkReExpStructure
139 type(netwkReExpStructure) :: downNetwkStr, upNetwkStr
141 !========================
142 ! Track gages from NWIS not in our param file and not in the intersection or
143 ! the Route_Link gages and the gages in the parameter file.
144 integer, parameter :: maxNwisNotRLAndParamsCount=20000
145 character(len=15), dimension(maxNwisNotRLAndParamsCount) :: nwisNotRLAndParams
146 integer :: nwisNotRLAndParamsCount
148 !========================
149 ! Random data
150 real, allocatable, dimension(:) :: t0Discharge
152 !========================
153 ! Control book keeping
154 logical :: nudgeThisLsmTimeStep
156 !! JLM: these options are hardcoded/hardwired for now, evnt go in namelist
157 integer, parameter :: obsResolutionInt = 15 ! minutes
158 real, parameter :: obsCheckFreq = 90 ! minutes
159 logical, parameter :: filterObsToDom = .TRUE.
160 character(len=15), parameter :: missingGage = ' '
161 logical, parameter :: nudgeWAdvance = .FALSE.
162 character(len=4), parameter :: nudgingClockType = 'wall'
163 logical, parameter :: sanityQcDischarge = .true.
164 logical, parameter :: readTimesliceErrFatal = .false.
165 logical, parameter :: futureLastObsFatal = .true.
166 logical, parameter :: futureLastObsZero = .false.
167 real, parameter :: invDistTimeWeightExp = 5.000
168 real, parameter :: noConstInterfCoeff = 1.0 !0.500
170 ! hydro.namelist: &NUDGING_nlist variables.
171 character(len=256) :: nudgingParamFile
172 character(len=256) :: netwkReExFile
173 character(len=256) :: nudgingLastObsFile !! passed by namelist
174 character(len=256) :: nudgingLastObsFileTry !! either the passed or, if not passed, a default
175 logical :: readTimesliceParallel
176 logical :: temporalPersistence
177 logical :: persistBias
178 logical :: biasWindowBeforeT0
179 integer :: nTimesLastObs !! aka nlst_rt(did)%nLastObs
180 integer :: minNumPairsBiasPersist
181 integer :: maxAgePairsBiasPersist
182 logical :: invDistTimeWeightBias
183 logical :: noConstInterfBias
185 !========================
186 ! Space book keeping
187 logical :: nudgeSpatial = .true. !! is spatial interpolation of nudging active?
188 integer, parameter :: did=1 !! 1 for WRF-uncoupled runs...
190 !========================
191 ! Time book keeping
192 real :: maxTau
193 integer :: nObsTimes
194 character(len=19) :: lsmTime, initTime
195 integer :: lsmDt !1234567890123456789
196 character(len=19), parameter :: missingLastObsTime='9999999999999999999'
197 character(len=2) :: obsResolution
198 logical :: gotT0Discharge
200 #ifdef HYDRO_D
201 integer, parameter :: flushUnit=6
202 logical, parameter :: flushAll=.true.
203 #endif
205 #ifdef MPP_LAND
206 !========================
207 ! Parallel book keeping
208 real, allocatable, dimension(:) :: chanlen_image0 !A global version kept on image 0
209 #endif
212 contains
214 !===================================================================================================
215 ! Program Names:
216 ! init_stream_nudging_clock
217 ! Author(s)/Contact(s):
218 ! James L McCreight <jamesmcc><ucar><edu>
219 ! Abstract:
220 ! One-time initialization of diagnostic stream nuding clock.
221 ! History Log:
222 ! 10/08/15 -Created, JLM.
223 ! Usage:
224 ! Parameters:
225 ! Input Files:
226 ! Output Files: None.
227 ! Condition codes: this only gets called if #ifdef HYDRO_D?
228 ! User controllable options: None.
229 ! Notes:
231 subroutine init_stream_nudging_clock
232 use module_nudging_utils, only: totalNudgeTime, &
233 sysClockCountRate, &
234 sysClockCountMax, &
235 clockType
237 implicit none
238 totalNudgeTime = 0. ! Nudging time accumulation init.
239 call system_clock(count_rate=sysClockCountRate, count_max=sysClockCountMax)
240 clockType = trim(nudgingClockType)
241 !!$#ifdef HYDRO_D /* un-ifdef to get timing results */
242 !!$print*,'Ndg: totalNudgeTime: ', totalNudgeTime
243 !!$print*,'Ndg: sysClockCountRate: ', sysClockCountRate
244 !!$print*,'Ndg: sysClockCountMax: ', sysClockCountMax
245 !!$print*,'Ndg: clockType: ', trim(nudgingClockType)
246 !!$if(flushAll) flush(flushUnit)
247 !!$#endif /* HYDRO_D un-ifdef to get timing results */
248 end subroutine init_stream_nudging_clock
251 !===================================================================================================
252 ! Program Names:
253 ! init_stream_nudging
254 ! Author(s)/Contact(s):
255 ! James L McCreight <jamesmcc><ucar><edu>
256 ! Abstract:
257 ! One-time initialization of certain stream nuding information. Some of this infomation
258 ! may be updated later in the run, but probably not frequently.
259 ! History Log:
260 ! 7/23/15 -Created, JLM.
261 ! Usage:
262 ! Parameters:
263 ! Input Files: currently hardwired module variables
264 ! nudgingParamFile
265 ! gageGageDistFile
266 ! netwkReExFile
267 ! Output Files: None.
268 ! Condition codes:
269 ! User controllable options: None.
270 ! Notes:
272 subroutine init_stream_nudging
274 use module_RT_data, only: rt_domain
275 use module_nudging_utils, only: whichInLoop2, &
276 nudging_timer, &
277 accum_nudging_time
278 use module_nudging_io, only: get_netcdf_dim, &
279 read_gridded_nudging_frxst_gage_csv, &
280 read_reach_gage_collocation, &
281 read_nudging_param_file, &
282 read_network_reexpression, &
283 find_nudging_last_obs_file, &
284 read_nudging_last_obs, &
285 read_nudging_last_obs_soa
286 use module_mpp_reachls, only: ReachLS_decomp
288 implicit none
290 integer :: nLinks, nLinksL
291 !integer, dimension(nlinks) :: strmFrxstPts
292 integer :: did=1
293 integer :: ii, kk, ll, tt, count, toSize, fromSize
294 integer :: nParamGages, nParamMonth, nParamThresh, nParamThreshCat, nGgDists, nGgDistsKeep
295 integer :: nWhGageDists1, nWhGageDists2, nGagesWParamsDom
296 integer, allocatable, dimension(:) :: whParamsDom
297 real, allocatable, dimension(:) :: g_nudge
298 integer :: downSize, upSize, baseSize, nStnLastObs
299 character(len=256) :: lastObsFile !! confirms existence of a file
300 logical :: lastObsFileFound
301 !!$#ifdef HYDRO_D /* un-ifdef to get timing results */
302 !!$real :: startCodeTimeAcc, endCodeTimeAcc
303 !!$
304 !!$#ifdef MPP_LAND
305 !!$if(my_id .eq. io_id) &
306 !!$#endif
307 !!$ call nudging_timer(startCodeTimeAcc)
308 !!$#endif /* HYDRO_D */ /* un-ifdef to get timing results */
310 !!$#ifdef HYDRO_D
311 !!$print*, "Ndg: Start init_stream_nudging"
312 !!$#ifdef MPP_LAND
313 !!$print*, 'Ndg: PARALLEL NUDGING!'
314 !!$#endif
315 !!$if(flushAll) flush(flushUnit)
316 !!$#endif /* HYDRO_D */
318 !! this ifdef is bizarre problem with pgf compiler
319 !#ifdef WRF_HYDRO_NUDGING !JLM
320 nudgingParamFile = nlst(did)%nudgingParamFile
321 netwkReExFile = nlst(did)%netwkReExFile
322 readTimesliceParallel = nlst(did)%readTimesliceParallel
323 temporalPersistence = nlst(did)%temporalPersistence
324 persistBias = nlst(did)%persistBias
325 biasWindowBeforeT0 = nlst(did)%biasWindowBeforeT0
326 nudgingLastObsFile = nlst(did)%nudgingLastObsFile
327 nTimesLastObs = nlst(did)%nLastObs
328 minNumPairsBiasPersist = nlst(did)%minNumPairsBiasPersist
329 maxAgePairsBiasPersist = nlst(did)%maxAgePairsBiasPersist
330 invDistTimeWeightBias = nlst(did)%invDistTimeWeightBias
331 noConstInterfBias = nlst(did)%noConstInterfBias
332 !#endif
334 nLinks = RT_DOMAIN(did)%NLINKS ! For gridded channel routing
335 #ifdef MPP_LAND
336 nLinksL = RT_DOMAIN(did)%gNLINKSL ! For reach-based routing in parallel, no decomp for nudging
337 #else
338 nLinksL = RT_DOMAIN(did)%NLINKSL ! For reach-based routing
339 #endif
341 !Variable init
342 nwisNotRLAndParams = missingGage
343 nwisNotRLAndParamsCount=0
344 gotT0Discharge = .false.
346 !=================================================
347 ! 0. This routine is called with condition that chanrtswcrt.ne.0 in module_HYDRO_drv
349 !=================================================
350 ! 1. Gage link/node collocation.
351 ! which gages are actually in the domain?
352 ! Musk routines: use routelink csv! New column for associated gage Id for each link.
353 ! Grid channel : frxst_pts layer in Fulldom, and a new Nudge_frxst_gage.csv file, which
354 ! simply contains frxst_pts index and the associated gage ID, *** may be blank?***
355 !========================
356 !!$! Gridded channel routing is option 3
357 !!$if (nlst_rt(did)%channel_option .eq. 3) then
358 !!$ strmfrxstpts = RT_DOMAIN(did)%STRMFRXSTPTS
359 !!$ ! For gridded channel routing, setup the relationship between frxst_pts and gage IDs via
360 !!$ ! the Nudging_frxst_gage.csv.
361 !!$ ! For now this is a csv, but it should be netcdf in the long run.
362 !!$#ifdef HYDRO_D
363 !!$ print*, 'Ndg: Start initializing Nudging_frxst_gage.csv'
364 !!$#endif
365 !!$ !! allocate the maximum number of gages, say 8000.
366 !!$ allocate(nodeGageTmp%nodeId(maxGages), nodeGageTmp%usgsId(maxGages))
367 !!$ ! This actually returns frxst point from the file...
368 !!$ call read_gridded_nudging_frxst_gage_csv(nodeGageTmp%nodeId, nodeGageTmp%usgsId, nGagesDomain)
369 !!$ ! ... we'll convert this to index on the stream network.
370 !!$ allocate(nodeGageStr%nodeId(nGagesDomain), nodeGageStr%usgsId(nGagesDomain))
371 !!$
372 !!$ !! JLM: need to handle desired frxst points which are not gages.
373 !!$ !! Need to make sure these are a complete set 1:nGages.
374 !!$ do ll=1,nLinks
375 !!$ if(strmFrxstPts(ll) .ne. -9999) then
376 !!$ nodeGageStr%nodeId(strmFrxstPts(ll)) = ll
377 !!$ nodeGageStr%usgsId(strmFrxstPts(ll)) = adjustr(nodeGageTmp%usgsId(strmFrxstPts(ll)))
378 !!$ end if
379 !!$ end do
380 !!$
381 !!$ deallocate(nodeGageTmp%nodeId, nodeGageTmp%usgsId)
382 !!$
383 !!$#ifdef HYDRO_D
384 !!$ print*,nGagesDomain
385 !!$ print*,nodeGageStr%nodeId
386 !!$ print*,nodeGageStr%usgsId
387 !!$ print*, 'Ndg: Finish initializing Nudging_frxst_gage.csv'
388 !!$#endif
389 !!$end if ! gridded channel models
391 !========================
392 ! Muskingum routines are channel_options 1 and 2
393 if (nlst(did)%channel_option .eq. 1 .or. &
394 nlst(did)%channel_option .eq. 2) then
396 ! For reach-based/muskingum routing methods, we are currently requiring
397 ! the netcdf file for input of gages.
399 #ifdef MPP_LAND
400 if(my_id .eq. io_id) then
401 #endif
403 #ifdef HYDRO_D
404 print*, 'Ndg: Start initializing reach gages (netcdf)'
405 if(flushAll) flush(flushUnit)
406 #endif
408 allocate(nodeGageTmp%usgsId(nLinksL))
409 call read_reach_gage_collocation(nodeGageTmp%usgsId)
411 nGagesDomain=0
412 !check: nLinksL .eq. size(nodeGageTmp%usgsId)
413 !do ll=1,nLinksL
414 do ll=1,size(nodeGageTmp%usgsId)
415 if(nodeGageTmp%usgsId(ll) .ne. missingGage) nGagesDomain=nGagesDomain+1
416 end do
418 if(nGagesDomain .gt. 0) then
419 allocate(nodeGageStr%nodeId(nGagesDomain), nodeGageStr%usgsId(nGagesDomain))
420 nodeGageStr%usgsId = pack(nodeGageTmp%usgsId, mask=nodeGageTmp%usgsId .ne. missingGage)
421 ! This just index, we are NOT using comIds in nudging.
422 ! nodeGageStr%nodeId = pack(rt_domain(did)%linkId, mask=nodeGageTmp%usgsId .ne. '')
423 nodeGageStr%nodeId = pack((/(ii, ii=1,nLinksL)/), mask=nodeGageTmp%usgsId .ne. missingGage)
424 end if
425 deallocate(nodeGageTmp%usgsId)
427 #ifdef HYDRO_D
428 print*,'Ndg: nGagesDomain:',nGagesDomain
429 print*,'Ndg: nLinksL', nLinksL
430 print*,"Ndg: size(nodeGageStr%nodeId):", size(nodeGageStr%nodeId)
431 print*,&
432 'Ndg: nodeGageStr%usgsId((size(nodeGageStr%nodeId)-nGagesDomain+1):(size(nodeGageStr%nodeId))):',&
433 nodeGageStr%usgsId((size(nodeGageStr%nodeId)-nGagesDomain+1):(size(nodeGageStr%nodeId)))
434 print*,'Ndg: Finish initializing reach gages (netcdf)'
435 if(flushAll) flush(flushUnit)
436 #endif
438 #ifdef MPP_LAND
439 end if ! my_id .eq. io_id
440 !! Broadcast
441 call mpp_land_bcast_int1(nGagesDomain)
442 if(my_id .ne. io_id) then
443 allocate(nodeGageStr%nodeId(nGagesDomain))
444 allocate(nodeGageStr%usgsId(nGagesDomain))
445 endif
446 call mpp_land_bcast_char1d(nodeGageStr%usgsId)
447 call mpp_land_bcast_int1d(nodeGageStr%nodeId)
448 #endif
450 end if ! muskingum channel models
452 !=================================================
453 ! 3. Read nudging parameter files and reduce to the gages in the domain (nGagesDomain, etc above).
454 #ifdef MPP_LAND
455 if(my_id .eq. IO_id) then
456 #endif
457 nParamGages = get_netcdf_dim(nudgingParamFile, 'stationIdInd', 'init_stream_nudging')
459 nParamMonth=0
460 nParamThresh=0
461 nParamThreshCat=0
462 if(temporalPersistence) then
463 nParamMonth = get_netcdf_dim(nudgingParamFile, 'monthInd', 'init_stream_nudging')
464 nParamThresh = get_netcdf_dim(nudgingParamFile, 'threshInd', 'init_stream_nudging')
465 nParamThreshCat = get_netcdf_dim(nudgingParamFile, 'threshCatInd', 'init_stream_nudging')
467 #ifdef HYDRO_D
468 print*,'Ndg: nParamGages: ', nParamGages
469 print*,'Ndg: nParamGages', nParamGages
470 print*,'Ndg: nParamMonth', nParamMonth
471 print*,'Ndg: nParamThresh', nParamThresh
472 print*,'Ndg: nParamThreshCat', nParamThreshCat
473 if(flushAll) flush(flushUnit)
474 #endif
476 if(nParamMonth.eq.0 .or. nParamThresh.eq.0 .or. nParamThreshCat.eq.0) then
477 temporalPersistence = .false.
478 nParamThresh=0
479 nParamMonth=0
480 else if((nParamThresh+1) .ne. nParamThreshCat) then
481 temporalPersistence = .false.
482 nParamThresh=0
483 nParamMonth=0
484 end if
485 endif !temporal persistence
487 #ifdef HYDRO_D
488 print*,'Ndg: nParamGages: ', nParamGages
489 print*,'Ndg: nParamMonth', nParamMonth
490 print*,'Ndg: nParamThresh', nParamThresh
491 print*,'Ndg: nParamThreshCat', nParamThreshCat
492 print*,'Ndg: temporalPersistence: ', temporalPersistence
493 if(flushAll) flush(flushUnit)
494 #endif
496 allocate(nudgingParamsTmp%usgsId( nParamGages))
497 allocate(nudgingParamsTmp%R( nParamGages))
498 allocate(nudgingParamsTmp%G( nParamGages))
499 allocate(nudgingParamsTmp%tau( nParamGages))
500 allocate(nudgingParamsTmp%qThresh( nParamThresh, nParamMonth, nParamGages))
501 allocate(nudgingParamsTmp%expCoeff(nParamThresh+1, nParamMonth, nParamGages))
503 call read_nudging_param_file(nudgingParamFile, &
504 nudgingParamsTmp%usgsId, &
505 nudgingParamsTmp%R, &
506 nudgingParamsTmp%G, &
507 nudgingParamsTmp%tau, &
508 nudgingParamsTmp%qThresh, &
509 nudgingParamsTmp%expCoeff )
511 #ifdef HYDRO_D
512 do ii=1,nParamThresh
513 print*,'Ndg: minval( nudgingParamsTmp%qThresh(ii,:,:)), ii=', ii,': ', &
514 minval( nudgingParamsTmp%qThresh(ii,:,:))
515 end do
516 do ii=1,nParamThreshCat
517 print*,'Ndg: minval( nudgingParamsTmp%expCoeff(ii,:,:)), ii=', ii,': ', &
518 minval( nudgingParamsTmp%expCoeff(ii,:,:))
519 end do
520 if(flushAll) flush(flushUnit)
521 #endif
523 ! Reduce the parameters to just the gages in the domain
524 allocate(whParamsDom(nParamGages))
526 call whichInLoop2(nudgingParamsTmp%usgsId, nodeGageStr%usgsId, whParamsDom, nGagesWParamsDom)
528 #ifdef HYDRO_D
529 if(nGagesWParamsDom .ne. nGagesDomain) then
530 print*,'Ndg: WARNING Gages are apparently missing from the nudgingParams.nc file'
531 print*,'Ndg: WARNING nGagesWParamsDom: ', nGagesWParamsDom
532 print*,'Ndg: WARNING nGagesDomain: ', nGagesDomain
533 if(flushAll) flush(flushUnit)
534 end if
535 #endif
537 allocate(nudgingParamsStr%usgsId(nGagesWParamsDom))
538 allocate(nudgingParamsStr%R( nGagesWParamsDom))
539 allocate(nudgingParamsStr%G( nGagesWParamsDom))
540 allocate(nudgingParamsStr%tau( nGagesWParamsDom))
541 if(temporalPersistence) then
542 allocate(nudgingParamsStr%qThresh( nParamThresh, nParamMonth, nGagesWParamsDom))
543 allocate(nudgingParamsStr%expCoeff(nParamThresh+1, nParamMonth, nGagesWParamsDom))
544 end if
546 count=1
547 do kk=1,nParamGages
548 if(whParamsDom(kk) .gt. 0) then
549 nudgingParamsStr%usgsId(count) = nudgingParamsTmp%usgsId(kk)
550 nudgingParamsStr%R(count) = nudgingParamsTmp%R(kk)
551 nudgingParamsStr%G(count) = nudgingParamsTmp%G(kk)
552 nudgingParamsStr%tau(count) = nudgingParamsTmp%tau(kk)
553 if(temporalPersistence) then
554 nudgingParamsStr%qThresh( :,:,count) = nudgingParamsTmp%qThresh( :,:,kk)
555 nudgingParamsStr%expCoeff(:,:,count) = nudgingParamsTmp%expCoeff(:,:,kk)
556 endif
557 count=count+1
558 end if
559 end do
561 deallocate(whParamsDom)
562 deallocate(nudgingParamsTmp%usgsId, nudgingParamsTmp%R )
563 deallocate(nudgingParamsTmp%G, nudgingParamsTmp%tau)
564 if(temporalPersistence) then
565 deallocate(nudgingParamsTmp%qThresh)
566 deallocate(nudgingParamsTmp%expCoeff)
567 end if
569 #ifdef HYDRO_D
570 print*,'Ndg: nudgingParamsStr%usgsId', nudgingParamsStr%usgsId(size(nudgingParamsStr%usgsId))
571 print*,'Ndg: nudgingParamsStr%R', nudgingParamsStr%R(size(nudgingParamsStr%R))
572 print*,'Ndg: nudgingParamsStr%G', nudgingParamsStr%G(size(nudgingParamsStr%G))
573 print*,'Ndg: nudgingParamsStr%tau', nudgingParamsStr%tau(size(nudgingParamsStr%tau))
574 if(temporalPersistence) then
575 print*,'Ndg: nudgingParamsStr%qThresh', nudgingParamsStr%qThresh(1,1,size(nudgingParamsStr%tau))
576 print*,'Ndg: nudgingParamsStr%expCoeff',nudgingParamsStr%expCoeff(1,1,size(nudgingParamsStr%tau))
577 end if
578 if(flushAll) flush(flushUnit)
579 #endif /* HYDRO_D */
580 #ifdef MPP_LAND
581 endif ! my_id .eq. io_id
583 !! Broadcast
584 call mpp_land_bcast_int1(nGagesWParamsDom)
585 if(my_id .ne. io_id) then
586 allocate(nudgingParamsStr%usgsId(nGagesWParamsDom))
587 allocate(nudgingParamsStr%R(nGagesWParamsDom))
588 allocate(nudgingParamsStr%G(nGagesWParamsDom))
589 allocate(nudgingParamsStr%tau(nGagesWParamsDom))
590 endif
591 call mpp_land_bcast_char1d(nudgingParamsStr%usgsId)
592 call mpp_land_bcast_real_1d(nudgingParamsStr%R)
593 call mpp_land_bcast_real_1d(nudgingParamsStr%G)
594 call mpp_land_bcast_real_1d(nudgingParamsStr%tau)
595 call mpp_land_bcast_logical(temporalPersistence)
596 if(temporalPersistence) then
597 call mpp_land_bcast_int1(nParamThresh)
598 call mpp_land_bcast_int1(nParamMonth)
599 if(my_id .ne. io_id) then
600 allocate(nudgingParamsStr%qThresh( nParamThresh, nParamMonth, nGagesWParamsDom))
601 allocate(nudgingParamsStr%expCoeff(nParamThresh+1, nParamMonth, nGagesWParamsDom))
602 endif
604 call mpp_land_bcast_real3d(nudgingParamsStr%qThresh)
605 call mpp_land_bcast_real3d(nudgingParamsStr%expCoeff)
606 end if ! temporalPersistence
607 #endif /* MPP_LAND */
609 allocate(obsStaticStr(size(nudgingParamsStr%usgsId)))
611 allocate(obsStaticStr_SoA%usgsId(size(nudgingParamsStr%usgsId)))
612 allocate(obsStaticStr_SoA%R(size(nudgingParamsStr%usgsId)))
613 allocate(obsStaticStr_SoA%tau(size(nudgingParamsStr%usgsId)))
614 allocate(obsStaticStr_SoA%G(size(nudgingParamsStr%usgsId)))
615 allocate(obsStaticStr_SoA%obsCellInd(size(nudgingParamsStr%usgsId)))
616 ! This seems silly now that I'm solving the whole thing at init.
617 obsStaticStr_SoA%R(:) = nudgingParamsStr%R(:)
618 obsStaticStr_SoA%G(:) = nudgingParamsStr%G(:)
619 obsStaticStr_SoA%tau(:) = nudgingParamsStr%tau(:)
620 obsStaticStr_SoA%usgsId(:) = nudgingParamsStr%usgsId(:)
621 obsStaticStr(:)%usgsId = nudgingParamsStr%usgsId
623 allocate(obsStaticStr_SoA%lastObsDischarge(nlst(did)%nLastObs, size(nudgingParamsStr%usgsId)))
624 allocate(obsStaticStr_SoA%lastObsModelDischarge(nlst(did)%nLastObs, size(nudgingParamsStr%usgsId)))
625 allocate(obsStaticStr_SoA%lastObsTime(nlst(did)%nLastObs, size(nudgingParamsStr%usgsId)))
626 allocate(obsStaticStr_SoA%lastObsQuality(nlst(did)%nLastObs, size(nudgingParamsStr%usgsId)))
628 !=================================================
629 ! 4. Read in 'nudgingLastObs.nc' file (initialization and broadcasting come later)
630 #ifdef MPP_LAND
631 if(my_id .eq. io_id) then
632 if(temporalPersistence) then
634 !! Initalizing the structure is not contingent upon there
635 !! being a file on disk
636 obsStaticStr_SoA%lastObsDischarge = real(-9999)
637 obsStaticStr_SoA%lastObsModelDischarge = real(-9999)
638 obsStaticStr_SoA%lastObsTime = missingLastObsTime
639 obsStaticStr_SoA%lastObsQuality = real(0) !! keep this zero
641 !! if blank, look for file at the current/init time
642 if(trim(nudgingLastObsFile) .eq. '') then
643 nudgingLastObsFileTry = 'nudgingLastObs.' // nlst(did)%olddate // '.nc'
644 else
645 nudgingLastObsFileTry = nudgingLastObsFile
646 end if
647 lastObsFile = find_nudging_last_obs_file(nudgingLastObsFileTry)
648 if(trim(lastObsFile) .ne. '') then
649 nStnLastObs = get_netcdf_dim(nudgingLastObsFileTry, 'stationIdInd', 'init_stream_nudging')
650 if(nStnLastObs .gt. 0) then
651 print*,'Reading in nudgingLastObsFileTry: ', trim(nudgingLastObsFileTry)
652 allocate(lastObsStr_SoA%usgsId(nStnLastObs))
653 allocate(lastObsStr_SoA%lastObsDischarge(nlst(did)%nLastObs,nStnLastObs))
654 allocate(lastObsStr_SoA%lastObsModelDischarge(nlst(did)%nLastObs,nStnLastObs))
655 allocate(lastObsStr_SoA%lastObsTime(nlst(did)%nLastObs,nStnLastObs))
656 allocate(lastObsStr_SoA%lastObsQuality(nlst(did)%nLastObs,nStnLastObs))
658 allocate(lastObsStr(nStnLastObs))
660 allocate(g_nudge(rt_domain(1)%gnlinksl))
661 call read_nudging_last_obs_soa(nudgingLastObsFileTry, lastObsStr_SoA, g_nudge)
662 end if ! nStnLastObs .gt. 0 -> read the file
663 endif! trim(lastObsFile).ne.''
664 end if ! temporalPersistence
665 else
666 allocate(g_nudge(1))
667 endif
669 if(my_id .eq. io_id) lastObsFileFound=allocated(g_nudge)
670 call mpp_land_bcast_logical(lastObsFileFound)
671 if(lastObsFileFound) then
672 call ReachLS_decomp(g_nudge, RT_DOMAIN(1)%nudge)
673 if(my_id .eq. io_id) deallocate(g_nudge)
674 endif
675 if(my_id .NE. io_id) deallocate(g_nudge)
676 #endif
678 !=================================================
679 ! 5. Sort out which gages are actually in the domain.
680 ! JLM: there can be an inconsistency between the gages with parameters
681 ! JLM: and the gages with distances.
682 ! JLM: How to handle locations without parameters?
683 ! JLM: How to determine the functional list of gages in the domain?
684 !
686 !=================================================
687 ! 6. Set up the global time parameters
688 #ifdef MPP_LAND
689 ! MPP need to broadcast the initial time and setup the dt
690 if(my_id .eq. io_id) then
691 lsmDt = nlst(did)%dt
692 endif
693 call mpp_land_bcast_int1(lsmDt)
694 #else
695 lsmDt = nlst(did)%dt
696 #endif
698 !S - -+- - !
699 !E ! - -+- -
700 ! |- - - -L* * * *|- - - -|- - - -|
701 ! t t w w w w t t
702 ! - time chunks of obsResolution
703 ! | separators denoting lsmDt chunks
704 ! * obsRresolution points in the current lsmDt/hydro time advance
705 ! + is the center of the assim window size = tau*2, at that the hydro time
706 ! ! are the bounds of all assim windows
707 ! L is LSM time
708 ! S-E The start and end times for hydro advance
709 ! t denotes times needed by size of tau
710 ! w denotes times needed by size of lsmDt
711 ! Regardless of tau, need window of size (2*tau+lsmDt)
712 ! which will have (2*tau+lsmDt)/obsResolution+1 observation times in it
713 ! the +1 is for the zeroth time.
715 maxTau = maxval(nudgingParamsStr%tau) ! This is fixed as tau does not change for a cycle.
716 write(obsResolution, '(i0.2)') obsResolutionInt
717 nObsTimes = 2*ceiling(maxTau/obsResolutionInt) + ceiling(real(lsmDt)/(60.*obsResolutionInt)) + 1
719 #ifdef HYDRO_D
720 #ifdef MPP_LAND
721 if(my_id .eq. io_id) then
722 #endif
723 print*,'Ndg: obsResolution:',obsResolution
724 print*,'Ndg: maxTau: ', maxTau
725 print*,'Ndg: nObsTimes: ', nObsTimes
726 if(flushAll) flush(flushUnit)
727 #ifdef MPP_LAND
728 end if
729 #endif
730 #endif
732 allocate(obsTimeStr(nObsTimes))
733 !! initialize these as blanks and 'none'
734 obsTimeStr(:)%time = ''
735 obsTimeStr(:)%updateTime = 'none'
737 #ifdef MPP_LAND
738 !=================================================
739 ! 8. MPP: Keep a copy of the full chanlen on image 0
740 !! only necessary if doing spatial nudging interpolation
741 if(my_id .eq. io_id) allocate(chanlen_image0(nLinksL))
742 if(my_id .ne. io_id) allocate(chanlen_image0(1)) !! memory fix, the next call modifies it's second arg
743 call ReachLS_write_io(rt_domain(did)%chanlen, chanlen_image0)
744 #endif
747 !=================================================
748 ! 9. Solve obsStaticStr once and for all on image 0 (could parallelize)
749 ! This also solves nudgeSpatial: is spatial nudging active? Needed in 10.
750 ! Remove lastObsStr (only needed for ingest of nudingLastObs/restart file).
751 call obs_static_to_struct()
752 #ifdef MPP_LAND
753 if(my_id .eq. io_id) then
754 #endif
755 if(allocated(lastObsStr)) deallocate(lastObsStr)
756 #ifdef MPP_LAND
757 endif
758 #endif
760 !=================================================
761 !10. Solve the bias terms for when biasWindowBeforeT0 = .TRUE. (in the forecast)
763 !if(persistBias .and. biasWindowBeforeT0) then
764 ! do gg=1,nGages.........
765 ! ! solve the bias
766 ! end do ! gg=1,nGages
767 !end if ! if(persistBias) then
769 !=================================================
770 #ifdef HYDRO_D
771 #ifdef MPP_LAND
772 if(my_id .eq. io_id) &
773 #endif
774 print*, "Ndg: Finish init_stream_nudging"
775 if(flushAll) flush(flushUnit)
776 #endif
778 !!$#ifdef HYDRO_D /* un-ifdef to get timing results */
779 !!$#ifdef MPP_LAND
780 !!$if(my_id .eq. io_id) then
781 !!$#endif
782 !!$ call nudging_timer(endCodeTimeAcc)
783 !!$ call accum_nudging_time(startCodeTimeAcc, endCodeTimeAcc, 'init_stream_nudging', .true.)
784 !!$#ifdef MPP_LAND
785 !!$end if
786 !!$#endif
787 !!$if(flushAll) flush(flushUnit)
788 !!$#endif /* HYDRO_D */ /* un-ifdef to get timing results */
790 end subroutine init_stream_nudging
792 !===================================================================================================
793 ! Program Names:
794 ! subroutine setup_stream_nudging
795 ! Author(s)/Contact(s):
796 ! James L McCreight <jamesmcc><ucar><edu>
797 ! Abstract:
798 ! Setup the nudging for the current hydroTime, only establishes the
799 ! shared obsTimeStr above.
800 ! History Log:
801 ! 6/04/15 -Created, JLM.
802 ! Usage:
803 ! Parameters:
804 ! Input Files:
805 ! Output Files:
806 ! Condition codes:
807 ! User controllable options:
808 ! Notes:
810 subroutine setup_stream_nudging(hydroDT)
812 use module_RT_data, only: rt_domain
813 use module_nudging_utils, only: whichLoop, &
814 whUniLoop, &
815 accum_nudging_time, &
816 nudging_timer
818 use module_date_utils_nudging, only: geth_newdate, &
819 round_resolution_minute, &
820 geth_idts
821 implicit none
823 integer, intent(in) :: hydroDT ! the number of seconds of hydro advance from lsmTime
824 !integer :: ff
825 integer :: ii, tt, oo
826 character(len=19) :: hydroTime, obsHydroTime ! hydro model time and corresponding observation
827 character(len=19), dimension(nObsTimes) :: obsTimes ! obs times in the current window
828 integer :: oldDiff, nShiftLeft
830 logical, allocatable, dimension(:) :: theMask
831 integer, allocatable, dimension(:) :: whObsMiss
832 logical, allocatable, dimension(:) :: obsTimeStrAllocated
834 integer :: nObsMiss
835 integer :: did=1 !! jlm: assuming did=1
836 integer :: nWhObsMiss
838 #ifdef MPP_LAND
839 integer :: nGages, nCellsInR, cc, nLinkAff, nStatic, iiImage
840 #endif
842 !!$#ifdef HYDRO_D /* un-ifdef to get timing results */
843 !!$real :: startCodeTime, endCodeTime
844 !!$real :: startCodeTimeOuter
845 !!$
846 !!$#ifdef MPP_LAND
847 !!$if(my_id .eq. io_id) &
848 !!$#endif
849 !!$ call nudging_timer(startCodeTimeOuter)
850 !!$if(flushAll) flush(flushUnit)
851 !!$#endif /* HYDRO_D :: un-ifdef to get timing results */
853 #ifdef HYDRO_D
854 #ifdef MPP_LAND
855 if(my_id .eq. io_id) &
856 #endif
857 print*,'Ndg: start setup_stream_nudging'
858 if(flushAll) flush(flushUnit)
859 #endif /* HYDRO_D */
861 !!$#ifdef HYDRO_D
862 !!$#ifdef MPP_LAND
863 !!$if(my_id .eq. io_id) &
864 !!$#endif
865 !!$ call nudging_timer(startCodeTime)
866 !!$if(flushAll) flush(flushUnit)
867 !!$#endif /* HYDRO_D */
869 !#ifdef MPP_LAND - just do this section on all images. There is no IO.
870 ! The hydro model time as a string.
871 #ifdef MPP_LAND
872 !update and broadcast the lsm time
873 if(my_id .eq. io_id) lsmTime = nlst(did)%olddate
874 if(my_id .eq. io_id) initTime = nlst(did)%startdate
875 call mpp_land_bcast_char(19, lsmTime )
876 call mpp_land_bcast_char(19, initTime)
877 # else
878 lsmTime = nlst(did)%olddate
879 initTime = nlst(did)%startdate
880 #endif
882 ! This is apparently fixed so that olddate is lsmTime and not
883 ! behind by 1 lsm timestep when the hydro model is run.
884 !call geth_newdate(hydroTime, nlst(1)%olddate, hydroDT+nint(nlst(1)%dt))
885 !call geth_newdate(hydroTime, lsmTimeMinusDt, hydroDT+lsmDt)
886 call geth_newdate(hydroTime, lsmTime, hydroDT)
888 ! Calculate the closest multiple of obsResolution to the hydroTime
889 obsHydroTime = round_resolution_minute(hydroTime, obsResolutionInt)
890 #ifdef HYDRO_D
891 #ifdef MPP_LAND
892 if(my_id .eq. io_id) &
893 #endif
894 print*,'Ndg: obsHydroTime: ', obsHydroTime
895 if(flushAll) flush(flushUnit)
896 #endif
898 ! Now solve all of the observation times in the current nudging window.
899 ! nObsTimes is set at init from maxTau and obsResolution.
900 ! These times correspond to the timestamps on the observation files, which
901 ! mark the center of the time period they represent with width obsResolution.
902 do tt=1,nObsTimes
903 call geth_newdate(obsTimes(tt), obsHydroTime, &
904 !obsResolutionInt*(tt - (nObsTimes+1)/2 )*60)
905 (tt-1-ceiling(maxTau/obsResolutionInt))*60*obsResolutionInt )
906 end do
908 ! If this is the first setup, these are all blank. This is their init value.
909 if(all(obsTimeStr(:)%time .eq. '')) obsTimeStr(:)%time = obsTimes
911 !print*,'Ndg: hydroTime: ', hydroTime
912 !print*,'Ndg: obsHydroTime: ', obsHydroTime
913 !print*,'Ndg: obsTimes: ', obsTimes
914 !print*,'Ndg: obsTimeStr: before obsTimeStr(:)%time: ', obsTimeStr(:)%time
916 ! If there are existing observations older than the first obsTime for this
917 ! nudging window, shift them left. This should just shift one position left
918 ! with each advance in time resolution.
919 call geth_idts(obsTimeStr(1)%time, obsTimes(1), oldDiff)
920 !print*,"Ndg: olddiff:", oldDiff
921 if(oldDiff .eq. 0) then
922 nShiftLeft=0
923 else
924 ! nShiftLeft should not exceed nObsTimes
925 nShiftLeft = min( abs(oldDiff/obsResolutionInt/60), nObsTimes )
926 end if
927 !print*,"Ndg: nShiftLeft:", nShiftLeft
929 if(nShiftLeft .gt. 0 .and. nShiftLeft .lt. nObsTimes) then
930 do tt=1,nObsTimes-nShiftLeft
932 obsTimeStr(tt)%time = obsTimeStr(tt+nShiftLeft)%time
933 obsTimeStr(tt)%updateTime = obsTimeStr(tt+nShiftLeft)%updateTime
935 if(allocated(obsTimeStr(tt)%allCellInds)) deallocate(obsTimeStr(tt)%allCellInds)
936 if(allocated(obsTimeStr(tt)%nGageCell)) deallocate(obsTimeStr(tt)%nGageCell)
937 if(allocated(obsTimeStr(tt)%obsStr)) deallocate(obsTimeStr(tt)%obsStr)
939 if(allocated(obsTimeStr(tt+nShiftLeft)%allCellInds)) then
940 allocate(obsTimeStr(tt)%allCellInds(size(obsTimeStr(tt+nShiftLeft)%allCellInds)))
941 obsTimeStr(tt)%allCellInds = obsTimeStr(tt+nShiftLeft)%allCellInds
942 end if
944 if(allocated(obsTimeStr(tt+nShiftLeft)%nGageCell)) then
945 allocate(obsTimeStr(tt)%nGageCell(size(obsTimeStr(tt+nShiftLeft)%nGageCell)))
946 obsTimeStr(tt)%nGageCell = obsTimeStr(tt+nShiftLeft)%nGageCell
947 end if
949 if(allocated(obsTimeStr(tt+nShiftLeft)%obsStr)) then
950 allocate(obsTimeStr(tt)%obsStr(size(obsTimeStr(tt+nShiftLeft)%obsStr)))
951 obsTimeStr(tt)%obsStr = obsTimeStr(tt+nShiftLeft)%obsStr
952 end if
954 end do
955 endif
956 if(nShiftLeft .gt. 0) then
957 ! here tt= nObsTimes
958 do tt=nObsTimes-nShiftLeft+1,nObsTimes
959 obsTimeStr(tt)%time = obsTimes(tt) !'' !! JLM is this a fix? why?
960 obsTimeStr(tt)%updateTime = 'none'
961 if(allocated(obsTimeStr(tt)%allCellInds)) deallocate(obsTimeStr(tt)%allCellInds)
962 if(allocated(obsTimeStr(tt)%nGageCell)) deallocate(obsTimeStr(tt)%nGageCell)
963 if(allocated(obsTimeStr(tt)%obsStr)) deallocate(obsTimeStr(tt)%obsStr)
964 end do
965 end if
967 !if(.NOT. all(obsTimeStr(:)%time .EQ. obsTimes)) &
968 ! call hydro_stop("obsTimeStr(:)%times not what they should be. Please investigate.")
970 ! Updates obs already in memory if flag is set.
971 ! Not going to be used for IOC. This will happen later.
972 ! Should the frequency be in model time or in real time? Probably real time
973 ! and will require a clock time of 'obsLast' checked in obsTimeStr?
974 ! call check_for_new_obs()
976 ! Read in obs at times not already checked, eg. new obs this window
977 ! Obs not already checked have updateTime='none'
978 ! Obs already checked for with:
979 ! : a missing file have updateTime='no file'
980 ! : no observations in the domain have updateTime='no obs'.
981 #ifdef MPP_LAND
982 if(my_id .eq. io_id) then
983 #endif
984 nWhObsMiss = size(obsTimeStr(:)%updateTime)
985 allocate(theMask(nWhObsMiss), whObsMiss(nWhObsMiss) )
986 do ii=1,size(obsTimeStr(:)%updateTime)
987 theMask(ii) = trim(obsTimeStr(ii)%updateTime) .eq. 'none'
988 end do
989 call whichLoop(theMask, whObsMiss, nObsMiss)
990 deallocate(theMask)
991 #ifdef MPP_LAND
992 end if
994 ! Broadcast the basic info above.
995 call mpp_land_bcast_int1(nWhObsMiss)
996 if(my_id .ne. io_id) then
997 if(allocated(whObsMiss)) deallocate(whObsMiss)
998 allocate(whObsMiss(nWhObsMiss))
999 endif
1000 call mpp_land_bcast_int1d(whObsMiss)
1001 call mpp_land_bcast_int1(nObsMiss)
1002 #endif
1004 !!$#ifdef HYDRO_D
1005 !!$if(my_id .eq. io_id) then
1006 !!$ call nudging_timer(endCodeTime)
1007 !!$ call accum_nudging_time(startCodeTime, endCodeTime, &
1008 !!$ 'setup1: prelim', .false.)
1009 !!$ call nudging_timer(startCodeTime)
1010 !!$if(flushAll) flush(flushUnit)
1011 !!$endif
1012 !!$#endif
1014 !bring in timeslice files
1015 allocate(obsTimeStrAllocated(nObsMiss))
1017 do ii=1,nObsMiss ! if nObsMiss is zero, this loop is skipped?
1018 #ifdef MPP_LAND
1019 if(readTimesliceParallel) then
1020 iiImage = mod(ii-1,numprocs)
1021 else
1022 iiImage = 0
1023 end if
1024 if(my_id .eq. iiImage) then !! this would give parallel IO
1025 #endif
1026 tt = whObsMiss(ii)
1027 ! set/reset this obsTime
1028 obsTimeStr(tt)%time = obsTimes(tt)
1029 obsTimeStr(tt)%updateTime = 'none'
1030 !! This might be paranoid/overkill...
1031 if(allocated(obsTimeStr(tt)%allCellInds)) deallocate(obsTimeStr(tt)%allCellInds)
1032 if(allocated(obsTimeStr(tt)%nGageCell)) deallocate(obsTimeStr(tt)%nGageCell)
1033 if(allocated(obsTimeStr(tt)%obsStr)) deallocate(obsTimeStr(tt)%obsStr)
1034 call timeslice_file_to_struct(tt) ! uses obsTimeStr%time to get file
1035 obsTimeStrAllocated(ii) = allocated(obsTimeStr(tt)%obsStr)
1036 !obsTimeStrAllocated = allocated(obsTimeStr(tt)%obsStr)
1037 #ifdef MPP_LAND
1038 end if ! my_id .eq. iiImage
1039 #endif
1040 end do
1042 !!$#ifdef HYDRO_D
1043 !!$if(my_id .eq. io_id) then
1044 !!$ !print*,'Ndg: obsTimeStr(tt)%time: ',obsTimeStr(tt)%time
1045 !!$ !print*,'Ndg: obsTimeStrAllocated(ii), ii: ', obsTimeStrAllocated(ii), ii
1046 !!$ call nudging_timer(endCodeTime)
1047 !!$ call accum_nudging_time(startCodeTime, endCodeTime, &
1048 !!$ 'setup2.0: timeslice_file_to_struct before bcast', .false.)
1049 !!$ !call nudging_timer(startCodeTime) ! skip this b/c there's no mpi sync till after next section
1050 !!$ if(flushAll) flush(flushUnit)
1051 !!$endif
1052 !!$#endif
1054 #ifdef MPP_LAND
1055 ! broadcast the IO from above
1056 do ii=1,nObsMiss ! if nObsMiss is zero, this loop is skipped?
1057 if(readTimesliceParallel) then
1058 iiImage = mod(ii-1,numprocs)
1059 else
1060 iiImage = 0
1061 end if
1062 tt = whObsMiss(ii)
1063 ! broadcast
1064 ! Here have to expose some of the guts of timeslice_file_to_struct
1065 ! Note nGages depends on ii/tt
1066 call mpp_land_bcast_int1_root(tt, iiImage)
1067 call mpp_land_bcast_logical_root(obsTimeStrAllocated(ii), iiImage)
1068 if(obsTimeStrAllocated(ii)) then
1069 if(my_id .eq. iiImage) nGages=size(obsTimeStr(tt)%obsStr)
1070 call mpp_land_bcast_int1_root(nGages, iiImage)
1071 if(my_id .ne. iiImage) then
1072 if(allocated(obsTimeStr(tt)%obsStr)) deallocate(obsTimeStr(tt)%obsStr)
1073 allocate(obsTimeStr(tt)%obsStr(nGages))
1074 endif
1075 ! Variables in order assigned in the call in timeslice_file_to_struct
1076 call mpp_land_bcast_char_root(19,obsTimeStr(tt)%time, iiImage)
1077 call mpp_land_bcast_char_root(19,obsTimeStr(tt)%updateTime, iiImage)
1078 call mpp_land_bcast_char1d_root(obsTimeStr(tt)%obsStr(:)%usgsId, iiImage)
1079 call mpp_land_bcast_char1d_root(obsTimeStr(tt)%obsStr(:)%obsTime, iiImage)
1080 call mpp_land_bcast_real_1d_root(obsTimeStr(tt)%obsStr(:)%obsQC, iiImage)
1081 call mpp_land_bcast_real_1d_root(obsTimeStr(tt)%obsStr(:)%obsDischarge, iiImage)
1082 end if
1083 end do
1084 #endif
1086 !!$#ifdef HYDRO_D
1087 !!$if(my_id .eq. io_id) then
1088 !!$ call nudging_timer(endCodeTime)
1089 !!$ call accum_nudging_time(startCodeTime, endCodeTime, &
1090 !!$ 'setup2.1: timeslice_file_to_struct after bcast', .false.)
1091 !!$ call nudging_timer(startCodeTime) ! images are synched so reset timer
1092 !!$ if(flushAll) flush(flushUnit)
1093 !!$endif
1094 !!$#endif
1096 !! get index of static info in obsStaticStr
1097 do ii=1,nObsMiss ! if nObsMiss is zero, this loop is skipped?
1098 #ifdef MPP_LAND
1099 if(readTimesliceParallel) then
1100 iiImage = mod(ii-1,numprocs)
1101 else
1102 iiImage = 0
1103 end if
1104 if(my_id .eq. iiImage) then
1105 tt = whObsMiss(ii)
1106 #endif
1107 if(obsTimeStrAllocated(ii)) then
1108 allocate(theMask(size(nudgingParamsStr%usgsId)))
1109 do oo=1,size(obsTimeStr(tt)%obsStr(:)%obsStaticInd)
1110 ! If the gage is not in the parameter file, skip it
1111 if(.not. any(nudgingParamsStr%usgsId .eq. obsTimeStr(tt)%obsStr(oo)%usgsId)) then
1112 !! If you ended up here, then filterObsToDom should not be on...
1113 if(filterObsToDom) call hydro_stop('obs_static_to_struct: logical clash with filterObsToDom')
1114 call accumulate_nwis_not_in_RLAndParams(nwisNotRLAndParams, &
1115 nwisNotRLAndParamsCount, &
1116 obsTimeStr(tt)%obsStr(oo)%usgsId)
1117 cycle !! skip this observation.
1118 endif
1119 theMask = obsStaticStr_SoA%usgsId .eq. obsTimeStr(tt)%obsStr(oo)%usgsId
1120 !! This is a double (tripple?!) for loop... thanks to whUniLoop
1121 obsTimeStr(tt)%obsStr(oo)%obsStaticInd = whUniLoop(theMask)
1122 end do
1123 deallocate(theMask)
1124 end if ! oo
1125 #ifdef MPP_LAND
1126 end if ! ii
1127 #endif
1128 end do
1130 !!$#ifdef HYDRO_D
1131 !!$if(my_id .eq. io_id) then
1132 !!$ call nudging_timer(endCodeTime)
1133 !!$ call accum_nudging_time(startCodeTime, endCodeTime, &
1134 !!$ 'setup3: obs_static_to_struct', .false.)
1135 !!$ call nudging_timer(startCodeTime)
1136 !!$ if(flushAll) flush(flushUnit)
1137 !!$endif
1138 !!$#endif
1140 #ifdef MPP_LAND
1141 do ii=1,nObsMiss ! broadcast IO from above. (if nObsMiss is zero, this loop is skipped)
1142 if(readTimesliceParallel) then
1143 iiImage = mod(ii-1,numprocs)
1144 else
1145 iiImage = 0
1146 end if
1147 ! broadcast
1148 if(obsTimeStrAllocated(ii)) then
1149 tt = whObsMiss(ii)
1150 ! Treat obs_static_to_strutct
1151 call mpp_land_bcast_int1d_root(obsTimeStr(tt)%obsStr(:)%obsStaticInd, iiImage)
1152 call mpp_land_bcast_char1d_root(obsStaticStr_SoA%usgsId(:), iiImage)
1153 call mpp_land_bcast_int1d_root(obsStaticStr_SoA%obsCellInd(:), iiImage)
1154 call mpp_land_bcast_real_1d_root(obsStaticStr_SoA%R(:), iiImage)
1155 call mpp_land_bcast_real_1d_root(obsStaticStr_SoA%G(:), iiImage)
1156 call mpp_land_bcast_real_1d_root(obsStaticStr_SoA%tau(:), iiImage)
1158 nStatic = sum( (/ (1, cc=1,size(obsStaticStr_SoA%obsCellInd(:))) /), &
1159 mask=obsStaticStr_SoA%obsCellInd(:) .ne. 0 )
1160 do cc=1,nStatic
1161 if(my_id .eq. iiImage) nCellsInR = size(obsStaticStr(cc)%cellsAffected)
1162 call mpp_land_bcast_int1_root(nCellsInR, iiImage)
1163 if(my_id .ne. iiImage) then
1164 if(allocated(obsStaticStr(cc)%cellsAffected)) deallocate(obsStaticStr(cc)%cellsAffected)
1165 if(allocated(obsStaticStr(cc)%dist)) deallocate(obsStaticStr(cc)%dist)
1166 if(allocated(obsStaticStr(cc)%ws)) deallocate(obsStaticStr(cc)%ws)
1167 allocate(obsStaticStr(cc)%cellsAffected(nCellsInR))
1168 allocate(obsStaticStr(cc)%dist(nCellsInR))
1169 allocate(obsStaticStr(cc)%ws(nCellsInR))
1170 end if
1171 call mpp_land_bcast_int1d_root(obsStaticStr(cc)%cellsAffected, iiImage)
1172 call mpp_land_bcast_real_1d_root(obsStaticStr(cc)%dist, iiImage)
1173 call mpp_land_bcast_real_1d_root(obsStaticStr(cc)%ws, iiImage)
1174 end do
1175 end if
1176 end do
1177 #endif
1179 !!$#ifdef HYDRO_D
1180 !!$if(my_id .eq. io_id) then
1181 !!$ call nudging_timer(endCodeTime)
1182 !!$ call accum_nudging_time(startCodeTime, endCodeTime, &
1183 !!$ 'setup4: bcast obs_static_struct', .false.)
1184 !!$ call nudging_timer(startCodeTime)
1185 !!$ if(flushAll) flush(flushUnit)
1186 !!$endif
1187 !!$#endif
1189 do ii=1,nObsMiss ! if nObsMiss is zero, this loop is skipped?
1190 #ifdef MPP_LAND
1191 iiImage = 0! causes issues also refactor tally_affected_links? mod(ii-1,numprocs)
1192 if(my_id .eq. iiImage) then
1193 #endif
1194 tt = whObsMiss(ii)
1195 if(obsTimeStrAllocated(ii)) call tally_affected_links(tt)
1196 endif
1197 #ifdef MPP_LAND
1198 end do
1199 #endif
1201 !!$#ifdef HYDRO_D
1202 !!$if(my_id .eq. io_id) then
1203 !!$ call nudging_timer(endCodeTime)
1204 !!$ call accum_nudging_time(startCodeTime, endCodeTime, &
1205 !!$ 'setup5: tally_affected_links', .false.)
1206 !!$ call nudging_timer(startCodeTime)
1207 !!$ if(flushAll) flush(flushUnit)
1208 !!$endif
1209 !!$#endif
1211 #ifdef MPP_LAND
1212 do ii=1,nObsMiss ! if nObsMiss is zero, this loop is skipped?
1213 iiImage = 0 ! this must happen on the same images as the above two mod(ii-1,numprocs)
1214 ! broadcast
1215 if(obsTimeStrAllocated(ii)) then
1216 tt = whObsMiss(ii)
1217 ! Treat tally_affected_links
1218 if(my_id .eq. iiImage) nLinkAff = size(obsTimeStr(tt)%allCellInds)
1219 call mpp_land_bcast_int1_root(nLinkAff, iiImage)
1221 if(my_id .ne. iiImage) then
1222 if(allocated(obsTimeStr(tt)%allCellInds)) deallocate(obsTimeStr(tt)%allCellInds)
1223 if(allocated(obsTimeStr(tt)%nGageCell)) deallocate(obsTimeStr(tt)%nGageCell)
1224 allocate(obsTimeStr(tt)%allCellInds(nLinkAff))
1225 allocate(obsTimeStr(tt)%nGageCell(nLinkAff))
1226 end if
1227 call mpp_land_bcast_int1d_root(obsTimeStr(tt)%allCellInds, iiImage)
1228 call mpp_land_bcast_int1d_root(obsTimeStr(tt)%nGageCell, iiImage)
1229 end if
1230 end do
1231 #endif
1233 !!$#ifdef HYDRO_D
1234 !!$if(my_id .eq. io_id) then
1235 !!$ call nudging_timer(endCodeTime)
1236 !!$ call accum_nudging_time(startCodeTime, endCodeTime, &
1237 !!$ 'setup6: after bcast tally_aff', .false.)
1238 !!$ call nudging_timer(startCodeTime)
1239 !!$ if(flushAll) flush(flushUnit)
1240 !!$endif
1241 !!$#endif
1243 #ifdef HYDRO_D
1244 #ifdef MPP_LAND
1245 !if(my_id .eq. io_id) then
1246 if(my_id .eq. -5) then
1247 #endif
1248 print*,'Ndg: '
1249 print*,'Ndg: !-------------------------------------------------'
1250 print*,'Ndg: obsTimeStr(tt=',tt,')'
1251 print*,'Ndg: obsTimeStr(tt)%time:', obsTimeStr(tt)%time
1252 print*,'Ndg: obsTimeStr(tt)%updateTime:', obsTimeStr(tt)%updateTime
1253 print*,'Ndg: obsTimeStr(tt)%allCellInds:', obsTimeStr(tt)%allCellInds
1254 print*,'Ndg: obsTimeStr(tt)%nGageCell:', obsTimeStr(tt)%nGageCell
1255 if(trim(obsTimeStr(tt)%updateTime(1:7)) .ne. 'no file') then
1256 print*,'Ndg: nobs:', size(obsTimeStr(tt)%obsStr)
1257 print*,'Ndg: obsTimeStr(tt)%obsStr(1)%usgsId:', obsTimeStr(tt)%obsStr(1)%usgsId
1258 print*,'Ndg: obsTimeStr(tt)%obsStr(1)%obsTime:', obsTimeStr(tt)%obsStr(1)%obsTime
1259 print*,'Ndg: obsTimeStr(tt)%obsStr(1)%obsDischarge:', obsTimeStr(tt)%obsStr(1)%obsDischarge
1260 print*,'Ndg: obsTimeStr(tt)%obsStr(1)%obsQC:', obsTimeStr(tt)%obsStr(1)%obsQC
1261 print*,'Ndg: obsTimeStr(tt)%obsStr(1)%obsStaticInd:', obsTimeStr(tt)%obsStr(1)%obsStaticInd
1262 end if
1263 print*,'Ndg: !-------------------------------------------------'
1264 print*,'Ndg: '
1265 #ifdef MPP_LAND
1266 if(flushAll) flush(flushUnit)
1267 end if
1268 #endif
1269 #endif /* HYDRO_D */
1271 !!$#ifdef HYDRO_D
1272 !!$#ifdef MPP_LAND
1273 !!$if(my_id .eq. io_id) then
1274 !!$#endif
1275 !!$ call nudging_timer(endCodeTime)
1276 !!$ call accum_nudging_time(startCodeTime, endCodeTime, &
1277 !!$ 'setup6: diagnostics at end', .false.)
1278 !!$ if(flushAll) flush(flushUnit)
1279 !!$#ifdef MPP_LAND
1280 !!$endif
1281 !!$#endif
1282 !!$#endif /* HYDRO_D */
1284 deallocate(obsTimeStrAllocated, whObsMiss)
1286 #ifdef HYDRO_D
1287 #ifdef MPP_LAND
1288 if(my_id .eq. io_id) &
1289 #endif
1290 print*,'Ndg: finish setup_stream_nudging'
1292 !!$#ifdef MPP_LAND
1293 !!$if(my_id .eq. io_id) then
1294 !!$#endif
1295 !!$ call nudging_timer(endCodeTime)
1296 !!$ call accum_nudging_time(startCodeTimeOuter, endCodeTime, &
1297 !!$ 'setup7: finish/full setup_stream_nudging', .true. )
1298 !!$#ifdef MPP_LAND
1299 !!$endif
1300 !!$#endif
1302 if(flushAll) flush(flushUnit)
1303 #endif /* HYDRO_D */
1305 end subroutine setup_stream_nudging
1307 !===================================================================================================
1308 ! Program Names:
1309 ! timeslice_file_to_struct
1310 ! Author(s)/Contact(s):
1311 ! James L McCreight <jamesmcc><ucar><edu>
1312 ! Abstract:
1313 ! Read a timeslice file, subset to observations in the domain, put into the
1314 ! obsTimeStr(timeIndex)%obsStr where obsTimeStr(timeIndex)%time impiles the
1315 ! file being read in and timeIndex is the only argument to this subroutine.
1316 ! History Log:
1317 ! 7/23/15 -Created, JLM.
1318 ! Usage:
1319 ! Parameters:
1320 ! timeIndex: the index corresponding to time in the obsTimeStr
1321 ! Input Files: Specified argument.
1322 ! Output Files: None.
1323 ! Condition codes:
1324 ! User controllable options:
1325 ! Namelist option, logical "filterObsToDom". This removes observations which are not in the
1326 ! domain. May be useful for running small domains using large data files.
1327 ! Notes:
1329 subroutine timeslice_file_to_struct(structIndex)
1330 use module_nudging_utils, only: whichInLoop2
1332 use module_nudging_io, only: find_timeslice_file, &
1333 read_timeslice_file, &
1334 get_netcdf_dim
1336 implicit none
1337 integer, intent(in) :: structIndex
1339 character(len=19) :: thisTime !! the of this time/slice of the structure
1340 character(len=256) :: fileName !! the corresponding obs file.
1342 integer :: nGages
1343 integer, allocatable, dimension(:) :: whSliceInDom
1344 integer :: nWhSliceInDom
1345 integer :: count, ww, errStatus
1347 character(len=19) :: timeIn, updateTimeIn !! the time of file/slice and when it was updated
1348 character(len=2) :: sliceResoIn !! the temporal resolution of the slice
1349 character(len=19), allocatable, dimension(:) :: gageTimeIn !!
1350 real, allocatable, dimension(:) :: gageQCIn !!
1351 character(len=15), allocatable, dimension(:) :: usgsIdIn !! USGS ID
1352 real, allocatable, dimension(:) :: gageDischargeIn !! m3/s
1354 ! Is there a timeslice file?
1355 thisTime = obsTimeStr(structIndex)%time
1356 fileName = find_timeslice_file(thisTime, obsResolution)
1358 ! If no file, note in updateTime and get out!
1359 if(fileName .eq. '') then
1360 obsTimeStr(structIndex)%updateTime='no file'
1361 #ifdef HYDRO_D
1362 print*,'Ndg: no timeSliceFile at this time: ', thisTime
1363 if(flushAll) flush(flushUnit)
1364 #endif
1365 return
1366 end if
1368 #ifdef HYDRO_D
1369 print*,'Ndg: Found file:', fileName
1370 print*,'Ndg: timeSlice: ',thisTime
1371 print*,'Ndg: timeSliceFile:',trim(fileName), my_id
1372 if(flushAll) flush(flushUnit)
1373 #endif
1375 nGages=get_netcdf_dim(fileName, 'stationIdInd', &
1376 'timeslice_file_to_struct', &
1377 readTimesliceErrFatal )
1378 ! If the dimension comes back zero, there's an issue with the file.
1379 ! If the error is not fatal, then handle the file as missing & print an
1380 ! additional message
1381 if(nGages .eq. 0) then
1382 obsTimeStr(structIndex)%updateTime='no file'
1383 print*,'Ndg: WARNING: issues with skipped timeSliceFile : ', thisTime
1384 return
1385 end if
1387 !! Reduce to just the observations in the domain or not?
1388 if(filterObsToDom) then
1390 ! Bring in the full file to intermediate local variables.
1391 !! JLM:: it would probably be more efficient to do this when readingin the netcdf files.
1392 !! This capability is probably not needed for IOC, so I'm not doing it now.
1393 allocate(usgsIdIn(nGages))
1394 allocate(gageTimeIn(nGages))
1395 allocate(gageQCIn(nGages))
1396 allocate(gageDischargeIn(nGages))
1398 call read_timeslice_file(fileName, &
1399 sanityQcDischarge, &
1400 timeIn, &
1401 updateTimeIn, &
1402 sliceResoIn, &
1403 usgsIdIn, &
1404 gageTimeIn, &
1405 gageQCIn, &
1406 gageDischargeIn, &
1407 readTimesliceErrFatal,&
1408 errStatus )
1410 if(errStatus .ne. 0) then
1411 obsTimeStr(structIndex)%updateTime='no file'
1412 print*,'Ndg: WARNING: issues with skipped timeSliceFile : ', thisTime
1413 deallocate(usgsIdIn, gageTimeIn, gageQCIn, gageDischargeIn)
1414 return
1415 end if
1417 if(timeIn .NE. obsTimeStr(structIndex)%time) &
1418 call hydro_stop('timeslice_file_to_struct: file time does not match structure')
1419 if(obsResolution .NE. sliceResoIn) &
1420 call hydro_stop('timeslice_file_to_struct: model and file timeslice resolution do not match.')
1421 !save when the file was updated
1422 obsTimeStr(structIndex)%updateTime = updateTimeIn
1424 allocate(whSliceInDom(size(usgsIdIn)))
1425 call whichInLoop2(usgsIdIn, nodeGageStr%usgsId, whSliceInDom, nWhSliceInDom)
1427 #ifdef HYDRO_D
1428 ! print*,'Ndg: usgsIdIn: ', usgsIdIn
1429 ! print*,'Ndg: nodeGageStr%usgsId', nodeGageStr%usgsId
1430 print*,'Ndg: nWhSliceInDom:', nWhSliceInDom
1431 if(flushAll) flush(flushUnit)
1432 #endif
1433 allocate(obsTimeStr(structIndex)%obsStr(nWhSliceInDom))
1434 !! because these dont get set here, give them default values.
1435 obsTimeStr(structIndex)%obsStr%obsStaticInd = 0
1436 obsTimeStr(structIndex)%obsStr%innov = -9999
1438 if(nWhSliceInDom .gt. 0) then
1439 count=1
1440 do ww=1,size(whSliceInDom)
1441 if(whSliceInDom(ww) .gt. 0) then
1442 obsTimeStr(structIndex)%obsStr(count)%usgsId = usgsIdIn(ww)
1443 obsTimeStr(structIndex)%obsStr(count)%obsTime = gageTimeIn(ww)
1444 obsTimeStr(structIndex)%obsStr(count)%obsQC = gageQCIn(ww)
1445 obsTimeStr(structIndex)%obsStr(count)%obsDischarge = gageDischargeIn(ww)
1446 count=count+1
1447 else
1448 !! the gage is not in the domain/parameter file, record that this
1449 !! gage is available but unable to be assimilated.
1450 call accumulate_nwis_not_in_RLAndParams(nwisNotRLAndParams, &
1451 nwisNotRLAndParamsCount, usgsIdIn(ww))
1452 end if
1453 end do
1454 end if
1456 deallocate(whSliceInDom)
1457 deallocate(usgsIdIn, gageTimeIn, gageQCIn, gageDischargeIn)
1459 else ! dont filterObsToDom
1461 ! not reducing the obs file to the domain, can
1462 ! bring in the full file directly to the structure.
1463 allocate(obsTimeStr(structIndex)%obsStr(nGages))
1464 obsTimeStr(structIndex)%obsStr%obsStaticInd = 0
1465 obsTimeStr(structIndex)%obsStr%innov = -9999
1467 call read_timeslice_file(fileName, &
1468 sanityQcDischarge, &
1469 timeIn, &
1470 obsTimeStr(structIndex)%updateTime, &
1471 sliceResoIn, &
1472 obsTimeStr(structIndex)%obsStr(:)%usgsId, &
1473 obsTimeStr(structIndex)%obsStr(:)%obsTime, &
1474 obsTimeStr(structIndex)%obsStr(:)%obsQC, &
1475 obsTimeStr(structIndex)%obsStr(:)%obsDischarge, &
1476 readTimesliceErrFatal, &
1477 errStatus )
1479 if(errStatus .ne. 0) then
1480 obsTimeStr(structIndex)%updateTime='no file'
1481 print*,'Ndg: WARNING: issues with skipped timeSliceFile : ', thisTime
1482 deallocate(obsTimeStr(structIndex)%obsStr)
1483 return
1484 end if
1486 if(timeIn .NE. obsTimeStr(structIndex)%time) &
1487 call hydro_stop('timeslice_file_to_struct: file time does not match that in obsTimeStr')
1488 if(obsResolution .NE. sliceResoIn) &
1489 call hydro_stop('timeslice_file_to_struct: model and file timeslice resolution do not match.')
1491 endif !filterObsToDom
1493 #ifdef HYDRO_D
1494 #ifdef MPP_LAND
1495 if(my_id .eq. io_id) &
1496 #endif
1497 print*,'Ndg: finish timeslice_file_to_struct'
1498 if(flushAll) flush(flushUnit)
1499 #endif
1501 end subroutine timeslice_file_to_struct
1503 !===================================================================================================
1504 ! Program Name:
1505 ! obs_static_to_struct
1506 ! Author(s)/Contact(s):
1507 ! James L McCreight <jamesmcc><ucar><edu>
1508 ! Abstract:
1509 ! For new obs read in from file, determine the associate static data:
1510 ! link index, parameters, cellsAffected, distances, and weights.
1511 ! History Log:
1512 ! 8/5/15 -Created, JLM.
1513 ! Usage:
1514 ! Parameters:
1515 ! Input Files: Specified argument.
1516 ! Output Files: None.
1517 ! Condition codes:
1518 ! User controllable options:
1519 ! Notes:
1521 subroutine obs_static_to_struct()
1522 use module_nudging_utils, only: whUniLoop
1524 implicit none
1525 logical, allocatable, dimension(:) :: theMask
1526 integer, dimension(maxNeighLinksDim) :: upAllInds, downAllInds ! the collected inds/links
1527 real, dimension(maxNeighLinksDim) :: upAllDists, downAllDists ! distance at each collected ind
1528 integer :: upLastInd, downLastInd ! # of inds/links collected so far
1529 integer :: ll, whLastObsStr, tt
1530 integer, allocatable, dimension(:) :: nCellsInR
1531 logical :: setLastObsInfo
1533 #ifdef HYDRO_D
1534 #ifdef MPP_LAND
1535 if(my_id .eq. io_id) &
1536 #endif
1537 print*,'Ndg: start obs_static_to_struct'
1538 if(flushAll) flush(flushUnit)
1539 #endif
1541 ! nudgingParamsStr was already reduced to the intersection of all the
1542 ! gages in the parameter file and all those defined in nodeGageStr (Route_Link)
1543 ! so this is minimal by construction.
1547 #ifdef MPP_LAND
1548 if(my_id .eq. io_id) then
1549 #endif
1550 ! have to search for the corresponding nodeId/obsCellInd. seems like this
1551 ! could have been done with making nudgingParamsStr, but it's not (necessary).
1552 allocate(theMask(size(nodeGageStr%usgsId)))
1553 ! loop over the entire structure and solve %obsCellind
1554 nudgeSpatial=.false.
1555 do ll=1,size(obsStaticStr_SoA%usgsId(:))
1556 !! this will always have a match b/c all
1557 !! obsStaticStr%usgsId=nudgingParamsStr$usgsId
1558 !! are in nodeGageStr by construction
1559 !! (if we were looping on nodeGageStr%usgsId this would not be assured)
1560 theMask = nodeGageStr%usgsId .eq. obsStaticStr_SoA%usgsId(ll)
1561 !! This is a double for loop... thanks to whUniLoop
1562 obsStaticStr_SoA%obsCellInd(ll) = nodeGageStr%nodeId(whUniLoop(theMask))
1563 if(obsStaticStr_SoA%R(ll) .ge. chanlen_image0(obsStaticStr_SoA%obsCellInd(ll))/2) &
1564 nudgeSpatial=.true.
1565 end do
1566 deallocate(theMask)
1567 #ifdef MPP_LAND
1568 end if
1569 !broadcast
1570 call mpp_land_bcast_logical(nudgeSpatial)
1571 call mpp_land_bcast_int1d(obsStaticStr_SoA%obsCellInd(:))
1572 #endif
1575 #ifdef MPP_LAND
1576 if(my_id .eq. io_id) &
1577 #endif
1578 print*,'Ndg: nudgeSpatial = ', nudgeSpatial
1580 !! get the requsite files if necessary
1581 !! this solves/allocates upNetwkStr and downNetwkStr
1582 if(nudgeSpatial) call get_netwk_reexpression()
1584 ! loop over the entire structure and solve
1585 ! 1) the spatial info: %cellsAffected, %dist, %ws
1586 ! 2) the lastObs info: %lastObsDischarge, %lastObsModelDischarge, %lastObsTime, %lastObsQuality
1587 ! 3) print diagnostics with or witout spatial info
1589 allocate(nCellsInR(size(obsStaticStr_SoA%usgsId(:))))
1590 #ifdef MPP_LAND
1591 if(my_id .eq. io_id) then
1592 #endif
1593 do ll=1,size(obsStaticStr_SoA%usgsId(:))
1595 !! spatial static info
1596 if(nudgeSpatial) then
1598 !! Now solve the spaital part of obsStaticStr
1599 !---------------------
1600 ! Calculate cells affected by this gage
1601 ! upstream to R keep distances
1602 upLastInd = 0
1603 call distance_along_channel( &
1604 upNetwkStr, & ! traversal structure in up/down direction
1605 obsStaticStr_SoA%obsCellInd(ll), & ! the starting link
1606 0.0000000000, & ! distance (m) at the starting node (this iter.)
1607 obsStaticStr_SoA%R(ll), & ! to traverse, in meters.
1608 upAllInds, & ! collected links/inds
1609 upAllDists, & ! distance to collected links
1610 upLastInd ) ! index of last collected link in above 2 arrays.
1612 ! downstream to R keep distances
1613 downLastInd = 0
1614 call distance_along_channel( &
1615 downNetwkStr, & ! traversal structure in up/down direction
1616 obsStaticStr_SoA%obsCellInd(ll), & ! the starting link
1617 0.0000000000, & ! distance (m) at the starting node (this iter.)
1618 obsStaticStr_SoA%R(ll), & ! to traverse, in meters.
1619 downAllInds, & ! collected links/inds
1620 downAllDists, & ! distance to collected links
1621 downLastInd ) ! index of last collected link in above 2 arrays.
1623 ! Collect the up and down stream cells
1624 nCellsInR(ll) = 1 + upLastInd + downLastInd
1625 allocate(obsStaticStr(ll)%cellsAffected(nCellsInR(ll)))
1626 allocate(obsStaticStr(ll)%dist(nCellsInR(ll)))
1627 allocate(obsStaticStr(ll)%ws(nCellsInR(ll)))
1629 obsStaticStr(ll)%cellsAffected(1) = obsStaticStr_SoA%obsCellInd(ll)
1630 obsStaticStr(ll)%cellsAffected(2:(upLastInd+1)) = upAllInds(1:upLastInd)
1631 obsStaticStr(ll)%cellsAffected((upLastInd+2):nCellsInR(ll)) = downAllInds(1:downLastInd)
1633 ! Distance is optional in obsStaticStr. Right now keeping it for diagnostic potential.
1634 ! May replace with a local variable later when we dont want to keep it.
1635 obsStaticStr(ll)%dist(1) = 0
1636 obsStaticStr(ll)%dist(2:(upLastInd+1)) = -1. * upAllDists(1:upLastInd)
1637 obsStaticStr(ll)%dist((upLastInd+2):nCellsInR(ll)) = downAllDists(1:downLastInd)
1639 ! Calculate the cressman spatial weights
1640 obsStaticStr(ll)%ws = &
1641 ( obsStaticStr_SoA%R(ll)**2 - obsStaticStr(ll)%dist**2 ) / &
1642 ( obsStaticStr_SoA%R(ll)**2 + obsStaticStr(ll)%dist**2 )
1644 else
1646 ! Collect the up and down stream cells
1647 nCellsInR(ll) = 1
1648 allocate(obsStaticStr(ll)%cellsAffected(nCellsInR(ll)))
1649 allocate(obsStaticStr(ll)%dist(nCellsInR(ll)))
1650 allocate(obsStaticStr(ll)%ws(nCellsInR(ll)))
1652 obsStaticStr(ll)%cellsAffected(1) = obsStaticStr_SoA%obsCellInd(ll)
1653 obsStaticStr(ll)%dist(1) = 0
1654 obsStaticStr(ll)%ws = 1
1656 end if ! nudgeSpatial
1658 if(temporalPersistence) then
1659 ! the last obs "static" info
1660 ! this puts in dummy values if none are found
1661 setLastObsInfo = .false.
1662 if(allocated(lastObsStr_SoA%usgsId)) then
1663 if(any(lastObsStr_SoA%usgsId(:) .eq. obsStaticStr_SoA%usgsId(ll))) then
1664 allocate(theMask(size(lastObsStr_SoA%usgsId(:))))
1665 theMask = lastObsStr_SoA%usgsId(:) .eq. obsStaticStr_SoA%usgsId(ll)
1666 whLastObsStr = whUniLoop(theMask)
1667 deallocate(theMask)
1668 !! Drop in the entire derived type for each array index! Very nice.
1669 obsStaticStr_SoA%lastObsStructure_SoA%usgsId(ll) = lastObsStr_SoA%usgsId(whLastObsStr)
1670 obsStaticStr_SoA%lastObsStructure_SoA%lastObsDischarge(:,ll) = lastObsStr_SoA%lastObsDischarge(:,whLastObsStr)
1671 obsStaticStr_SoA%lastObsStructure_SoA%lastObsTime(:,ll) = lastObsStr_SoA%lastObsTime(:,whLastObsStr)
1672 obsStaticStr_SoA%lastObsStructure_SoA%lastObsQuality(:,ll) = lastObsStr_SoA%lastObsQuality(:,whLastObsStr)
1673 obsStaticStr_SoA%lastObsStructure_SoA%lastObsModelDischarge(:,ll) = lastObsStr_SoA%lastObsModelDischarge(:,whLastObsStr)
1675 setLastObsInfo = .true.
1676 endif
1677 endif
1679 if(.not. setLastObsInfo) then ! the gage did not have a last observation
1680 obsStaticStr_SoA%lastObsDischarge(:,ll) = real(-9999)
1681 obsStaticStr_SoA%lastObsModelDischarge(:,ll) = real(-9999)
1682 do tt=1,nTimesLastObs
1683 obsStaticStr_SoA%lastObsTime(tt,ll) = missingLastObsTime
1684 end do
1685 obsStaticStr_SoA%lastObsQuality(:,ll) = real(0) !! keep this zero
1686 end if
1687 end if !temporalPersistence
1689 #ifdef HYDRO_D
1690 ! diagnostics
1691 if (my_id .eq. -5) then
1692 ! 1 2 3 4 5 6 7 8 9 10
1693 !if( any( (/ 136, 981, 2242, 2845, 2946, 3014, 3066, 3068, 3072, 3158, 3228, &
1694 ! 3325, 3328, 3330, 3353, 3374, 3398, 3416, 3445, 3487, 3536, &
1695 ! 3621, 3667, 13554, 13651, 15062, 15311, 17450, 17394 /) &
1696 ! .eq. obsStaticStr(ll)%obsCellInd ) )then
1697 print*,'Ndg: '
1698 print*,'Ndg: !----------------------------------'
1699 print*,'Ndg: ! obsStaticStr(',ll,'):'
1700 print*,'Ndg: obsStaticStr(ll)%usgsId: ', obsStaticStr_SoA%usgsId(ll)
1701 print*,'Ndg: obsStaticStr(ll)%obsCellInd: ', obsStaticStr_SoA%obsCellInd(ll)
1702 print*,'Ndg: obsStaticStr(ll)%R: ', obsStaticStr_SoA%R(ll)
1703 print*,'Ndg: obsStaticStr(ll)%G: ', obsStaticStr_SoA%G(ll)
1704 print*,'Ndg: obsStaticStr(ll)%tau: ', obsStaticStr_SoA%tau(ll)
1705 print*,'Ndg: obsStaticStr(ll)%cellsAffected: ', obsStaticStr(ll)%cellsAffected
1706 print*,'Ndg: obsStaticStr(ll)%dist: ', obsStaticStr(ll)%dist
1707 print*,'Ndg: obsStaticStr(ll)%ws: ', obsStaticStr(ll)%ws
1708 if(temporalPersistence) then
1709 print*,'Ndg: obsStaticStr(ll)%lastObsDischarge: ', obsStaticStr_SoA%lastObsDischarge(:,ll)
1710 print*,'Ndg: obsStaticStr(ll)%lastObsDischarge: ', obsStaticStr_SoA%lastObsModelDischarge(:,ll)
1711 print*,'Ndg: obsStaticStr(ll)%lastObsTime: ', obsStaticStr_SoA%lastObsTime(:,ll)
1712 print*,'Ndg: obsStaticStr(ll)%lastObsQuality: ', obsStaticStr_SoA%lastObsQuality(:,ll)
1713 endif
1714 print*,'Ndg: !----------------------------------'
1715 print*,'Ndg: '
1716 endif ! my_id .eq. -5 (off) or any of the diagnostic inds - toggle the comments.
1717 if(flushAll) flush(flushUnit)
1718 #endif /* HYDRO_D */
1720 end do ! ll
1721 #ifdef MPP_LAND
1722 end if !! my_id .eq. io_id
1724 !broadcast %cellsAffected, %dist, %ws
1725 do ll=1,size(obsStaticStr_SoA%usgsId)
1726 call mpp_land_bcast_int1(nCellsInR(ll))
1727 if(my_id .ne. io_id) then
1728 allocate(obsStaticStr(ll)%cellsAffected(nCellsInR(ll)))
1729 allocate(obsStaticStr(ll)%dist(nCellsInR(ll)))
1730 allocate(obsStaticStr(ll)%ws(nCellsInR(ll)))
1731 endif
1732 call mpp_land_bcast_int1d(obsStaticStr(ll)%cellsAffected)
1733 call mpp_land_bcast_real_1d(obsStaticStr(ll)%dist)
1734 call mpp_land_bcast_real_1d(obsStaticStr(ll)%ws)
1736 if(temporalPersistence) then
1737 call mpp_land_bcast_real_1d(obsStaticStr_SoA%lastObsDischarge(:,ll))
1738 call mpp_land_bcast_real_1d(obsStaticStr_SoA%lastObsModelDischarge(:,ll))
1739 call mpp_land_bcast_char1d( obsStaticStr_SoA%lastObsTime(:,ll))
1740 call mpp_land_bcast_real_1d(obsStaticStr_SoA%lastObsQuality(:,ll))
1741 endif
1742 end do
1743 #endif /* MPP_LAND */
1745 deallocate(nCellsInR)
1747 #ifdef HYDRO_D
1748 print*,'Finnish obs_static_to_struct'
1749 if(flushAll) flush(flushUnit)
1750 #endif
1752 end subroutine obs_static_to_struct
1755 !===================================================================================================
1756 ! Subroutine Name:
1757 ! subroutine output_nudging_last_obs
1758 ! Author(s)/Contact(s):
1759 ! James L McCreight <jamesmcc><ucar><edu>
1760 ! Abstract:
1761 ! Collect and Write out the last observations collected over time.
1762 ! History Log:
1763 ! 02/09/16 -Created, JLM.
1764 ! Usage:
1765 ! Parameters:
1766 ! Input Files:
1767 ! Output Files:
1768 ! Condition codes:
1769 ! User controllable options: None.
1770 ! Notes: Needs better error handling...
1772 subroutine output_nudging_last_obs()
1773 use module_RT_data, only: rt_domain
1774 use module_nudging_utils, only: whUniLoop
1775 use module_nudging_io, only: write_nudging_last_obs, write_nudging_last_obs_soa
1776 #ifdef MPP_LAND
1777 use MODULE_mpp_ReachLS, only: linkls_s, linkls_e
1778 implicit none
1779 real, allocatable, dimension(:) :: g_nudge
1780 logical, allocatable, dimension(:) :: theMask
1781 integer :: whImage, oo
1782 #endif
1784 if(.not. temporalPersistence) return
1786 #ifdef MPP_LAND
1787 !! if MPP: last obs are being written on different processors,
1788 !! get them back to image0 for output.
1789 !! 1. loop over obsStaticStr.
1790 !! 2. find the image holding obsCellInd
1791 !! 3. communicate it's last ob back to image0.
1792 allocate(theMask(size(linkls_s))) ! number of processors
1794 do oo=1, size(obsStaticStr_SoA%obsCellInd(:)) ! loop over the gages
1796 ! Is this gage in the range for this image?
1797 theMask = ( obsStaticStr_SoA%obsCellInd(oo) .ge. linkls_s .and. &
1798 obsStaticStr_SoA%obsCellInd(oo) .le. linkls_e )
1799 ! The image where the gage info is located
1800 whImage = whUniLoop(theMask) - 1
1802 ! if the data is on this image or this image is the receiver, then communicate.
1803 if(my_id .ne. whImage .and. my_id .ne. io_id) cycle
1804 ! no communication necessary on the io image
1805 if(whImage .eq. io_id) cycle
1807 call mpp_comm_1d_real(obsStaticStr_SoA%lastObsDischarge(:,oo), whImage, io_id)
1808 call mpp_comm_1d_real(obsStaticStr_SoA%lastObsModelDischarge(:,oo), whImage, io_id)
1809 call mpp_comm_1d_char(obsStaticStr_SoA%lastObsTime(:,oo), whImage, io_id)
1810 call mpp_comm_1d_real(obsStaticStr_SoA%lastObsQuality(:,oo), whImage, io_id)
1812 end do
1814 deallocate(theMask)
1816 if(my_id .eq. io_id) then
1817 allocate(g_nudge(rt_domain(1)%gnlinksl))
1818 else
1819 allocate(g_nudge(1))
1820 end if
1821 call ReachLS_write_io(RT_DOMAIN(1)%nudge,g_nudge)
1823 if(my_id .eq. io_id) &
1824 #endif /* MPP_LAND */
1825 !! write out the last obs to file
1826 call write_nudging_last_obs_soa(obsStaticStr_SoA%lastObsStructure_SoA, &
1827 nlst(did)%olddate, &
1828 g_nudge )
1830 deallocate(g_nudge)
1832 end subroutine output_nudging_last_obs
1835 !===================================================================================================
1836 ! Program Name:
1837 ! accumulate_nwis_not_in_RLAndParams
1838 ! Author(s)/Contact(s):
1839 ! James L McCreight <jamesmcc><ucar><edu>
1840 ! Abstract:
1841 ! Keep a running, non-redundant log of gages seen in the forecast cycle
1842 ! (from nwis) which were not found in the intersection of the Route_Link
1843 ! gages and the parameter file gages.
1844 ! History Log:
1845 ! 11/4/15 - Created, JLM.
1846 ! Usage:
1847 ! Parameters:
1848 ! Input Files:
1849 ! Output Files:
1850 ! Condition codes:
1851 ! User controllable options:
1852 ! Notes:
1854 subroutine accumulate_nwis_not_in_RLAndParams(nwisNotRLAndParams_local, &
1855 nwisNotRLAndParamsCount_local, &
1856 gageId )
1857 implicit none
1858 character(len=15), dimension(:), intent(inout) :: nwisNotRLAndParams_local
1859 integer, intent(inout) :: nwisNotRLAndParamsCount_local
1860 character(len=15), intent(in) :: gageId
1861 if(.not. any(nwisNotRLAndParams_local .eq. gageId)) then
1862 nwisNotRLAndParamsCount_local=nwisNotRLAndParamsCount_local+1
1863 if(nwisNotRLAndParamsCount .gt. size(nwisNotRLAndParams)) &
1864 call hydro_stop('accumulate_nwis_not_in_RLAndParams: coutn of gages from NWIS not found ' // &
1865 'in intersection of RL and Param gages are exceeding the hardcoded limit:,'// &
1866 ' maxNwisNotRLAndParamsCount')
1867 nwisNotRLAndParams_local(nwisNotRLAndParamsCount_local)=gageId
1868 end if
1869 end subroutine accumulate_nwis_not_in_RLAndParams
1872 !===================================================================================================
1873 ! Program Name:
1874 ! distance_along_channel
1875 ! Author(s)/Contact(s):
1876 ! James L McCreight <jamesmcc><ucar><edu>
1877 ! Abstract:
1878 ! Traverse the gridded channel network up or down stream, return the cumulative distance
1879 ! and indices of points visited. Going from one link to the next counts half the length
1880 ! of each as the distance. Conceptually that's midpoint to midpoint though there's really
1881 ! no midpoint. If half the length to the first midpoint exceeds R, no index is returned.
1882 ! That is, lastInd is zero.
1884 ! History Log:
1885 ! 6/4/15 - Created, JLM.
1886 ! 8/7/15 - Heavily refactored to remove searching, JLM.
1887 ! Usage:
1888 ! Parameters:
1889 ! See formal arguments and their declarations
1890 ! Input Files: Specified argument.
1891 ! Output Files: None.
1892 ! Condition codes:
1893 ! User controllable options:
1894 ! Notes:
1895 ! The total number of links gathered is hardwired to 10,000 which is for both directions. In NHD+
1896 ! the shortest link is 1m. If you managed to traverse just the shortest 10,000 links then you only
1897 ! go [R language: > sum(head(sort(reInd$length),10000))/1000 = 60.0745] 60km. This is the minimum
1898 ! upper bound on R implied by the choice of 10,000. If lastInd ever becomes 10,001 a fatal error
1899 ! is issued.
1901 recursive subroutine distance_along_channel(direction, & ! traversal structure in up/down direction
1902 startInd, & ! the starting link
1903 startDist, & ! distance at the starting node (this iter.)
1904 radius, & ! to traverse, in meters.
1905 allInds, & ! collected links/inds
1906 allDists, & ! distance to collected links
1907 lastInd ) ! index of last collected link in above.
1908 use module_RT_data, only: rt_domain
1910 implicit none
1911 type(netwkReExpStructure), intent(in) :: direction ! up/down NtwkStr with tarversal inds, pass by ref!
1912 integer, intent(in) :: startInd ! the starting link
1913 real, intent(in) :: startDist ! distance at the starting node (for this iteration)
1914 real, intent(in) :: radius ! in meteres
1915 integer, intent(inout), dimension(maxNeighLinksDim) :: allInds ! the collected inds/links
1916 real, intent(inout), dimension(maxNeighLinksDim) :: allDists ! the distance at each collected ind
1917 integer, intent(inout) :: lastInd ! the number of inds/links collected so far.
1919 ! whGo is only > 1 upstream with a current max of 17 in NHD+
1920 integer :: go, nGo
1921 integer :: gg
1922 real :: newDist
1923 integer, parameter :: did=1
1925 !! this routine is only called on io_id
1927 #ifdef HYDRO_D
1928 #ifdef MPP_LAND
1929 if(my_id .eq. io_id) &
1930 #endif
1931 print*,'Ndg: start distance_along_channel'
1932 if(flushAll) flush(flushUnit)
1933 #endif
1935 if(direction%end(startInd) .eq. 0) return ! a pour point (downstream) or a 1st order link (upstream)
1937 nGo = direction%end(startInd) - direction%start(startInd)
1938 nGo = nGo + 1 ! end-start+1 if end-start > 0
1939 do gg=0,nGo-1
1940 go = direction%go( direction%start(startInd) + gg )
1941 allInds(lastInd+1) = go
1942 #ifdef MPP_LAND
1943 newDist = ( chanlen_image0(startInd) + chanlen_image0(go) ) / 2.
1944 #else
1945 newDist = ( rt_domain(did)%chanlen(startInd) + rt_domain(did)%chanlen(go) ) / 2.
1946 #endif
1947 if(startDist + newDist .gt. radius) return ! strictly greater than.
1948 allDists(lastInd+1) = startDist + newDist
1949 lastInd = lastInd+1
1950 if(lastInd .ge. 10001) &
1951 call hydro_stop('distance_along_channel: hardwired 10000 variable size exceeded. FIX.')
1952 call distance_along_channel( &
1953 direction, & ! the traversal structure, pass by reference.
1954 go, & ! the new startInd is where we go from the old startInd
1955 startDist+newDist, & ! a little bit further is where we start the next call.
1956 radius, & ! static
1957 allInds, & ! collected inds
1958 allDists, & ! collected dists
1959 lastInd ) ! the number collected so far
1960 end do
1962 #ifdef HYDRO_D
1963 print*,'Ndg: end distance_along_channel'
1964 if(flushAll) flush(flushUnit)
1965 #endif
1967 end subroutine distance_along_channel
1969 !===================================================================================================
1970 ! Program Name:
1971 ! tally_affected_links
1972 ! Author(s)/Contact(s):
1973 ! James L McCreight <jamesmcc><ucar><edu>
1974 ! Abstract:
1975 ! History Log:
1976 ! 8/11/15
1977 ! Usage: call tally_affected_links(tt)
1978 ! Parameters:
1979 ! Input Files: None.
1980 ! Output Files: None.
1981 ! Condition codes:
1982 ! User controllable options:
1983 ! Notes:
1985 subroutine tally_affected_links(timeIndex)
1986 use module_RT_data, only: rt_domain
1988 implicit none
1989 integer, intent(in) :: timeIndex
1990 integer, allocatable, dimension(:) :: affectedInds, nGageAffect
1991 integer :: nlinks, ii, oo, cc, nLinkAff, indAff, staticInd
1992 integer, parameter :: did=1
1994 if (nlst(did)%channel_option .eq. 3) nLinks = RT_DOMAIN(did)%NLINKS ! For gridded channel routing
1995 if (nlst(did)%channel_option .eq. 1 .or. &
1996 nlst(did)%channel_option .eq. 2 ) &
1997 nLinks = &
1998 #ifdef MPP_LAND
1999 RT_DOMAIN(did)%gNLINKSL ! For reach-based routing in parallel
2000 #else
2001 RT_DOMAIN(did)%NLINKSL ! For reach-based routing
2002 #endif
2004 allocate(affectedInds(nLinks), nGageAffect(nLinks))
2006 affectedInds = 0
2007 nGageAffect = 0
2009 ! parallelize this ||||||||||||||||||||||||||||||||||||||||
2010 do oo=1,size(obsTimeStr(timeIndex)%obsStr)
2011 staticInd = obsTimeStr(timeIndex)%obsStr(oo)%obsStaticInd
2012 if(staticInd .eq. 0) cycle
2013 nLinkAff = size(obsStaticStr(staticInd)%cellsAffected)
2014 do cc=1, nLinkAff
2015 indAff = obsStaticStr(staticInd)%cellsAffected(cc)
2016 affectedInds(indAff) = indAff
2017 nGageAffect(indAff) = nGageAffect(indAff) + 1
2018 end do
2019 end do
2021 ! how many total cells affected?
2022 nLinkAff = sum((/(1, ii=1,nLinks)/), mask=affectedInds .ne. 0)
2024 allocate(obsTimeStr(timeIndex)%allCellInds(nLinkAff))
2025 allocate(obsTimeStr(timeIndex)%nGageCell(nLinkAff))
2027 obsTimeStr(timeIndex)%allCellInds = pack(affectedInds, mask=affectedInds .ne. 0)
2028 obsTimeStr(timeIndex)%nGageCell = pack(nGageAffect, mask=nGageAffect .ne. 0)
2030 deallocate(affectedInds, nGageAffect)
2032 #ifdef HYDRO_D
2033 print*,'Ndg: end tally_affected_links'
2034 if(flushAll) flush(flushUnit)
2035 #endif
2037 end subroutine tally_affected_links
2040 !===================================================================================================
2041 ! program name:
2042 ! time_wt
2043 ! author(s)/contact(s):
2044 ! James McCreight
2045 ! abstract:
2046 ! compute the temporal weight factor for an observation
2047 ! history log:
2048 ! 6/4/15 - created,
2049 ! 12/1/15 - moved to inverse distance
2050 ! usage:
2051 ! parameters:
2052 ! q
2053 ! input files:
2054 ! output files:
2055 ! condition codes:
2056 ! user controllable options:
2057 ! notes:
2059 real function time_wt(modelTime, tau, obsTime)
2060 use module_date_utils_nudging, only: geth_idts
2061 implicit none
2062 character(len=19), intent(in) :: modelTime ! model time string
2063 real, intent(in) :: tau ! tau half window (minutes)
2064 character(len=19), intent(in) :: obsTime ! observation time string
2065 ! local variables
2066 integer :: timeDiff
2067 real, parameter :: ten = 10.000000000
2068 real, parameter :: one = 1.000000000
2069 real, parameter :: zero = 0.000000000
2070 integer, parameter :: zeroInt = 0
2071 !! returned timeDiff is in seconds
2072 call geth_idts(obsTime, modelTime, timeDiff)
2073 timeDiff = timeDiff / 60. ! minutes to match tau
2075 !! this is the old, ramped, wrf-style time weighting.
2076 !time_wt = 0. ! = if(abs(timeDiff) .gt. tau)
2077 !if(abs(timeDiff) .lt. tau/2.) time_wt = 1.
2078 !if(abs(timeDiff) .ge. tau/2. .and. &
2079 ! abs(timeDiff) .le. tau ) time_wt = (tau - abs(timeDiff)) / (tau/2.)
2081 !this is the new, inverse distance weighting
2082 if(abs(timeDiff) .gt. tau) then
2083 time_wt = zero
2084 return
2085 end if
2086 !! this function has a range of [1, 10^10] so that
2087 !! the variable "weighting" in the function 'nudge_term_link' is
2088 !! assuredly greater than 1e-8 if weights are calculated.
2089 time_wt = (ten ** ten) * ((one/ten) ** (abs(timeDiff)/(tau/ten)))
2091 end function time_wt
2093 !!$!===================================================================================================
2094 !!$! Program Name:
2095 !!$! x
2096 !!$! Author(s)/Contact(s):
2097 !!$! y
2098 !!$! Abstract:
2099 !!$! z
2100 !!$! History Log:
2101 !!$! 6/4/15 - Created,
2102 !!$! Usage:
2103 !!$! Parameters:
2104 !!$! q
2105 !!$! Input Files:
2106 !!$! Output Files:
2107 !!$! Condition codes:
2108 !!$! User controllable options:
2109 !!$! Notes:
2110 !!$real function cress(rr,rrmax)
2111 !!$implicit none
2112 !!$real,intent(in):: rr,rrmax
2113 !!$if(rr .ge. rrmax)then
2114 !!$ cress = 0.
2115 !!$else
2116 !!$ cress = (rrmax - rr)/(rrmax + rr)
2117 !!$endif
2118 !!$end function cress
2121 !===================================================================================================
2122 ! Program Name:
2123 ! x
2124 ! Author(s)/Contact(s):
2125 ! Wu YH
2126 ! James McCreight >jamesmcc<>at<>ucar<>dot<>edu<
2127 ! Abstract:
2128 ! Calculate the innovations of obs.
2129 ! History Log:
2130 ! 2015.07.15 - Created.
2131 ! 2015.09.19
2132 ! Usage:
2133 ! Parameters:
2134 ! q
2135 ! Input Files:
2136 ! Output Files:
2137 ! Condition codes:
2138 ! User controllable options:
2139 ! Notes:
2141 subroutine nudge_innov(discharge)
2142 implicit none
2143 real, dimension(:,:), intent(in) :: discharge !! modeled discharge (m3/s)
2144 integer :: tt, oo, nObsTt, linkInd, staticInd
2145 if(.not. allocated(obsTimeStr)) return
2146 do tt=1,size(obsTimeStr)
2147 if(.not. allocated(obsTimeStr(tt)%obsStr)) cycle
2148 nObsTt = size(obsTimeStr(tt)%obsStr)
2149 do oo=1,nObsTt
2150 staticInd = obsTimeStr(tt)%obsStr(oo)%obsStaticInd
2151 linkInd = obsStaticStr_SoA%obsCellInd(staticInd)
2152 obsTimeStr(tt)%obsStr(oo)%innov = &
2153 ( obsTimeStr(tt)%obsStr(oo)%obsDischarge - discharge(linkInd,2) ) &
2154 *obsTimeStr(tt)%obsStr(oo)%obsQC
2155 enddo
2156 end do
2157 end subroutine nudge_innov
2159 !!$
2160 !===================================================================================================
2161 ! Program Name:
2162 ! nudge_term_all
2163 ! Author(s)/Contact(s):
2164 ! Wu YH, James McCreight
2165 ! Abstract:
2166 ! Calculate the nudging term for one
2167 ! History Log:
2168 ! 2015.07.21 - Created,
2169 ! Usage:
2170 ! Parameters:
2171 ! q
2172 ! Input Files:
2173 ! Output Files:
2174 ! Condition codes:
2175 ! User controllable options:
2176 ! Notes:
2178 subroutine nudge_term_all(discharge, nudgeAdj, hydroAdv)
2179 use module_RT_data, only: rt_domain
2180 use module_nudging_utils, only: nudging_timer, &
2181 accum_nudging_time
2182 #ifdef MPP_LAND
2183 use MODULE_mpp_ReachLS, only: linkls_s, linkls_e, gNLinksL, ReachLS_write_io
2184 #endif
2186 implicit none
2187 real, dimension(:,:), intent(inout) :: discharge !! modeled discharge (m3/s)
2188 real, dimension(:), intent(out) :: nudgeAdj !! nudge to modeled discharge (m3/s)
2189 integer, intent(in) :: hydroAdv !! number of seconds the channel model has advanced
2191 real, allocatable, dimension(:) :: global_discharge !! modeled discharge (m3/s)
2193 integer :: ll, startInd, endInd, checkInd
2194 real :: theNudge
2196 !!$#ifdef HYDRO_D /* un-ifdef to get timing results */
2197 !!$real :: startCodeTimeAcc, endCodeTimeAcc
2198 !!$#ifdef MPP_LAND
2199 !!$if(my_id .eq. io_id) &
2200 !!$#endif /* MPP_LAND */
2201 !!$ call nudging_timer(startCodeTimeAcc)
2202 !!$if(flushAll) flush(flushUnit)
2203 !!$#endif /* HYDRO_D */ /* un-ifdef to get timing results */
2205 #ifdef MPP_LAND
2206 !! get global discharge. This is only needed if spatial interpolation.
2207 !! eventually logic: if spatialNudging determines use of global_discharge
2208 allocate(global_discharge(gNLinksL))
2209 call ReachLS_write_io(discharge(:,2), global_discharge)
2210 call mpp_land_bcast_real_1d(global_discharge)
2211 !! need to transform passed index to appropriate index for image
2212 startInd=linkls_s(my_id+1)
2213 endInd =linkls_e(my_id+1)
2214 #else
2215 startInd=1
2216 endInd =RT_DOMAIN(did)%NLINKSL
2217 #endif
2219 if(.not. gotT0Discharge) then
2220 allocate(t0Discharge(size(discharge(:,1))))
2221 t0Discharge = discharge(:,1)
2222 gotT0Discharge = .true.
2223 end if
2225 do ll=startInd, endInd ! ll is in the index of the global Route_Link
2227 #ifdef HYDRO_D
2228 checkInd = -9999 ! 136 !2569347 !2139306 !213095 ! 2211014 ! see below as well
2229 if(ll .eq. checkInd) then
2230 print*,'Ndg: checkInd: -------------------------'
2231 print*,"Ndg: checkInd: checkInd: ",checkInd
2232 print*,'Ndg: checkInd: discharge(ll,2) before:',discharge(ll-startInd+1,2)
2233 if(flushAll) flush(flushUnit)
2234 end if
2235 #endif
2236 theNudge = nudge_term_link(ll-startInd+1, hydroAdv, global_discharge)
2237 discharge(ll-startInd+1,2) = discharge(ll-startInd+1,2) + theNudge
2238 #ifdef HYDRO_D
2239 if(ll .eq. checkInd) then
2240 print*,'Ndg: checkInd: theNudge:',theNudge
2241 print*,'Ndg: checkInd: discharge(ll,2) after:',discharge(ll-startInd+1,2)
2242 if(flushAll) flush(flushUnit)
2243 endif
2244 #endif
2245 end do
2246 !! the following only valid since nudge_term_all is applied after the time step.
2247 nudgeAdj=discharge(:,2)-discharge(:,1)
2248 discharge(:,1)=discharge(:,2)
2249 deallocate(global_discharge)
2251 #ifdef HYDRO_D /* un-ifdef to get timing results */
2252 !!$#ifdef MPP_LAND
2253 !!$if(my_id .eq. io_id) then
2254 !!$#endif
2255 !!$ call nudging_timer(endCodeTimeAcc)
2256 !!$ call accum_nudging_time(startCodeTimeAcc, endCodeTimeAcc, 'nudge_term_all', .true.)
2257 !!$#ifdef MPP_LAND
2258 !!$endif
2259 !!$#endif
2260 #ifdef MPP_LAND
2261 if(my_id .eq. io_id) &
2262 #endif
2263 print*,'Ndg: finish nudge_term_all'
2264 if(flushAll) flush(flushUnit)
2265 #endif /* HYDRO_D */ /* un-ifdef to get timing results */
2267 end subroutine nudge_term_all
2270 !===================================================================================================
2271 ! Program Name:
2272 ! nudge_term_link
2273 ! Author(s)/Contact(s):
2274 ! Wu YH,
2275 ! James L McCreight jamesmcc><at><ucar><dot><edu
2276 ! Abstract:
2277 ! Calculate the nudging term for one model cell
2278 ! History Log:
2279 ! 2015.07.21 - Created,
2280 ! Usage:
2281 ! Parameters:
2282 ! q
2283 ! Input Files:
2284 ! Output Files:
2285 ! Condition codes:
2286 ! User controllable options:
2287 ! Notes:
2288 ! note well that the weighting term has to sum to more than 1.E-8
2289 ! or the weighting is ignored.
2291 function nudge_term_link(linkIndIn, hydroAdv, discharge)
2292 use module_nudging_utils, only: whUniLoop
2293 use module_date_utils_nudging, only: geth_newdate, geth_idts
2294 #ifdef MPP_LAND
2295 use MODULE_mpp_ReachLS, only: linkls_s, linkls_e
2296 #endif
2298 implicit none
2299 integer, intent(in) :: linkIndIn !! stream cell index (local)
2300 integer, intent(in) :: hydroAdv !! number of seconds the channel model has advanced
2301 real, dimension(:), intent(in) :: discharge !! modeled discharge (m3/s)
2302 real :: nudge_term_link
2304 logical :: persistBias_local
2305 integer :: linkInd
2306 character(len=19) :: hydroTime !! actual time of the routing model
2307 real :: theInnov, weighting, theBias
2308 integer :: tt, oo, ll, jj, nObsTt, staticInd, obsInd
2309 real, parameter :: smallestWeight=1.e-8
2310 !! persistence related variables in following block
2311 integer :: whGageId, flowMonth, whThresh, lastObsDt, prstDtErrInt, obsDtInt, pairCount
2312 logical, dimension(nTimesLastObs) :: biasMask, whObsTooOld
2313 real, dimension(nTimesLastObs) :: biasWeights
2314 real :: prstB, prstDtErr, prstErr, biasWeightSum, coefVarX, coefVarY
2315 character(len=15) :: prstGageId
2316 integer :: nCorr
2317 real :: corr, sumx, sumx2, sumxy, sumy, sumy2, xCorr, yCorr, deltaT0Discharge
2318 integer, parameter :: sixty=60
2319 real, parameter :: smallestNudge=1.e-10
2320 real, parameter :: zero=1.e-37
2321 real, parameter :: zeroInt=0
2322 real, parameter :: one= 1.00000000000000000000000
2324 integer :: whPrstParams
2325 logical, allocatable, dimension(:) :: theMask
2326 #ifdef MPP_LAND
2327 integer :: whImage, checkInd
2329 !! need to transform passed index to appropriate index for image
2330 linkInd = linkIndIn + linkls_s(my_id+1) - 1
2331 #else
2332 linkInd = linkIndIn
2333 #endif
2334 checkInd = -9999 ! 136 !2569347 !2139306 !2565959 !213095!2211014 ! see above
2336 call geth_newdate(hydroTime, lsmTime, hydroAdv)
2338 nudge_term_link = zero
2339 weighting = zero
2341 #ifdef HYDRO_D
2342 if(linkInd .eq. checkInd) then
2343 print*,'Ndg: checkInd: hydroTime: ',hydroTime
2344 print*,'Ndg: checkInd: linkInd', linkInd
2345 if(flushAll) flush(flushUnit)
2346 end if
2347 #endif
2349 !! JLM: document what is this loop doing?
2350 !! 1) Moving the observation structure in to nudginLastObs struct.
2351 !! 2) Calculating the nudging terms & weights at different positions along the assim window.
2352 do tt=1,size(obsTimeStr)
2354 if(.not. allocated(obsTimeStr(tt)%allCellInds)) cycle
2355 if(.not. any(obsTimeStr(tt)%allCellInds .eq. linkInd)) cycle
2356 if(.not. allocated(obsTimeStr(tt)%obsStr)) cycle
2358 do oo=1,size(obsTimeStr(tt)%obsStr)
2360 !! if no spatial interp... this could be sped up, there's at most one match.
2361 !! this might just be an exit statement at the end of the loop for non-spatial nudges.
2362 staticInd = obsTimeStr(tt)%obsStr(oo)%obsStaticInd
2363 if(.not. any(obsStaticStr(staticInd)%cellsAffected .eq. linkInd)) cycle
2365 if(temporalPersistence) then
2366 ! If the model time is greater than or at the obs time AND
2367 ! the obs time is greater than the last written obs: then write to last obs.
2368 call geth_idts(hydroTime, obsTimeStr(tt)%obsStr(oo)%obsTime, prstDtErrInt)
2370 if(missingLastObsTime .eq. &
2371 obsStaticStr_SoA%lastObsTime(nTimesLastObs,staticInd)) then
2372 obsDtInt=1
2373 else
2374 call geth_idts(obsTimeStr(tt)%obsStr(oo)%obsTime, &
2375 obsStaticStr_SoA%lastObsTime(nTimesLastObs,staticInd), obsDtInt)
2376 end if
2378 if(prstDtErrInt .ge. 0 .and. &
2379 obsDtInt .gt. 0 .and. &
2380 obsTimeStr(tt)%obsStr(oo)%obsQC .eq. 1) then
2381 !! shift (note cshift does not work for all these)
2382 obsStaticStr_SoA%lastObsTime(1:(nTimesLastObs-1),staticInd) = &
2383 obsStaticStr_SoA%lastObsTime(2:nTimesLastObs,staticInd)
2384 obsStaticStr_SoA%lastObsDischarge(1:(nTimesLastObs-1),staticInd) = &
2385 obsStaticStr_SoA%lastObsDischarge(2:nTimesLastObs,staticInd)
2386 obsStaticStr_SoA%lastObsModelDischarge(1:(nTimesLastObs-1),staticInd) = &
2387 obsStaticStr_SoA%lastObsModelDischarge(2:nTimesLastObs,staticInd)
2388 obsStaticStr_SoA%lastObsQuality(1:(nTimesLastObs-1),staticInd) = &
2389 obsStaticStr_SoA%lastObsQuality(2:nTimesLastObs,staticInd)
2391 obsStaticStr_SoA%lastObsTime(nTimesLastObs,staticInd) = &
2392 obsTimeStr(tt)%obsStr(oo)%obsTime
2393 obsStaticStr_SoA%lastObsDischarge(nTimesLastObs,staticInd) = &
2394 obsTimeStr(tt)%obsStr(oo)%obsDischarge
2395 obsStaticStr_SoA%lastObsModelDischarge(nTimesLastObs,staticInd) = &
2396 discharge(obsStaticStr_SoA%obsCellInd(staticInd))
2397 obsStaticStr_SoA%lastObsQuality(nTimesLastObs,staticInd) = &
2398 obsTimeStr(tt)%obsStr(oo)%obsQC
2399 end if
2400 endif !temporalPersistence
2402 allocate(theMask(size(obsStaticStr(staticInd)%cellsAffected)))
2403 theMask = obsStaticStr(staticInd)%cellsAffected .eq. linkInd
2404 ll = whUniLoop(theMask)
2405 deallocate(theMask)
2407 obsInd = obsStaticStr_SoA%obsCellInd(staticInd)
2409 theInnov = &
2410 ( obsTimeStr(tt)%obsStr(oo)%obsDischarge - discharge(obsInd) ) &
2411 *obsTimeStr(tt)%obsStr(oo)%obsQC
2413 nudge_term_link = nudge_term_link &
2414 +theInnov &
2415 *obsStaticStr_SoA%G(staticInd) &
2416 *time_wt(hydroTime, &
2417 obsStaticStr_SoA%tau(staticInd), &
2418 obsTimeStr(tt)%obsStr(oo)%obsTime) &
2419 *time_wt(hydroTime, &
2420 obsStaticStr_SoA%tau(staticInd), &
2421 obsTimeStr(tt)%obsStr(oo)%obsTime) &
2422 *obsStaticStr(staticInd)%ws(ll) &
2423 *obsStaticStr(staticInd)%ws(ll)
2425 !! note well that the weighting has to sum to more than 1.E-8
2426 !! or the whole nudge is ignored.
2427 weighting = weighting &
2428 +time_wt(hydroTime, &
2429 obsStaticStr_SoA%tau(staticInd), &
2430 obsTimeStr(tt)%obsStr(oo)%obsTime) &
2431 *time_wt(hydroTime, &
2432 obsStaticStr_SoA%tau(staticInd), &
2433 obsTimeStr(tt)%obsStr(oo)%obsTime) &
2434 *obsStaticStr(staticInd)%ws(ll) &
2435 *obsStaticStr(staticInd)%ws(ll)
2437 #ifdef HYDRO_D
2438 if(linkInd .eq. checkInd) then
2439 !if( (obsTimeStr(tt)%time .eq. obsTimeStr(tt)%obsStr(oo)%obsTime) .and. &
2440 ! ( abs(discharge(obsInd)+theInnov - obsTimeStr(tt)%obsStr(oo)%obsDischarge) .gt. .01) ) then
2441 print*,'Ndg: checkInd: -*-*-*-*-*-*-*-*-*-*-*-*'
2442 print*,'Ndg: checkInd: linkInd: ', linkInd
2443 print*,'Ndg: checkInd: tt: ', tt
2444 print*,'Ndg: checkInd: oo: ', oo
2445 print*,'Ndg: checkInd: obsCellInd: ', obsInd
2447 print*,'Ndg: checkInd: discharge: ', discharge(obsInd)
2448 print*,'Ndg: checkInd: theInnov: ',theInnov
2450 print*,'Ndg: checkInd: obsDischarge: ', obsTimeStr(tt)%obsStr(oo)%obsDischarge
2452 print*,'Ndg: checkInd: obsQC: ', obsTimeStr(tt)%obsStr(oo)%obsQC
2453 print*,'Ndg: checkInd: usgsId: ', obsTimeStr(tt)%obsStr(oo)%usgsId
2455 print*,'Ndg: checkInd: time_wt: ', time_wt(hydroTime,obsStaticStr_SoA%tau(staticInd), &
2456 obsTimeStr(tt)%obsStr(oo)%obsTime )
2457 print*,'Ndg: checkInd: hydroTime : ', hydroTime
2458 print*,'Ndg: checkInd: obsTimeStr(tt)%time: ', obsTimeStr(tt)%time
2459 print*,'Ndg: checkInd: obsTime: ', obsTimeStr(tt)%obsStr(oo)%obsTime
2460 print*,'Ndg: checkInd: ws: ', obsStaticStr(staticInd)%ws(ll)
2461 print*,'Ndg: checkInd: nudge_term_link: ', nudge_term_link
2462 print*,'Ndg: checkInd: weighting: ', weighting
2463 if(flushAll) flush(flushUnit)
2464 endif
2465 #endif /* HYDRO_D */
2467 enddo ! oo
2469 enddo ! tt
2471 !!---------------------------------------------
2472 !! Was there a nudge in +-tau or do we apply persistence?
2473 if(abs(weighting) .gt. smallestWeight) then !! 1.e-10
2475 !!$#ifdef HYDRO_D
2476 !!$ if(linkInd .eq. 4) then
2477 !!$ print*,'Ndg: nudge_term_link: ', nudge_term_link
2478 !!$ print*,'Ndg: weighting: ', weighting
2479 !!$ print*,'Ndg: nudge_term_link/weighting: ',nudge_term_link/weighting
2480 !!$ end if
2481 !!$#endif
2482 ! normalize the nudge
2483 nudge_term_link = nudge_term_link/weighting
2485 else if (temporalPersistence) then
2487 !! PERSISTENCE
2489 !! -1) Make a local copy of the global/namelist setting persistBias so that
2490 !! if there are not enough observations available to persistBias, the
2491 !! standard observation persistence is applied.
2492 persistBias_local = persistBias
2494 !! if no nudge was applied then fall back on to persistence of last observation
2495 !! 0) assume we dont find a solution so we can freely bail out.
2496 nudge_term_link = zero
2498 !! 1) what is the current link's gageId?
2499 if(.not. any(obsStaticStr_SoA%obsCellInd(:) .eq. linkInd)) return
2501 allocate(theMask(size(obsStaticStr_SoA%obsCellInd(:))))
2503 theMask = obsStaticStr_SoA%obsCellInd(:) .eq. linkInd
2504 whGageId = whUniLoop(theMask)
2506 deallocate(theMask)
2507 prstGageId = obsStaticStr_SoA%usgsId(whGageId)
2508 obsInd = obsStaticStr_SoA%obsCellInd(whGageId)
2510 !! 2) minimum number of observations
2511 if(persistBias_local) then
2512 !! This section creates a mask of usable values to use in the bias calculation
2513 !! subject to two parameters of the method.
2515 !!-----------------------------------
2516 !! maxAgePairsBiasPersist
2517 !! if greater than 0: apply an age-based filter
2518 !! if zero : apply no additional use all available obs.
2519 !! if less than zero: apply an count-based filter
2522 if(maxAgePairsBiasPersist .gt. 0) then
2523 !! assume the obs are all too old, find the not so old ones
2524 whObsTooOld = .true.
2525 !! assume all weights are bad
2526 if(persistBias_local .and. invDistTimeWeightBias) biasWeights = -9999
2527 do tt=nTimesLastObs,1,-1
2528 !! missing obs are not mixed with good obs. The first missing obs means the rest are missing.
2529 if(obsStaticStr_SoA%lastObsTime(tt,whGageId) .eq. missingLastObsTime) exit
2530 if(biasWindowBeforeT0) then
2531 !! use the init time to define the end of the bias window
2532 call geth_idts(initTime, obsStaticStr_SoA%lastObsTime(tt,whGageId), obsDtInt)
2533 else
2534 !! use the hydro model time to define the end of the bias window
2535 call geth_idts(hydroTime, obsStaticStr_SoA%lastObsTime(tt,whGageId), obsDtInt)
2536 endif
2537 if(persistBias_local .and. invDistTimeWeightBias) &
2538 !JLM OLD
2539 biasWeights(tt) = real(obsDtInt + nlst(did)%dtrt_ch)
2540 !JLM NEW
2541 !biasWeights(tt) = real(obsDtInt + obsResolutionInt)
2543 if(obsDtInt .gt. (maxAgePairsBiasPersist*60*60)) exit
2544 whObsTooOld(tt) = .false.
2545 end do
2546 biasMask= obsStaticStr_SoA%lastObsQuality(:,whGageId) .eq. one .and. &
2547 obsStaticStr_SoA%lastObsTime(:,whGageId) .ne. missingLastObsTime .and. &
2548 obsStaticStr_SoA%lastObsModelDischarge(:,whGageId) .gt. zero .and. &
2549 (.not. whObsTooOld)
2550 end if
2552 if(maxAgePairsBiasPersist .eq. 0) then
2553 biasMask= obsStaticStr_SoA%lastObsQuality(:,whGageId) .eq. one .and. &
2554 obsStaticStr_SoA%lastObsTime(:,whGageId) .ne. missingLastObsTime .and. &
2555 obsStaticStr_SoA%lastObsModelDischarge(:,whGageId) .gt. zero
2556 endif
2558 if(maxAgePairsBiasPersist .lt. 0) then
2559 biasMask= obsStaticStr_SoA%lastObsQuality(:,whGageId) .eq. one .and. &
2560 obsStaticStr_SoA%lastObsTime(:,whGageId) .ne. missingLastObsTime .and. &
2561 obsStaticStr_SoA%lastObsModelDischarge(:,whGageId) .gt. zero
2563 !if((-1*maxAgePairs) .ge. nTimesLastObs) then there's nothing to do.
2564 if((-1*maxAgePairsBiasPersist) .lt. nTimesLastObs) then
2565 pairCount=0
2566 do tt=nTimesLastObs,1,-1
2567 if(biasMask(tt)) pairCount = pairCount + 1
2568 if(pairCount .gt. (-1*maxAgePairsBiasPersist)) biasMask(tt) = .false.
2569 end do
2570 end if
2571 end if
2574 !!-----------------------------------
2575 !! minNumPairsBiasPersist
2576 !! If the pair count starts to dwindle during the forecast period, how does the transition
2577 !! happen from bias correction to open loop?
2578 !! Looks like it would get switched instantaneously. Avoid this problem with the
2579 !! negative maxAgePairsBiasPersist which always use the last maxAgePairsBiasPersist count
2580 !! of obs available. If there enough obs for bias correction at a site, this will ensure
2581 !! that there's no transition away from bias correction in the forecast. If there are not
2582 !! enough obs for bias correction, then observation persistence+decay will be applied as
2583 !! in v1.0.
2584 if(count(biasMask) .lt. minNumPairsBiasPersist) then
2585 persistBias_local = .false.
2586 !! JLM move this to hydro_d
2587 !print*,'JLM: Ndg: not enough observations for bias persistence, using obs persistence', &
2588 ! prstGageId
2589 end if
2590 end if ! if(persistBias_local)
2592 !! Must check this for both cases obs and persistance nudging
2593 !! JLM: currently not also checking the quality. May be ok, but not ideal.
2594 if(obsStaticStr_SoA%lastObsTime(nTimesLastObs,whGageId) .eq. missingLastObsTime) return
2596 !! 3) are there parameters for this link? no -> return
2597 if(.not. any(nudgingParamsStr%usgsId .eq. prstGageId)) return
2599 #ifdef HYDRO_D
2600 print*,'Ndg: Persistence nudging execution'
2601 if(flushAll) flush(flushUnit)
2602 #endif
2604 !! 4) what is the index of this link/gage in nudgingParamsStr?
2605 allocate(theMask(size(nudgingParamsStr%usgsId)))
2606 theMask = nudgingParamsStr%usgsId .eq. prstGageId
2607 whPrstParams = whUniLoop(theMask)
2608 deallocate(theMask)
2610 !! --- Deal with the parameters ---
2611 !! 5) what month is the flow in [1,12]
2612 read(hydroTime(6:7),*) flowMonth
2614 !! 6) what thereshold is the flow?
2615 whThresh=1 !! init value
2616 do tt=1,size(nudgingParamsStr%qThresh,1)
2617 if(discharge(obsInd) .gt. nudgingParamsStr%qThresh(tt, flowMonth, whPrstParams)) &
2618 whThresh=tt+1
2619 end do
2621 !!7) get the parameter of the exponential decay
2622 prstB=nudgingParamsStr%expCoeff(whThresh, flowMonth, whPrstParams)
2624 !! --- Deal with the last ob ---
2625 !! 9) what is the time difference with the last ob/obs?
2626 prstDtErrInt = -9999*60 ! JLM: magic value? set as a parameter above
2627 do tt=nTimesLastObs,1,-1
2628 ! Conditons for rejecting obs
2629 if(persistBias_local) then
2630 if(.not. biasMask(tt)) cycle
2631 else
2632 if(obsStaticStr_SoA%lastObsQuality(tt,whGageId) .eq. 0) cycle
2633 end if
2634 call geth_idts(hydroTime, obsStaticStr_SoA%lastObsTime(tt,whGageId), prstDtErrInt)
2635 if(prstDtErrInt .lt. zeroInt) then
2636 print*,'initTime:', hydroTime
2637 print*,'obsStaticStr(whGageId)%lastObsTime(tt):', obsStaticStr_SoA%lastObsTime(tt,whGageId)
2638 print*,'prstDtErrInt:',prstDtErrInt
2639 print*,'Ndg: WARNING: lastObs times are in the future of the model time.'
2640 if(futureLastObsFatal) then
2641 call hydro_stop('nudge_term_link: negative prstDtErrInt: last obs in future')
2642 else
2643 prstDtErrInt=abs(prstDtErrInt)
2644 end if
2645 end if
2646 exit
2647 end do
2648 ! By design, should not get here
2649 if(prstDtErrInt .eq. -9999*60) then
2650 persistBias_local=.false.
2651 call hydro_stop('nudge_term_link: negative prstDtErrInt: still default value=-9999*60')
2652 end if
2653 prstDtErr = real(prstDtErrInt)/real(sixty) !! second -> minutes
2655 !! 10) what is the error or bias ?
2657 !! Always calculate the observation nudge error.
2658 prstErr = obsStaticStr_SoA%lastObsDischarge(nTimesLastObs,whGageId) - discharge(obsInd)
2660 !! Bias term.
2661 if(persistBias_local) then
2663 if(invDistTimeWeightBias) then
2665 if(count(biasMask .and. (biasWeights.gt.zero)) .lt. minNumPairsBiasPersist) then
2666 persistBias_local=.false.
2667 else
2669 !! normalize to the latest (least) time difference and invert.
2670 biasWeights=1./(biasWeights/minval(biasWeights, &
2671 mask=biasMask .and. (biasWeights.gt.zero) ))
2672 biasWeights=biasWeights**invDistTimeWeightExp !! = 5.000, hard coded above
2674 !! The less than equal to one prevents division by huge/infinite
2675 !! numbers which may result from zero time differences.
2677 ! JLM OLD
2678 biasWeightSum=sum(biasWeights, mask=biasMask .and. (biasWeights.gt.zero) )
2679 ! JLM NEW
2680 !biasWeightSum=sum(biasWeights, mask=biasMask .and. &
2681 ! (biasWeights.gt.zero) .and. &
2682 ! (biasWeights.le.one) )
2684 theBias = sum( biasWeights * &
2685 (obsStaticStr_SoA%lastObsDischarge(:,whGageId) - &
2686 obsStaticStr_SoA%lastObsModelDischarge(:,whGageId)), &
2687 mask=biasMask .and. (biasWeights.gt.zero)) / &
2688 biasWeightSum
2690 endif
2692 else
2694 !! Simple, straight average / unweighted, way of calculating bias
2695 theBias = sum( obsStaticStr_SoA%lastObsDischarge(:,whGageId) - &
2696 obsStaticStr_SoA%lastObsModelDischarge(:,whGageId), mask=biasMask ) / &
2697 max( 1, count(biasMask) )
2699 end if
2701 if(noConstInterfBias) then
2702 !! Game to reduce bias issues.
2703 !! Maybe make this an option later if it works.
2704 !! do I need to stash the t0 flow? It is not the same as the last ob or the lastObsModelDischarge
2705 deltaT0Discharge = discharge(obsInd) - t0Discharge(linkIndIn)
2707 if(theBias .lt. zero .and. &
2708 deltaT0Discharge .lt. zero ) then
2710 if(t0Discharge(linkIndIn) .gt. .01) then
2711 !! Tweaks
2712 !1!theBias=theBias - noConstInterfCoeff*(deltaT0Discharge) ! negative - negative
2713 !2!theBias=theBias*(1-(noConstInterfCoeff*abs(deltaT0Discharge)/t0Discharge(linkIndIn)))
2714 !3
2715 theBias = theBias*&
2716 (1.0-(noConstInterfCoeff* &
2717 ( (abs(deltaT0Discharge)/t0Discharge(linkIndIn)+.25)**2.5 - .25**2.5) ) )
2718 else
2719 !! if the toDischarge is really small, zero out the bias.
2720 theBias = zero
2721 end if
2722 theBias = min( theBias, zero ) ! do not let the bias go positive
2723 theBias = max( theBias, -1*discharge(obsInd)/2. ) ! do not let the bias to remove be more than half the total flow.
2725 end if
2728 if(theBias .gt. zero .and. &
2729 deltaT0Discharge .gt. zero ) then
2731 if(t0Discharge(linkIndIn) .gt. .01) then
2732 !! Tweaks
2733 !1 !theBias = theBias - noConstInterfCoeff*(deltaT0Discharge) ! positive - positive
2734 !2!theBias=theBias*(1-(noConstInterfCoeff*abs(deltaT0Discharge)/t0Discharge(linkIndIn)))
2735 !3
2736 theBias = theBias*&
2737 (1.0-(noConstInterfCoeff* &
2738 ( (abs(deltaT0Discharge)/t0Discharge(linkIndIn)+.25)**2.5 - .25**2.5) ) )
2739 else
2740 !! if the toDischarge is really small, use a naminal t0Discharge of .01
2741 theBias = theBias*&
2742 (1.0-(noConstInterfCoeff* &
2743 ( (abs(deltaT0Discharge)/.01+.25)**2.5 - .25**2.5) ) )
2744 end if
2745 theBias = max( theBias, zero ) ! do not let the bias go negative
2747 end if
2748 end if
2750 endif ! if (persistBias_local)
2753 !! 11) solution
2754 if(.not. persistBias_local) then
2755 nudge_term_link = &
2756 obsStaticStr_SoA%lastObsQuality(nTimesLastObs,whGageId) * &
2757 prstErr * exp(real(-1)*prstDtErr/prstB)
2758 else
2759 nudge_term_link = &
2760 obsStaticStr_SoA%lastObsQuality(nTimesLastObs,whGageId) * &
2761 prstErr * exp(real(-1)*prstDtErr/prstB) + &
2762 theBias * (1.0 - exp(real(-1)*prstDtErr/prstB) )
2764 endif !! if(.not. persistBias_local) then; else;
2766 !! Currently cant get here...
2767 if((prstDtErr .lt. zero) .and. futureLastObsZero) nudge_term_link = zero
2769 !! Ensure the nudge does not make the streamflow less than zero.
2770 !! This is particularly important if persistBias
2771 !! Is this OK as ZERO?? Musk-cunge allows it (in theory).
2772 if (discharge(obsInd)+nudge_term_link .lt. zero) then
2773 nudge_term_link = zero - discharge(obsInd) + smallestNudge
2774 !print*,'JLM setting nudge_term_link + observed discharge to zero (or very small)', prstGageId
2775 endif
2778 #ifdef HYDRO_D
2779 if(linkInd .eq. checkInd) then
2780 print*,'Ndg: checkInd pst: -------------------'
2781 print*,'Ndg: checkInd pst: usgsId: ', prstGageId
2782 print*,'Ndg: checkInd pst: -------------------'
2783 print*,'Ndg: checkInd pst: prstGageId: ', prstGageId
2784 print*,'Ndg: checkInd pst: obsInd: ', obsInd
2785 print*,'Ndg: checkInd pst: initTime:', initTime
2786 print*,'Ndg: checkInd pst: hydroTime: ', hydroTime
2787 print*,'Ndg: checkInd pst: flowmonth: ', flowMonth
2788 print*,'Ndg: checkInd pst: lastObsTime: ', obsStaticStr_SoA%lastObsTime(nTimesLastObs,whGageId)
2789 print*,'Ndg: checkInd pst: time difference mins: ', prstDtErr
2790 print*,'Ndg: checkInd pst: obsStaticStr(whGageId)%lastObsQuality: ', &
2791 obsStaticStr_SoA%lastObsQuality(nTimesLastObs,whGageId)
2792 print*,'Ndg: checkInd pst: obsStaticStr(whGageId)%lastObsDischarge: ', &
2793 obsStaticStr_SoA%lastObsDischarge(nTimesLastObs,whGageId)
2794 print*,'Ndg: checkInd pst: modeled discharge(obsInd): ', discharge(obsInd)
2795 print*,'Ndg: checkInd pst: qThresh: ', nudgingParamsStr%qThresh(:, flowMonth, whPrstParams)
2796 print*,'Ndg: checkInd pst: whThresh: ', whThresh
2797 print*,'Ndg: checkInd pst: discharge difference:', prstErr
2798 print*,'Ndg: checkInd pst: prstB:', prstB
2799 print*,'Ndg: checkInd pst: time difference mins prstDtErr: ', prstDtErr
2800 print*,'Ndg: checkInd pst: exponent term: ', exp(real(-1)*prstDtErr/prstB)
2802 if(persistBias_local) then
2803 print*,'Ndg: checkInd pst: BIAS PERSISTENCE INFO'
2804 do tt=1,nTimesLastObs
2805 if(biasMask(tt)) then
2806 print*,'checkInd: ---- tt:', tt
2807 print*,'checkInd: Diff: ', obsStaticStr_SoA%lastObsDischarge(tt,whGageId) - &
2808 obsStaticStr_SoA%lastObsModelDischarge(tt,whGageId)
2809 print*,'checkInd: o - : ', obsStaticStr_SoA%lastObsDischarge(tt,whGageId)
2810 print*,'checkInd: m: ', obsStaticStr_Soa%lastObsModelDischarge(tt,whGageId)
2811 print*,'checkInd: biasWeights: ', biasWeights(tt)
2812 end if
2813 end do
2814 print*,'Ndg: checkInd pst: max( 1, count(biasMask) ): ',max( 1, count(biasMask) )
2815 print*,'Ndg: checkInd pst: biasWeightSum: ', biasWeightSum
2816 print*,'Ndg: checkInd pst: corr: ', corr
2817 print*,'Ndg: checkInd pst: nCorr: ', nCorr
2818 print*,'Ndg: checkInd pst: theBias:', theBias
2819 print*,'Ndg: checkInd pst: theBias/corr:', theBias/corr
2820 print*,'Ndg: checkInd pst: obsStaticStr(whGageId)%lastObsQuality(nTimesLastObs) *'
2821 print*,'Ndg: checkInd pst: prstErr * exp(real(-1)*prstDtErr/prstB)', &
2822 obsStaticStr_SoA%lastObsQuality(nTimesLastObs,whGageId) * &
2823 prstErr * exp(real(-1)*prstDtErr/prstB)
2824 print*,'Ndg: checkInd pst: theBias * (1-exp(real(-1)*prstDtErr/prstB)):', &
2825 theBias * (1-exp(real(-1)*prstDtErr/prstB))
2826 endif
2828 print*,'Ndg: checkInd pst: nudge_term_link: ', nudge_term_link
2830 if(flushAll) flush(flushUnit)
2831 end if
2832 #endif /* HYDRO_D */
2833 endif !! if(abs(weighting) .gt. smallestWeight) then !! 1.e-10
2835 #ifdef HYDRO_D
2836 #ifdef MPP_LAND
2837 if(my_id .eq. io_id) &
2838 #endif
2839 print*,'Ndg: finish nudge_term_link'
2840 if(flushAll) flush(flushUnit)
2841 #endif
2844 end function nudge_term_link
2847 !===================================================================================================
2848 ! Program Name:
2849 ! nudge_apply_upstream_muskingumCunge
2850 ! Author(s)/Contact(s):
2851 ! James L McCreight, <jamesmcc><ucar><edu>
2852 ! Abstract:
2853 ! Gets the previous nudge into the collected upstream current and previous fluxes.
2854 ! History Log:
2855 ! 4/5/17 - Created, JLM
2856 ! Usage:
2857 ! Parameters:
2858 ! Input Files:
2859 ! Output Files:
2860 ! Condition codes:
2861 ! User controllable options:
2862 ! Notes:
2864 subroutine nudge_apply_upstream_muskingumCunge( qup, quc, np, k )
2865 use module_RT_data, only: rt_domain
2866 implicit none
2867 real, intent(inout) :: qup ! previous upstream inflows
2868 real, intent(inout) :: quc ! current upstream inflows
2869 real, intent(in) :: np ! previous nudge
2870 integer, intent(in) :: k ! index of flow/routlink on local image
2872 !! If not on a gage... get out.
2873 if(rt_domain(1)%gages(k) .eq. rt_domain(1)%gageMiss) return
2875 qup = qup + np
2876 quc = quc + np
2878 end subroutine nudge_apply_upstream_muskingumCunge
2880 !===================================================================================================
2881 ! Program Name:
2882 ! get_netwk_reexpression
2883 ! Author(s)/Contact(s):
2884 ! James L McCreight, <jamesmcc><ucar><edu>
2885 ! Abstract:
2886 ! Bring in the network re-expression for indexed traversal
2887 ! of the stream network.
2888 ! History Log:
2889 ! 8/20/15 - Created, JLM
2890 ! Usage:
2891 ! Parameters:
2892 ! Input Files:
2893 ! Output Files:
2894 ! Condition codes:
2895 ! User controllable options:
2896 ! Notes: Sets module derived type strs: upNetwkStr and downNetwkStr
2897 subroutine get_netwk_reexpression
2898 use module_RT_data, only: rt_domain
2899 use module_nudging_io,only: read_network_reexpression, get_netcdf_dim
2901 implicit none
2902 integer :: downSize, upSize, baseSize, nLinksL
2904 #ifdef MPP_LAND
2905 nLinksL = RT_DOMAIN(did)%gNLINKSL ! For reach-based routing in parallel, no decomp for nudging
2906 #else
2907 nLinksL = RT_DOMAIN(did)%NLINKSL ! For reach-based routing
2908 #endif
2910 #ifdef MPP_LAND
2911 if(my_id .eq. IO_id) then
2912 #endif
2913 downSize = get_netcdf_dim(netwkReExFile, 'downDim', 'init_stream_nudging')
2914 upSize = get_netcdf_dim(netwkReExFile, 'upDim', 'init_stream_nudging')
2915 baseSize = get_netcdf_dim(netwkReExFile, 'baseDim', 'init_stream_nudging')
2916 !if(baseSize .ne. nLinksL) call hydro_stop('init_stream_nudging: baseSize .ne. nLinksL')
2918 allocate(upNetwkStr%go(upSize))
2919 allocate(upNetwkStr%start(nLinksL))
2920 allocate(upNetwkStr%end(nLinksL))
2921 allocate(downNetwkStr%go(downSize))
2922 allocate(downNetwkStr%start(nLinksL))
2923 allocate(downNetwkStr%end(nLinksL))
2925 call read_network_reexpression( &
2926 netwkReExFile, & ! file with dims of the stream netwk
2927 upNetwkStr%go, & ! where each ind came from, upstream
2928 upNetwkStr%start, & ! where each ind's upstream links start in upGo
2929 upNetwkStr%end, & ! where each ind's upstream links end in upGo
2930 downNetwkStr%go, & ! where each ind goes, downstream
2931 downNetwkStr%start, & ! where each ind's downstream links start in downGo
2932 downNetwkStr%end ) ! where each ind's downstream links end in downGo
2934 #ifdef MPP_LAND
2935 endif ! my_id .eq. io_id
2937 ! Broadcast
2938 call mpp_land_bcast_int1(upSize)
2939 call mpp_land_bcast_int1(downSize)
2940 if(my_id .ne. io_id) then
2941 allocate(upNetwkStr%go(upSize))
2942 allocate(upNetwkStr%start(nLinksL))
2943 allocate(upNetwkStr%end(nLinksL))
2944 allocate(downNetwkStr%go(downSize))
2945 allocate(downNetwkStr%start(nLinksL))
2946 allocate(downNetwkStr%end(nLinksL))
2947 endif
2948 call mpp_land_bcast_int1d(upNetwkStr%go)
2949 call mpp_land_bcast_int1d(upNetwkStr%start)
2950 call mpp_land_bcast_int1d(upNetwkStr%end)
2951 call mpp_land_bcast_int1d(downNetwkStr%go)
2952 call mpp_land_bcast_int1d(downNetwkStr%start)
2953 call mpp_land_bcast_int1d(downNetwkStr%end)
2954 #endif
2955 end subroutine get_netwk_reexpression
2958 !===================================================================================================
2959 ! Program Name:
2960 ! finish_stream_nudging
2961 ! Author(s)/Contact(s):
2962 ! James L McCreight, <jamesmcc><ucar><edu>
2963 ! Abstract:
2964 ! Finish off the nudging routines, memory and timing.
2965 ! History Log:
2966 ! 8/20/15 - Created, JLM
2967 ! Usage:
2968 ! Parameters:
2969 ! Input Files:
2970 ! Output Files:
2971 ! Condition codes:
2972 ! User controllable options:
2973 ! Notes:
2975 subroutine finish_stream_nudging
2976 use module_nudging_utils, only: totalNudgeTime, &
2977 nudging_timer, &
2978 accum_nudging_time, &
2979 whUniLoop
2980 use module_nudging_io, only: write_nwis_not_in_RLAndParams
2981 implicit none
2983 character(len=15), allocatable, dimension(:) :: nwisNotRLAndParamsGathered
2984 character(len=15), allocatable, dimension(:) :: nwisNotRLAndParamsConsolidated
2985 integer, dimension(numprocs) :: nwisNotRLAndParamsCountVector
2986 integer :: nwisNotRLAndParamsCountConsolidated, gg, ii
2988 #ifdef HYDRO_D /* un-ifdef to get timing results */
2989 !!$real :: startCodeTimeAcc, endCodeTimeAcc
2990 #ifdef MPP_LAND
2991 if(my_id .eq. io_id) &
2992 #endif
2993 print*,'Ndg: start finish_stream_nudging'
2994 if(flushAll) flush(flushUnit)
2995 !!$#ifdef MPP_LAND
2996 !!$if(my_id .eq. io_id) &
2997 !!$#endif
2998 !!$ call nudging_timer(startCodeTimeAcc)
2999 #endif /* HYDRO_D un-ifdef to get timing results */
3001 !! accumulate_nwis_not_in_RLAndParams
3002 !! because of parallel IO, accumulation is happening on
3003 !! multiple images and these need to be consolidated before outputting.
3004 ! get the total of nwisNotRLAndParam gages on all images
3005 nwisNotRLAndParamsCountVector(my_id+1) = nwisNotRLAndParamsCount
3006 do ii=0,numprocs-1
3007 call mpp_land_bcast_int1_root(nwisNotRLAndParamsCountVector(ii+1), ii)
3008 end do
3010 allocate(nwisNotRLAndParamsGathered(sum(nwisNotRLAndParamsCountVector)))
3012 #ifdef MPP_LAND
3013 if(my_id .eq. io_id) &
3014 #endif
3015 allocate(nwisNotRLAndParamsConsolidated(sum(nwisNotRLAndParamsCountVector)))
3017 call write_IO_char_head(nwisNotRLAndParams, nwisNotRLAndParamsGathered, nwisNotRLAndParamsCountVector)
3019 if(my_id .eq. io_id) then
3020 ! now consolidate these again...
3021 nwisNotRLAndParamsCountConsolidated=0
3022 do gg=1,size(nwisNotRLAndParamsGathered)
3023 call accumulate_nwis_not_in_RLAndParams(nwisNotRLAndParamsConsolidated, &
3024 nwisNotRLAndParamsCountConsolidated, &
3025 nwisNotRLAndParamsGathered(gg) )
3026 end do
3027 if(nwisNotRLAndParamsCountConsolidated .gt. 0) &
3028 call write_nwis_not_in_RLAndParams(nwisNotRLAndParamsConsolidated, &
3029 nwisNotRLAndParamsCountConsolidated )
3030 endif
3032 !! deallocate local variables
3033 if(allocated(nwisNotRLAndParamsGathered)) deallocate(nwisNotRLAndParamsGathered)
3034 if(allocated(nwisNotRLAndParamsConsolidated)) deallocate(nwisNotRLAndParamsConsolidated)
3036 !! deallocate module variables here and time this
3037 if(allocated(obsTimeStr)) deallocate(obsTimeStr)
3038 if(allocated(obsStaticStr)) deallocate(obsStaticStr)
3039 !if(allocated(nodeGageTmp%nodeId)) deallocate(nodeGageTmp%nodeId)
3040 if(allocated(nodeGageTmp%usgsId)) deallocate(nodeGageTmp%usgsId) !redundant
3041 if(allocated(nodeGageStr%nodeId)) deallocate(nodeGageStr%nodeId)
3042 if(allocated(nodeGageStr%usgsId)) deallocate(nodeGageStr%usgsId)
3043 if(allocated(nudgingParamsTmp%usgsId)) deallocate(nudgingParamsTmp%usgsId) !redundant
3044 if(allocated(nudgingParamsTmp%R)) deallocate(nudgingParamsTmp%R) !redundant
3045 if(allocated(nudgingParamsTmp%G)) deallocate(nudgingParamsTmp%G) !redundant
3046 if(allocated(nudgingParamsTmp%tau)) deallocate(nudgingParamsTmp%tau) !redundant
3047 if(allocated(nudgingParamsStr%usgsId)) deallocate(nudgingParamsStr%usgsId)
3048 if(allocated(nudgingParamsStr%R)) deallocate(nudgingParamsStr%R)
3049 if(allocated(nudgingParamsStr%G)) deallocate(nudgingParamsStr%G)
3050 if(allocated(nudgingParamsStr%tau)) deallocate(nudgingParamsStr%tau)
3051 if(allocated(nudgingParamsStr%qThresh)) deallocate(nudgingParamsStr%qThresh)
3052 if(allocated(nudgingParamsStr%expCoeff)) deallocate(nudgingParamsStr%expCoeff)
3053 if(allocated(upNetwkStr%go)) deallocate(upNetwkStr%go)
3054 if(allocated(upNetwkStr%start)) deallocate(upNetwkStr%start)
3055 if(allocated(upNetwkStr%end)) deallocate(upNetwkStr%end)
3056 if(allocated(downNetwkStr%go)) deallocate(downNetwkStr%go)
3057 if(allocated(downNetwkStr%start)) deallocate(downNetwkStr%start)
3058 if(allocated(downNetwkStr%end)) deallocate(downNetwkStr%end)
3059 if(allocated(chanlen_image0)) deallocate(chanlen_image0)
3061 !! print ouf the full timing results here.
3062 !! dont i need to make another accumultion call here?
3063 #ifdef HYDRO_D
3064 #ifdef MPP_LAND
3065 if(my_id .eq. io_id) then
3066 #endif
3067 !!$ call nudging_timer(endCodeTimeAcc)
3068 !!$ call accum_nudging_time(startCodeTimeAcc, endCodeTimeAcc, 'very end of nudging', .true.)
3069 !!$ print*,'Ndg: *nudge timing* TOTAL NUDGING TIME (seconds): ', totalNudgeTime
3070 print*,'Ndg: end of finish_stream_nudging'
3071 if(flushAll) flush(flushUnit)
3072 #ifdef MPP_LAND
3073 end if
3074 #endif
3075 #endif /* HYDRO_D */
3077 end subroutine finish_stream_nudging
3081 end module module_stream_nudging