1module sstBias_mod
2 ! MODULE sstBias_mod (prefix='sstb' category='1. High-level functionality')
3 !
4 !:Purpose: Compute SST satellite data bias estimation and correction.
5 !
6 use obsSpaceData_mod
7 use horizontalCoord_mod
8 use verticalCoord_mod
9 use timeCoord_mod
10 use kdTree2_mod
11 use codePrecision_mod
12 use mathPhysConstants_mod
13 use utilities_mod
14 use midasMpi_mod
15 use codtyp_mod
16 use gridStateVector_mod
17 use gridStateVectorFileIO_mod
18 use oceanMask_mod
19 use localizationFunction_mod
20 use columnData_mod
21 use statetocolumn_mod
22
23 implicit none
24 save
25 private
26
27 ! public subroutines
28 public :: sstb_computeBias, sstb_applySatelliteSSTBiasCorrection
29
30 ! external
31 integer, external :: fnom, fclos
32
33 ! mpi topology
34 integer :: myLatBeg, myLatEnd
35 integer :: myLonBeg, myLonEnd
36 integer :: latPerPE, latPerPEmax, lonPerPE, lonPerPEmax
37
38 integer, parameter :: maxNumberSensors = 10
39
40 ! namelist variables
41 real(8) :: searchRadius ! horizontal search radius, in km, for obs gridding
42 real(4) :: maxBias ! max acceptable difference (insitu - satellite)
43 real(4) :: iceFractionThreshold ! consider no ice condition below this threshold
44 integer :: numberPointsBG ! parameter, number of matchups of the background bias estimation
45 character(len=20) :: timeInterpType_nl ! 'NEAREST' or 'LINEAR'
46 integer :: numObsBatches ! number of batches for calling interp setup
47 integer :: numberSensors ! MUST NOT BE INCLUDED IN NAMELIST!
48 character(len=10) :: sensorList(maxNumberSensors) ! list of sensors
49 logical :: saveAuxFields ! to store or not auxiliary fields: nobs and weight
50 real(4) :: weightMin ! minimum value of weight for the current day bias
51 real(4) :: weightMax ! maximum value of weight for the current day bias
52 real(4) :: bgTermZeroBias ! background term of zero bias estimate
53
54 contains
55
56 !--------------------------------------------------------------------------
57 ! sstb_computeBias
58 !--------------------------------------------------------------------------
59 subroutine sstb_computeBias(obsData, hco, vco)
60 !
61 !:Purpose: compute bias for SST satellite data with respect to insitu data
62 !
63 implicit none
64
65 ! Arguments:
66 type(struct_obs), intent(inout) :: obsData ! obsSpaceData
67 type(struct_hco), intent(inout), pointer :: hco ! horizontal grid structure
68 type(struct_vco), intent(in) , pointer :: vco ! vertical grid structure
69
70 ! Locals:
71 integer :: sensorIndex, productIndex
72 real(8) :: insituGrid (hco%ni, hco%nj)
73 real(8) :: satelliteGrid(hco%ni, hco%nj)
74 logical :: mask(hco%ni, hco%nj), openWater(hco%ni, hco%nj)
75 type(struct_ocm) :: oceanMask
76 integer :: numberOpenWaterPoints, lonIndex, latIndex
77 integer :: nobsFoundInsituGlob, nobsFoundInsituLoc
78 integer :: nobsFoundSatGlob, nobsFoundSatLoc
79 type(struct_gsv) :: stateVector_ice
80 real(4), pointer :: seaice_ptr(:,:,:)
81 integer , parameter :: numberProducts = 2 ! day and night
82 character(len=*), parameter :: listProducts(numberProducts)= (/'day', 'night'/)
83
84 write(*,*) 'sstb_computeBias: Starting...'
85 write(*,*) 'sstb_computeBias: Sea-ice Fraction threshold: ', iceFractionThreshold
86
87 ! read the namelist
88 call readNml()
89
90 ! get mpi topology
91 call mmpi_setup_lonbands(hco%ni, lonPerPE, lonPerPEmax, myLonBeg, myLonEnd)
92 call mmpi_setup_latbands(hco%nj, latPerPE, latPerPEmax, myLatBeg, myLatEnd)
93
94 ! get latest sea-ice analysis
95 call gsv_allocate(stateVector_ice, 1, hco, vco, dataKind_opt = 4, &
96 datestamp_opt = -1, mpi_local_opt = .false., &
97 varNames_opt = (/'LG'/), hInterpolateDegree_opt ='LINEAR')
98 call gio_readFromFile(stateVector_ice, './seaice_analysis', ' ','A', &
99 unitConversion_opt=.false., containsFullField_opt=.true.)
100 call gsv_getField(stateVector_ice, seaice_ptr)
101
102 ! Get land mask from analysisgrid file (1=water, 0=land)
103 ! and the number of open water points
104 call ocm_readMaskFromFile(oceanMask, hco, vco, './analysisgrid')
105
106 numberOpenWaterPoints = 0
107 mask(:, :) = .false.
108 openWater(:, :) = .false.
109 do latIndex = 1, hco%nj
110 do lonIndex = 1, hco%ni
111 if (oceanMask%mask(lonIndex, latIndex, 1)) then
112 mask(lonIndex, latIndex) = .true.
113 if (seaice_ptr(lonIndex, latIndex, 1) <= iceFractionThreshold) then
114 openWater(lonIndex, latIndex) = .true.
115 numberOpenWaterPoints = numberOpenWaterPoints + 1
116 end if
117 end if
118 end do
119 end do
120 call ocm_deallocate(oceanMask)
121 call gsv_deallocate(stateVector_ice)
122 write(*,*) 'sstb_computeBias: computing bias for ', numberOpenWaterPoints, ' open water points'
123
124 insituGrid(:, :) = MPC_missingValue_R8
125
126 call sstb_getGriddedObs(obsData, insituGrid, nobsFoundInsituGlob, &
127 nobsFoundInsituLoc, hco, openWater, 'insitu')
128
129 do sensorIndex = 1, numberSensors
130 do productIndex = 1, numberProducts
131 if (nobsFoundInsituGlob > 0) then
132 satelliteGrid(:, :) = MPC_missingValue_R8
133 call sstb_getGriddedObs(obsData, satelliteGrid(:, :), nobsFoundSatGlob, &
134 nobsFoundSatLoc, hco, openWater, sensorList(sensorIndex), &
135 dayOrNight_opt = listProducts(productIndex))
136 if (nobsFoundSatGlob > 0) then
137 call sstb_getGriddedBias(satelliteGrid(:, :), insituGrid, nobsFoundSatGlob, &
138 nobsFoundSatLoc, hco, vco, mask, openWater, &
139 sensorList(sensorIndex), numberOpenWaterPoints, &
140 listProducts(productIndex))
141 else
142 write(*,*) 'sstb_computeBias: WARNING: missing ', trim(sensorList(sensorIndex)), ' ', &
143 trim(listProducts(productIndex)),' data.'
144 write(*,*) 'Bias estimate will be read from the previous state...'
145 call sstb_getBiasFromPreviousState(hco, vco, &
146 sensorList(sensorIndex), &
147 listProducts(productIndex))
148 end if
149 else
150 write(*,*) 'sstb_computeBias: WARNING: missing insitu data.'
151 write(*,*) 'sstb_computeBias: bias estimates for all sensors will be read from the previous state...'
152 call sstb_getBiasFromPreviousState(hco, vco, &
153 sensorList(sensorIndex), &
154 listProducts(productIndex))
155 end if
156 end do
157 end do
158
159 write(*,*) 'sstb_computeBias: done.'
160
161 end subroutine sstb_computeBias
162
163 !--------------------------------------------------------------------------
164 ! sstb_getGriddedObs
165 !--------------------------------------------------------------------------
166 subroutine sstb_getGriddedObs(obsData, obsGrid, countObsGlob, countObsLoc, &
167 hco, openWater, instrument, dayOrNight_opt)
168 !
169 !:Purpose: put observations of a given family on the grid
170 !
171 implicit none
172
173 ! Arguments:
174 type(struct_obs), intent(inout) :: obsData ! obsSpaceData
175 real(8) , intent(inout) :: obsGrid(:,:) ! observations on the grid
176 integer , intent(out) :: countObsGlob ! global number of data found (all procs)
177 integer , intent(out) :: countObsLoc ! number of data found (current MPI proc)
178 type(struct_hco), intent(in) , pointer :: hco ! horizontal grid structure
179 logical , intent(in) :: openWater(:,:) ! open water points (.true.)
180 character(len=*), intent(in) :: instrument ! name of instrument
181 character(len=*), intent(in), optional :: dayOrNight_opt ! look for daytime or nighttime obs
182
183 ! Locals:
184 integer, parameter :: maxPointsSearch = 200000
185 real(4), parameter :: solarZenithThreshold = 90.0 ! to distinguish day and night
186 type(kdtree2), pointer :: tree => null()
187 real(kdkind), allocatable :: positionArray(:,:)
188 type(kdtree2_result) :: searchResults(maxPointsSearch)
189 real(kdkind) :: refPosition(3)
190 real(kdkind) :: lon_grd, lat_grd
191 real(pre_obsReal) :: lon_obs, lat_obs
192 integer :: lonIndex, latIndex
193 integer :: bodyIndex, headerIndex, ierr, headerCounter, codtyp
194 integer :: localObsIndex
195 integer :: ndataFoundGridLoc(hco%ni, hco%nj) ! kd-tree output: number of data found within the search radius for every grid point
196 integer :: ndataFoundGridGlob(hco%ni, hco%nj) ! to compute mpi_allreduce of ndataFoundGridLoc
197 real(kdkind) :: searchRadiusSquared
198 integer, allocatable :: headerIndexes(:)
199 real(8) :: currentObs
200 character(len=50) :: instrumentString
201
202 countObsLoc = 0
203
204 ! count local observations for the given instrument
205 if (trim(instrument) == 'insitu') then
206 instrumentString = trim(instrument)
207 do headerIndex = 1, obs_numheader(obsData)
208 codtyp = obs_headElem_i(obsData, obs_ity, headerIndex)
209 if (codtyp == codtyp_get_codtyp('shipnonauto') .or. &
210 codtyp == codtyp_get_codtyp('drifter') .or. &
211 codtyp == codtyp_get_codtyp('ashipauto')) then
212 countObsLoc = countObsLoc + 1
213 end if
214 end do
215 else
216 instrumentString = trim(instrument)//' '//dayOrNight_opt
217 do headerIndex = 1, obs_numheader(obsData)
218 if (obs_elem_c(obsData, 'STID' , headerIndex) == trim(instrument)) then
219 if (trim(dayOrNight_opt) == 'day') then
220 if (obs_headElem_r(obsData, obs_sun, headerIndex) < solarZenithThreshold) &
221 countObsLoc = countObsLoc + 1
222 else if (trim(dayOrNight_opt) == 'night') then
223 if (obs_headElem_r(obsData, obs_sun, headerIndex) >= solarZenithThreshold) &
224 countObsLoc = countObsLoc + 1
225 end if
226 end if
227 end do
228 end if
229 write(*,*) ''
230 write(*,"(a, i10, a)") 'sstb_getGriddedObs: found ', countObsLoc, ' '//trim(instrumentString)//' data'
231
232 call rpn_comm_allreduce(countObsLoc, countObsGlob, 1, "mpi_integer", "mpi_sum", "grid", ierr)
233
234 obsGrid(:, :) = 0.0d0
235 ndataFoundGridLoc(:,:) = 0
236 ndataFoundGridGlob(:,:) = 0
237
238 POSITIVECOUNTOBSLOC: if (countObsLoc > 0) then
239
240 write(*,*)'sstb_getGriddedObs: define kd-tree using data positions...'
241 allocate(positionArray(3, countObsLoc))
242 allocate(headerIndexes(countObsLoc))
243
244 headerCounter = 0
245 do headerIndex = 1, obs_numheader(obsData)
246 if (trim(instrument) == 'insitu') then
247 codtyp = obs_headElem_i(obsData, obs_ity, headerIndex)
248 if (codtyp == codtyp_get_codtyp('shipnonauto') .or. &
249 codtyp == codtyp_get_codtyp('drifter') .or. &
250 codtyp == codtyp_get_codtyp('ashipauto')) then
251 lon_obs = obs_headElem_r(obsData, obs_lon, headerIndex)
252 lat_obs = obs_headElem_r(obsData, obs_lat, headerIndex)
253 headerCounter = headerCounter + 1
254 positionArray(:, headerCounter) = kdtree2_3dPosition(lon_obs, lat_obs)
255 headerIndexes(headerCounter) = headerIndex
256 end if
257 else
258 if (obs_elem_c(obsData, 'STID' , headerIndex) == trim(instrument)) then
259 if (trim(dayOrNight_opt) == 'day') then
260 if (obs_headElem_r(obsData, obs_sun, headerIndex) < solarZenithThreshold) then
261 lon_obs = obs_headElem_r(obsData, obs_lon, headerIndex)
262 lat_obs = obs_headElem_r(obsData, obs_lat, headerIndex)
263 headerCounter = headerCounter + 1
264 positionArray(:, headerCounter) = kdtree2_3dPosition(lon_obs, lat_obs)
265 headerIndexes(headerCounter) = headerIndex
266 end if
267 else if (trim(dayOrNight_opt) == 'night') then
268 if (obs_headElem_r(obsData, obs_sun, headerIndex) >= solarZenithThreshold) then
269 lon_obs = obs_headElem_r(obsData, obs_lon, headerIndex)
270 lat_obs = obs_headElem_r(obsData, obs_lat, headerIndex)
271 headerCounter = headerCounter + 1
272 positionArray(:, headerCounter) = kdtree2_3dPosition(lon_obs, lat_obs)
273 headerIndexes(headerCounter) = headerIndex
274 end if
275 end if
276 end if
277 end if
278 end do
279
280 nullify(tree)
281 tree => kdtree2_create(positionArray, sort=.true., rearrange=.true.)
282
283 ! do the search
284 write(*, '(a,f5.1,a)') 'sstb_getGriddedObs: Collocation radius: ', searchRadius, ' km'
285 searchRadiusSquared = (1.1d0 * searchRadius * 1000.d0)**2 ! convert from km to m2
286
287 write(*,*) 'sstb_getGriddedObs: computing the sum of data values and '//&
288 'their number within the collocation radius for every grid point...'
289 do latIndex = 1, hco%nj
290 do lonIndex = 1, hco%ni
291
292 ! compute gridded obs for every open water point
293 OPENWATERPTS: if (openWater(lonIndex, latIndex)) then
294
295 lon_grd = real(hco % lon2d_4 (lonIndex, latIndex), 8)
296 lat_grd = real(hco % lat2d_4 (lonIndex, latIndex), 8)
297 refPosition(:) = kdtree2_3dPosition(lon_grd, lat_grd)
298 call kdtree2_r_nearest(tp = tree, qv = refPosition, r2 = searchRadiusSquared, &
299 nfound = ndataFoundGridLoc(lonIndex, latIndex), &
300 nalloc = maxPointsSearch, results = searchResults)
301 if (ndataFoundGridLoc(lonIndex, latIndex) > maxPointsSearch) &
302 call utl_abort('sstb_getGriddedObs: the parameter maxPointsSearch must be increased')
303
304 if (ndataFoundGridLoc(lonIndex, latIndex) > 0) then
305 do localObsIndex = 1, ndataFoundGridLoc(lonIndex, latIndex)
306 bodyIndex = obs_headElem_i(obsData, obs_rln, headerIndexes(searchResults(localObsIndex)%idx))
307 currentObs = obs_bodyElem_r(obsData, obs_var, bodyIndex)
308 obsGrid(lonIndex, latIndex) = obsGrid(lonIndex, latIndex) + currentObs
309 end do
310 end if
311
312 end if OPENWATERPTS
313 end do
314 end do
315
316 deallocate(headerIndexes)
317 deallocate(positionArray)
318
319 end if POSITIVECOUNTOBSLOC
320
321 POSITIVECOUNTOBSGLOB: if (countObsGlob > 0) then
322 write(*,*) 'sstb_getGriddedObs: computing average values for every grid point...'
323 do latIndex = 1, hco%nj
324 do lonIndex = 1, hco%ni
325 ! summing the values over all mpi tasks and sending them back to all tasks preserving the order of summation
326 call mmpi_allreduce_sumreal8scalar(obsGrid(lonIndex, latIndex), "grid")
327 call rpn_comm_allreduce(ndataFoundGridLoc(lonIndex, latIndex), &
328 ndataFoundGridGlob(lonIndex, latIndex), 1, &
329 'mpi_integer', 'mpi_sum', 'grid', ierr)
330 if (ndataFoundGridGlob(lonIndex, latIndex) > 0) then
331 obsGrid(lonIndex, latIndex) = obsGrid(lonIndex, latIndex) / &
332 real(ndataFoundGridGlob(lonIndex, latIndex))
333 else
334 obsGrid(lonIndex, latIndex) = MPC_missingValue_R8
335 end if
336 if (.not.openWater(lonIndex, latIndex)) obsGrid(lonIndex, latIndex) = MPC_missingValue_R8
337 end do
338 end do
339 call rpn_comm_barrier('GRID', ierr)
340 write(*,*) 'sstb_getGriddedObs: gridding for '//trim(instrumentString)//' data completed'
341 end if POSITIVECOUNTOBSGLOB
342
343 end subroutine sstb_getGriddedObs
344
345 !--------------------------------------------------------------------------
346 ! sstb_getGriddedBias
347 !--------------------------------------------------------------------------
348 subroutine sstb_getGriddedBias(satelliteGrid, insituGrid, nobsGlob, nobsLoc, &
349 hco, vco, mask, openWater, sensor, &
350 numberOpenWaterPoints, dayOrNight)
351 !
352 !:Purpose: compute the satellite SST data bias estimation field on a grid
353 !
354 implicit none
355
356 ! Arguments:
357 real(8) , intent(inout) :: satelliteGrid(:,:) ! gridded satellite data
358 real(8) , intent(inout) :: insituGrid(:,:) ! gridded insitu data
359 integer , intent(in) :: nobsGlob ! number of data on all procs
360 integer , intent(in) :: nobsLoc ! number of data on the current MPI proc
361 type(struct_hco), intent(in) , pointer :: hco ! horizontal grid structure
362 type(struct_vco), intent(in) , pointer :: vco ! vertical grid structure
363 logical , intent(in) :: mask(:,:) ! land-ocean mask
364 logical , intent(in) :: openWater(:,:) ! open water points (.true.)
365 character(len=*), intent(in) :: sensor ! current sensor name
366 integer , intent(in) :: numberOpenWaterPoints! number of open water points to allocate kd-tree work arrays
367 character(len=*), intent(in) :: dayOrNight ! look for daytime or nighttime obs
368
369 ! Locals:
370 real(4), parameter :: solarZenithThreshold = 90.0 ! to distinguish day and night
371 type(kdtree2), pointer :: tree => null()
372 real(kdkind), allocatable :: positionArray(:,:)
373 integer, parameter :: maxPointsSearch = 200000
374 type(kdtree2_result) :: searchResults(maxPointsSearch)
375 real(kdkind) :: refPosition(3)
376 real(kdkind) :: lon_grd, lat_grd
377 integer :: lonIndex, latIndex
378 integer :: ierr
379 integer :: localIndex, indexCounter, localLonIndex, localLatIndex
380 integer :: numPointsFound
381 real(kdkind) :: searchRadiusSquared
382 type(struct_gsv) :: stateVector ! state vector containing bias estimation field
383 type(struct_gsv) :: stateVector_searchRadius, stateVector_previous
384 type(struct_gsv) :: stateVectorNobs, stateVectorWeight
385 real(4), pointer :: griddedBias_r4_ptr(:,:,:), searchRadius_ptr(:,:,:)
386 real(4), pointer :: nobsField_r4_ptr(:,:,:), weightField_r4_ptr(:,:,:)
387 real(4), pointer :: griddedBias_r4_previous_ptr(:,:,:)
388 integer, allocatable :: gridPointIndexes(:,:)
389 real(8) :: weight, distance, correlation, lengthscale, difference, numberPoints
390 character(len=1) :: extension
391 character(len=*), parameter :: outputFileName = './satellite_bias.fst'
392 character(len=*), parameter :: outputAuxFileName = './auxOutput.fst' ! auxiliary file name to store nobs^a and weight
393
394 write(*,*) ''
395 write(*,*) 'sstb_getGriddedBias: computing bias for: '//trim(sensor)//' '//trim(dayOrNight)
396 write(*,*) 'sstb_getGriddedBias: the current processor contains ', nobsLoc, ' data out of ', nobsGlob
397
398 if (dayOrNight == 'day') then
399 extension = 'D'
400 else if (dayOrNight == 'night') then
401 extension = 'N'
402 else
403 call utl_abort('sstb_getGriddedBias: wrong extension: '//trim(extension))
404 end if
405
406 POSITIVEOBSNUMBER: if (nobsLoc > 0) then
407
408 write(*,*) 'sstb_getGriddedBias: compute kd-tree with all the grid points...'
409 allocate(positionArray(3, numberOpenWaterPoints))
410 allocate(gridPointIndexes(2, numberOpenWaterPoints))
411
412 call lfn_setup('FifthOrder')
413
414 indexCounter = 0
415 do latIndex = 1, hco%nj
416 do lonIndex = 1, hco%ni
417 if (openWater(lonIndex, latIndex)) then
418 indexCounter = indexCounter + 1
419 lon_grd = real(hco%lon2d_4(lonIndex, latIndex), 8)
420 lat_grd = real(hco%lat2d_4(lonIndex, latIndex), 8)
421 positionArray(:, indexCounter) = kdtree2_3dPosition(lon_grd, lat_grd)
422 gridPointIndexes(1, indexCounter) = lonIndex
423 gridPointIndexes(2, indexCounter) = latIndex
424 end if
425 end do
426 end do
427
428 nullify(tree)
429 tree => kdtree2_create(positionArray, sort=.true., rearrange=.true.)
430
431 end if POSITIVEOBSNUMBER
432
433 ! get search radius field
434 call gsv_allocate(stateVector_searchRadius, 1, hco, vco, dataKind_opt = 4, &
435 datestamp_opt = -1, mpi_local_opt = .true., varNames_opt = (/'TM'/))
436 call gio_readFromFile(stateVector_searchRadius, './searchRadius', 'RADIUS','A', &
437 unitConversion_opt=.false., containsFullField_opt=.true.)
438 call gsv_getField(stateVector_searchRadius, searchRadius_ptr)
439
440 if (saveAuxFields) then
441 ! output nobs state vector
442 call gsv_allocate(stateVectorNobs, 1, hco, vco, dataKind_opt = 4, &
443 datestamp_opt = tim_getDateStamp(), mpi_local_opt = .true., &
444 varNames_opt = (/'TM'/))
445 ! pointer for nobs stateVector
446 call gsv_getField(stateVectorNobs, nobsField_r4_ptr)
447 ! output weight state vector
448 call gsv_allocate(stateVectorWeight, 1, hco, vco, dataKind_opt = 4, &
449 datestamp_opt = tim_getDateStamp(), mpi_local_opt = .true., &
450 varNames_opt = (/'TM'/))
451 ! pointer for weight stateVector
452 call gsv_getField(stateVectorWeight, weightField_r4_ptr)
453 nobsField_r4_ptr(:,:,1) = 0.0d0
454 weightField_r4_ptr(:,:,1) = 0.0d0
455 end if
456
457 ! previous bias estimation
458 call gsv_allocate(stateVector_previous, 1, hco, vco, dataKind_opt = 4, &
459 datestamp_opt = -1, mpi_local_opt = .true., varNames_opt = (/'TM'/))
460 call gio_readFromFile(stateVector_previous, './trlm_01', 'B_'//trim(sensor)//'_'//trim(extension), &
461 ' ', unitConversion_opt=.false., containsFullField_opt=.true.)
462 call gsv_getField(stateVector_previous, griddedBias_r4_previous_ptr)
463
464 ! resulting bias estimation state vector
465 call gsv_allocate(stateVector, 1, hco, vco, dataKind_opt = 4, &
466 datestamp_opt = tim_getDateStamp(), mpi_local_opt = .true., varNames_opt = (/'TM'/))
467 ! pointer for bias estimation stateVector
468 call gsv_getField(stateVector, griddedBias_r4_ptr)
469
470 griddedBias_r4_ptr(:,:,1) = 0.0d0
471
472 if (nobsLoc > 0) then
473 write(*,*) 'sstb_getGriddedBias: do the search for '//trim(sensor)//' '//trim(dayOrNight)//'...'
474 else
475 write(*,*) 'sstb_getGriddedBias: no '//trim(sensor)//' '&
476 //trim(dayOrNight)//' data on the current processor'
477 write(*,*) 'sstb_getGriddedBias: previous estimation state will be used.'
478 end if
479
480 do latIndex = myLatBeg, myLatEnd
481 do lonIndex = myLonBeg, myLonEnd
482
483 POSITIVEOBSNUMBER2: if (nobsLoc > 0) then
484
485 numberPoints = 0.0d0
486 LANDSEAMASK: if (mask(lonIndex, latIndex)) then
487 lon_grd = real(hco%lon2d_4(lonIndex, latIndex), 8)
488 lat_grd = real(hco%lat2d_4(lonIndex, latIndex), 8)
489 refPosition(:) = kdtree2_3dPosition(lon_grd, lat_grd)
490 ! convert from km to m2
491 searchRadiusSquared = (1.1d0 * searchRadius_ptr(lonIndex, latIndex, 1) * 1000.d0)**2
492 call kdtree2_r_nearest(tp = tree, qv = refPosition, &
493 r2 = searchRadiusSquared, &
494 nfound = numPointsFound, &
495 nalloc = maxPointsSearch, &
496 results = searchResults)
497 POSITIVENPTSFOUND: if (numPointsFound > 0) then
498 do localIndex = 1, numPointsFound
499 localLonIndex = gridPointIndexes(1, searchResults(localIndex)%idx)
500 localLatIndex = gridPointIndexes(2, searchResults(localIndex)%idx)
501 difference = satelliteGrid(localLonIndex, localLatIndex) - &
502 insituGrid(localLonIndex, localLatIndex)
503 if (insituGrid (localLonIndex, localLatIndex) /= MPC_missingValue_R8 .and. &
504 satelliteGrid(localLonIndex, localLatIndex) /= MPC_missingValue_R8 .and. &
505 abs(difference) < maxBias) then
506 distance = sqrt(searchResults(localIndex)%dis)
507 lengthscale = 1000.d0 * searchRadius_ptr(lonIndex, latIndex, 1)
508 correlation = lfn_response(distance, lengthscale)
509 griddedBias_r4_ptr(lonIndex, latIndex, 1) = griddedBias_r4_ptr(lonIndex, latIndex, 1) + &
510 correlation * difference
511 numberPoints = numberPoints + correlation
512 end if
513 end do
514
515 if (numberPoints > 0.0d0) &
516 griddedBias_r4_ptr(lonIndex, latIndex, 1) = griddedBias_r4_ptr(lonIndex, latIndex, 1) / numberPoints
517 end if POSITIVENPTSFOUND
518 end if LANDSEAMASK
519
520 weight = numberPoints / (numberPoints + numberPointsBG)
521 if (weight < weightMin) weight = weightMin
522 if (weight > weightMax) weight = weightMax
523
524 if (saveAuxFields) then
525 weightField_r4_ptr(lonIndex, latIndex, 1) = weight
526 nobsField_r4_ptr(lonIndex, latIndex, 1) = numberPoints
527 end if
528
529 ! computation of the bias:
530 griddedBias_r4_ptr(lonIndex, latIndex, 1) = (1.0d0 - weight) * bgTermZeroBias * &
531 griddedBias_r4_previous_ptr(lonIndex, latIndex, 1) + &
532 weight * griddedBias_r4_ptr(lonIndex, latIndex, 1)
533 else ! no data on the current processor
534
535 ! the bias estimation on the current processor is the estimation from previous state:
536 griddedBias_r4_ptr(lonIndex, latIndex, 1) = griddedBias_r4_previous_ptr(lonIndex, latIndex, 1)
537
538 end if POSITIVEOBSNUMBER2
539 end do
540 end do
541
542 call rpn_comm_barrier('GRID', ierr)
543 call gio_writeToFile(stateVector, outputFileName, 'B_'//trim(sensor)//'_'//trim(extension))
544
545 if (nobsLoc > 0) then
546 deallocate(gridPointIndexes)
547 deallocate(positionArray)
548 end if
549
550 call gsv_deallocate(stateVector_searchRadius)
551 call gsv_deallocate(stateVector_previous)
552 call gsv_deallocate(stateVector)
553 if (saveAuxFields) then
554 call gio_writeToFile(stateVectorNobs, outputAuxFileName, 'N_'//trim(sensor)//'_'//trim(extension))
555 call gio_writeToFile(stateVectorWeight, outputAuxFileName, 'W_'//trim(sensor)//'_'//trim(extension))
556 call gsv_deallocate(stateVectorNobs)
557 call gsv_deallocate(stateVectorWeight)
558 end if
559
560 write(*,*) 'sstb_getGriddedBias: completed for: '//trim(sensor)//' '//trim(dayOrNight)
561
562 end subroutine sstb_getGriddedBias
563
564 !--------------------------------------------------------------------------
565 ! sstb_getBiasCorrection
566 !--------------------------------------------------------------------------
567 subroutine sstb_getBiasCorrection(stateVector, column, obsData, hco, sensor, dayOrNight)
568 !
569 !:Purpose: To compute bias correction and put it into obsSpace data.
570 ! Columns from input field are interpolated to obs location
571 !
572 implicit none
573
574 ! Arguments:
575 type(struct_gsv) , intent(inout) :: stateVector ! state vector containing bias estimation field
576 type(struct_columnData), intent(inout) :: column ! column data
577 type(struct_obs) , intent(inout) :: obsData ! obsSpaceData
578 type(struct_hco) , intent(in), pointer :: hco ! horizontal grid
579 character(len=*) , intent(in) :: sensor ! current sensor name
580 character(len=*) , intent(in) :: dayOrNight ! look for daytime or nighttime obs
581
582 ! Locals:
583 real(4), parameter :: solarZenithThreshold = 90.0 ! to distinguish day and night
584 integer :: bodyIndex, headerIndex
585 real(8) :: currentObs
586
587 write(*,*)
588 write(*,*) 'sstb_getBiasCorrection: computing bias correction for ', trim(sensor), ' ', trim(dayOrNight), 'time ************'
589
590 call s2c_nl(stateVector, obsData, column, hco, timeInterpType = timeInterpType_nl, &
591 moveObsAtPole_opt = .true., numObsBatches_opt = numObsBatches, dealloc_opt = .false.)
592
593 do headerIndex = 1, obs_numheader(obsData)
594
595 if (obs_elem_c(obsData, 'STID' , headerIndex) == trim(sensor)) then
596
597 bodyIndex = obs_headElem_i(obsData, obs_rln, headerIndex)
598 currentObs = obs_bodyElem_r(obsData, obs_var, bodyIndex)
599
600 if (trim(dayOrNight) == 'day') then
601 if (obs_headElem_r(obsData, obs_sun, headerIndex) < solarZenithThreshold) then
602 call obs_bodySet_r(obsData, obs_bcor, bodyIndex, col_getElem(column, 1, headerIndex, 'TM'))
603 call obs_bodySet_r(obsData, obs_var , bodyIndex, currentObs - col_getElem(column, 1, headerIndex, 'TM'))
604 end if
605 else if (trim(dayOrNight) == 'night') then
606 if (obs_headElem_r(obsData, obs_sun, headerIndex) >= solarZenithThreshold) then
607 call obs_bodySet_r(obsData, obs_bcor, bodyIndex, col_getElem(column, 1, headerIndex, 'TM'))
608 call obs_bodySet_r(obsData, obs_var , bodyIndex, currentObs - col_getElem(column, 1, headerIndex, 'TM'))
609 end if
610 end if
611 end if
612
613 end do
614
615 write(*,*) 'sstb_getBiasCorrection: END'
616
617 end subroutine sstb_getBiasCorrection
618
619 !--------------------------------------------------------------------------
620 ! readNml (private subroutine)
621 !--------------------------------------------------------------------------
622 subroutine readNml()
623 !
624 !:Purpose: Read the namelist `namSSTbiasEstimate`
625 !
626 implicit none
627
628 ! Locals:
629 integer :: ierr, nulnam, sensorIndex
630
631 namelist /namSSTbiasEstimate/ searchRadius, maxBias, iceFractionThreshold, numberPointsBG, &
632 timeInterpType_nl, numObsBatches, numberSensors, sensorList, &
633 weightMin, weightMax, saveAuxFields, bgTermZeroBias
634
635 ! Setting default namelist variable values
636 searchRadius = 10.
637 maxBias = 1.
638 iceFractionThreshold = 0.6
639 numberSensors = MPC_missingValue_INT
640 numberPointsBG = 0
641 timeInterpType_nl = 'NEAREST'
642 numObsBatches = 20
643 sensorList(:) = ''
644 weightMin = 0.0
645 weightMax = 1.0
646 saveAuxFields = .False.
647 bgTermZeroBias = 1.0
648
649 ! Read the namelist
650 nulnam = 0
651 ierr = fnom( nulnam, './flnml', 'FTN+SEQ+R/O', 0 )
652 read(nulnam, nml = namSSTbiasEstimate, iostat = ierr )
653 if (ierr /= 0) call utl_abort('readNml (sstb): Error reading namelist')
654 if (mmpi_myid == 0) write(*, nml = namSSTbiasEstimate )
655 ierr = fclos( nulnam )
656
657 if (numberSensors /= MPC_missingValue_INT) then
658 call utl_abort('readNml (sstb): check namSSTbiasEstimate namelist section: numberSensors should be removed')
659 end if
660
661 numberSensors = 0
662 do sensorIndex = 1, maxNumberSensors
663 if (trim(sensorList(sensorIndex))=='') exit
664 numberSensors = numberSensors + 1
665 end do
666 if (numberSensors == 0) call utl_abort('readNml (sstb): check namSSTbiasEstimate namelist section: empty sensorList')
667
668 if (mmpi_myid == 0) then
669 write(*,*)
670 write(*,*) 'readNml (sstb): sensors to treat: '
671 do sensorIndex = 1, numberSensors
672 write(*,*) 'readNml (sstb): sensor index: ', sensorIndex, ', sensor: ', sensorList( sensorIndex )
673 end do
674 write(*,*) 'readNml (sstb): interpolation type: ', timeInterpType_nl
675 write(*,*) 'readNml (sstb): number obs batches: ', numObsBatches
676 write(*,*) 'readNml (sstb): weight limits for current bias estimate: ', weightMin, weightMax
677 end if
678
679 end subroutine readNml
680
681 !--------------------------------------------------------------------------
682 ! sstb_applySatelliteSSTBiasCorrection
683 !--------------------------------------------------------------------------
684 subroutine sstb_applySatelliteSSTBiasCorrection(obsData, hco, vco, column)
685 !
686 !:Purpose: To apply bias satellite SST data bias correction and put it into obsSpace data.
687 ! Columns from input field are interpolated to obs location
688 !
689 implicit none
690
691 ! Arguments:
692 type(struct_obs) , intent(inout) :: obsData ! obsSpaceData
693 type(struct_hco) , intent(inout), pointer :: hco ! horizontal grid structure
694 type(struct_vco) , intent(in) , pointer :: vco ! vertical grid structure
695 type(struct_columnData), intent(inout) :: column ! column data
696
697 ! Locals:
698 type(struct_gsv) :: stateVector
699 integer , parameter :: numberProducts = 2 ! day and night
700 character(len=*), parameter :: listProducts(numberProducts)= (/'day', 'night'/) ! day and night biases
701 integer :: sensorIndex, productIndex
702 character(len=1) :: extension
703 character(len=*), parameter :: biasFileName = './satellite_bias.fst'
704
705 ! read the namelist
706 call readNml()
707
708 ! allocate state vector for bias estimation field
709 call gsv_allocate(stateVector, 1, hco, vco, dataKind_opt = 4, hInterpolateDegree_opt = 'LINEAR', &
710 datestamp_opt = -1, mpi_local_opt = .true., varNames_opt = (/'TM'/))
711
712 do sensorIndex = 1, numberSensors
713 do productIndex = 1, numberProducts
714
715 if (trim(listProducts(productIndex)) == 'day') then
716 extension = 'D'
717 else if (trim(listProducts(productIndex)) == 'night') then
718 extension = 'N'
719 else
720 call utl_abort('sstb_applySatelliteSSTBiasCorrection: wrong extension: '//trim(extension))
721 end if
722
723 call gio_readFromFile(stateVector, biasFileName, 'B_'//trim(sensorList(sensorIndex))//'_'//trim(extension), &
724 ' ', unitConversion_opt = .false., containsFullField_opt = .true.)
725 call sstb_getBiasCorrection(stateVector, column, obsData, hco, trim(sensorList(sensorIndex)), &
726 trim(listProducts(productIndex)))
727 end do
728 end do
729
730 call gsv_deallocate(stateVector)
731
732 end subroutine sstb_applySatelliteSSTBiasCorrection
733
734 !--------------------------------------------------------------------------
735 ! sstb_getBiasFromPreviousState
736 !--------------------------------------------------------------------------
737 subroutine sstb_getBiasFromPreviousState(hco, vco, sensor, dayOrNight)
738 !
739 !:Purpose: to get a satellite SST data bias estimate from the previous state if data is missing.
740 ! or there are no insitu data for the current dateStamp,
741 ! hence, unable to compute bias estimates for the satellite data.
742 !
743 implicit none
744
745 ! Arguments:
746 type(struct_hco), intent(inout), pointer :: hco ! horizontal grid structure
747 type(struct_vco), intent(in) , pointer :: vco ! vertical grid structure
748 character(len=*), intent(in) :: sensor ! sensor name
749 character(len=*), intent(in) :: dayOrNight ! look for daytime or nighttime bias estimation
750
751 ! Locals:
752 type(struct_gsv) :: stateVector ! state vector containing current bias estimation field
753 type(struct_gsv) :: stateVector_previous ! state vector containing previous bias estimation field
754 real(4), pointer :: griddedBias_r4_ptr(:, :, :)
755 real(4), pointer :: griddedBias_r4_previous_ptr(:, :, :)
756 character(len=1) :: extension
757 character(len=*), parameter :: outputFileName = './satellite_bias.fst'
758
759 write(*,*) ''
760 write(*,*) 'sstb_getBiasFromPreviousState: for ', trim(sensor), ' ', trim(dayOrNight), ' data...'
761
762 if (trim(dayOrNight) == 'day') then
763 extension = 'D'
764 else if (trim(dayOrNight) == 'night') then
765 extension = 'N'
766 else
767 call utl_abort('sstb_getBiasFromPreviousState: wrong extension: '//trim(extension))
768 end if
769
770 ! read previous bias estimation
771 call gsv_allocate(stateVector_previous, 1, hco, vco, dataKind_opt = 4, &
772 datestamp_opt = -1, mpi_local_opt = .true., varNames_opt = (/'TM'/))
773 call gio_readFromFile(stateVector_previous, './trlm_01', 'B_'//trim(sensor)//'_'//trim(extension), &
774 ' ', unitConversion_opt=.false., containsFullField_opt=.true.)
775 call gsv_getField(stateVector_previous, griddedBias_r4_previous_ptr)
776
777 ! resulting bias estimation state vector
778 call gsv_allocate(stateVector, 1, hco, vco, dataKind_opt = 4, &
779 datestamp_opt = tim_getDateStamp(), mpi_local_opt = .true., varNames_opt = (/'TM'/))
780 ! pointer for bias estimation stateVector
781 call gsv_getField(stateVector, griddedBias_r4_ptr)
782
783 griddedBias_r4_ptr(:, :, :) = griddedBias_r4_previous_ptr(:, :, :)
784 call gio_writeToFile(stateVector, outputFileName, 'B_'//trim(sensor)//'_'//trim(extension))
785
786 call gsv_deallocate(stateVector)
787 call gsv_deallocate(stateVector_previous)
788
789 write(*,*) 'sstb_getBiasFromPreviousState: done ', trim(sensor), ' ', trim(dayOrNight), ' data.'
790
791 end subroutine sstb_getBiasFromPreviousState
792
793end module SSTbias_mod
